C++ RegionType::SetIndex方法代码示例

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;
    }
}

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;
}

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;
}

  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);
	
}

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();
}

void
NrrdPlugin::readSlice(int idx[3], int sz[3],
		       int nbytes, uchar *slice)
{
  typedef itk::Image<T, 3> ImageType;

  typedef itk::ImageFileReader<ImageType> ReaderType;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(m_fileName[0].toAscii().data());
  typedef itk::NrrdImageIO NrrdIOType;
  NrrdIOType::Pointer nrrdIO = NrrdIOType::New();
  reader->SetImageIO(nrrdIO);

  typedef itk::RegionOfInterestImageFilter<ImageType,ImageType> RegionExtractor;

  ImageType::RegionType region;
  ImageType::SizeType size;
  ImageType::IndexType index;
  index[2] = idx[2];
  index[1] = idx[1];
  index[0] = idx[0];
  size[2] = sz[2];
  size[1] = sz[1];
  size[0] = sz[0];
  region.SetIndex(index);
  region.SetSize(size);
  
  // Extract the relevant sub-region.
  RegionExtractor::Pointer extractor = RegionExtractor::New();
  extractor->SetInput(reader->GetOutput());
  extractor->SetRegionOfInterest(region);
  extractor->Update();
  ImageType *dimg = extractor->GetOutput();
  char *tdata = (char*)(dimg->GetBufferPointer());
  memcpy(slice, tdata, nbytes);
}

ITKImplicit::ITKImplicit(void)
{
	new SurfParamString(this,&dirFile,"","fileDir", "Input file or dicom directory");
	new SurfParamButton(this,new LoadButtonCallback(this),"load","","load a dicom or imagetype file (extension .vtk suggested)");
	new SurfParamButton(this,new SaveButtonCallback(this),"save","","save a dicom or imagetype file (extension .vtk suggested)");
	new SurfParamDouble(this,&threshold,0,"Threshold", "CT Number");
	new SurfParamDouble(this,&scale,1.0,"Scale", "Volume Scaling Factor");
	new SurfParamgmVector3(this,&tr,gmVector3(0,0,0),"Translation","Volume Translation");
	new SurfParamButton(this, new ApplyParameterChangesCallback(this),"Apply","Apply parameter changes",  "apply parameter changes: Namely translation, scaling,"
																											"and spline order");

	new SurfSurfRefParam(this, &surface,"<invalid/optional>","Surf(opt.)","SurfaceReference(optional)",	
																			"The ITK Volume can also be initialized using a" 
																			"surface. In the case of an implicit the bounding box should be meaningful."
																			"If no bounding box is provided a unit cube centered at the origin is assumed."
																			"The implicit is sampled inside of this bounding box."
																			"In the case of a surface mesh, the behavior is special:"
																			"a distance field is calculated inside of ITS bounding box");


	new SurfParamInt(this,&(voxelRes[0]),64,"sizeX", "NbX", "Nb of gridcells along x-Axis");
	new SurfParamInt(this,&(voxelRes[1]),64,"sizeY","nbY", "Nb of gridcells along y-Axis");
	new SurfParamInt(this,&(voxelRes[2]),64,"sizeZ","nbZ", "Nb of gridcells along z-Axis");
	new SurfParamDouble(this, &boundingBoxScaling, 1.3,"scale","scalebb", "Scaling of the bounding box to make it slightly bigger");

	new SurfAttrRefParam(this,(ParticleAttribute **)&bbox,"invalid:ParticleBoundingBox","bbox","box",
		"Particle bounding box attribute (it is advised to provide one for implicit surfaces).");

	new SurfParamButton(this, new SetITKImplicit(this),"initImp","Initialize ITK implicit","set ITK imp");
	new SurfParamInt(this,&order,3,"Order","Spline Order");


	//DICOM reader
	dicomIO = ImageIOType::New();
	reader = ReaderType::New();
	reader->SetImageIO( dicomIO );
	nameGenerator = NamesGeneratorType::New();


	//now we initialize a simple default image
	//and all the rest of this constructor is just to do this...
	//...I love the ITK interface ... lol
	//init default image, this can still be overwritten
	myImage = ImageType::New();

	ImageType::SpacingType spacing;
	spacing[0]=0.1;
	spacing[1]=0.1;;
	spacing[2]=0.1;;
	myImage->SetSpacing(spacing);
	
	//we set the origin such that center lies at (0,0,0)
	ImageType::PointType origin;
	origin[0]=-0.05;
	origin[1]=-0.05;
	origin[2]=-0.05;
	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]=1;
	size[1]=1;
	size[2]=1;
	ImageType::RegionType region;
	region.SetIndex(start);
	region.SetSize(size);
	myImage->SetRegions(region);
	
	//allocate memory
	myImage->Allocate();


	//interpolating function
	myInterpFunc = InterpFunc::New();	
	myInterpFunc->SetSplineOrder((unsigned int) 0);
	myInterpFunc->SetInputImage(myImage);

	
}

