C++ imagetype::RegionType类代码示例

MriWatcherGUI::ImageType::Pointer MriWatcherGUI::CreateNewImage(int sizex, int sizey, int sizez, float dimx, float dimy,
                                                                float dimz)
{
  int imagesize[3];

  imagesize[0] = sizex;
  imagesize[1] = sizey;
  imagesize[2] = sizez;

  ImageType::Pointer m_outputimage = ImageType::New();
  float              values[3];

  values[0] = dimx;
  values[1] = dimy;
  values[2] = dimz;

  float origin_x = ( (imagesize[0] / 2) * values[0] * (-1) );
  float origin_y = ( (imagesize[1] / 2) * values[1] * (-1) );
  float origin_z = ( (imagesize[0] / 2) * values[2] * (-1) );
  float origin[3] = {origin_x, origin_y, origin_z};

  ImageType::RegionType region;
  ImageType::SizeType   size;
  size[0] = imagesize[0];
  size[1] = imagesize[1];
  size[2] = imagesize[2];
  region.SetSize( size );
  m_outputimage->SetRegions( region );
  m_outputimage->Allocate();
  m_outputimage->SetOrigin(origin);
  m_outputimage->SetSpacing(values);

  return m_outputimage;
}

bool VolumeProcess::RunBinaryForDistanceMapUsingGraphCuts()
{
	ImageType::RegionType region = m_outputImage->GetBufferedRegion();
	int numColumns = region.GetSize(0);
	int numRows = region.GetSize(1);
	int numStacks = region.GetSize(2);
	long numPix = numStacks*numColumns*numRows;

	unsigned short * binImagePtr = new unsigned short[numPix];
	PixelType * dataImagePtr = m_outputImage->GetBufferPointer();

	//int ok = Cell_Binarization_3D(dataImagePtr, binImagePtr, numRows, numColumns, numStacks, 0, 1);	//Do Binarization
    int ok = Neuron_Binarization_3D(dataImagePtr, binImagePtr, numRows, numColumns, numStacks, 0, 1);
	if(!ok)
		return false;

	// background is bright and forground is dark:
	for(long i=0; i<numPix; ++i)
	{
		if( binImagePtr[i] == 0 )
		{
			dataImagePtr[i] = 255;
		}
		else 
            dataImagePtr[i] = 0;
	}

	delete [] binImagePtr;

	return true;
}

void CreateKernel(ImageType::Pointer kernel, unsigned int width)
{
  ImageType::IndexType start;
  start.Fill(0);
 
  ImageType::SizeType size;
  size.Fill(width);
 
  ImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);
 
  kernel->SetRegions(region);
  kernel->Allocate();
 
  itk::ImageRegionIterator<ImageType> imageIterator(kernel, region);
 
   while(!imageIterator.IsAtEnd())
    {
    //imageIterator.Set(255);
    imageIterator.Set(1);
 
    ++imageIterator;
    }
}

ImageType::Pointer TDimImage::toItkImage( const unsigned int &channel ) const {

  int sample = std::min<int>( this->samples()-1, channel); 
  int bitdepth = sizeof(PixelType);
  D_DataFormat pf = D_FMT_UNSIGNED;
  if (!std::numeric_limits<PixelType>::is_integer) 
    pf = D_FMT_FLOAT;
  else
  if (std::numeric_limits<PixelType>::is_signed) 
    pf = D_FMT_SIGNED;
  
  TDimImage img;
  if (this->depth()!=bitdepth || this->pixelType()!=pf) {
    img = this->ensureTypedDepth();
    img = img.convertToDepth( bitdepth, DimLut::ltLinearDataRange, pf );
  } else
    img = *this;

  //------------------------------------------------------------------
  // Create Itk Image and copy data
  //------------------------------------------------------------------
  ImageType::Pointer image = ImageType::New();
  ImageType::SizeType size;
  size[0] = nx;
  size[1] = ny;

  ImageType::IndexType start;
  start[0] = 0;
  start[1] = 0;  

  ImageType::RegionType region;
  region.SetSize( size );
  region.SetIndex( start );
  image->SetRegions( region );

  image->Allocate();

  double spacing[2];
  spacing[0] = 1.0;
  spacing[1] = 1.0;
  image->SetSpacing( spacing );

  double origin[2];
  origin[0] = 0.0;
  origin[1] = 0.0;
  image->SetOrigin( origin );

  typedef itk::ImageRegionIterator< ImageType >  IteratorType;
  IteratorType it( image, region );
  it.GoToBegin();

  // copy data
  PixelType *data = (PixelType *) img.sampleBits(sample);
  while( ! it.IsAtEnd() ) {
    it.Set( *data );
    ++it;
    ++data;
  }
  return image;
}

