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

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

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

        virtual ReadResult readImage(const std::string& file, const osgDB::ReaderWriter::Options* options) const
        {
            std::string ext = osgDB::getLowerCaseFileExtension(file);
            if (!acceptsExtension(ext)) return ReadResult::FILE_NOT_HANDLED;

            std::string fileName = osgDB::findDataFile( file, options );
            if (fileName.empty()) return ReadResult::FILE_NOT_FOUND;

            notice()<<"Reading DICOM file "<<fileName<<std::endl;

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


            typedef itk::ImageFileReader< ImageType > ReaderType;
            ReaderType::Pointer reader = ReaderType::New();
            reader->SetFileName( fileName.c_str() );

            typedef itk::GDCMImageIO           ImageIOType;
            ImageIOType::Pointer gdcmImageIO = ImageIOType::New();
            reader->SetImageIO( gdcmImageIO );

            try
            {
                reader->Update();
            }
            catch (itk::ExceptionObject & e)
            {
                std::cerr << "exception in file reader " << std::endl;
                std::cerr << e.GetDescription() << std::endl;
                std::cerr << e.GetLocation() << std::endl;
                return ReadResult::ERROR_IN_READING_FILE;
            }
            
            ImageType::Pointer inputImage = reader->GetOutput();
            
            ImageType::RegionType region = inputImage->GetBufferedRegion();
            ImageType::SizeType size = region.GetSize();
            ImageType::IndexType start = region.GetIndex();
            
            //inputImage->GetSpacing();
            //inputImage->GetOrigin();
            
            unsigned int width = size[0];
            unsigned int height = size[1];
            unsigned int depth = size[2];
            
            osg::RefMatrix* matrix = new osg::RefMatrix;
            
            notice()<<"width = "<<width<<" height = "<<height<<" depth = "<<depth<<std::endl;
            for(unsigned int i=0; i<Dimension; ++i)
            {
                (*matrix)(i,i) = inputImage->GetSpacing()[i];
                (*matrix)(3,i) = inputImage->GetOrigin()[i];
            }
            
            osg::Image* image = new osg::Image;
            image->allocateImage(width, height, depth, GL_LUMINANCE, GL_UNSIGNED_BYTE, 1);

            unsigned char* data = image->data();            
            typedef itk::ImageRegionConstIterator< ImageType > IteratorType;
            IteratorType it(inputImage, region);
            
            it.GoToBegin();
            while (!it.IsAtEnd())
            {
                *data = it.Get();
                ++data;
                ++it;
            }
            
            image->setUserData(matrix);
            
            matrix->preMult(osg::Matrix::scale(double(image->s()), double(image->t()), double(image->r())));

            return image;
        }