/**
 * draws the tree from the TreeDrawer into a volumetric 3D
 * image as a series of 2D png slices
 */
void drawImage(TreeDrawer * td, const char* rootName){
	typedef unsigned char PixelType;
	const unsigned int Dimension = 3;
	typedef itk::Image< PixelType, Dimension > ImageType;

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

	ImageType::SizeType size;
	size[0] = td->dim[0]; // size along X
	size[1] = td->dim[1]; // size along Y
	size[2] = td->dim[2]; // size along Z

	ImageType::IndexType start;
	start[0] = 0; // first index on X
	start[1] = 0; // first index on Y
	start[2] = 0; // first index on Z

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

	image->SetRegions( region );
	image->Allocate();
	
	ImageType::IndexType pixelIndex;
	pixelIndex[0] = 0; // x position
	pixelIndex[1] = 0; // y position
	pixelIndex[2] = 0; // z position

	for(int i = 0; i < td->dim[0]; i++){
		for(int j = 0; j < td->dim[1]; j++){
			for(int k = 0 ; k < td->dim[2]; k++){
				pixelIndex[0] = i;
				pixelIndex[1] = j;
				pixelIndex[2] = k;
				
				image->SetPixel(pixelIndex, td->imageAt(i, j, k));
			}
		}
	}


	typedef itk::Image< unsigned char, 2 > Image2DType;
	typedef itk::ImageSeriesWriter< ImageType, Image2DType > WriterType;
	WriterType::Pointer writer = WriterType::New();
	writer->SetInput( image);

	typedef itk::NumericSeriesFileNames NameGeneratorType;
	NameGeneratorType::Pointer nameGenerator = NameGeneratorType::New();

	std::string format = rootName;
	format += "%03d";
	format += ".jpg";
	nameGenerator->SetSeriesFormat( format.c_str() );

	const unsigned int firstSlice = start[2];
	const unsigned int lastSlice = start[2] + size[2] - 1;
	nameGenerator->SetStartIndex( firstSlice );
	nameGenerator->SetEndIndex( lastSlice );
	nameGenerator->SetIncrementIndex( 1 );

	writer->SetFileNames( nameGenerator->GetFileNames() );

	try{
		writer->Update();
	}catch( itk::ExceptionObject & excp ){
        
        throw "Exception thrown while reading the image";
        
	}

	return;
}

