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

void Binarize(InputImageType::Pointer im)
{
	IteratorType iter(im,im->GetLargestPossibleRegion());
	for(iter.GoToBegin(); !iter.IsAtEnd(); ++iter)
	{
		iter.Set((iter.Get()!=0)?255:0);
	}
}

Input2DImageType::Pointer Curvelet::getSlice(InputImageType::Pointer im, int slice)
{
	Input2DImageType::Pointer out = Input2DImageType::New();
	Input2DImageType::SizeType size;
	size[0] = im->GetLargestPossibleRegion().GetSize()[0];
	size[1] = im->GetLargestPossibleRegion().GetSize()[1];
	Input2DImageType::IndexType index;
	index.Fill(0);
	Input2DImageType::RegionType region;
	region.SetSize(size);
	region.SetIndex(index);
	out->SetRegions(region);
	out->Allocate();
	if(out->GetBufferPointer()==NULL)
		printf("Could not allocate memory -1 ... I'm going to crash any moment now.. \n");
	memcpy(out->GetBufferPointer(),im->GetBufferPointer()+slice*size[0]*size[1],size[0]*size[1]*sizeof(unsigned char));
	return out;
}

LabelImageType::Pointer getYousefSegmented(InputImageType::Pointer im_input,std::list<Seed> &seed_list,char *filename)
{
	// copy the image into a unsigned char *
	char configfile[1024];
	strcpy(configfile,filename);

	printf("Entering YousefSeg\n");
	InputImageType::SizeType size = im_input->GetLargestPossibleRegion().GetSize();
	unsigned char * in_Image;
	in_Image = (unsigned char*) malloc( size[0]*size[1]*(size[2]+1)*sizeof(unsigned char));
	if(in_Image == NULL)
	{
		printf("Couldn't allocate memory\n");
	}

	memset(in_Image,0,size[0]*size[1]*(size[2]+1)*sizeof(unsigned char));

	ConstIteratorType pix_buf(im_input,im_input->GetLargestPossibleRegion());
	int ind = 0;
	for ( pix_buf.GoToBegin(); !pix_buf.IsAtEnd(); ++pix_buf, ++ind )
		in_Image[ind]=(pix_buf.Get());

	printf("Copied input data\n");

	yousef_nucleus_seg *NucleusSeg = new yousef_nucleus_seg();
	NucleusSeg->readParametersFromFile(configfile);
	NucleusSeg->setDataImage(in_Image,size[0],size[1],size[2]+1,"null");

	unsigned short * output_img;
	//	int *bounds_img;
	NucleusSeg->runBinarization();
	output_img = NucleusSeg->getBinImage();
	//	getITKImage(output_img);
	//	getProcessedBinaryImage(
	NucleusSeg->runSeedDetection();
	std::vector<Seed> seeds = NucleusSeg->getSeeds();
	printf("In yousef_seg Seed size = %d\n", (int)seeds.size());
	std::vector<Seed>::iterator iter = seeds.begin();
	for(;iter!=seeds.end();iter++)
	{
		seed_list.push_back(*iter);
	}
	NucleusSeg->runClustering();
	printf("Finished Clustering\n");
	if(NucleusSeg->isSegmentationFinEnabled())
	{
		NucleusSeg->runAlphaExpansion3D();		
		output_img=NucleusSeg->getSegImage();
	}
	else
	{
		output_img=NucleusSeg->getClustImage();
	}

	//	bounds_img = NucleusSeg->getBoundsImage();

	printf("Finished segmentation\n");

	LabelImageType::Pointer label = LabelImageType::New();
	label->SetRegions(im_input->GetLargestPossibleRegion());
	label->Allocate();

	LabelIteratorType liter(label,label->GetLargestPossibleRegion());
	ind = 0;
	for(liter.GoToBegin();!liter.IsAtEnd();++liter,++ind)
	{
		liter.Set(output_img[ind]);
	}
	delete NucleusSeg;
	free(in_Image);
	return label;
}

