C++ Pointer::GetBufferPointer方法代码示例

void LabelsVolumeGenerator::postProcessTargetITK()
{
    typedef unsigned char       PixelType;
    typedef Image<PixelType, 3> ImageType;

    // create a new image from the target data

    ImageType::Pointer input = ImageType::New();

    ImageType::IndexType start;
    start.Fill(0);

    ImageType::SizeType size;
    for (int i = 0; i < 3; ++i) size[i] = m_targetVolume->dimensions[i];

    ImageType::RegionType region;
    region.SetSize(size);
    region.SetIndex(start);

    input->SetRegions(region);
    input->SetSpacing(m_targetVolume->spacing.data());
    input->Allocate();

    memcpy(input->GetBufferPointer(), m_targetVolume->data, m_bufferSize);

    // create a grayscale dilation filter and a structuring element

    typedef BinaryBallStructuringElement<PixelType, 3> StructureType;
    typedef GrayscaleDilateImageFilter<ImageType, ImageType, StructureType> DilateFilterType;

    DilateFilterType::Pointer filter = DilateFilterType::New();

    StructureType structure;
    structure.SetRadius(1);

    filter->SetKernel(structure);

    // set up progress reporting
    if (m_progressReporter)
    {
        CStyleCommand::Pointer command = CStyleCommand::New();
        command->SetClientData(m_progressReporter);
        command->SetCallback(ProgressCallback);
        filter->AddObserver(ProgressEvent(), command);
        m_progressReporter->start("Post-Processing Label volume",
                                  "Dilating label field...");
    }

    // hook up the filters and run

    filter->SetInput(input);
    filter->Update();

    if (m_progressReporter)
        m_progressReporter->finish();

    // copy back into the target volume data
    memcpy(m_targetVolume->data, filter->GetOutput()->GetBufferPointer(), m_bufferSize);

    // threshold and gaussian blur to put a smooth version into the alpha channel

    typedef BinaryThresholdImageFilter<ImageType, ImageType> ThresholderType;
//    typedef DiscreteGaussianImageFilter<ImageType, ImageType> SmoothFilterType;
    typedef SmoothingRecursiveGaussianImageFilter<ImageType, ImageType> SmoothFilterType;

    ThresholderType::Pointer thresholder = ThresholderType::New();
    thresholder->SetLowerThreshold(1);
    thresholder->SetUpperThreshold(255);
    thresholder->SetInsideValue(255);
    thresholder->SetOutsideValue(0);

    SmoothFilterType::Pointer smoother = SmoothFilterType::New();
//    smoother->SetVariance(0.05); // in physical units
//    smoother->SetMaximumKernelWidth(5);
    smoother->SetSigma(0.2);

    // set up progress reporting (again)
    if (m_progressReporter)
    {
        CStyleCommand::Pointer command = CStyleCommand::New();
        command->SetClientData(m_progressReporter);
        command->SetCallback(ProgressCallback);
        smoother->AddObserver(ProgressEvent(), command);
        m_progressReporter->start("Post-Processing Label volume",
                                  "Smoothing alpha mask...");
    }

    // hook up the filters and run

    thresholder->SetInput(input);
    smoother->SetInput(thresholder->GetOutput());
    smoother->Update();

    // copy back into the target volume data
    memcpy(m_targetMask->data, smoother->GetOutput()->GetBufferPointer(), m_bufferSize);

    if (m_progressReporter)
        m_progressReporter->finish();
}