int main()
{
	std::string votesfilename = "../../half_vessel_votes2.txt";
	std::string outputfilename = "half_vessel_votes2.tif";
	FILE * fp = fopen(votesfilename.c_str(), "r");


	int x,y,z,v;
	int maxx = -1,maxy = -1, maxz = -1;
	while(fscanf(fp,"%d %d %d %d",&x,&y,&z,&v)>0)
	{
		maxx = MAX(maxx, x);
		maxy = MAX(maxy, y);
		maxz = MAX(maxz, z);
	}

	printf("maxx = %d maxy = %d maxz = %d\n",maxx,maxy,maxz);
	fclose(fp);


	fp = fopen(votesfilename.c_str(), "r");


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

	ImageType::SizeType size;
	ImageType::IndexType index;
	ImageType::RegionType region;

	size[0] = maxx+1;
	size[1] = maxy+1;
	size[2] = maxz+1;

	index[0] = index[1] = index[2] = 0;

	region.SetIndex(index);
	region.SetSize(size);
	im->SetRegions(region);
	im->Allocate();

	while(fscanf(fp,"%d %d %d %d",&x,&y,&z,&v)>0)
	{
		index[0] = x;
		index[1] = y;
		index[2] = z;
		im->SetPixel(index,v);
	}

	//	scanf("%*d");

	fclose(fp);

	FileWriterType::Pointer writer = FileWriterType::New();
	writer->SetFileName(outputfilename.c_str());
	writer->SetInput(im);
	writer->Update();


	return 0;
}

