C++ transformtype::Pointer类代码示例

int Invert( int argc , char* argv[] )
{
  //check arguments
  if( argc != 4 )
  {
    std::cout<< argv[ 0 ] << " " << argv[ 1 ] << " inputTransform outputTransform" << std::endl ;
    return 1 ;
  }
  typedef itk::MatrixOffsetTransformBase< double , 3 , 3 > TransformType ;
  itk::TransformFactory< TransformType >::RegisterTransform();
  itk::TransformFileReader::Pointer transformFile = itk::TransformFileReader::New() ;
  transformFile->SetFileName( argv[ 2 ] ) ;
  transformFile->Update() ;
  if( transformFile->GetTransformList()->size() != 1 )
  {
    std::cerr << "Please give a transform file containing only one transformation" << std::endl ;
    return 1 ;
  }
  TransformType::Pointer transform ;
  transform = dynamic_cast< TransformType* >
        ( transformFile->GetTransformList()->front().GetPointer() ) ;
  if( !transform )
  {
    std::cerr << "Transform type is not handled. Please convert your transform first" << std::endl ;
    return 1 ;
  }
  TransformType::Pointer inverse = TransformType::New() ;
  transform->GetInverse( inverse ) ;
  itk::TransformFileWriter::Pointer transformWriter = itk::TransformFileWriter::New() ;
  transformWriter->SetFileName( argv[ 3 ] ) ;
  transformWriter->AddTransform( inverse ) ;
  transformWriter->Update() ;
  return 0 ;
}

mitk::Vector3D
  mitk::SlicedGeometry3D::AdjustNormal( const mitk::Vector3D &normal ) const
{
  TransformType::Pointer inverse = TransformType::New();
  m_ReferenceGeometry->GetIndexToWorldTransform()->GetInverse( inverse );

  Vector3D transformedNormal = inverse->TransformVector( normal );

  transformedNormal.Normalize();
  return transformedNormal;
}

int main(int argc, char ** argv)
{
	// load the image and the bounding box
	BoundingBox::Pointer boundingBox = BoundingBox::New();
	boundingBox->SetInfation(atof(argv[3]));
	boundingBox->Load(argv[1]);
	
	// load hte images and compute the reference coordinates
	CMRFileExtractor::Pointer extractor = CMRFileExtractor::New();
	extractor->SetFolderName(argv[2]);
	extractor->Extract();

	ValveOriginFinder::Pointer originFinder = ValveOriginFinder::New();
	originFinder->Set2CImage(extractor->Get2CImage(0));
	originFinder->Set3CImage(extractor->Get3CImage(0));
	originFinder->SetImageStack(extractor->GetStackImage(0));
	originFinder->Compute();


	// apply the transform to the bounding box
	typedef itk::Similarity3DTransform<double> TransformType;
	TransformType::Pointer transform = TransformType::New();
	transform->SetMatrix(originFinder->GetRotation());
	transform->SetTranslation(originFinder->GetTranslation());
	boundingBox->TransformBoundingBox(transform);


	BoundingBox::MaskType::Pointer mask = BoundingBox::MaskType::New();
	boundingBox->ComputeImageMask(extractor->Get2CImage(0), 1, mask);


	utils::LabelVolumeIO::Write("mask.nrrd", mask);
	utils::ImageVolumeIO::Write("image.nrrd", extractor->Get2CImage(0));


	SimpleMRFSegmenter::Pointer segmenter = SimpleMRFSegmenter::New();
	segmenter->SetImage(extractor->Get2CImage(0));
	segmenter->SetSmoothnessCost(atof(argv[4]));
	segmenter->SetMask(mask);
	segmenter->Segment();


	utils::LabelVolumeIO::Write("seg.nrrd", segmenter->GetOutput());


	return 0;
}

void ImageScaleTransform::process()
    {
    foreach (const ElementBase *source, mSourceElementsReadySet)
        for (int i = 0; i < source->getFramesNo(); ++i)
            {
            const FrameBase *frame = source->getFrame(i);
            if (frame->getMaxDimension() == ColorImageFrame::Dimensions)
                {
                mImageFrame.setSourceName(frame->getSourceName());
                mSrcFrame.resizeAndCopyFrame(*frame);
                ColorImageFrame::ImageType::Pointer srcImg = mSrcFrame;

                typedef ScaleTransform<double, 2> TransformType;
                TransformType::Pointer scaleTransform = TransformType::New();
                FixedArray<float, 2> scale;
                scale[0] = property("widthScale").toDouble();
                scale[1] = property("heightScale").toDouble();
                scaleTransform->SetScale(scale);
                Point<float, 2> center;
                center[0] = srcImg->GetLargestPossibleRegion().GetSize()[0]/2;
                center[1] = srcImg->GetLargestPossibleRegion().GetSize()[1]/2;
                scaleTransform->SetCenter(center);

                typedef ResampleImageFilter<ColorImageFrame::ImageType, ColorImageFrame::ImageType> ResampleImageFilterType;
                ResampleImageFilterType::Pointer resampleFilter = ResampleImageFilterType::New();
                resampleFilter->SetTransform(scaleTransform);
                resampleFilter->SetInput(srcImg);
                resampleFilter->SetSize(srcImg->GetLargestPossibleRegion().GetSize());
                resampleFilter->Update();
                mImageFrame = resampleFilter->GetOutput();

                emit framesReady();
                break;
                }
            }
    }