//.........这里部分代码省略.........
		else
		  v = ((ushort*)tmp)[j*m_nZ + k];
		ok = checkPathBlend(po, v, lut);
	      }
	    
	    if (ok)
	      cropped[jk] = mop;
	    else
	      cropped[jk] = 0;
	    
	    jk ++;
	  }
      
      if (m_voxelType == 0)
	{
	  for(int j=0; j<m_nY*m_nZ; j++)
	    {
	      if (cropped[j] == 0)
		tmp[j] = 0;
	    }
	}
      else
	{
	  for(int j=0; j<m_nY*m_nZ; j++)
	    {
	      if (cropped[j] == 0)
		((ushort*)tmp)[j] = 0;
	    }
	}
      
      applyOpacity(iv, cropped, lut, tmp);
      memcpy(opacityVol + i0*m_nY*m_nZ, tmp, m_nY*m_nZ);
      
      i0++;
    }
  delete [] tmp;
  delete [] cropped;
  m_meshProgress->setValue(100);
  qApp->processEvents();

  //------------
  if (m_tearPresent)
    {
      uchar *data0 = new uchar[m_nX*m_nY*m_nZ];
      memcpy(data0, opacityVol, m_nX*m_nY*m_nZ);
      applyTear(d0, d1, 0,
		data0, opacityVol, false);
      
      delete [] data0;
    }


  typedef uchar PixelType;
  const unsigned int Dimension = 3;
  typedef itk::Image< PixelType, Dimension > ImageType;

  ImageType::IndexType start;
  start.Fill(0);

  ImageType::SizeType size;
  size[0] = m_nZ;
  size[1] = m_nY;
  size[2] = m_nX;

  ImageType::RegionType region(start, size);

  ImageType::Pointer image = ImageType::New();
  image->SetRegions(region);
  image->Allocate();
  image->FillBuffer(0);
  uchar *iptr = (uchar*)image->GetBufferPointer();
  memcpy(iptr, opacityVol, m_nX*m_nY*m_nZ);

  typedef itk::Image< float, 3 > OutputImageType;
  typedef itk::SignedMaurerDistanceMapImageFilter<ImageType, OutputImageType> DistanceMapFilter;
  DistanceMapFilter::Pointer filter = DistanceMapFilter::New();
  filter->SetInput( image );
  filter->SetSquaredDistance(0);
  filter->SetUseImageSpacing(0);
  filter->SetInsideIsPositive(1);
  filter->Update();

  QFile fp;
  fp.setFileName(flnm);
  fp.open(QFile::WriteOnly);
  uchar vt = 8;
  fp.write((char*)&vt, 1);
  fp.write((char*)&m_nX, 4);
  fp.write((char*)&m_nY, 4);
  fp.write((char*)&m_nZ, 4);
  OutputImageType *dimg = filter->GetOutput();
  char *tdata = (char*)(dimg->GetBufferPointer());
  fp.write(tdata, 4*m_nX*m_nY*m_nZ);
  fp.close();

  m_meshLog->moveCursor(QTextCursor::End);
  m_meshLog->insertPlainText("Signed Distance Map data saved in "+flnm);

  QMessageBox::information(0, "", QString("Signed Distance Map saved in "+flnm));
}