void Cell::GetMask(const std::string & soma_filename)
{
	typedef itk::ImageFileReader< SomaImageType > ReaderType;
	ReaderType::Pointer reader = ReaderType::New();
	reader->SetFileName(soma_filename);

	//**WARNING** PURPOSELY COMMENTED OUT SO YOU DO NOT TRY TO WRITE THIS CODE, THIS KILLS PERFORMANCE BECAUSE IT ATTEMPTS TO READ THE ENTIRE IMAGE FOR EACH CELL INSTEAD OF JUST THE ROI
	//try
	//{
	//	reader->Update();
	//}
	//catch (itk::ExceptionObject &err)
	//{
	//	std::cerr << "reader Exception: " << err << std::endl;
	//}

	typedef itk::RegionOfInterestImageFilter< SomaImageType, MaskImageType > ROIFilterType;
	ROIFilterType::Pointer roi_filter = ROIFilterType::New();

	ImageType::IndexType start;
	start[0] = this->roi_origin[0];
	start[1] = this->roi_origin[1];
	start[2] = this->roi_origin[2];

	ImageType::SizeType size = this->roi_size;

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

	roi_filter->SetRegionOfInterest(desiredRegion);
	roi_filter->SetInput(reader->GetOutput());

	try
	{
		//roi_filter->Update();
		ftk::TimeStampOverflowSafeUpdate(roi_filter.GetPointer());
	}
	catch (itk::ExceptionObject &err)
	{
		std::cout << "roi_filter Exception: " << err << std::endl;
	}

	this->mask = roi_filter->GetOutput();
	this->mask->DisconnectPipeline();	//Disconnect pipeline so we don't propagate...

	ImageType::PointType origin;
	origin[0] = 0;
	origin[1] = 0;
	origin[2] = 0;
	this->mask->SetOrigin(origin);

	//Make the file name of the mask image
	std::stringstream mask_filename_stream;
	mask_filename_stream << cell_x << "_" << cell_y << "_" << cell_z << "_mask.TIF";	//X_Y_Z_masked.TIF

	//Write the masked cell image
	//WriteImage(mask_filename_stream.str(), this->mask);

	//Get the label image from the binary image
	typedef itk::BinaryImageToLabelMapFilter< MaskImageType > BinaryToLabelFilterType;
	BinaryToLabelFilterType::Pointer labelMapFilter = BinaryToLabelFilterType::New();
	labelMapFilter->SetInput(this->mask);

	try
	{
		ftk::TimeStampOverflowSafeUpdate( labelMapFilter.GetPointer() );
		//labelMapFilter->Update();
	}
	catch (itk::ExceptionObject &err)
	{
		std::cerr << "labelMapFilter exception: " << err << std::endl;
		std::cerr << "Mask image: " << this->mask << std::endl;
		std::cerr << this->mask << std::endl;
		std::cerr << labelMapFilter << std::endl;
	}

	BinaryToLabelFilterType::OutputImageType::Pointer label_map_image = labelMapFilter->GetOutput();
	label_map_image->DisconnectPipeline();

	typedef itk::LabelMapToLabelImageFilter< BinaryToLabelFilterType::OutputImageType, LabelImageType > LabelMapToLabelImageFilterType;
	LabelMapToLabelImageFilterType::Pointer labelImageFilter = LabelMapToLabelImageFilterType::New();
	labelImageFilter->SetInput(label_map_image);
	try
	{
		ftk::TimeStampOverflowSafeUpdate( labelImageFilter.GetPointer() );
		//labelImageFilter->Update();
	}
	catch (itk::ExceptionObject &err)
	{
		std::cerr << "labelImageFilter exception: " << err << std::endl;
	}

	this->soma_label_image = labelImageFilter->GetOutput();
	this->soma_label_image->DisconnectPipeline();	//Disconnect pipeline so we don't propagate...
}

