+ m_ArgsInfo=a;
+ if (m_ArgsInfo.reference_given) AddInputFilename(m_ArgsInfo.reference_arg);
+
+ if (m_ArgsInfo.target_given) {
+ AddInputFilename(m_ArgsInfo.target_arg);
+ }
+
+ if (m_ArgsInfo.output_given) SetOutputFilename(m_ArgsInfo.output_arg);
+}
+//==============================================================================
+// Update with the number of dimensions and pixeltype
+//==============================================================================
+template<class InputImageType>
+void AffineRegistrationGenericFilter::UpdateWithInputImageType()
+{
+ //=============================================================================
+ //Input
+ //=============================================================================
+
+ typedef typename InputImageType::PixelType PixelType;
+//typedef typename InputImageType::ImageDimension Dimension;
+
+
+#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
+ bool threadsGiven=m_ArgsInfo.threads_given;
+ int threads=m_ArgsInfo.threads_arg;
+#endif
+
+ //Coordinate Representation
+ typedef double TCoordRep;
+
+
+ typename InputImageType::Pointer fixedImage = this->template GetInput<InputImageType>(0);
+
+ typename InputImageType::Pointer inputFixedImage = this->template GetInput<InputImageType>(0);
+
+ // typedef input2
+ typename InputImageType::Pointer movingImage = this->template GetInput<InputImageType>(1);
+
+ typename InputImageType::Pointer inputMovingImage = this->template GetInput<InputImageType>(1);
+
+
+
+ //The pixeltype of the fixed image will be used for output
+ typedef itk::Image< PixelType, InputImageType::ImageDimension > FixedImageType;
+
+ //Whatever the pixel type, internally we work with an image represented in float
+ typedef typename InputImageType::PixelType InternalPixelType;
+ typedef itk::Image< PixelType, InputImageType::ImageDimension > InternalImageType;
+
+
+ //Read in the reference/fixed image
+// typedef itk::ImageFileReader< InternalImageType > ReaderType;
+// typename ReaderType::Pointer fixedImageReader = ReaderType::New();
+// fixedImageReader->SetFileName( m_ArgsInfo.reference_arg);
+
+
+ //Read in the object/moving image
+// typename ReaderType::Pointer movingImageReader = ReaderType::New();
+// movingImageReader->SetFileName( m_ArgsInfo.target_arg );
+ if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
+// fixedImageReader->Update();
+// movingImageReader->Update();
+
+ if (m_Verbose) std::cout << "Reading images... " << std::endl;
+
+ //we connect pointers to these internal images
+ // typedef typename fixedImageReader fixedImage;
+ // typedef typename movingImageReader movingImage;
+
+ //We keep the images used for input for possible output
+// typedef typename fixedImageReader inputFixedImage;
+// typedef typename movingImageReader inputMovingImage;
+
+
+ //============================================================================
+ // Preprocessing
+ //============================================================================
+
+ //If given, the intensities of both images are first normalized to a zero mean and SD of 1
+ // (usefull for MI, not necessary for Mattes' MI but performed anyway for the ouput)
+ if ( m_ArgsInfo.normalize_flag ) {
+ typedef itk::NormalizeImageFilter< InternalImageType,InternalImageType > NormalizeFilterType;
+
+ typename NormalizeFilterType::Pointer fixedNormalizeFilter = NormalizeFilterType::New();
+ typename NormalizeFilterType::Pointer movingNormalizeFilter = NormalizeFilterType::New();
+
+ fixedNormalizeFilter->SetInput( fixedImage );
+ movingNormalizeFilter->SetInput( movingImage );
+
+ fixedNormalizeFilter->Update();
+ movingNormalizeFilter->Update();
+
+ //We keep the images used for input for possible output
+ inputFixedImage= fixedNormalizeFilter->GetOutput();
+ inputMovingImage= movingNormalizeFilter->GetOutput();
+
+ //the pointers are reconnected for further output
+ fixedImage=fixedNormalizeFilter->GetOutput();
+ movingImage=movingNormalizeFilter->GetOutput();
+
+ if (m_Verbose) std::cout << "Normalizing image intensities to zero mean and SD of 1..." << std::endl;
+ }
+
+
+ //If given, the images are blurred before processing
+ if ( m_ArgsInfo.blur_arg!= 0) {
+ typedef itk::DiscreteGaussianImageFilter<InternalImageType,InternalImageType> GaussianFilterType;
+ typename GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
+ typename GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
+ fixedSmoother->SetVariance( m_ArgsInfo.blur_arg );
+ movingSmoother->SetVariance(m_ArgsInfo.blur_arg );
+
+ fixedSmoother->SetInput( fixedImage );
+ movingSmoother->SetInput( movingImage );
+
+ fixedSmoother->Update();
+ movingSmoother->Update();
+
+ fixedImage=fixedSmoother->GetOutput();
+ movingImage=movingSmoother->GetOutput();
+
+ if (m_Verbose) std::cout << "Blurring images with a Gaussian with standard deviation of " << m_ArgsInfo.blur_arg <<"..." << std::endl;
+ }
+
+
+ //============================================================================
+ // Setting up the moving image in a reference system
+ //============================================================================
+ const itk::Vector<double, InputImageType::ImageDimension> movingResolution = movingImage->GetSpacing();
+ typename InternalImageType::RegionType movingRegion = movingImage->GetLargestPossibleRegion();
+ typename InternalImageType::RegionType::SizeType movingSize = movingRegion.GetSize();
+
+ // Print the parameters of the moving image
+ if (m_Verbose) {
+ std::cout << "Object or Moving image:"<<std::endl;
+ std::cout << "Size: " << movingSize[0] << ", " << movingSize[1];
+ if (InputImageType::ImageDimension==3) std::cout<<", " << movingSize[2];
+ std::cout << std::endl;
+
+ std::cout<< "Resolution: "<< movingResolution[0] << ", " << movingResolution[1];
+ if (InputImageType::ImageDimension==3) std::cout<< ", " << movingResolution[2];
+ std::cout << std::endl;
+ }
+
+
+ //============================================================================
+ // Setting up the fixed image in a reference system
+ //============================================================================
+ const itk::Vector<double, InputImageType::ImageDimension> fixedResolution = fixedImage->GetSpacing();
+ typename InternalImageType::RegionType fixedRegion = fixedImage->GetLargestPossibleRegion();
+ typename InternalImageType::RegionType::SizeType fixedSize = fixedRegion.GetSize();
+
+ // Print the parameters of the moving image and the transform
+ if (m_Verbose) {
+ std::cout << "Target or Moving image:"<<std::endl;
+ std::cout << "Size: " << fixedSize[0] << ", " << fixedSize[1];
+ if (InputImageType::ImageDimension==3) std::cout<<", " << fixedSize[2];
+ std::cout << std::endl;
+
+ std::cout<< "Resolution: "<< fixedResolution[0] << ", " << fixedResolution[1];
+ if (InputImageType::ImageDimension==3) std::cout<< ", " << fixedResolution[2];
+ std::cout << std::endl;
+ }
+
+
+
+ //===========================================================================
+ // If given, we connect a mask to reference or target
+ //============================================================================
+ typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
+ typename MaskType::Pointer fixedMask=NULL;
+ if (m_ArgsInfo.referenceMask_given) {
+ fixedMask= MaskType::New();
+ typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageMaskType;
+ typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
+ typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
+ maskReader->SetFileName(m_ArgsInfo.referenceMask_arg);
+ try {
+ maskReader->Update();
+ } catch ( itk::ExceptionObject & err ) {
+ std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
+ std::cerr << err << std::endl;
+ return;
+ }
+ if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
+
+ // Set the image to the spatialObject
+ fixedMask->SetImage( maskReader->GetOutput() );
+ }
+
+ typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
+ typename MaskType::Pointer movingMask=NULL;
+ if (m_ArgsInfo.targetMask_given) {
+ movingMask= MaskType::New();
+ typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageMaskType;
+ typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
+ typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
+ maskReader->SetFileName(m_ArgsInfo.targetMask_arg);
+ try {
+ maskReader->Update();
+ } catch ( itk::ExceptionObject & err ) {
+ std::cerr << "ExceptionObject caught !" << std::endl;
+ std::cerr << err << std::endl;
+ }
+ if (m_Verbose)std::cout <<"Target image mask was read..." <<std::endl;
+
+ movingMask->SetImage( maskReader->GetOutput() );
+ }
+
+
+ //============================================================================
+ // The image pyramids
+ //============================================================================
+ typedef itk::RecursiveMultiResolutionPyramidImageFilter<InternalImageType,InternalImageType > FixedImagePyramidType;
+ typedef itk::RecursiveMultiResolutionPyramidImageFilter<InternalImageType,InternalImageType > MovingImagePyramidType;
+ typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
+ typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
+ fixedImagePyramid->SetUseShrinkImageFilter(false);
+ fixedImagePyramid->SetInput(fixedImage);
+ fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
+ movingImagePyramid->SetUseShrinkImageFilter(false);
+ movingImagePyramid->SetInput(movingImage);
+ movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
+ if (m_Verbose) std::cout<<"Creating the image pyramid..."<<std::endl;
+
+ fixedImagePyramid->Update();
+ movingImagePyramid->Update();
+
+
+
+ //============================================================================
+ // We retrieve the type of metric from the command line
+ //============================================================================
+ typedef clitk::GenericMetric<args_info_clitkAffineRegistration, InternalImageType, InternalImageType> GenericMetricType;
+ typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
+ genericMetric->SetArgsInfo(m_ArgsInfo);
+ genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
+ if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
+ typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
+ typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
+ if (movingMask) metric->SetMovingImageMask(movingMask);
+
+#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
+ if (threadsGiven) metric->SetNumberOfThreads( threads );
+#else
+ if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
+#endif
+
+ //============================================================================
+ // Initialize using image moments.
+ //============================================================================
+ if (m_ArgsInfo.moment_flag) {
+ typedef itk::ImageMomentsCalculator< InternalImageType > CalculatorType;
+ typename CalculatorType::Pointer fixedCalculator= CalculatorType::New();
+
+ typename InternalImageType::Pointer fixedThresh;
+ if (m_ArgsInfo.intThreshold_given) {
+ typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
+ typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New();
+ thresholder->SetInput(fixedImage);
+ thresholder->SetLower(m_ArgsInfo.intThreshold_arg);
+ thresholder->Update();
+ fixedThresh=thresholder->GetOutput();
+ } else fixedThresh=fixedImage;
+
+ fixedCalculator->SetImage(fixedThresh);
+ fixedCalculator->Compute();
+ Vector<double, InputImageType::ImageDimension> fixedCenter=fixedCalculator->GetCenterOfGravity();
+ if (m_Verbose)std::cout<<"The fixed center of gravity is "<<fixedCenter<<"..."<<std::endl;
+
+ typedef itk::ImageMomentsCalculator< InternalImageType > CalculatorType;
+ typename CalculatorType::Pointer movingCalculator= CalculatorType::New();
+
+ typename InternalImageType::Pointer movingThresh;
+ if (m_ArgsInfo.intThreshold_given) {
+ typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
+ typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New();
+ thresholder->SetInput(movingImage);
+ thresholder->SetLower(m_ArgsInfo.intThreshold_arg);
+ thresholder->Update();
+ movingThresh=thresholder->GetOutput();
+ } else movingThresh=movingImage;
+
+ movingCalculator->SetImage(movingThresh);
+ movingCalculator->Compute();
+ Vector<double, InputImageType::ImageDimension> movingCenter=movingCalculator->GetCenterOfGravity();
+ if (m_Verbose)std::cout<<"The moving center of gravity is "<<movingCenter<<"..."<<std::endl;
+
+ Vector<double, InputImageType::ImageDimension> shift= movingCenter-fixedCenter;
+ if (m_Verbose)std::cout<<"The initial shift applied is "<<shift<<"..."<<std::endl;
+
+ m_ArgsInfo.transX_arg= shift [0];
+ m_ArgsInfo.transY_arg= shift [1];
+ if (InputImageType::ImageDimension==3) m_ArgsInfo.transZ_arg=shift [2];
+ }
+
+ //============================================================================
+ // Transform
+ //============================================================================
+ typedef clitk::GenericAffineTransform<args_info_clitkAffineRegistration, TCoordRep, InputImageType::ImageDimension > GenericAffineTransformType;
+ typename GenericAffineTransformType::Pointer genericAffineTransform = GenericAffineTransformType::New();
+ genericAffineTransform->SetArgsInfo(m_ArgsInfo);
+ typedef itk::Transform< double, InputImageType::ImageDimension, InputImageType::ImageDimension > TransformType;
+ typename TransformType::Pointer transform = genericAffineTransform->GetTransform();
+ std::cout<<m_ArgsInfo.transform_arg<<std::endl;
+
+
+ //=======================================================
+ // Interpolator
+ //=======================================================
+ std::cout<<"setting Interpolator..."<<std::endl;
+ typedef clitk::GenericInterpolator<args_info_clitkAffineRegistration, InternalImageType,TCoordRep > GenericInterpolatorType;
+ typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
+ genericInterpolator->SetArgsInfo(m_ArgsInfo);
+ typedef itk::InterpolateImageFunction< InternalImageType, TCoordRep > InterpolatorType;
+ typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
+ std::cout<<"end of interpolator"<<std::endl;
+
+ //============================================================================
+ // Optimizer
+ //============================================================================
+ typedef clitk::GenericOptimizer<args_info_clitkAffineRegistration> GenericOptimizerType;
+ unsigned int nParam = transform->GetNumberOfParameters();
+ typename GenericOptimizerType::Pointer genericOptimizer=GenericOptimizerType::New();
+ genericOptimizer->SetArgsInfo(m_ArgsInfo);
+ genericOptimizer->SetOutputIteration(m_Verbose);
+ genericOptimizer->SetOutputPosition(m_Verbose);
+ genericOptimizer->SetOutputValue(m_Verbose);
+ genericOptimizer->SetOutputGradient(m_ArgsInfo.gradient_flag);
+ genericOptimizer->SetMaximize(genericMetric->GetMaximize());
+ genericOptimizer->SetNumberOfParameters(nParam);
+ typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
+ OptimizerType::Pointer optimizer=genericOptimizer->GetOptimizerPointer();
+
+ // Scales
+ itk::Optimizer::ScalesType scales( nParam );
+ for (unsigned int i=nParam-InputImageType::ImageDimension; i<nParam; i++) //Translations
+ scales[i] = m_ArgsInfo.tWeight_arg;
+ for (unsigned int i=0; i<nParam-InputImageType::ImageDimension; i++) //Rest
+ scales[i] = m_ArgsInfo.rWeight_arg*180./M_PI;
+ optimizer->SetScales(scales);
+ //============================================================================
+ // Multiresolution registration
+ //============================================================================
+ std::cout<<"start MultiResolution..."<<std::endl;
+ typedef itk::MultiResolutionImageRegistrationMethod< InternalImageType,InternalImageType > RegistrationType;
+ typename RegistrationType::Pointer registration = RegistrationType::New();
+ registration->SetFixedImage( fixedImage );
+ registration->SetFixedImageRegion(fixedImage->GetLargestPossibleRegion());
+ registration->SetMovingImage( movingImage );
+ registration->SetFixedImagePyramid( fixedImagePyramid );
+ registration->SetMovingImagePyramid( movingImagePyramid );
+ registration->SetTransform( transform );
+ registration->SetInitialTransformParameters( transform->GetParameters() );
+ registration->SetInterpolator( interpolator );
+ registration->SetMetric(metric);
+ registration->SetOptimizer(optimizer);
+ registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
+ if (m_Verbose) std::cout << "Setting "<< m_ArgsInfo.levels_arg <<" resolution levels..." << std::endl;
+ if (m_Verbose) std::cout << "Initial Transform: "<< registration->GetInitialTransformParameters()<<std::endl;
+
+ //============================================================================
+ // Connecting the commander to the registration to monitor it
+ //============================================================================
+ if (m_Verbose) {
+
+ // Output iteration info
+ CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
+ observer->SetOptimizer(genericOptimizer);
+ optimizer->AddObserver( itk::IterationEvent(), observer );
+
+
+ // Output level info
+ typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
+ typename CommandType::Pointer command = CommandType::New();
+ command->SetArgsInfo(m_ArgsInfo);
+ registration->AddObserver( itk::IterationEvent(), command );
+
+ }
+
+
+ //============================================================================
+ // Finally we can start the registration with the given amount of multiresolution levels
+ //============================================================================
+ if (m_Verbose) std::cout << "Starting the registration now..." << std::endl;
+
+ try {
+ registration->StartRegistration();
+ } catch ( itk::ExceptionObject & err ) {
+ std::cerr << "ExceptionObject caught !" << std::endl;
+ std::cerr << err << std::endl;