void QAngioSubstractionExtension::computeAutomateSingleImage()
{
    QApplication::setOverrideCursor(Qt::WaitCursor);
    const    unsigned int          Dimension = 2;
    typedef  Volume::ItkPixelType  PixelType;

    typedef itk::Image< PixelType, Dimension >  FixedImageType;
    typedef itk::Image< PixelType, Dimension >  MovingImageType;
    typedef   float     InternalPixelType;
    typedef itk::Image< InternalPixelType, Dimension > InternalImageType;

    typedef itk::TranslationTransform< double, Dimension > TransformType;
    typedef itk::GradientDescentOptimizer                  OptimizerType;
    typedef itk::LinearInterpolateImageFunction< 
                                    InternalImageType,
                                    double             > InterpolatorType;
    typedef itk::ImageRegistrationMethod< 
                                    InternalImageType, 
                                    InternalImageType >  RegistrationType;
    typedef itk::MutualInformationImageToImageMetric< 
                                          InternalImageType, 
                                          InternalImageType >    MetricType;

    TransformType::Pointer      transform     = TransformType::New();
    OptimizerType::Pointer      optimizer     = OptimizerType::New();
    InterpolatorType::Pointer   interpolator  = InterpolatorType::New();
    RegistrationType::Pointer   registration  = RegistrationType::New();

    registration->SetOptimizer(optimizer);
    registration->SetTransform(transform);
    registration->SetInterpolator(interpolator);

    MetricType::Pointer         metric        = MetricType::New();
    registration->SetMetric(metric);
    metric->SetFixedImageStandardDeviation(0.4);
    metric->SetMovingImageStandardDeviation(0.4);
    metric->SetNumberOfSpatialSamples(50);

    typedef itk::ExtractImageFilter< Volume::ItkImageType, FixedImageType > FilterType;
    
    FilterType::Pointer extractFixedImageFilter = FilterType::New();
    Volume::ItkImageType::RegionType inputRegion = m_mainVolume->getItkData()->GetLargestPossibleRegion();
    Volume::ItkImageType::SizeType size = inputRegion.GetSize();
    //Dividim la mida per dos per tal de quedar-nos només amb la part central
    // ja que si no ens registre el background
    size[0] = size[0] / 2;
    size[1] = size[1] / 2;
    size[2] = 0;
    Volume::ItkImageType::IndexType start = inputRegion.GetIndex();
    const unsigned int sliceReference = m_imageSelectorSpinBox->value();
    //comencem a un quart de la imatge
    start[0] = size[0] / 2;
    start[1] = size[1] / 2;
    start[2] = sliceReference;
    Volume::ItkImageType::RegionType desiredRegion;
    desiredRegion.SetSize(size);
    desiredRegion.SetIndex(start);
    extractFixedImageFilter->SetExtractionRegion(desiredRegion);
    extractFixedImageFilter->SetInput(m_mainVolume->getItkData());
    extractFixedImageFilter->Update();

    FilterType::Pointer extractMovingImageFilter = FilterType::New();
    Volume::ItkImageType::IndexType startMoving = inputRegion.GetIndex();
    const unsigned int sliceNumber = m_2DView_1->getViewer()->getCurrentSlice();
    startMoving[0] = size[0] / 2;
    startMoving[1] = size[1] / 2;
    startMoving[2] = sliceNumber;
    Volume::ItkImageType::RegionType desiredMovingRegion;
    desiredMovingRegion.SetSize(size);
    desiredMovingRegion.SetIndex(startMoving);
    extractMovingImageFilter->SetExtractionRegion(desiredMovingRegion);
    extractMovingImageFilter->SetInput(m_mainVolume->getItkData());
    extractMovingImageFilter->Update();

    typedef itk::NormalizeImageFilter< 
                                FixedImageType, 
                                InternalImageType 
                                        > FixedNormalizeFilterType;

    typedef itk::NormalizeImageFilter< 
                                MovingImageType, 
                                InternalImageType 
                                              > MovingNormalizeFilterType;

    FixedNormalizeFilterType::Pointer fixedNormalizer = 
                                            FixedNormalizeFilterType::New();

    MovingNormalizeFilterType::Pointer movingNormalizer =
                                            MovingNormalizeFilterType::New();
    typedef itk::DiscreteGaussianImageFilter<
                                      InternalImageType, 
                                      InternalImageType
                                                    > GaussianFilterType;

    GaussianFilterType::Pointer fixedSmoother  = GaussianFilterType::New();
    GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();

    fixedSmoother->SetVariance(2.0);
    movingSmoother->SetVariance(2.0);
    fixedNormalizer->SetInput(extractFixedImageFilter->GetOutput());
//.........这里部分代码省略.........

bool ShowSegmentationAsSmoothedSurface::ThreadedUpdateFunction()
{
  Image::Pointer image;
  GetPointerParameter("Input", image);

  float smoothing;
  GetParameter("Smoothing", smoothing);

  float decimation;
  GetParameter("Decimation", decimation);

  float closing;
  GetParameter("Closing", closing);

  int timeNr = 0;
  GetParameter("TimeNr", timeNr);

  if (image->GetDimension() == 4)
    MITK_INFO << "CREATING SMOOTHED POLYGON MODEL (t = " << timeNr << ')';
  else
    MITK_INFO << "CREATING SMOOTHED POLYGON MODEL";

  MITK_INFO << "  Smoothing  = " << smoothing;
  MITK_INFO << "  Decimation = " << decimation;
  MITK_INFO << "  Closing    = " << closing;

  Geometry3D::Pointer geometry = dynamic_cast<Geometry3D *>(image->GetGeometry()->Clone().GetPointer());

  // Make ITK image out of MITK image

  typedef itk::Image<unsigned char, 3> CharImageType;
  typedef itk::Image<unsigned short, 3> ShortImageType;
  typedef itk::Image<float, 3> FloatImageType;

  if (image->GetDimension() == 4)
  {
    ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New();
    imageTimeSelector->SetInput(image);
    imageTimeSelector->SetTimeNr(timeNr);
    imageTimeSelector->UpdateLargestPossibleRegion();
    image = imageTimeSelector->GetOutput(0);
  }

  ImageToItk<CharImageType>::Pointer imageToItkFilter = ImageToItk<CharImageType>::New();

  try
  {
    imageToItkFilter->SetInput(image);
  }
  catch (const itk::ExceptionObject &e)
  {
    // Most probably the input image type is wrong. Binary images are expected to be
    // >unsigned< char images.
    MITK_ERROR << e.GetDescription() << endl;
    return false;
  }

  imageToItkFilter->Update();

  CharImageType::Pointer itkImage = imageToItkFilter->GetOutput();

  // Get bounding box and relabel

  MITK_INFO << "Extracting VOI...";

  int imageLabel = 1;
  bool roiFound = false;

  CharImageType::IndexType minIndex;
  minIndex.Fill(numeric_limits<CharImageType::IndexValueType>::max());

  CharImageType::IndexType maxIndex;
  maxIndex.Fill(numeric_limits<CharImageType::IndexValueType>::min());

  itk::ImageRegionIteratorWithIndex<CharImageType> iter(itkImage, itkImage->GetLargestPossibleRegion());

  for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter)
  {
    if (iter.Get() == imageLabel)
    {
      roiFound = true;
      iter.Set(1);

      CharImageType::IndexType currentIndex = iter.GetIndex();

      for (unsigned int dim = 0; dim < 3; ++dim)
      {
        minIndex[dim] = min(currentIndex[dim], minIndex[dim]);
        maxIndex[dim] = max(currentIndex[dim], maxIndex[dim]);
      }
    }
    else
    {
      iter.Set(0);
    }
  }

  if (!roiFound)
  {
    ProgressBar::GetInstance()->Progress(8);
//.........这里部分代码省略.........

// perform B-spline registration for 2D image
void runBspline2D(StringVector& args) {
    typedef itk::BSplineTransform<double, 2, 3> TransformType;
    typedef itk::LBFGSOptimizer OptimizerType;
    typedef itk::MeanSquaresImageToImageMetric<RealImage2, RealImage2> MetricType;
    typedef itk:: LinearInterpolateImageFunction<RealImage2, double> InterpolatorType;
    typedef itk::ImageRegistrationMethod<RealImage2, RealImage2> RegistrationType;

    MetricType::Pointer         metric        = MetricType::New();
    OptimizerType::Pointer      optimizer     = OptimizerType::New();
    InterpolatorType::Pointer   interpolator  = InterpolatorType::New();
    RegistrationType::Pointer   registration  = RegistrationType::New();

    // The old registration framework has problems with multi-threading
    // For now, we set the number of threads to 1
    registration->SetNumberOfThreads(1);

    registration->SetMetric(        metric        );
    registration->SetOptimizer(     optimizer     );
    registration->SetInterpolator(  interpolator  );

    TransformType::Pointer  transform = TransformType::New();
    registration->SetTransform( transform );


    ImageIO<RealImage2> io;

    // Create the synthetic images
    RealImage2::Pointer  fixedImage  = io.ReadImage(args[0]);
    RealImage2::Pointer  movingImage  = io.ReadImage(args[1]);

    // Setup the registration
    registration->SetFixedImage(  fixedImage   );
    registration->SetMovingImage(   movingImage);

    RealImage2::RegionType fixedRegion = fixedImage->GetBufferedRegion();
    registration->SetFixedImageRegion( fixedRegion );

    TransformType::PhysicalDimensionsType   fixedPhysicalDimensions;
    TransformType::MeshSizeType             meshSize;
    for( unsigned int i=0; i < 2; i++ )
    {
        fixedPhysicalDimensions[i] = fixedImage->GetSpacing()[i] *
        static_cast<double>(
                            fixedImage->GetLargestPossibleRegion().GetSize()[i] - 1 );
    }
    unsigned int numberOfGridNodesInOneDimension = 18;
    meshSize.Fill( numberOfGridNodesInOneDimension - 3 );
    transform->SetTransformDomainOrigin( fixedImage->GetOrigin() );
    transform->SetTransformDomainPhysicalDimensions( fixedPhysicalDimensions );
    transform->SetTransformDomainMeshSize( meshSize );
    transform->SetTransformDomainDirection( fixedImage->GetDirection() );

    typedef TransformType::ParametersType     ParametersType;

    const unsigned int numberOfParameters =
    transform->GetNumberOfParameters();

    ParametersType parameters( numberOfParameters );

    parameters.Fill( 0.0 );

    transform->SetParameters( parameters );

    //  We now pass the parameters of the current transform as the initial
    //  parameters to be used when the registration process starts.

    registration->SetInitialTransformParameters( transform->GetParameters() );

    std::cout << "Intial Parameters = " << std::endl;
    std::cout << transform->GetParameters() << std::endl;

    //  Next we set the parameters of the LBFGS Optimizer.

    optimizer->SetGradientConvergenceTolerance( 0.005 );
    optimizer->SetLineSearchAccuracy( 0.9 );
    optimizer->SetDefaultStepLength( .1 );
    optimizer->TraceOn();
    optimizer->SetMaximumNumberOfFunctionEvaluations( 1000 );

    std::cout << std::endl << "Starting Registration" << std::endl;

    try
    {
        registration->Update();
        std::cout << "Optimizer stop condition = "
        << registration->GetOptimizer()->GetStopConditionDescription()
        << std::endl;
    }
    catch( itk::ExceptionObject & err )
    {
        std::cerr << "ExceptionObject caught !" << std::endl;
        std::cerr << err << std::endl;
        return;
    }

    OptimizerType::ParametersType finalParameters =
    registration->GetLastTransformParameters();

    std::cout << "Last Transform Parameters" << std::endl;
//.........这里部分代码省略.........

/**
 * @brief      3D resample data to new grid size
 *
 * @param  M   Incoming data
 * @param  f   Resampling factor in all 3 dimensions
 * @param  im  Interpolation method (LINEAR|BSPLINE)
 *
 * @return     Resampled data
 */
template<class T> static Matrix<T> 
resample (const Matrix<T>& M, const Matrix<double>& f, const InterpMethod& im) {
	

	Matrix <T> res = M;
	
#ifdef HAVE_INSIGHT
	
	typedef typename itk::OrientedImage< T, 3 > InputImageType;
	typedef typename itk::OrientedImage< T, 3 > OutputImageType;
	typedef typename itk::IdentityTransform< double, 3 > TransformType;
	typedef typename itk::LinearInterpolateImageFunction< InputImageType, double > InterpolatorType;
	typedef typename itk::ResampleImageFilter< InputImageType, InputImageType > ResampleFilterType;
	
	typename InterpolatorType::Pointer linterp = InterpolatorType::New();
	
	TransformType::Pointer trafo = TransformType::New();
	trafo->SetIdentity();
	
	typename InputImageType::SpacingType space;
	space[0] = 1.0/f[0];
	space[1] = 1.0/f[1];
	space[2] = 1.0/f[2];
	
	typedef typename InputImageType::SizeType::SizeValueType SizeValueType;
	typename InputImageType::SizeType size; 
	size[0] = static_cast<SizeValueType>(res.Dim(0));
	size[1] = static_cast<SizeValueType>(res.Dim(1));
	size[2] = static_cast<SizeValueType>(res.Dim(2));
	
	typename itk::OrientedImage< T, 3 >::Pointer input = itk::OrientedImage< T, 3 >::New();
	typename itk::OrientedImage< T, 3 >::Pointer output = itk::OrientedImage< T, 3 >::New();
	
	typename itk::Image< T, 3 >::IndexType ipos;
	ipos[0] = 0; ipos[1] = 0; ipos[2] = 0;
	typename itk::Image< T, 3 >::IndexType opos;
	opos[0] = 0; opos[1] = 0; opos[2] = 0;
	
	typename itk::Image< T, 3 >::RegionType ireg;
	ireg.SetSize(size);
	ireg.SetIndex(ipos);
	input->SetRegions(ireg);
	input->Allocate();
	
	typename itk::Image< T, 3 >::RegionType oreg;
	oreg.SetSize(size);
	ireg.SetIndex(opos);
	output->SetRegions(oreg);
	output->Allocate();
	
	for (size_t z = 0; z < res.Dim(2); z++)
		for (size_t y = 0; y < res.Dim(1); y++)
			for (size_t x = 0; x < res.Dim(0); x++) {
				ipos[0] = x; ipos[1] = y; ipos[2] = z;
				input->SetPixel (ipos, res.At(x,y,z));
			}
	
	typename ResampleFilterType::Pointer rs = ResampleFilterType::New();
	rs->SetInput( input );
	rs->SetTransform( trafo );
	rs->SetInterpolator( linterp );
	rs->SetOutputOrigin ( input->GetOrigin());
	rs->SetOutputSpacing ( space );
	rs->SetOutputDirection ( input->GetDirection());
	rs->SetSize ( size );
	rs->Update ();
	
	output = rs->GetOutput();
	
	res = Matrix<T> (res.Dim(0)*f[0], res.Dim(1)*f[1], res.Dim(2)*f[2]);
	res.Res(0) = res.Res(0)/f[0];
	res.Res(1) = res.Dim(1)/f[1];
	res.Res(2) = res.Dim(2)/f[2];
	
	for (size_t z = 0; z < res.Dim(2); z++)
		for (size_t y = 0; y < res.Dim(1); y++)
			for (size_t x = 0; x < res.Dim(0); x++) {
				opos[0] = x; opos[1] = y; opos[2] = z;
				res.At(x,y,z) = output->GetPixel (opos);
			}
	
#else 
	
	printf ("ITK ERROR - Resampling not performed without ITK!\n");
	
#endif
	
	return res;
	
}

int main (int argc, char **argv)
{
  int verbose=0, clobber=0,skip_grid=0;
  int order=2;
  std::string like_f,xfm_f,output_f,input_f;
  double uniformize=0.0;
  int invert=0;
  char *history = time_stamp(argc, argv); 
  
  static struct option long_options[] = {
		{"verbose", no_argument,       &verbose, 1},
		{"quiet",   no_argument,       &verbose, 0},
		{"clobber", no_argument,       &clobber, 1},
		{"like",    required_argument, 0, 'l'},
		{"transform",    required_argument, 0, 't'},
    {"order",    required_argument, 0, 'o'},
    {"uniformize",    required_argument, 0, 'u'},
    {"invert_transform", no_argument, &invert, 1},
		{0, 0, 0, 0}
		};
  
  for (;;) {
      /* getopt_long stores the option index here. */
      int option_index = 0;

      int c = getopt_long (argc, argv, "vqcl:t:o:u:", long_options, &option_index);

      /* Detect the end of the options. */
      if (c == -1) break;

      switch (c)
			{
			case 0:
				break;
			case 'v':
				cout << "Version: 0.1" << endl;
				return 0;
      case 'l':
        like_f=optarg; break;
      case 't':
        xfm_f=optarg; break;
      case 'o':
        order=atoi(optarg);break;
      case 'u':
        uniformize=atof(optarg);break;
			case '?':
				/* getopt_long already printed an error message. */
			default:
				show_usage (argv[0]);
				return 1;
			}
    }

	if ((argc - optind) < 2) {
		show_usage(argv[0]);
		return 1;
	}
  input_f=argv[optind];
  output_f=argv[optind+1];
  
  if (!clobber && !access (output_f.c_str (), F_OK))
  {
    std::cerr << output_f.c_str () << " Exists!" << std::endl;
    return 1;
  }
  
	try
  {
    itk::ObjectFactoryBase::RegisterFactory(itk::MincImageIOFactory::New());
    itk::ImageFileReader<minc::image3d >::Pointer reader = itk::ImageFileReader<minc::image3d >::New();
    
    //initializing the reader
    reader->SetFileName(input_f.c_str());
    reader->Update();
    
		minc::image3d::Pointer in=reader->GetOutput();

		FilterType::Pointer filter  = FilterType::New();
    
    //creating coordinate transformation objects
		TransformType::Pointer transform = TransformType::New();
    if(!xfm_f.empty())
    {
      //reading a minc style xfm file
      transform->OpenXfm(xfm_f.c_str());
      if(!invert) transform->Invert(); //should be inverted by default to walk through target space
      filter->SetTransform( transform );
    }

    //creating the interpolator
		InterpolatorType::Pointer interpolator = InterpolatorType::New();
		interpolator->SetSplineOrder(order);
		filter->SetInterpolator( interpolator );
		filter->SetDefaultPixelValue( 0 );
    
    //this is for processing using batch system
    filter->SetNumberOfThreads(1);
    
    if(!like_f.empty())
    {
//.........这里部分代码省略.........

int main( int argc, char *argv[] )
{
  string input_name;
  string output_dir;
  if (argc == 3) {
    input_name = argv[1];
    output_dir = argv[2];
  }

  const     unsigned int   Dimension = 3;
  const     unsigned int   OutDimension = 2;
  typedef short InputPixelType;
  typedef int FilterPixelType;
  typedef itk::Image< InputPixelType,  Dimension >   InputImageType;
  typedef itk::Image< FilterPixelType, Dimension >   FilterImageType;
  typedef itk::Image< FilterPixelType, OutDimension >   OutFilterImageType;

  InputImageType::Pointer image;
  itk::MetaDataDictionary dict;


  if (input_name.size() && output_dir.size()) 
    {
      if (boost::filesystem::is_regular_file( input_name )) {
	typedef itk::ImageFileReader< InputImageType >  ReaderType;
	ReaderType::Pointer reader = ReaderType::New();
	reader->SetFileName( input_name );
	try 
	  { 
	  reader->Update();
	  } 
	catch( itk::ExceptionObject & err ) 
	  { 
	  std::cerr << "ERROR: ExceptionObject caught !" << std::endl; 
	  std::cerr << err << std::endl; 
	  return EXIT_FAILURE;
	  } 
	image = reader->GetOutput();
	dict = reader->GetMetaDataDictionary();
      } else if (boost::filesystem::is_directory( input_name )) {
        itkBasic::SeriesReader sreader( input_name );
	sreader.readSeriesData( 2 );
	try 
	{
	    itkBasic::ReaderType::Pointer imageReader = itkBasic::ReaderType::New();
	    itkBasic::FileNamesContainer fc;
	    sreader.getSeriesFileNames(0, fc);
	    image = itkBasic::getDicomSerie( fc, imageReader, 1 ); 
	    dict = *((*imageReader->GetMetaDataDictionaryArray())[0]);
	}
	catch( itk::ExceptionObject & err ) 
	  { 
	  std::cerr << "ERROR: ExceptionObject caught !" << std::endl; 
	  std::cerr << err << std::endl; 
	  return EXIT_FAILURE;
	  } 
      }
    }
    
    if (!image) {
	std::cerr << argv[0] << ": input output" << std::endl;
	exit(1);
    }
  

  typedef itk::SigmoidImageFilter< InputImageType, FilterImageType > SigmoidCasterType;
  SigmoidCasterType::Pointer sigmoidcaster = SigmoidCasterType::New();
  
  sigmoidcaster->SetInput( image );
  sigmoidcaster->SetOutputMaximum( 4000 );
  sigmoidcaster->SetOutputMinimum( 1000 );

  
  typedef itk::AccumulateImageFilter< FilterImageType, FilterImageType > AccumulateFilter;
  AccumulateFilter::Pointer accumulator = AccumulateFilter::New();
  accumulator->SetAccumulateDimension(1);
  accumulator->SetInput( sigmoidcaster->GetOutput() );

  typedef itk::ExtractImageFilter< FilterImageType, OutFilterImageType > ExtractFilter;
  ExtractFilter::Pointer extractor = ExtractFilter::New();
  extractor->SetInput( accumulator->GetOutput() );
  FilterImageType::Pointer accuOut = accumulator->GetOutput();
  accuOut->UpdateOutputInformation();
  FilterImageType::RegionType extractRegion = accuOut->GetLargestPossibleRegion();
  
  extractRegion.SetSize(1,0);
  
  extractor->SetExtractionRegion( extractRegion );

  typedef itk::ResampleImageFilter<OutFilterImageType, OutFilterImageType > ResampleFilter;
  ResampleFilter::Pointer resampler = ResampleFilter::New();
  resampler->SetInput( extractor->GetOutput() );
  
  typedef itk::BSplineInterpolateImageFunction< OutFilterImageType > InterpolatorType;
  InterpolatorType::Pointer interpolator = InterpolatorType::New();
  interpolator->SetSplineOrder(3);
  
  resampler->SetInterpolator( interpolator );
  OutFilterImageType::Pointer exOut = extractor->GetOutput();
  exOut->UpdateOutputInformation();
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
	std::string inputFilenamesFilename = argv[1];
	double keyPointIntensityThreshold = atof(argv[2]);
	double dogSplitsPerOctave = atof(argv[3]);
	double statingScale = atof(argv[4]);
	double eLocation = atof(argv[5]);
	double eScale = std::log(atof(argv[6]));
	double eOrientation = atof(argv[7]);
	double gammaValue = atof(argv[8]);

	std::string distanceMapFilenamesFilename = argv[9];
	double extractionDistanceThreshold = atof(argv[10]);




	// load up the set of aligned images
	FilenamesReader::FilenamesType inputFilenames = FilenamesReader::Read(inputFilenamesFilename);
	FilenamesReader::FilenamesType distanceMapFilenames   = FilenamesReader::Read(distanceMapFilenamesFilename);
	ImageVolumeList images;
	RealVolumeList distanceMaps;
	for(unsigned int i = 0; i < inputFilenames.size(); i++)
	{
		ImageVolume::Pointer image = ImageVolumeIO::Read(inputFilenames[i]);
		images.push_back(image);
		RealVolume::Pointer distMap = RealVolumeIO::Read(distanceMapFilenames[i]);
		distanceMaps.push_back(distMap);
	}


	unsigned int sliceToTest = 7;

	// for each slice we want to learn the features
	const unsigned int sliceNum = images.front()->GetLargestPossibleRegion().GetSize()[2];
	for(unsigned int slice = sliceToTest; slice < sliceNum; slice++)
	{

		// get the set of slices that have some image data in them
		ImageSliceList validImages;
		RealSliceList validDistanceMaps;
		
		for(unsigned int im = 0; im < images.size(); im++)
		{
			ImageSlice::Pointer extractedSlice = ImageSlice::New();
			RealSlice::Pointer distanceSlice = RealSlice::New();
			ExtractSlice<ImageVolume, ImageSlice>(images[im], slice, extractedSlice);
			ExtractSlice<RealVolume, RealSlice>(distanceMaps[im], slice, distanceSlice);

			if(ImageContainsData(extractedSlice))
			{
				validDistanceMaps.push_back(distanceSlice);
				validImages.push_back(extractedSlice);
			}
		}

		/*
		if(validImages.size() < 3)
			continue;
		*/

		std::cout << "Slice Num: " << slice << " Image Number: " << validImages.size() << std::endl;



		typedef itk::Vector<double, 2> VectorType;
		typedef itk::Image<VectorType, 2> GradientType;
		typedef filter::HistogramOfGradeintsFeatureExtractor<GradientType> FeatureBuilderType;
		typedef FeatureBuilderType::FeatureType HoGFeatureType;
		std::vector<HoGFeatureType> allFeatures;
		std::vector<HoGFeatureType> allFeatures1;
		std::vector<HoGFeatureType> allFeatures2;

		

		unsigned int featureCount = 0;
		for(unsigned int im = 0; im < validImages.size(); im++)
		{
			ImageSlice::Pointer extractedSlice = validImages[im];

			// first we extract all of the keypoints points
			typedef filter::DoGKeyPointExtractor<utils::ImageSlice> ExtractorType;
			ExtractorType::Pointer extractor = ExtractorType::New();
			extractor->SetInput(extractedSlice);
			extractor->SetKeypointThreshold(keyPointIntensityThreshold);
			extractor->SetSplitsPerOctave(dogSplitsPerOctave);
			extractor->SetStartingSigma(statingScale);
			extractor->SetDistanceMap(validDistanceMaps[im]);
			extractor->SetDistanceThreshold(extractionDistanceThreshold);
			extractor->Update();

			// orientate the feature points
			typedef filter::KeyPointOrientator<utils::ImageSlice> Orientator;
			Orientator::Pointer orientator  = Orientator::New();
			orientator->SetInput(extractedSlice);
			orientator->SetKeyPoints(extractor->GetOutput());
			orientator->SetHistogramBins(32);
			orientator->SetSigmaScale(2);
			orientator->SetSampleRadius(5);
			orientator->Update();
//.........这里部分代码省略.........

int main( int argc, char *argv[] )
{
if( argc < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " outputImagefile [differenceBeforeRegistration] ";
std::cerr << " [differenceAfterRegistration] ";
std::cerr << " [sliceBeforeRegistration] ";
std::cerr << " [sliceDifferenceBeforeRegistration] ";
std::cerr << " [sliceDifferenceAfterRegistration] ";
std::cerr << " [sliceAfterRegistration] " << std::endl;
return EXIT_FAILURE;
}
const unsigned int Dimension = 3;
typedef float PixelType;
typedef itk::Image< PixelType, Dimension > FixedImageType;
typedef itk::Image< PixelType, Dimension > MovingImageType;
// Software Guide : BeginLatex
//
// The Transform class is instantiated using the code below. The only
// template parameter to this class is the representation type of the
// space coordinates.
//
// \index{itk::Versor\-Rigid3D\-Transform!Instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet

// Software Guide : EndCodeSnippet


typedef itk:: LinearInterpolateImageFunction< MovingImageType, double > InterpolatorType;
typedef itk::ImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType;

MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetInterpolator( interpolator );
// Software Guide : BeginLatex
//
// The transform object is constructed below and passed to the registration
// method.
//
// \index{itk::Versor\-Rigid3D\-Transform!New()}
// \index{itk::Versor\-Rigid3D\-Transform!Pointer}
// \index{itk::Registration\-Method!SetTransform()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
TransformType::Pointer transform = TransformType::New();
registration->SetTransform( transform );
// Software Guide : EndCodeSnippet
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
fixedImageReader->SetFileName( argv[1] );
movingImageReader->SetFileName( argv[2] );
registration->SetFixedImage( fixedImageReader->GetOutput() );
registration->SetMovingImage( movingImageReader->GetOutput() );
fixedImageReader->Update();
registration->SetFixedImageRegion(
fixedImageReader->GetOutput()->GetBufferedRegion() );
// Software Guide : BeginLatex
//
// The input images are taken from readers. It is not necessary here to
// explicitly call \code{Update()} on the readers since the
// \doxygen{CenteredTransformInitializer} will do it as part of its
// computations. The following code instantiates the type of the
// initializer. This class is templated over the fixed and moving image type
// as well as the transform type. An initializer is then constructed by
// calling the \code{New()} method and assigning the result to a smart
// pointer.
//
// \index{itk::Centered\-Transform\-Initializer!Instantiation}
// \index{itk::Centered\-Transform\-Initializer!New()}
// \index{itk::Centered\-Transform\-Initializer!SmartPointer}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
// Software Guide : BeginLatex
//
// Let's execute this example over some of the images available in the ftp
// site
//
// \url{ftp://public.kitware.com/pub/itk/Data/BrainWeb}
//
// Note that the images in the ftp site are compressed in \code{.tgz} files.
// You should download these files an uncompress them in your local system.
// After decompressing and extracting the files you could take a pair of
// volumes, for example the pair:
//
// \begin{itemize}
// \item \code{brainweb1e1a10f20.mha}
// \item \code{brainweb1e1a10f20Rot10Tx15.mha}
//.........这里部分代码省略.........

bool mitk::NavigationDataLandmarkTransformFilter::FindCorrespondentLandmarks(LandmarkPointContainer& sources, const LandmarkPointContainer& targets) const
{
  if (sources.size() < 6 || targets.size() < 6)
    return false;
  //throw std::invalid_argument("ICP correspondence finding needs at least 6 landmarks");

  /* lots of type definitions */
  typedef itk::PointSet<mitk::ScalarType, 3> PointSetType;
  //typedef itk::BoundingBox<PointSetType::PointIdentifier, PointSetType::PointDimension> BoundingBoxType;

  typedef itk::EuclideanDistancePointMetric< PointSetType, PointSetType> MetricType;
  //typedef MetricType::TransformType TransformBaseType;
  //typedef MetricType::TransformType::ParametersType ParametersType;
  //typedef TransformBaseType::JacobianType JacobianType;
  //typedef itk::Euler3DTransform< double > TransformType;
  typedef itk::VersorRigid3DTransform< double > TransformType;
  typedef TransformType ParametersType;
  typedef itk::PointSetToPointSetRegistrationMethod< PointSetType, PointSetType > RegistrationType;

  /* copy landmarks to itk pointsets for registration */
  PointSetType::Pointer sourcePointSet = PointSetType::New();
  unsigned int i = 0;
  for (LandmarkPointContainer::const_iterator it = sources.begin(); it != sources.end(); ++it)
  {
    PointSetType::PointType doublePoint;
    mitk::itk2vtk(*it, doublePoint); // copy mitk::ScalarType point into double point as workaround to ITK 3.10 bug
    sourcePointSet->SetPoint(i++, doublePoint /**it*/);
  }

  i = 0;
  PointSetType::Pointer targetPointSet = PointSetType::New();
  for (LandmarkPointContainer::const_iterator it = targets.begin(); it != targets.end(); ++it)
  {
    PointSetType::PointType doublePoint;
    mitk::itk2vtk(*it, doublePoint); // copy mitk::ScalarType point into double point as workaround to ITK 3.10 bug
    targetPointSet->SetPoint(i++, doublePoint /**it*/);
  }

  /* get centroid and extends of our pointsets */
  //BoundingBoxType::Pointer sourceBoundingBox = BoundingBoxType::New();
  //sourceBoundingBox->SetPoints(sourcePointSet->GetPoints());
  //sourceBoundingBox->ComputeBoundingBox();
  //BoundingBoxType::Pointer targetBoundingBox = BoundingBoxType::New();
  //targetBoundingBox->SetPoints(targetPointSet->GetPoints());
  //targetBoundingBox->ComputeBoundingBox();


  TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();
  //transform->SetTranslation(targetBoundingBox->GetCenter() - sourceBoundingBox->GetCenter());

  itk::LevenbergMarquardtOptimizer::Pointer optimizer = itk::LevenbergMarquardtOptimizer::New();
  optimizer->SetUseCostFunctionGradient(false);

  RegistrationType::Pointer registration = RegistrationType::New();

  // Scale the translation components of the Transform in the Optimizer
  itk::LevenbergMarquardtOptimizer::ScalesType scales(transform->GetNumberOfParameters());
  const double translationScale = 5000; //sqrtf(targetBoundingBox->GetDiagonalLength2())  * 1000; // dynamic range of translations
  const double rotationScale = 1.0; // dynamic range of rotations
  scales[0] = 1.0 / rotationScale;
  scales[1] = 1.0 / rotationScale;
  scales[2] = 1.0 / rotationScale;
  scales[3] = 1.0 / translationScale;
  scales[4] = 1.0 / translationScale;
  scales[5] = 1.0 / translationScale;
  //scales.Fill(0.01);
  unsigned long numberOfIterations = 80000;
  double gradientTolerance = 1e-10; // convergence criterion
  double valueTolerance = 1e-10; // convergence criterion
  double epsilonFunction = 1e-10; // convergence criterion
  optimizer->SetScales( scales );
  optimizer->SetNumberOfIterations( numberOfIterations );
  optimizer->SetValueTolerance( valueTolerance );
  optimizer->SetGradientTolerance( gradientTolerance );
  optimizer->SetEpsilonFunction( epsilonFunction );


  registration->SetInitialTransformParameters( transform->GetParameters() );
  //------------------------------------------------------
  // Connect all the components required for Registration
  //------------------------------------------------------
  MetricType::Pointer metric = MetricType::New();

  registration->SetMetric( metric );
  registration->SetOptimizer( optimizer );
  registration->SetTransform( transform );
  registration->SetFixedPointSet( targetPointSet );
  registration->SetMovingPointSet( sourcePointSet );

  try
  {
    //registration->StartRegistration();
    registration->Update();
  }
  catch( itk::ExceptionObject & e )
  {
    MITK_INFO << "Exception caught during ICP optimization: " << e;
    return false;
    //throw e;
//.........这里部分代码省略.........

int mitkPyramidImageRegistrationMethodTest( int argc, char* argv[] )
{
  if( argc < 4 )
  {
    MITK_ERROR << "Not enough input \n Usage: <TEST_NAME> fixed moving type [output_image [output_transform]]"
               << "\n \t fixed : the path to the fixed image \n"
               << " \t moving : path to the image to be registered"
               << " \t type : Affine or Rigid defining the type of the transformation"
               << " \t output_image : output file optional, (full) path, and optionally output_transform : also (full)path to file";
    return EXIT_FAILURE;
  }

  MITK_TEST_BEGIN("PyramidImageRegistrationMethodTest");

  mitk::Image::Pointer fixedImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load( argv[1] )[0].GetPointer());
  mitk::Image::Pointer movingImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load( argv[2] )[0].GetPointer());

  std::string type_flag( argv[3] );

  mitk::PyramidImageRegistrationMethod::Pointer registrationMethod = mitk::PyramidImageRegistrationMethod::New();
  registrationMethod->SetFixedImage( fixedImage );
  registrationMethod->SetMovingImage( movingImage );

  if( type_flag == "Rigid" )
  {
    registrationMethod->SetTransformToRigid();
  }
  else if( type_flag == "Affine" )
  {
    registrationMethod->SetTransformToAffine();
  }
  else
  {
    MITK_WARN << " No type specified, using 'Affine' .";
  }

  registrationMethod->Update();

  bool imageOutput = false;
  bool transformOutput = false;

  std::string image_out_filename, transform_out_filename;

  std::string first_output( argv[4] );
  // check for txt, otherwise suppose it is an image
  if( first_output.find(".txt") != std::string::npos )
  {
    transformOutput = true;
    transform_out_filename = first_output;
  }
  else
  {
    imageOutput = true;
    image_out_filename = first_output;
  }

  if( argc > 4 )
  {
    std::string second_output( argv[5] );
    if( second_output.find(".txt") != std::string::npos )
    {
      transformOutput = true;
      transform_out_filename = second_output;
    }
  }

  MITK_INFO << " Selected output: " << transform_out_filename  << " " << image_out_filename;

  try{

    unsigned int paramCount = registrationMethod->GetNumberOfParameters();
    double* params = new double[ paramCount ];
    registrationMethod->GetParameters( &params[0] );

    std::cout << "Parameters: ";
    for( unsigned int i=0; i< paramCount; i++)
    {
      std::cout << params[ i ] << " ";
    }
    std::cout << std::endl;

    if( imageOutput )
    {
      mitk::IOUtil::Save( registrationMethod->GetResampledMovingImage(), image_out_filename.c_str() );
    }


    if( transformOutput )
    {

      itk::TransformFileWriter::Pointer writer = itk::TransformFileWriter::New();

      // Get transform parameter for resampling / saving
      // Affine
      if( paramCount == 12 )
      {
        typedef itk::AffineTransform< double > TransformType;
        TransformType::Pointer transform = TransformType::New();

        TransformType::ParametersType affine_params( paramCount );
//.........这里部分代码省略.........

int main (int argc, char **argv)
{
  int verbose=0, clobber=0,skip_grid=0;
  double max=5.0;
  double extent=300;

  static struct option long_options[] = {
		{"verbose", no_argument,       &verbose, 1},
		{"quiet",   no_argument,       &verbose, 0},
		{"clobber", no_argument,       &clobber, 1},
		{"spacing", required_argument, 0, 's'},
		{"max",     required_argument, 0, 'm'},
		{"version", no_argument,       0, 'v'},
    {"extent", required_argument,  0, 'e'},
		{0, 0, 0, 0}
		};
  
  double spacing=4.0;
  for (;;) {
      /* getopt_long stores the option index here. */
      int option_index = 0;

      int c = getopt_long (argc, argv, "s:m:v", long_options, &option_index);

      /* Detect the end of the options. */
      if (c == -1) break;

      switch (c)
			{
			case 0:
				break;
			case 's':
				spacing=atof(optarg);
				break;
			case 'v':
				cout << "Version: 1.0" << endl;
				return 0;
      case 'm':
        max=atof(optarg); break;
      case 'e':
        extent=atof(optarg); break;
			case '?':
				/* getopt_long already printed an error message. */
			default:
				show_usage ();
				return 1;
			}
    }

	if ((argc - optind) < 2) {
		show_usage ();
		return 1;
	}
  std::string input=argv[optind];
  std::string output=argv[optind+1];
	try
  {
    gsl_rng_env_setup();
    
		typedef minc::SphericalHarmonicsTransform TransformType;
    
		TransformType::ParametersType finalParameters;
		load_parameters(input.c_str(),finalParameters);
		TransformType::Pointer finalTransform = TransformType::New();
    cout<<"Loaded parameters:"<<finalParameters<<endl;
		finalTransform->ImportParameters( finalParameters , true);
    
    cout<<"Imported!"<<endl;
		minc::def3d::Pointer grid(minc::def3d::New());
    
    allocate_image3d<minc::def3d>(grid, 
                      fixed_vec<3, unsigned int>(extent/spacing), 
                      fixed_vec<3, double>(spacing), 
                      fixed_vec<3, double>(-extent/2));
		
		if(verbose) 
		{
			cout<<"Generating a grid file, ";
			cout<<"extent: "<<extent<<" spacing: "<<spacing<<" ..."<<std::flush;
		}
		
    def3d_iterator it(grid,grid->GetLargestPossibleRegion());
		for(it.GoToBegin();!it.IsAtEnd();++it) 
    {
      tag_point p,p2;
      grid->TransformIndexToPhysicalPoint(it.GetIndex(),p);
			p2=finalTransform->TransformPointUnCached(p);
      def_vector moved;
			moved[0]=p2[0]-p[0];
			moved[1]=p2[1]-p[1];
			moved[2]=p2[2]-p[2];
      if(fabs(moved[0])>max || fabs(moved[1])>max ||fabs(moved[2])>max)
        moved[0]=moved[1]=moved[2]=0.0;
			
      it.Value()=moved;
    }
		if(verbose)
			cout<<"Done!"<<endl;
    save_minc<def3d>(output.c_str(), grid);
		
//.........这里部分代码省略.........