bool CCropDialog::GetCropX(QString strSrcImage, quint32 &nStartX)
{
	if(QFile::exists(strSrcImage) != true)
	{
		return false; // error
	}

	QString strTemplateFilename(tr("%1/%2")
								.arg(ui.lineEditSourceFolder->text())
								.arg(sTEMPLATE_IMAGE_FILENAME));
	///////////////////////////////
	// read the template image file
	
	ReaderType::Pointer pTemplateReader = ReaderType::New();
	pTemplateReader->SetFileName(strTemplateFilename.toStdString().c_str());
	try
	{
		pTemplateReader->Update();
	}
	catch( itk::ExceptionObject & err ) 
	{ 
		std::cerr << "ExceptionObject caught !" << std::endl; 
		std::cerr << err << std::endl; 
		return false; // EXIT_FAILURE;
	} 

	quint32 nX, nY;
	quint32 nMinSAD, nSAD, nMinSADX;
	quint32 nRightEdgeX;
	
	nMinSADX = nDEF_START_X; // just in case
	nMinSAD = 0xFFFFFFFF;
	
	ImageType::SizeType size;
	size[0] = nTEMPLATE_WIDTH;
	size[1] = nTEMPLATE_LENGHT;
	ImageType::IndexType start;

	ReaderType::Pointer pReader = ReaderType::New();
	pReader->SetFileName( strSrcImage.toStdString().c_str() );
	//////////////////////
	// loop start
	
	for(nY = 0, nX = nSTART_X; nX < nEND_X; nX++)
	{
		///////////////////////////////////////
		// read a region from the search image
		start[0] = nX, start[1] = nY;
	
		ImageType::RegionType DesiredRegion;
		DesiredRegion.SetSize(  size  );
		DesiredRegion.SetIndex( start );
		FilterType::Pointer pFilter = FilterType::New();
		pFilter->SetRegionOfInterest( DesiredRegion );
		pFilter->SetInput( pReader->GetOutput() );
		
		//bool bbErrorInSourceImage = false;
		try 
		{ 
			pFilter->Update(); 
		} 
		catch( itk::ExceptionObject & err ) 
		{ 
			std::cerr << "ExceptionObject caught !" << std::endl; 
			std::cerr << err << std::endl; 
			return false; // EXIT_FAILURE;
		} 
		
		ConstIteratorType SearchIt(pFilter->GetOutput(), pFilter->GetOutput()->GetRequestedRegion());
		ConstIteratorType TemplateIt(pTemplateReader->GetOutput(), pTemplateReader->GetOutput()->GetRequestedRegion());
	
		nSAD = SumOfAbsoluteDifferences(SearchIt, TemplateIt);
		if(nSAD < nMinSAD)
		{
			nMinSAD = nSAD;
			nMinSADX = nX;
		}
	}
	////////////////
	// loop end
	nRightEdgeX = nMinSADX +  nTEMPLATE_WIDTH/2;
	nStartX = nRightEdgeX - m_TissueBlock.m_nWidth;
	return true;
}

void Initialisation::savePointAsBinaryImage(ImageType::Pointer initialImage, string filename, OrientationType orientation)
{
    if (points_.size() > 0)
    {
        typedef itk::Image< unsigned char, 3 > BinaryImageType;
        BinaryImageType::Pointer binary = BinaryImageType::New();
        ImageType::RegionType region;
        ImageType::IndexType start;
        start[0] = 0; start[1] = 0; start[2] = 0;
        ImageType::SizeType size, imSize = initialImage->GetLargestPossibleRegion().GetSize();
        size[0] = imSize[0]; size[1] = imSize[1]; size[2] = imSize[2];
        region.SetSize(size);
        region.SetIndex(start);
        binary->CopyInformation(initialImage);
        binary->SetRegions(region);
        binary->Allocate();
        binary->FillBuffer(false);
        
        typedef ImageType::IndexType IndexType;
        ContinuousIndex ind;
        IndexType ind2;
        unsigned int pSize = points_.size();
        unsigned int indexMiddle = 0;
        for (unsigned int i=0; i<pSize; i++) {
            if (points_[i][1] == startSlice_)
                indexMiddle = i;
        }
        ind[0] = points_[indexMiddle][0]; ind[1] = points_[indexMiddle][1]; ind[2] = points_[indexMiddle][2];
        PointType pt;
        inputImage_->TransformContinuousIndexToPhysicalPoint(ind, pt);
        initialImage->TransformPhysicalPointToIndex(pt, ind2);
        binary->SetPixel(ind2,true);
        
        OrientImage<BinaryImageType> orientationFilter;
        orientationFilter.setInputImage(binary);
        orientationFilter.orientation(orientation);
        binary = orientationFilter.getOutputImage();
        
        ImageIterator it( binary, binary->GetRequestedRegion() );
        it.GoToBegin();
        while(!it.IsAtEnd())
        {
            if (it.Get()==true)
            {
                ind2 = it.GetIndex();
                break;
            }
            ++it;
        }
        if (verbose_) cout << "Center of spinal cord saved on pixel : " << ind2 << endl;
        
        WriterBinaryType::Pointer writer = WriterBinaryType::New();
        itk::NiftiImageIO::Pointer io = itk::NiftiImageIO::New();
        writer->SetImageIO(io);
        writer->SetFileName(filename);
        writer->SetInput(binary);
        try {
            writer->Write();
        } catch( itk::ExceptionObject & e ) {
            std::cerr << "Exception caught while writing image " << std::endl;
            std::cerr << e << std::endl;
        }
    }
    else cout << "Error: Spinal cord center not detected" << endl;
}