// ------------------------------------------------------------------------
void buildOutput(const SeriesTransform::List &series, ImageType::Pointer &outputImage,
		ImageType::Pointer &outputLabel, const unsigned int & timestep)
{
	// get the output parameters
	unsigned int slices = series.size();
	unsigned int timeSteps = series.front().images.size();
	ImageType::Pointer ref = series.front().images.front();




	
	ImageType::SpacingType spacing = ref->GetSpacing();
	spacing[2] = series.front().sliceThickness;
	ImageType::DirectionType direction = ref->GetDirection();
	ImageType::PointType origin = ref->GetOrigin();
	ImageType::RegionType region = ref->GetLargestPossibleRegion();
	region.SetSize(2,slices);


	// create the outputs
	outputImage->SetSpacing(spacing);
	outputImage->SetDirection(direction);
	outputImage->SetOrigin(origin);
	outputImage->SetRegions(region);
	outputImage->Allocate();
	

	outputLabel->SetSpacing(spacing);
	outputLabel->SetDirection(direction);
	outputLabel->SetOrigin(origin);
	outputLabel->SetRegions(region);
	outputLabel->Allocate();

	itk::ImageRegionIterator<ImageType> outLIt(outputLabel, outputLabel->GetLargestPossibleRegion());
	itk::ImageRegionIterator<ImageType> outImIt(outputImage, outputImage->GetLargestPossibleRegion());



	// loop through the slices
	for(unsigned int i = 0; i < slices; i++)
	{
		ImageType::Pointer im = series[i].images[timestep];
		ImageType::Pointer label = series[i].labelImages[timestep];

		itk::ImageRegionConstIterator<ImageType> imIt(im, im->GetLargestPossibleRegion());
		itk::ImageRegionConstIterator<ImageType> lIt(label, label->GetLargestPossibleRegion());

		while(!imIt.IsAtEnd())
		{
			outLIt.Set(lIt.Get());
			outImIt.Set(imIt.Get());

			++imIt; ++lIt;
			++outLIt; ++outImIt;
		}
	}
}

  void setUp()
  {
    typedef itk::Image<double, 3> ImageType;
    typedef itk::VectorImage<double, 3> VectorImageType;
    typedef itk::ImageRegionIterator<ImageType> ImageIteratorType;
    typedef itk::ImageDuplicator<ImageType> DuplicatorType;
    typedef itk::ComposeImageFilter<ImageType> CompositeFilterType;

    // generate two images with one component
    ImageType::Pointer imageComponent1 = itk::Image<double, 3>::New();
    ImageType::IndexType start;
    start.Fill(0);
    ImageType::SizeType size;
    size.Fill(5);
    ImageType::RegionType region;
    region.SetSize(size);
    region.SetIndex(start);
    imageComponent1->SetRegions(region);
    imageComponent1->Allocate();

    DuplicatorType::Pointer duplicator = DuplicatorType::New();
    duplicator->SetInputImage(imageComponent1);
    duplicator->Update();
    ImageType::Pointer imageComponent2 = duplicator->GetOutput();

    // give them differing data
    ImageIteratorType iterator1(imageComponent1, imageComponent1->GetLargestPossibleRegion());
    iterator1.GoToBegin();
    int i = 0;
    while (!iterator1.IsAtEnd())
    {
      iterator1.Set((double)i);
      ++iterator1;
      ++i;
    }

    ImageIteratorType iterator2(imageComponent2, imageComponent2->GetLargestPossibleRegion());
    iterator2.GoToBegin();
    i = 2000;
    while (!iterator2.IsAtEnd())
    {
      iterator2.Set((double)i);
      ++iterator2;
      ++i;
    }

    // copy into single VectorImage
    CompositeFilterType::Pointer compositeFilter = CompositeFilterType::New();
    compositeFilter->SetInput(0, imageComponent1);
    compositeFilter->SetInput(1, imageComponent2);
    compositeFilter->Update();
    itk::VectorImage<double, 3>::Pointer multiComponentImage = compositeFilter->GetOutput();

    // cast images to mitk
    mitk::CastToMitkImage(multiComponentImage, m_mitkMultiComponentImage);
    mitk::CastToMitkImage(imageComponent1, m_mitkImageComponent1);
    mitk::CastToMitkImage(imageComponent2, m_mitkImageComponent2);
  }

  /*
   * random a voxel. define plane through this voxel. reslice at the plane. compare the pixel vaues of the voxel
   * in the volume with the pixel value in the resliced image.
   * there are some indice shifting problems which causes the test to fail for oblique planes. seems like the chosen
   * worldcoordinate is not corrresponding to the index in the 2D image. and so the pixel values are not the same as
   * expected.
   */
  static void PixelvalueBasedTest()
  {
    /* setup itk image */
    typedef itk::Image<unsigned short, 3> ImageType;

    typedef itk::ImageRegionConstIterator< ImageType > ImageIterator;

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

    ImageType::IndexType start;
    start[0] = start[1] = start[2] = 0;

    ImageType::SizeType size;
    size[0] = size[1] = size[2] = 32;

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

    image->SetRegions(imgRegion);
    image->SetSpacing(1.0);
    image->Allocate();


    ImageIterator imageIterator( image, image->GetLargestPossibleRegion() );
    imageIterator.GoToBegin();


    unsigned short pixelValue = 0;

    //fill the image with distinct values
    while ( !imageIterator.IsAtEnd() )
    {
      image->SetPixel(imageIterator.GetIndex(), pixelValue);
      ++imageIterator;
      ++pixelValue;
    }
    /* end setup itk image */



    mitk::Image::Pointer imageInMitk;
    CastToMitkImage(image, imageInMitk);



    /*mitk::ImageWriter::Pointer writer = mitk::ImageWriter::New();
    writer->SetInput(imageInMitk);
    std::string file = "C:\\Users\\schroedt\\Desktop\\cube.nrrd";
    writer->SetFileName(file);
    writer->Update();*/

                PixelvalueBasedTestByPlane(imageInMitk, mitk::PlaneGeometry::Frontal);
                PixelvalueBasedTestByPlane(imageInMitk, mitk::PlaneGeometry::Sagittal);
                PixelvalueBasedTestByPlane(imageInMitk, mitk::PlaneGeometry::Axial);

  }

// ------------------------------------------------------------------------
void buildShortAxisVolume(const SeriesTransform::Map &transforms, 
		const unsigned int seriesNumber, ImageType::Pointer &saVolume)
{
	// get the short axis transforms
	SeriesTransform::Map::const_iterator mapIt = transforms.begin();
	std::vector<SeriesTransform> saSlices;
	while(mapIt != transforms.end())
	{
		if(mapIt->second.series == seriesNumber)
		{
			unsigned int sliceNum = mapIt->second.slice;

			if(saSlices.size() < (sliceNum+1))
				saSlices.resize(sliceNum+1);

			saSlices[sliceNum] = mapIt->second;
		}
		
		++mapIt;
	}


	// get the dimensions of the output image
	ImageType::Pointer reference = saSlices[0].images[0];
	unsigned int x,y,z;
	x = reference->GetLargestPossibleRegion().GetSize()[0];
	y = reference->GetLargestPossibleRegion().GetSize()[1];
	z = saSlices.size();

	ImageType::RegionType region;
	ImageType::SizeType size;
	ImageType::IndexType index;
	size[0] = x;
	size[1] = y;
	size[2] = z;

	index.Fill(0);

	region.SetSize(size);
	region.SetIndex(index);
	
	// get the other parts
	ImageType::SpacingType spacing = reference->GetSpacing();
	spacing[2] = saSlices[0].sliceThickness;
	ImageType::DirectionType direction = reference->GetDirection();
	ImageType::PointType origin = reference->GetOrigin();


	saVolume->SetOrigin(origin);
	saVolume->SetDirection(direction);
	saVolume->SetSpacing(spacing);
	saVolume->SetRegions(region);
	saVolume->Allocate();
	saVolume->FillBuffer(0);
	
}