int main(int argc, char ** argv)
{
	std::string folderName = argv[1];
	std::string filter = argv[2];

	typedef utils::Directory Directory;
	Directory::FilenamesType filenames = Directory::GetFiles(folderName, ".nrrd");

	Directory::FilenamesType goodFilenames;
	filterFilenames(filenames, filter, goodFilenames);

	// load the images
	std::vector<ImageType::Pointer> imageList;
	for(unsigned int i = 0; i < goodFilenames.size(); i++)
	{
		std::cout << goodFilenames[i] << std::endl;

		typedef itk::ImageFileReader<ImageType> ReaderType;

		ReaderType::Pointer reader = ReaderType::New();
		reader->SetFileName(goodFilenames[i]);
		reader->SetImageIO(itk::NrrdImageIO::New());
		reader->Update();

		imageList.push_back(reader->GetOutput());

	}

	std::sort(imageList.begin(), imageList.end(), imageSort);

	ImageType::Pointer outputImage = ImageType::New();
	ImageType::RegionType outputRegion;
	ImageType::SizeType outputSize = imageList.front()->GetLargestPossibleRegion().GetSize();
	outputSize[2] = imageList.size(); 
	ImageType::IndexType outputIndex;
	outputIndex.Fill(0);
	outputRegion.SetSize(outputSize);
	outputRegion.SetIndex(outputIndex);

	// compute the spacing 
	double meanSpacing = 0.0;
	for(unsigned int i = 1; i < imageList.size(); i++)
	{
		ImageType::Pointer im1 = imageList[i-1];
		ImageType::Pointer im2 = imageList[i];

		meanSpacing += (im1->GetOrigin()-im2->GetOrigin()).GetNorm();
	}

	meanSpacing /= (double) (imageList.size()-1);
	ImageType::SpacingType spacing = imageList.front()->GetSpacing();
	spacing[2] = meanSpacing;

	outputImage->SetRegions(outputRegion);
	outputImage->SetSpacing(spacing);
	outputImage->SetOrigin(imageList.front()->GetOrigin());
	outputImage->SetDirection(imageList.front()->GetDirection());
	outputImage->Allocate();

	itk::ImageRegionIterator<ImageType> outIt(outputImage, outputImage->GetLargestPossibleRegion());
	outIt.GoToBegin();

	while(!outIt.IsAtEnd())
	{
		ImageType::IndexType index = outIt.GetIndex();
		ImageType::IndexType testIndex = index;
		testIndex[2] = 0;

		outIt.Set(imageList[index[2]]->GetPixel(testIndex));


		++outIt;
	}

	std::string stackFilename = folderName + "/stack.nrrd";


	typedef itk::ImageFileWriter<ImageType> WriterType;
	WriterType::Pointer writer = WriterType::New();
	writer->SetInput(outputImage);
	writer->SetFileName(stackFilename);
	writer->SetImageIO(itk::NrrdImageIO::New());
	writer->Update();

	return 0;

}

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

  //check that the user specified the right number of command line arguments
  if(argc < 3)
    {
    cerr << argv[0] << " <input file> <output file>" << endl;
    cerr << argv[0] << " <raw input file> <sizeX> <sizeY> <sizeZ> <output file>"
         << endl;
    return EXIT_FAILURE;
    }

  //check if the input image is .raw
  bool rawInput = false;
  string inputFileName = argv[1];
  const char *outputFileName;
  //double space[3]={1,1,1};
  if(inputFileName.rfind(".raw") != string::npos)
    {
    //if so, the user is also required to pass in image dimensions
    if(argc < 6)
      {
      cerr << "Usage: <raw input file> <sizeX> <sizeY> <sizeZ> <output file>" << endl;
      return EXIT_FAILURE;
      }
    rawInput = true;
    sizeX = atoi(argv[2]);
    sizeY = atoi(argv[3]);
    sizeZ = atoi(argv[4]);
    outputFileName = argv[5];
    }
  else
    {
    outputFileName = argv[2];
    }

  FILE *infile = 0;
  FILE *outfile;
  int i,j,k;
  int ii, jj, kk;
  DATATYPEOUT *volout;
  long idx;

 
  int sizeExpand = 0;
  DATATYPEOUT blockMax;
  int timesDilate;
  int border;
  unsigned short *GraphcutResults;

  cout << "Volume Processing..." << endl;
  //make sure we can write to the output file
  if((outfile=fopen(outputFileName, "wb")) == NULL)
    {
    cerr << "Output file open error!" << endl;
    return EXIT_FAILURE;
    }

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

   if(rawInput)
    {
    if((infile=fopen(inputFileName.c_str(), "rb"))==NULL)
      {
      cout << "Input file open error!" << endl;
      return EXIT_FAILURE;
      }

    volin = (DATATYPEIN*)malloc(sizeX*sizeY*(sizeZ+sizeExpand*2)*sizeof(DATATYPEIN));

    if (fread(volin, sizeof(DATATYPEIN), sizeX*sizeY*sizeZ, infile) < (unsigned int)(sizeX*sizeY*sizeZ))
      {
      cerr << "File size is not the same as volume size" << endl;
      return EXIT_FAILURE;
      }
    
	inputImage = ImageType::New();

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

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

	ImageType::SizeType  size;
    size[0]  = sizeX;  
    size[1]  = sizeY;  
	size[2]  = sizeZ;  

    ImageType::RegionType region;
    region.SetSize( size );
//.........这里部分代码省略.........