int main(int, char ** argv)
{
	std::string inputFolder = argv[1];
	std::string outputFolder = argv[2];
	std::string lookupFilename = argv[3];


	// load the lookup file and read the image filenames
	LookupMap lookup;
	FilenamesType imageFilenames;
	readFilenames(inputFolder, imageFilenames);
	loadLookupFile(lookupFilename, lookup);

	// looop through the list of images
	for(unsigned int i = 0; i < imageFilenames.size(); i++)
	{
		std::string fname = imageFilenames[i];
		QFileInfo finfo(QString::fromStdString(fname));
		std::string basename = finfo.completeBaseName().toStdString() + ".nrrd";
		FilenamesType dicomFilenames = lookup[basename];

		// slice up the input image to recreate the output
		typedef utils::ImageVolume ImageType;
		ImageType::Pointer input = utils::ImageVolumeIO::Read(fname);

		const unsigned int slices = input->GetLargestPossibleRegion().GetSize()[2];

		for(unsigned int slice = 0; slice < slices; slice++)
		{
			typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> ROIFilter;
			ROIFilter::Pointer roiFilter = ROIFilter::New();
			roiFilter->SetInput(input);

			ImageType::RegionType roi = input->GetLargestPossibleRegion();
			ImageType::IndexType roiIndex = roi.GetIndex();
			ImageType::SizeType roiSize = roi.GetSize();

			roiIndex[2] = slice;
			roiSize[2] = 1;
			roi.SetSize(roiSize);
			roi.SetIndex(roiIndex);

			roiFilter->SetRegionOfInterest(roi);
			roiFilter->Update();

			ImageType::Pointer imSlice = roiFilter->GetOutput();


			// load the dicom file
			std::string dicomFilename = dicomFilenames[slice];
			DcmFileFormat fileFormat;
			fileFormat.loadFile(dicomFilename.c_str());

			unsigned int numberOfPixels = imSlice->GetLargestPossibleRegion().GetNumberOfPixels();
			unsigned short * buffer = new unsigned short[numberOfPixels];
			fileFormat.getDataset()->putAndInsertUint16Array(DCM_PixelData, buffer, imSlice->GetLargestPossibleRegion().GetNumberOfPixels());

			// reset the pixel values
			itk::ImageRegionConstIterator<ImageType> it(imSlice, imSlice->GetLargestPossibleRegion());
			unsigned int count = 0;
			while(!it.IsAtEnd())
			{
				buffer[count] = it.Get();
				++it; ++count;
			}
			
			fileFormat.getDataset()->putAndInsertUint16Array(DCM_PixelData, buffer, imSlice->GetLargestPossibleRegion().GetNumberOfPixels());
			// create the output filename 
			std::stringstream ss;
			ss << outputFolder << "/dicom_" << i << "_" << slice << ".dcm";


			fileFormat.saveFile(ss.str().c_str());
		}
	}

	return 0;
}