int main(int argc,char ** argv){

	typedef itk::Image<float,3> ImageType;
	typedef itk::Image<float,2> SliceType;

	typedef itk::ImageFileReader<ImageType> ReaderType;
	typedef itk::ImageFileWriter<SliceType> WriterType;

	typedef itk::BoundedReciprocalImageFilter<ImageType,ImageType> BoundedReciprocalType;
	BoundedReciprocalType::Pointer boundedReciprocal = BoundedReciprocalType::New();


	typedef ttt::AdvectiveDiffusion2DIterationImageFilter<SliceType,SliceType> AdvectionDiffusion2DIterationType;

	ReaderType::Pointer reader = ReaderType::New();
	reader->SetFileName(argv[1]);
	reader->UpdateOutputInformation();


	typedef itk::ExtractImageFilter<ImageType,SliceType> ExtractorType;

	ExtractorType::Pointer extractor = ExtractorType::New();

	ImageType::RegionType extractionRegion = reader->GetOutput()->GetLargestPossibleRegion();
	extractionRegion.SetSize(2,0);

	boundedReciprocal->SetInput(reader->GetOutput());
	extractor->SetInput(boundedReciprocal->GetOutput());
	extractor->SetExtractionRegion(extractionRegion);
	extractor->SetDirectionCollapseToIdentity();

	AdvectionDiffusion2DIterationType::Pointer advectionDiffusionIteration =AdvectionDiffusion2DIterationType::New();

	advectionDiffusionIteration->SetInput(extractor->GetOutput());
	advectionDiffusionIteration->SetNumberOfThreads(1.0);


	WriterType::Pointer sliceWriter = WriterType::New();
	sliceWriter->SetInput(extractor->GetOutput());
	sliceWriter->SetFileName(argv[2]);
	sliceWriter->Update();


	WriterType::Pointer writer = WriterType::New();
	writer->SetFileName(argv[3]);

	writer->SetInput(advectionDiffusionIteration->GetOutput());
	writer->Update();

}

//get number of slices
int CTImageTreeItem::getNumberOfSlices() const {
	ImageType::Pointer itkImage = peekITKImage();
	if (itkImage.IsNotNull()) {
		ImageType::RegionType imageRegion = itkImage->GetLargestPossibleRegion();
		return static_cast<uint>(imageRegion.GetSize(2));
	} else {
		int num = m_fnList.size();
		if (num > 1) 
			return num;
		else {
			std::string t;
			itk::ExposeMetaData( m_dict, getNumberOfFramesTag(), t );
			std::istringstream buffer(t);
			buffer >> num;
			return num;
		}  
	}
}

bool VolumeProcess:: RunFillingZeroOnBouandary(int Bx,int By,int Bz)
{
	ImageType::Pointer tempImg = ImageType::New();
	tempImg->SetRegions( m_outputImage->GetLargestPossibleRegion() );
	tempImg->Allocate();
	tempImg->FillBuffer(0);

	ImageType::RegionType fullRegion = m_outputImage->GetBufferedRegion();

	int numColumns = fullRegion.GetSize(0);
	int numRows = fullRegion.GetSize(1);
	int numStacks = fullRegion.GetSize(2);
	int i, j,k; 

	itk::ImageRegionIteratorWithIndex< ImageType > itr( m_outputImage, fullRegion );

	for(itr.GoToBegin(); !itr.IsAtEnd(); ++itr)
	{
		
		ImageType::IndexType index = itr.GetIndex();
	    ImageType::PixelType pix = m_outputImage->GetPixel(index);
        i = index[0];
		j = index[1];
		k = index[2];
		if (i < Bx || i > numColumns -Bx || j < By || j > numRows-By ||k <Bz ||k > numStacks-Bz )
		    tempImg->SetPixel(itr.GetIndex(), 0);
		else 
			tempImg->SetPixel(itr.GetIndex(), pix);
     }
		//Copy temp img back to image 
	itk::ImageRegionIterator< ImageType > itr1( tempImg, tempImg->GetLargestPossibleRegion());
	itk::ImageRegionIterator< ImageType > itr2( m_outputImage, m_outputImage->GetLargestPossibleRegion());
	for(itr1.GoToBegin(), itr2.GoToBegin() ; !itr1.IsAtEnd(); ++itr1, ++itr2)
	{
	  itr2.Set( itr1.Get() );
	}

	if(debug)
		std::cerr << "RunFillingZero Done" << std::endl;

	return true;
  
}

