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

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

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

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

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

  /*
   * 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();

}

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

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

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

void MakeDices(const char * montagefileName, const char * seedfileName, int diceWidth, 
			   double alfa, double beta, int timethreshold, double curvatureScaling, double rmsThres, int holeSize, int minObjSize)
{
	typedef SomaExtractor::ProbImageType ImageType;
	typedef SomaExtractor::OutputImageType OutputImageType;
	typedef itk::RegionOfInterestImageFilter< ImageType, ImageType> RegionOfInterestFilter;
	SomaExtractor *Somas = new SomaExtractor();

	std::cout<< "ReadMontage..."<<std::endl;
	ImageType::Pointer image = Somas->SetInputImage(montagefileName);
	int SZ = image->GetLargestPossibleRegion().GetSize()[2];

	std::vector< itk::Index<3> > seedVector;
	Somas->ReadSeedpoints( seedfileName, seedVector, 0);
	std::cout<< "Seed points size:"<< seedVector.size()<<std::endl;	

	#pragma omp parallel for
	for(int i = 0; i < seedVector.size(); i++)
	{
		SomaExtractor *somaExtractor = new SomaExtractor();
		ImageType::IndexType start;
		start[0] = seedVector[i][0] - diceWidth / 2;
		start[1] = seedVector[i][1] - diceWidth / 2;
		start[2] = 0;
		ImageType::SizeType size;
		size[0] = diceWidth;
		size[1] = diceWidth;
		size[2] = SZ;

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

		RegionOfInterestFilter::Pointer regionFilter = RegionOfInterestFilter::New();
		regionFilter->SetInput(image);
		regionFilter->SetRegionOfInterest(desiredRegion);
		regionFilter->Update();
		ImageType::Pointer diceImage = regionFilter->GetOutput();
		std::ostringstream oss1;
		oss1<< "Dice_"<< i<<".tif";
		std::string str1 = oss1.str();
		somaExtractor->writeImage(str1.c_str(), diceImage);
		
		std::vector< itk::Index<3> > seedsForDice;
		itk::Index<3> seedIndex = seedVector[i];
		seedIndex[0] = diceWidth / 2;
		seedIndex[1] = diceWidth / 2;
		seedsForDice.push_back(seedIndex);

		float Origion[3] = {0,0,0};
		diceImage->SetOrigin(Origion);

		//SomaExtractor::SegmentedImageType::Pointer segImage= somaExtractor->SegmentSoma(diceImage, seedsForDice, alfa, 
		//					beta, timethreshold, curvatureScaling, rmsThres, holeSize, minObjSize);
		//double threshold = 0;
//		//ImageType::Pointer segImage = Somas->EnhanceContrast(diceImage, seedIndex[2], alfa, beta, threshold);
//
//#pragma omp critical
//		std::cout<< "segment "<< i<<"\t"<<seedsForDice.size()<<std::endl;
//		std::ostringstream oss;
//		oss<< "Dice_Segment_"<< i<<".tif";
//		std::string str = oss.str();
//		somaExtractor->writeImage(str.c_str(), segImage);
		delete somaExtractor;
	}
	delete Somas;
}

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

	
}