int main()
{
	InputImageType::Pointer im = readImage<InputImageType>("C:/Users/arun/Research/Farsight/exe/bin/CF_1_inverted_bg_sub.tif");
	FILE *fp = fopen("C:/Users/arun/Research/Farsight/exe/bin/seeds.txt","r");

	//IteratorType initer(im,im->GetLargestPossibleRegion());
	//initer.GoToBegin();
	//
	//for(;!initer.IsAtEnd(); ++initer)
	//{
	//	initer.Set(transfer_function1(initer.Get()));
	//}
	//
	//writeImage<InputImageType>(im,"C:/Users/arun/Research/Farsight/exe/bin/hp2_cropped2_filtered.tif");
	//
	//return 0;
	typedef itk::SymmetricSecondRankTensor<double,3> HessianType;
	typedef itk::Hessian3DToVesselnessMeasureImageFilter<float> MeasureType;
	typedef itk::Image<HessianType,3> HessianImageType;
	typedef itk::MultiScaleHessianBasedMeasureImageFilter< InputImageType, HessianImageType, FloatImageType> VesselnessFilterType;
	
	std::vector<InputImageType::RegionType> in1,in2,out1,out2;

	get_tiles(im->GetLargestPossibleRegion().GetSize(),1500,1500,1500,100,100,10,in1,in2,out1,out2);
	
	InputImageType::Pointer om;
/*
	om = InputImageType::New();
	om->SetRegions(im->GetLargestPossibleRegion());
	om->Allocate();
	for(int counter = 0; counter < in1.size(); counter++)
	{
		InputImageType::Pointer imtile = InputImageType::New();//
		imtile->SetRegions(in2[counter]);
		imtile->Allocate();

		in1[counter].Print(std::cout);
		in2[counter].Print(std::cout);
		IteratorType iter1(im,in1[counter]);
		IteratorType iter2(imtile,in2[counter]);
		for(iter1.GoToBegin(),iter2.GoToBegin();!iter1.IsAtEnd(); ++iter1,++iter2)
		{
			iter2.Set(iter1.Get());
		}

		VesselnessFilterType::Pointer vfilt = VesselnessFilterType::New();
		MeasureType::Superclass::Pointer measure = MeasureType::New();
		
		vfilt->SetInput(imtile);
		vfilt->SetHessianToMeasureFilter((VesselnessFilterType::HessianToMeasureFilterType *)measure);
		vfilt->SetSigmaMinimum(3.0);
		vfilt->SetSigmaMaximum(5.0);
		vfilt->SetNumberOfSigmaSteps(3);
		vfilt->SetSigmaStepMethod(VesselnessFilterType::EquispacedSigmaSteps);
		vfilt->Update();
		FloatImageType::Pointer omtile = vfilt->GetOutput();
		
		typedef itk::ImageRegionIterator<FloatImageType> FloatIteratorType;
		FloatIteratorType iter3;
		iter1 = IteratorType(om,out1[counter]);
		iter3 = FloatIteratorType(omtile,out2[counter]);
		for(iter1.GoToBegin(),iter3.GoToBegin();!iter1.IsAtEnd();++iter1,++iter3)
		{
			iter1.Set(iter3.Get());
		}
	}
	writeImage<InputImageType>(om,"C:/Users/arun/Research/Farsight/exe/bin/vesselnesstest.tif");
*/
	om = readImage<InputImageType>("C:/Users/arun/Research/Farsight/exe/bin/vesselnesstest.tif");
	typedef itk::BinaryBallStructuringElement<InputImageType::PixelType,3> StructElementType;
	typedef itk::GrayscaleDilateImageFilter<InputImageType,InputImageType,StructElementType> FilterType1;
	FilterType1::Pointer minfilt = FilterType1::New();
	minfilt->SetInput(om);
	FilterType1::RadiusType radius;
	radius[0] = 1;
	radius[1] = 1;
	radius[2] = 1;
	StructElementType strel;
	strel.SetRadius(radius);
	minfilt->SetKernel(strel);
	minfilt->Update();

	InputImageType::Pointer seed_out  = InputImageType::New();
	seed_out->SetRegions(om->GetLargestPossibleRegion());
	seed_out->Allocate();
	seed_out->FillBuffer(0);

	int thresh_value = 6;
	int number_of_seeds = 200;
	int tnum_seeds = 0;

	typedef itk::ImageRegionIteratorWithIndex<InputImageType> IndexIteratorType;

	IndexIteratorType it1(minfilt->GetOutput(),minfilt->GetOutput()->GetLargestPossibleRegion());
	IteratorType it2(om,om->GetLargestPossibleRegion());

	for(it2.GoToBegin();!it2.IsAtEnd(); ++it2)
	{
		if(it2.Get()>thresh_value)
			tnum_seeds++;
//.........这里部分代码省略.........

static mitk::Image::Pointer ResampleBySpacing(mitk::Image *input, float *spacing, bool useLinInt = true, bool useNN = false)
{
  if (!useNN)
  {
    InputImageType::Pointer itkImage = InputImageType::New();
    CastToItkImage(input,itkImage);

    /**
   * 1) Resampling
   *
   */
    // Identity transform.
    // We don't want any transform on our image except rescaling which is not
    // specified by a transform but by the input/output spacing as we will see
    // later.
    // So no transform will be specified.
    typedef itk::IdentityTransform<double, 3> T_Transform;

    // The resampler type itself.
    typedef itk::ResampleImageFilter<InputImageType, InputImageType>  T_ResampleFilter;

    // Prepare the resampler.
    // Instantiate the transform and specify it should be the id transform.
    T_Transform::Pointer _pTransform = T_Transform::New();
    _pTransform->SetIdentity();

    // Instantiate the resampler. Wire in the transform and the interpolator.
    T_ResampleFilter::Pointer _pResizeFilter = T_ResampleFilter::New();

    // Specify the input.
    _pResizeFilter->SetInput(itkImage);

    _pResizeFilter->SetTransform(_pTransform);

    // Set the output origin.
    _pResizeFilter->SetOutputOrigin(itkImage->GetOrigin());

    // Compute the size of the output.
    // The size (# of pixels) in the output is recomputed using
    // the ratio of the input and output sizes.
    InputImageType::SpacingType inputSpacing = itkImage->GetSpacing();
    InputImageType::SpacingType outputSpacing;
    const InputImageType::RegionType& inputSize = itkImage->GetLargestPossibleRegion();

    InputImageType::SizeType outputSize;
    typedef InputImageType::SizeType::SizeValueType SizeValueType;

    // Set the output spacing.
    outputSpacing[0] = spacing[0];
    outputSpacing[1] = spacing[1];
    outputSpacing[2] = spacing[2];

    outputSize[0] = static_cast<SizeValueType>(inputSize.GetSize()[0] * inputSpacing[0] / outputSpacing[0] + .5);
    outputSize[1] = static_cast<SizeValueType>(inputSize.GetSize()[1] * inputSpacing[1] / outputSpacing[1] + .5);
    outputSize[2] = static_cast<SizeValueType>(inputSize.GetSize()[2] * inputSpacing[2] / outputSpacing[2] + .5);

    _pResizeFilter->SetOutputSpacing(outputSpacing);
    _pResizeFilter->SetSize(outputSize);

    typedef itk::LinearInterpolateImageFunction< InputImageType > LinearInterpolatorType;
    LinearInterpolatorType::Pointer lin_interpolator = LinearInterpolatorType::New();


    typedef itk::WindowedSincInterpolateImageFunction< InputImageType, 4> WindowedSincInterpolatorType;
    WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New();

    if (useLinInt)
      _pResizeFilter->SetInterpolator(lin_interpolator);
    else
      _pResizeFilter->SetInterpolator(sinc_interpolator);

    _pResizeFilter->Update();

    mitk::Image::Pointer image = mitk::Image::New();
    image->InitializeByItk(_pResizeFilter->GetOutput());
    mitk::GrabItkImageMemory( _pResizeFilter->GetOutput(), image);
    return image;
  }

  BinaryImageType::Pointer itkImage = BinaryImageType::New();
  CastToItkImage(input,itkImage);

  /**
 * 1) Resampling
 *
 */
  // Identity transform.
  // We don't want any transform on our image except rescaling which is not
  // specified by a transform but by the input/output spacing as we will see
  // later.
  // So no transform will be specified.
  typedef itk::IdentityTransform<double, 3> T_Transform;

  // The resampler type itself.
  typedef itk::ResampleImageFilter<BinaryImageType, BinaryImageType>  T_ResampleFilter;

  // Prepare the resampler.
  // Instantiate the transform and specify it should be the id transform.
  T_Transform::Pointer _pTransform = T_Transform::New();
  _pTransform->SetIdentity();
//.........这里部分代码省略.........

  static void TestOverlay(mitk::Image::Pointer original,
                          mitk::Image::Pointer truth,
                          const double lower,
                          const double upper)
  {
    mitk::Image::Pointer overlayImage;
    const double th[] = {lower, upper};
    typedef itk::Image<unsigned int, 3> ImageType;
    ImageType::Pointer itkOverlayImage = ImageType::New();

    AccessByItk_2(original, InternalThreshold, overlayImage, th);
    /*
        AccessFixedDimensionByItk_2( original, InternalThreshold2,
                                     3, itkOverlayImage,
                                     th );

        overlayImage = mitk::ImportItkImage( itkOverlayImage );
    */
    // mitk::IOUtil::Save(truth, "/tmp/truth_TestOverlay.nii");
    try
    {
      // mitk::Image::Pointer temp = overlayImage;
      mitk::IOUtil::Save(overlayImage, "/tmp/overlayImage_TestOverlay.nrrd");
    }
    catch (const itk::ExceptionObject &e)
    {
      MITK_ERROR << "Save image: exception : " << e.what();
    }

    typedef itk::Image<unsigned int, 3> InputImageType;

    InputImageType::Pointer overlayItk;
    try
    {
      mitk::CastToItkImage(overlayImage, overlayItk);
    }
    catch (const mitk::Exception &e)
    {
      MITK_ERROR << "(CAST) Catched exception while creating accessor " << e.what();
      // MITK_TEST_FAILED_MSG("Exception for ouverlay image");
    }

    /*
        typedef  itk::ImageFileWriter< InputImageType  > WriterType;
        WriterType::Pointer writer = WriterType::New();
        writer->SetFileName("/tmp/overlayITK_TestOverlay.nii");
        writer->SetInput(overlayItk);
        writer->Update();
    */
    InputImageType::Pointer truthItk;
    mitk::CastToItkImage(truth, truthItk);

    bool difference = false;
    /*
        try
        {
          typedef unsigned int TPixel;

          itk::ImageRegionConstIteratorWithIndex< InputImageType > iter( truthItk, truthItk->GetLargestPossibleRegion()
       );
          iter.GoToBegin();

          mitk::ImagePixelReadAccessor< TPixel, 3 > readAccessor( overlayImage, overlayImage->GetVolumeData(0),
       mitk::ImageAccessorBase::ExceptionIfLocked );

          while( !iter.IsAtEnd() )
          {
            TPixel ref = iter.Get();
            TPixel val = readAccessor.GetPixelByIndex( iter.GetIndex() );

            difference |= ( ref != val );
            //if( difference )
            //{
              std::cout << iter.GetIndex() << ":" << ref << " ? " << val << "\n";
            //}
            ++iter;
          }
        }
        catch( const mitk::Exception &e)
        {
          MITK_ERROR << "Catched exception while creating accessor "<< e.what();
          //MITK_TEST_FAILED_MSG("Exception for ouverlay image");
        }



    */

    /*
        typedef itk::Testing::ComparisonImageFilter
          <InputImageType, InputImageType> ComparisonImageFilterType;
        ComparisonImageFilterType::Pointer comp =
          ComparisonImageFilterType::New();
        comp->SetValidInput(truthItk);
        comp->SetTestInput(overlayItk);
        try
        {
          comp->Update();
        }
        catch( const itk::ExceptionObject& e)
//.........这里部分代码省略.........

void segmentation::on_retinalLayerSegButton_clicked()
{
    fstream fin1, fin7, fin11, fin12; 
    char* s1Txt = ("../data/ped/txt/15715_1.txt");//s1.txt
    char* s7Txt = ("../data/ped/txt/15715_7.txt"); 
    char* s11Txt = ("../data/ped/txt/15715_11.txt"); 
    char* s12Txt = ("../data/ped/txt/15715_12.txt");
    //这样开辟内存容易出现堆栈溢出问题, 选择new & delete
    //int s1[64][512], s7[64][512], s11[64][512], s12[64][512]; 
    int (*s1)[512] = new int [64][512];
    int (*s7)[512] = new int [64][512];
    int (*s11)[512] = new int [64][512];
    int (*s12)[512] = new int [64][512];
    fin1.open(s1Txt, ios::in);
    fin7.open(s7Txt, ios::in);
    fin11.open(s11Txt, ios::in);
    fin12.open(s12Txt, ios::in);
    for(int z=0; z<64; z++)
    {
        for(int x=0; x<512; x++)
        {
            fin1>>s1[z][x];
            fin7>>s7[z][x];
            fin11>>s11[z][x];
            fin12>>s12[z][x];
        }
    }
    fin1.close();
    fin7.close();
    fin11.close();
    fin12.close();

    typedef itk::Image<unsigned short,3>InputImageType; 
    typedef itk::Image<unsigned short,3>OutputImageType; 
    InputImageType::Pointer inputImage = InputImageType::New();
    InputImageType::Pointer outputImage = InputImageType::New();
    typedef itk::ImageFileReader<InputImageType>ReaderType;
    ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName(inputFileName); 
    reader->Update();
    inputImage = reader->GetOutput();
    InputImageType::IndexType voxelIndex;
    InputImageType::SizeType imgSize = inputImage->GetLargestPossibleRegion().GetSize();
    OutputImageType::IndexType index;
    index[0]=0;
    index[1]=0;
    index[2]=0;
    OutputImageType::SizeType size;
    size[0]=imgSize[0];
    size[1]=imgSize[1];
    size[2]=imgSize[2];  
    //create a region for enhance result
    OutputImageType::RegionType region;
    region.SetIndex(index);
    region.SetSize(size);
    outputImage->SetRegions(region);
    outputImage->Allocate();
    //make four surfaces into a MHD file
    for(int z = 0; z < imgSize[2]; z++)
        for(int y = 0; y < imgSize[1]; y++)
            for(int x = 0; x < imgSize[0]; x++)
            {
                voxelIndex[0] = x;
                voxelIndex[1] = y;
                voxelIndex[2] = z;
                if(y == s1[z][x])
                {
                    inputImage->SetPixel(voxelIndex, 65535);
                }
                if(y == s7[z][x])
                {
                    inputImage->SetPixel(voxelIndex, 65535);
                }
                if(y == s11[z][x])
                {
                    inputImage->SetPixel(voxelIndex, 65535);
                }
                if(y == s12[z][x])
                {
                    inputImage->SetPixel(voxelIndex, 65535);
                }
            }
            delete []s1;//释放内存
            delete []s7;
            delete []s11;
            delete []s12;

            //文件前缀名
            filePrefix = inputFileName;//char* to string
            filePrefix = filePrefix.substr(0, filePrefix.length() - 4);
            filePrefix = filePrefix + "_layerSeg.mhd";
            strcpy(outputFileName, filePrefix.c_str());//string to char*
            typedef itk::ImageFileWriter<OutputImageType>WriterType;
            WriterType::Pointer writer = WriterType::New();
            writer->SetFileName(outputFileName);
            writer->SetInput(inputImage);
            writer->Update();
            emit returnOutputFileName(outputFileName);//发出信号
}

//腐蚀操作
void segmentation::erosionOperation()//未用到
{
    const unsigned int Dimension = 3;
    typedef unsigned char   InputPixelType;
    typedef unsigned char   OutputPixelType;
    typedef itk::Image< InputPixelType,  Dimension >   InputImageType;
    typedef itk::Image< OutputPixelType, Dimension >   OutputImageType;
    typedef itk::ImageFileReader< InputImageType  >  ReaderType;
    typedef itk::ImageFileWriter< OutputImageType >  WriterType;
    typedef itk::BinaryThresholdImageFilter< InputImageType, InputImageType >  ThresholdFilterType;
    typedef itk::BinaryBallStructuringElement< 
        InputPixelType,
        Dimension  >             StructuringElementType;
    typedef itk::BinaryErodeImageFilter<
        InputImageType, 
        OutputImageType,
        StructuringElementType >  ErodeFilterType;
    ReaderType::Pointer reader = ReaderType::New();
    WriterType::Pointer writerErosion  = WriterType::New();
    ThresholdFilterType::Pointer thresholder = ThresholdFilterType::New();
    ErodeFilterType::Pointer  binaryErode  = ErodeFilterType::New();
    StructuringElementType  structuringElementErosion;
    reader->SetFileName( initFileName );//读入二值图
    reader->Update();
    InputImageType::Pointer inputImage = InputImageType::New();
    InputImageType::IndexType voxelIndex;
    inputImage = reader->GetOutput();
    InputImageType::SizeType imgSize = inputImage->GetLargestPossibleRegion().GetSize();
    long vol = 0;
    unsigned char temp;
    float r;
    for(int z = 0; z < imgSize[2]; z++)
        for(int y = 0; y < imgSize[1]; y++)
            for(int x = 0; x < imgSize[0]; x++)
            {
                voxelIndex[0] = x;
                voxelIndex[1] = y;
                voxelIndex[2] = z;
                temp = inputImage->GetPixel(voxelIndex);
                if(temp == 255)// 255 for PED
                    vol += 1;
            }
            r = pow((3 * vol) / (4 * PI), (1.0 / 3)) ;
            r = r / 20;//experiment data
            structuringElementErosion.SetRadius( r );  // 3x3 structuring element
            structuringElementErosion.CreateStructuringElement();
            binaryErode->SetKernel(  structuringElementErosion );
            //文件前缀名
            filePrefix = inputFileName;//char* to string
            filePrefix = filePrefix.substr(0, filePrefix.length() - 4);
            string erosionFileName;
            erosionFileName = filePrefix + "_erosionResult.mhd";
            strcpy(outputFileName, erosionFileName.c_str());//string to char*
            writerErosion->SetFileName(outputFileName);
            const InputPixelType lowerThreshold = 255;
            const InputPixelType upperThreshold = 255;
            thresholder->SetInput( reader->GetOutput() );
            InputPixelType background =   0;
            InputPixelType foreground = 255;
            thresholder->SetOutsideValue( background );
            thresholder->SetInsideValue(  foreground );
            thresholder->SetLowerThreshold( lowerThreshold );
            thresholder->SetUpperThreshold( upperThreshold );
            binaryErode->SetInput( thresholder->GetOutput() );
            binaryErode->SetErodeValue( foreground );
            writerErosion->SetInput( binaryErode->GetOutput() );
            writerErosion->Update();
            //binaryErode->GetOutput()->GetPixel(index);//获取像素值失败

            //腐蚀结果叠加到原图
            typedef itk::Image< unsigned short,  Dimension >   OriginalImageType;
            typedef itk::ImageFileReader<OriginalImageType>OriginalReaderType;
            OriginalReaderType::Pointer orignalImgreader = OriginalReaderType::New();
            OriginalImageType::Pointer originalImage = OriginalImageType::New();
            OriginalReaderType::IndexType originalImgVoxelIndex;
            reader->SetFileName(outputFileName);//读入腐蚀结果图像
            reader->Update();
            inputImage = reader->GetOutput();
            orignalImgreader->SetFileName(inputFileName);//读入原图像
            orignalImgreader->Update();
            originalImage = orignalImgreader->GetOutput();
            for(int z = 0; z < imgSize[2]; z++)
                for(int y = 0; y < imgSize[1]; y++)
                    for(int x = 0; x < imgSize[0]; x++)
                    {
                        voxelIndex[0] = x;
                        voxelIndex[1] = y;
                        voxelIndex[2] = z;
                        originalImgVoxelIndex[0] = x;
                        originalImgVoxelIndex[1] = y;
                        originalImgVoxelIndex[2] = z;
                        temp = inputImage->GetPixel(voxelIndex);
                        if(temp == 255)
                            originalImage->SetPixel(originalImgVoxelIndex, 65535);
                    }
                    //输出结果
                    typedef itk::ImageFileWriter<OriginalImageType>NewWriterType;
                    NewWriterType::Pointer writer = NewWriterType::New();
                    //文件前缀名
//.........这里部分代码省略.........

void segmentation::on_abnormalRegionSegButton_clicked()
{
    fstream fin11, fin12; 
    char* s11Txt = ("../data/ped/txt/15715_11.txt");//s11.txt
    char* s12Txt = ("../data/ped/txt/15715_12.txt");
    int (*s11)[512] = new int[64][512];
    int (*s12)[512] = new int[64][512];
    int (*p)[512] = new int[64][512];
    //read TXT
    fin11.open(s11Txt, ios::in); 
    fin12.open(s12Txt, ios::in);
    //get the p matrix
    for(int z = 0; z < 64; z++)
    {
        for(int x = 0; x < 512; x++)
        {
            fin11>>s11[z][x];
            fin12>>s12[z][x];
            if(s12[z][x] - s11[z][x] > 5)  //5 threshold  
                p[z][x] = 1;
            else 
                p[z][x] = 0;
        }
    }
    fin11.close();
    fin12.close();
    //read filtered MHD data
    typedef itk::Image<unsigned short, 3>InputImageType; 
    typedef itk::Image<unsigned short, 3>OutputImageType;
    typedef itk::Image<unsigned char, 3>InitOutputImageType;
    InputImageType::Pointer inputImage = InputImageType::New();
    OutputImageType::Pointer outputImage = OutputImageType::New();
    InitOutputImageType::Pointer initOutputImage = InitOutputImageType::New();
    typedef itk::ImageFileReader<InputImageType>ReaderType;
    ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName(inputFileName);//读入原图像
    reader->Update();
    inputImage = reader->GetOutput();
    InputImageType::IndexType voxelIndex;
    InitOutputImageType::IndexType initvoxelIndex;
    OutputImageType::IndexType newvoxelIndex;
    InputImageType::SizeType imgSize = inputImage->GetLargestPossibleRegion().GetSize();
    OutputImageType::IndexType index;
    index[0] = 0;
    index[1] = 0;
    index[2] = 0;
    OutputImageType::SizeType size;
    size[0] = imgSize[0];
    size[1] = imgSize[1];
    size[2] = imgSize[2];
    //create a region for initial result
    InitOutputImageType::RegionType initRegion;
    initRegion.SetIndex(index);
    initRegion.SetSize(size);
    initOutputImage->SetRegions( initRegion);
    initOutputImage->Allocate();
    //create a region for enhance result
    OutputImageType::RegionType region;
    region.SetIndex(index);
    region.SetSize(size);
    outputImage->SetRegions(region);
    outputImage->Allocate();
    //Initial result for PED segmentation  (a binary image)
    for(int z = 0; z < imgSize[2]; z++)
        for(int x = 0; x < imgSize[0]; x++)
        {
            initvoxelIndex[0] = x;
            initvoxelIndex[2] = z;
            for(int y = 0; y < imgSize[1]; y++)
            {
                //set all background a black region
                initvoxelIndex[1] = y;
                initOutputImage->SetPixel(initvoxelIndex, 0);
            }
            //set the same intensity for all PED region (empirical value)
            if(p[z][x] == 1)
            {
                for(int y = s11[z][x]; y <= s12[z][x]; y++) 
                {
                    initvoxelIndex[1] = y;
                    initOutputImage->SetPixel(initvoxelIndex, 255);//亮区域
                }
            }
        }
        //输出中间分割结果
        //文件前缀名
        filePrefix = inputFileName;//char* to string
        filePrefix = filePrefix.substr(0, filePrefix.length() - 4);
        string strInitFileName;
        strInitFileName = filePrefix + "_initBinaryImg.mhd";
        strcpy(initFileName, strInitFileName.c_str());//string to char*
        typedef itk::ImageFileWriter<InitOutputImageType>InitWriterType;
        InitWriterType::Pointer initWriter = InitWriterType::New();
        initWriter->SetFileName(initFileName);//生成二值图
        initWriter->SetInput(initOutputImage);
        initWriter->Update();
        //Enhance PED region and overlay it on the original image
        for(int z = 0; z < imgSize[2]; z++)
            for(int x = 0; x < imgSize[0]; x++)
            {
//.........这里部分代码省略.........

InputImageType::Pointer Curvelet::RunOnInputImage(InputImageType::Pointer InputImage)
{
	//InputImage = NewInputImage;
	slices = InputImage->GetLargestPossibleRegion().GetSize()[2];
	InputImageType::Pointer outputim = InputImageType::New();
	outputim->SetRegions(InputImage->GetLargestPossibleRegion());
	outputim->Allocate();
	FloatImageType::Pointer cosim = FloatImageType::New();
	cosim->SetRegions(InputImage->GetLargestPossibleRegion());
	cosim->Allocate();
	FloatImageType::Pointer sinim = FloatImageType::New();
	sinim->SetRegions(InputImage->GetLargestPossibleRegion());
	sinim->Allocate();

	if(outputim->GetBufferPointer() == NULL || cosim->GetBufferPointer() == NULL || sinim->GetBufferPointer() == NULL)
	{
		printf("Couldnt' allocate memory - 3.. going to crash now\n");
	}
	int max_dim = tile_size;


	int xsize = InputImage->GetLargestPossibleRegion().GetSize()[0];
	int ysize = InputImage->GetLargestPossibleRegion().GetSize()[1];

	int kx = 0;int ky = 0;


	kx = xsize /(max_dim-this->border);
	ky = ysize /(max_dim-this->border);

	int remx = xsize % (max_dim-this->border);
	int remy = ysize % (max_dim-this->border);

	if ( remx > 0 )
		kx ++;
	if ( remy > 0 )
		ky ++;

	for(int xco = 0; xco < kx; xco++)
	{
		for(int yco = 0; yco < ky; yco++)
		{

			InputImageType::SizeType imsize = InputImage->GetLargestPossibleRegion().GetSize();
			InputImageType::IndexType index;
			InputImageType::SizeType size;
			InputImageType::RegionType region;

			index.Fill(0);
			size[0] =  MIN((xco)*(max_dim-this->border)+max_dim-1,imsize[0]-1) -  xco * (max_dim-this->border) +1;
			size[1] =  MIN((yco)*(max_dim-this->border)+max_dim-1,imsize[1]-1) -  yco * (max_dim-this->border) +1;
			size[2] = imsize[2];

			InputImageType::Pointer imtile = InputImageType::New();
			region.SetIndex(index);
			region.SetSize(size);
			imtile->SetRegions(region);
			imtile->Allocate();
			if(imtile->GetBufferPointer()==NULL)
				printf("Couldn't allocate memory - 4 .. going to crash now\n");
			InputImageType::RegionType region1;
			index[0] = xco *(max_dim-this->border);
			index[1] = yco *(max_dim-this->border);
			index[2] = 0;
			region1.SetIndex(index);
			region1.SetSize(size);

			typedef itk::ImageRegionIterator<InputImageType> IteratorType;
			IteratorType iter1(InputImage,region1);
			IteratorType iter2(imtile,region);


			//printf("xco = %d yco = %d :\n",xco,yco);
			region1.Print(std::cout);
			region.Print(std::cout);

			iter1.GoToBegin();
			iter2.GoToBegin();
			for(;!iter1.IsAtEnd();++iter1,++iter2)
			{
				iter2.Set(iter1.Get());
			}


			InputImageType::Pointer outputtile = InputImageType::New();
			outputtile->SetRegions(imtile->GetLargestPossibleRegion());
			outputtile->Allocate();
			FloatImageType::Pointer cosimtile = FloatImageType::New();
			cosimtile->SetRegions(imtile->GetLargestPossibleRegion());
			cosimtile->Allocate();
			FloatImageType::Pointer sinimtile = FloatImageType::New();
			sinimtile->SetRegions(imtile->GetLargestPossibleRegion());
			sinimtile->Allocate();
			if(outputtile->GetBufferPointer() == NULL || cosimtile->GetBufferPointer()==NULL || sinimtile->GetBufferPointer() == NULL )
			{
				printf("Couldn't allocate memory - 5 .. going to crash now ..\n");
			}

			{
#pragma omp parallel for shared(cosimtile,imtile,sinimtile,outputtile)  num_threads(numt)
//.........这里部分代码省略.........