void ITKImplicit::initScalingAndOriginFromBoundingBox(const gmVector3 & minV, const gmVector3 & maxV)
{
	//we make the bounding box slightly bigger
	
	gmVector3 originalCenter=(minV+maxV)/2;
	gmVector3 newMinV=(minV-originalCenter)*boundingBoxScaling+originalCenter;
	gmVector3 newMaxV=(maxV-originalCenter)*boundingBoxScaling+originalCenter;

	//for the moment we will use a scaling, 
	//that leads to a unit cube
	ImageType::SpacingType spacing;
	gmVector3 lengthV=newMaxV-newMinV;
	spacing[0]=fabs(lengthV[0])/((double)voxelRes[0]);
	spacing[1]=fabs(lengthV[1])/((double)voxelRes[1]);
	spacing[2]=fabs(lengthV[2])/((double)voxelRes[2]);
	myImage->SetSpacing(spacing);
	
	//we set the origin such that center lies at (0,0,0)
	ImageType::PointType origin;
	origin[0]=newMinV[0];
	origin[1]=newMinV[1];
	origin[2]=newMinV[2];
	myImage->SetOrigin(origin);
	
	//we want an image with voxelRes voxels
	ImageType::IndexType start;
	start[0]=0;
	start[1]=0;
	start[2]=0;
	ImageType::SizeType size;
	size[0]=max(1,voxelRes[0]);
	size[1]=max(1,voxelRes[1]);
	size[2]=max(1,voxelRes[2]);
	ImageType::RegionType region;
	region.SetIndex(start);
	region.SetSize(size);
	myImage->SetRegions(region);
	
	//allocate memory
	myImage->Allocate();
}

bool VolumeProcess::RunDistanceTransform(void)
{
	ImageType::RegionType region = m_outputImage->GetBufferedRegion();
	int numColumns = region.GetSize(0);
	int numRows = region.GetSize(1);
	int numStacks = region.GetSize(2);
	long numPix = numStacks*numColumns*numRows;

	PixelType * dataImagePtr = m_outputImage->GetBufferPointer();
	unsigned char * binImagePtr = new unsigned char[numPix];

	itk::ImageRegionIterator< ImageType > itr( m_outputImage, m_outputImage->GetLargestPossibleRegion() );

    long idx = 0;
	for(itr.GoToBegin(); !itr.IsAtEnd(); ++itr)
	{
		binImagePtr[idx++]=itr.Get();
	}

	bool ok= distTransform(binImagePtr, numRows, numColumns, numStacks);	//Do distance transform
	if(!ok)
		return false;

	if(debug)
	{
		std::cerr << " Distance Transform is done and we will replace image with the value of DT!" << std::endl;
	}
	//repalce the image by the distance transform
	
	idx = 0;
	for(itr.GoToBegin(); !itr.IsAtEnd(); ++itr)
	{
		itr.Set( binImagePtr[idx++] );
	}

	delete [] binImagePtr;
	return true;
}

void TDimImage::fromItkImage( const ImageType *image ) {
  
  ImageType::RegionType region;
  region = image->GetBufferedRegion();
  ImageType::SizeType size = region.GetSize();
  ImageType::IndexType start = region.GetIndex();
  unsigned int nx = size[0];
  unsigned int ny = size[1];

  // allocate this image
  if (this->alloc( nx, ny, 1, sizeof(PixelType))<=0) return;

  // copy data
  typedef itk::ImageRegionConstIterator< ImageType >  IteratorType;
  IteratorType it( image, region );
  it.GoToBegin();
  const PixelType *data = (PixelType *) this->sampleBits(0);
  while( !it.IsAtEnd() ) {
    *data = it.Get();
    ++it;
    ++data;
  }
}