void CUnitVolumeMakerDlg::GetImage(int iWidth, int iLength, int iHeight, int iZ,
								 EZoomOut eZoomOutLevel)
{
	int nCol;
	int nX, nY, nZ;
	int nRibbonWidth;
	int nImageWidth;
	int nRightWidth;
	int nTotalWidth;
	QString strPathName, strImageFileName, strImageVolPathName, strImage4VolFileName;
	
	nRibbonWidth = qRound(float(m_TissueBlock.m_Ribbon.m_nWidth)/eZoomOutLevel);
	nCol = (int)((m_TissueBlock.m_UnitCube.m_nWidth*iWidth)/nRibbonWidth);
	nX = nRibbonWidth - (nRibbonWidth*(nCol+1)-(iWidth*m_TissueBlock.m_UnitCube.m_nWidth));
	nY = iLength*m_TissueBlock.m_UnitCube.m_nLength;
	nZ = iHeight*m_TissueBlock.m_UnitCube.m_nHeight;
	strPathName = QString(tr("%1%2%3/%4/")
			.arg(sINIT_DATA_DIR).arg(sZOOM_OUT_DIR_PREFIX).arg(eZoomOutLevel).arg(nCol));
	strImageFileName.sprintf("%05d%s", nZ+eZoomOutLevel*iZ, sDATA_FILE_FORMAT);
		
	nRightWidth = nRibbonWidth - nX;
	
	if(nRightWidth > m_TissueBlock.m_UnitCube.m_nWidth)
	{
		nImageWidth = qMin(m_TissueBlock.m_UnitCube.m_nWidth, nRibbonWidth);
	}
	else
	{
		nImageWidth = nRightWidth;
	}
	
	// extract nX, nY, nImageWidth, m_TissueBlock.m_UnitCube.m_nLength
	//...
	qDebug() << "nX: " << nX << ", nY: " << nY << ", nImageWidth: " << nImageWidth << "\n";
	ui.listWidgetResult->addItem(tr("nCol: %1, nX: %2, nY: %3, nImageWidth: %4").arg(nCol).arg(nX).arg(nY).arg(nImageWidth));
	ui.listWidgetResult->addItem(strPathName);
	ui.listWidgetResult->addItem(strImageFileName);
	
	// source image filename (full path)
	QString strFileName = strPathName + strImageFileName;
	
	/////////////////////////
	// read an image
	// typedef
	
	typedef unsigned char PixelType;
	const unsigned int nDimension = 2;
	typedef itk::Image<PixelType, nDimension> ImageType;
	typedef itk::ImageFileReader<ImageType> ReaderType;
	typedef itk::ImageFileWriter<ImageType> WriterType;
	//typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> ROIFilterType;
	typedef itk::ImageRegionConstIterator< ImageType > ConstIteratorType;
	typedef itk::ImageRegionIterator< ImageType>       IteratorType;

	/////////////////////////////////////
	// input: a region of a source image 
	ImageType::RegionType InputRegion;
	ImageType::RegionType::IndexType InputStart;
	ImageType::RegionType::SizeType  InputSize;
	
	InputStart[0] = nX;
	InputStart[1] = nY;
	
	InputSize[0]  = nImageWidth;
	InputSize[1]  = m_TissueBlock.m_UnitCube.m_nLength;
	
	InputRegion.SetSize( InputSize );
	InputRegion.SetIndex( InputStart );
	
	/////////////////////////////////////
	// output: a region of a output image 
	ImageType::RegionType OutputRegion, OutputSubRegion;
	ImageType::RegionType::IndexType OutputStart;
	ImageType::RegionType::SizeType  OutputSize;

	OutputStart[0] = 0;
	OutputStart[1] = 0;
	OutputSize[0]  = m_TissueBlock.m_UnitCube.m_nWidth; //InputSize[0];
	OutputSize[1]  = m_TissueBlock.m_UnitCube.m_nLength; //InputSize[1];

	OutputRegion.SetSize( OutputSize );
	OutputRegion.SetIndex( OutputStart );

	OutputSize[0]  = InputSize[0];
	OutputSize[1]  = InputSize[1];

	OutputSubRegion.SetSize( OutputSize );
	OutputSubRegion.SetIndex( OutputStart );

	bool bErrorInSourceImage = false;

	/////////////////////////////////////
	// read the source image
	ReaderType::Pointer pReader = ReaderType::New();
	pReader->SetFileName(strFileName.toStdString());
	try
	{
	    pReader->Update();
	}
	catch ( itk::ExceptionObject &err)
//.........这里部分代码省略.........