int main(int argc, char**argv)
{
	int nx, ny, nz, nbytes, direction, nx8, ny8, nz8, k;
//	int n, x, y, z, dx, dy;
	long long width, height, depth, xysize;
	char *dataFile, *tiffFile;
	bool use_compression = true, data2tiff, compressdata = true;
	FILE *fpdata;

	if (argc != 4) {
		printf("Usage: maketiff dataFile tiffFile direction\n");
		printf("       direction = 0 (tiff to data) = 1 (data to tiff)\n");
		return 1;
	}

	dataFile = argv[1];
	tiffFile = argv[2];
	sscanf(argv[3],"%d",&direction);
	if (direction == 0) {
		data2tiff = false;
		printf("Output data file: %s\n",dataFile);
	} else {
		data2tiff = true;
		printf("Output image file: %s\n",tiffFile);
	}

	if (data2tiff) {		// create tiff file
		fpdata = fopen(dataFile,"rb");
		if (fpdata == NULL) {
			return 2;
		}
		fread(&nx,4,1,fpdata);
		fread(&ny,4,1,fpdata);
		fread(&nz,4,1,fpdata);
//		fread(&nbytes,4,1,fpdata);

		//if (nbytes != nx*ny*nz) {
		//	printf("Error: this is a compressed data file\n");
		//	return 10;
		//}
		width = nx;
		height = ny;
		depth = nz;
		xysize = width*height;
		printf("Desired image dimensions: width, height, depth: %d %d %d\n",width,height,depth);

		ImageType::Pointer im = ImageType::New();
		ImageType::SizeType imsize; 
		imsize[0] = width;
		imsize[1] = height;
		imsize[2] = depth;
		ImageType::IndexType imstart; 
		imstart[0] = 0;
		imstart[1] = 0;
		imstart[2] = 0;
		ImageType::RegionType imregion; 
		imregion.SetSize(imsize);
		imregion.SetIndex(imstart);
		im->SetRegions(imregion);
		im->Allocate();
		p = (unsigned char *)(im->GetBufferPointer());

		nbytes = nx*ny*nz;
		fread(p,nbytes,1,fpdata);

		typedef itk::ImageFileWriter<ImageType> FileWriterType;
		FileWriterType::Pointer writer = FileWriterType::New();
		writer->SetFileName(tiffFile);
		writer->SetInput(im);
		if (use_compression) {
			writer->UseCompressionOn();
		}
		try
		{
			writer->Update();
		}
		catch (itk::ExceptionObject &e)
		{
			std::cout << e << std::endl;
			return 3;
		}
		if (use_compression) {
			printf("Created compressed image file: %s\n",tiffFile);
		} else {
			printf("Created uncompressed image file: %s\n",tiffFile);
		}
	} else {		// create data file
		typedef itk::ImageFileReader<ImageType> FileReaderType;
		FileReaderType::Pointer reader = FileReaderType::New();

		reader->SetFileName(tiffFile);
		try
		{
			reader->Update();
		}
		catch (itk::ExceptionObject &e)
		{
			std::cout << e << std::endl;
			printf("Read error on input file\n");
			return 2;	// Read error on input tiff file
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    itk::ImageRegionIterator<ImageType>
        it2(comp2Image, comp2Image->GetLargestPossibleRegion());

    it1.GoToBegin();
    it2.GoToBegin();
    itprob.GoToBegin();

    while( !itprob.IsAtEnd() )
    {
      if(itprob.Get() > 0.2 * maxProb)
      {
        MeasurementVectorType mv;
        mv[0] = ( MeasurementType ) it1.Get();
        mv[1] = ( MeasurementType ) it2.Get();
        listSample->PushBack(mv);
      }
      ++it1;
      ++it2;
      ++itprob;
    }

    // generate a histogram from the list sample
    typedef float HistogramMeasurementType;
    typedef itk::Statistics::Histogram< HistogramMeasurementType, itk::Statistics::DenseFrequencyContainer2 > HistogramType;
    typedef itk::Statistics::SampleToHistogramFilter< ListSampleType, HistogramType > GeneratorType;
    GeneratorType::Pointer generator = GeneratorType::New();

    GeneratorType::HistogramType::SizeType size(2);
    size.Fill(30);
    generator->SetHistogramSize( size );

    generator->SetInput( listSample );
    generator->SetMarginalScale( 10.0 );

    generator->Update();

    // look for frequency mode in the histogram
    GeneratorType::HistogramType::ConstPointer histogram = generator->GetOutput();
    GeneratorType::HistogramType::ConstIterator iter = histogram->Begin();
    float maxFreq = 0;
    MeasurementVectorType maxValue;
    maxValue.Fill(0);
    while ( iter != histogram->End() )
    {
      if(iter.GetFrequency() > maxFreq)
      {
        maxFreq = iter.GetFrequency();
        maxValue[0] = iter.GetMeasurementVector()[0];
        maxValue[1] = iter.GetMeasurementVector()[1];
      }
      ++iter;
    }

    // generate return image that contains the angular
    // error of the voxels to the histogram max measurement
    ImageType::Pointer returnImage = ImageType::New();
    returnImage->SetSpacing( comp1Image->GetSpacing() );   // Set the image spacing
    returnImage->SetOrigin( comp1Image->GetOrigin() );     // Set the image origin
    returnImage->SetDirection( comp1Image->GetDirection() );  // Set the image direction
    returnImage->SetRegions( comp1Image->GetLargestPossibleRegion() );
    returnImage->Allocate();

    itk::ImageRegionConstIterator<ImageType>
        cit1(comp1Image, comp1Image->GetLargestPossibleRegion());

    itk::ImageRegionConstIterator<ImageType>
        cit2(comp2Image, comp2Image->GetLargestPossibleRegion());

    itk::ImageRegionIterator<ImageType>
        itout(returnImage, returnImage->GetLargestPossibleRegion());

    cit1.GoToBegin();
    cit2.GoToBegin();
    itout.GoToBegin();

    vnl_vector<float> v(3);
    v[0] = cos( maxValue[0] ) * sin( maxValue[1] );
    v[1] = sin( maxValue[0] ) * sin( maxValue[1] );
    v[2] = cos( maxValue[1] );
//    MITK_INFO << "max vector: " << v;
    while( !cit1.IsAtEnd() )
    {
      vnl_vector<float> v1(3);
      v1[0] = cos( cit1.Get() ) * sin( cit2.Get() );
      v1[1] = sin( cit1.Get() ) * sin( cit2.Get() );
      v1[2] = cos( cit2.Get() );

      itout.Set(fabs(angle(v,v1)));
//      MITK_INFO << "ang_error " << v1 << ": " << fabs(angle(v,v1));

      ++cit1;
      ++cit2;
      ++itout;
    }

    mitk::Image::Pointer retval = mitk::Image::New();
    retval->InitializeByItk(returnImage.GetPointer());
    retval->SetVolume(returnImage->GetBufferPointer());
    return retval;
  }