+ 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() );
+
+ // Find the bounding box of the "inside" label
+ typedef itk::LabelStatisticsImageFilter<ImageMaskType, ImageMaskType> StatisticsImageFilterType;
+ typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
+ statisticsImageFilter->SetInput(maskReader->GetOutput());
+ statisticsImageFilter->SetLabelInput(maskReader->GetOutput());
+ statisticsImageFilter->Update();
+ typename StatisticsImageFilterType::BoundingBoxType boundingBox = statisticsImageFilter->GetBoundingBox(1);
+
+ // Limit the transform region to the mask
+ for (unsigned int i=0; i<InputImageType::ImageDimension; i++)
+ {
+ transformRegionIndex[i]=boundingBox[2*i];
+ transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
+ }
+ transformRegion.SetSize(transformRegionSize);
+ transformRegion.SetIndex(transformRegionIndex);
+ fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
+
+ // Crop the fixedImage to the bounding box to facilitate multi-resolution
+ typedef itk::ExtractImageFilter<FixedImageType,FixedImageType> ExtractImageFilterType;
+ typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New();
+ extractImageFilter->SetInput(fixedImage);
+ extractImageFilter->SetExtractionRegion(transformRegion);
+ extractImageFilter->Update();
+ croppedFixedImage=extractImageFilter->GetOutput();
+
+ // Update the metric region
+ metricRegion = croppedFixedImage->GetLargestPossibleRegion();
+ metricRegionIndex=metricRegion.GetIndex();
+ croppedFixedImage->TransformIndexToPhysicalPoint(metricRegionIndex, metricRegionOrigin);
+
+ // Set start index to zero (with respect to croppedFixedImage/transform region)
+ metricRegionIndex.Fill(0);
+ metricRegion.SetIndex(metricRegionIndex);
+ croppedFixedImage->SetRegions(metricRegion);
+ croppedFixedImage->SetOrigin(metricRegionOrigin);
+ }
+
+ 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() );
+ }
+
+
+ //=======================================================
+ // Output Regions
+ //=======================================================
+
+ if (m_Verbose)
+ {
+ // Fixed image region
+ std::cout<<"The fixed image has its origin at "<<fixedImage->GetOrigin()<<std::endl
+ <<"The fixed image region starts at index "<<fixedImageRegion.GetIndex()<<std::endl
+ <<"The fixed image region has size "<< fixedImageRegion.GetSize()<<std::endl;
+
+ // Transform region
+ std::cout<<"The transform has its origin at "<<transformRegionOrigin<<std::endl
+ <<"The transform region will start at index "<<transformRegion.GetIndex()<<std::endl
+ <<"The transform region has size "<< transformRegion.GetSize()<<std::endl;
+
+ // Metric region
+ std::cout<<"The metric region has its origin at "<<metricRegionOrigin<<std::endl
+ <<"The metric region will start at index "<<metricRegion.GetIndex()<<std::endl
+ <<"The metric region has size "<< metricRegion.GetSize()<<std::endl;
+
+ }
+
+
+ //=======================================================
+ // Pyramids (update them for transform initializer)
+ //=======================================================
+ typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, FixedImageType> FixedImagePyramidType;
+ typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType;
+ typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
+ typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
+ fixedImagePyramid->SetUseShrinkImageFilter(false);
+ fixedImagePyramid->SetInput(croppedFixedImage);
+ 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();
+ typedef clitk::MultiResolutionPyramidRegionFilter<FixedImageType> FixedImageRegionPyramidType;
+ typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
+ fixedImageRegionPyramid->SetRegion(metricRegion);
+ fixedImageRegionPyramid->SetSchedule(fixedImagePyramid->GetSchedule());
+
+
+ //=======================================================
+ // Rigid or Affine Transform
+ //=======================================================
+ typedef itk::AffineTransform <double,3> RigidTransformType;
+ RigidTransformType::Pointer rigidTransform=NULL;
+ if (m_ArgsInfo.initMatrix_given)
+ {
+ if(m_Verbose) std::cout<<"Reading the prior transform matrix "<< m_ArgsInfo.initMatrix_arg<<"..."<<std::endl;
+ rigidTransform=RigidTransformType::New();
+ itk::Matrix<double,4,4> rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.initMatrix_arg);
+
+ //Set the rotation
+ itk::Matrix<double,3,3> finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix);
+ rigidTransform->SetMatrix(finalRotation);
+
+ //Set the translation
+ itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
+ rigidTransform->SetTranslation(finalTranslation);
+ }
+
+
+ //=======================================================
+ // B-LUT FFD Transform
+ //=======================================================
+ typedef clitk::BSplineDeformableTransform<TCoordRep,InputImageType::ImageDimension, SpaceDimension > TransformType;
+ typename TransformType::Pointer transform= TransformType::New();
+ if (fixedMask) transform->SetMask( fixedMask );
+ if (rigidTransform) transform->SetBulkTransform( rigidTransform );
+
+ //-------------------------------------------------------------------------
+ // The transform initializer
+ //-------------------------------------------------------------------------
+ typedef clitk::BSplineDeformableTransformInitializer< TransformType,FixedImageType> InitializerType;
+ typename InitializerType::Pointer initializer = InitializerType::New();
+ initializer->SetVerbose(m_Verbose);
+ initializer->SetImage(fixedImagePyramid->GetOutput(0));
+ initializer->SetTransform(transform);
+
+ //-------------------------------------------------------------------------
+ // Order
+ //-------------------------------------------------------------------------
+ typename FixedImageType::RegionType::SizeType splineOrders ;
+ splineOrders.Fill(3);
+ if (m_ArgsInfo.order_given)
+ for(unsigned int i=0; i<InputImageType::ImageDimension;i++)
+ splineOrders[i]=m_ArgsInfo.order_arg[i];
+ if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
+ initializer->SetSplineOrders(splineOrders);
+
+ //-------------------------------------------------------------------------
+ // Levels
+ //-------------------------------------------------------------------------
+
+ // Spacing
+ if (m_ArgsInfo.spacing_given)
+ {
+ initializer->m_ControlPointSpacingArray.resize(m_ArgsInfo.levels_arg);
+ initializer->SetControlPointSpacing(m_ArgsInfo.spacing_arg);
+ initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]=initializer->m_ControlPointSpacing;
+ if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]
+ <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
+
+ for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
+ {
+ initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]=initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-i]*2;
+ if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]
+ <<" at level "<<m_ArgsInfo.levels_arg-i<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
+ }
+
+ }
+
+ // Control
+ if (m_ArgsInfo.control_given)
+ {
+ initializer->m_NumberOfControlPointsInsideTheImageArray.resize(m_ArgsInfo.levels_arg);
+ initializer->SetNumberOfControlPointsInsideTheImage(m_ArgsInfo.control_arg);
+ initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]=initializer->m_NumberOfControlPointsInsideTheImage;
+ if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]<<"control points inside the image"
+ <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
+
+ for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
+ {
+ for(unsigned int j=0;j<InputImageType::ImageDimension;j++)
+ initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i][j]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i][j]/2.);
+ // initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i]/2.);
+ if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]<<"control points inside the image"
+ <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
+
+ }
+ }
+
+ // Inialize on the first level
+ if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level 1 of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
+ if (m_ArgsInfo.spacing_given) initializer->SetControlPointSpacing( initializer->m_ControlPointSpacingArray[0]);
+ if (m_ArgsInfo.control_given) initializer->SetNumberOfControlPointsInsideTheImage(initializer->m_NumberOfControlPointsInsideTheImageArray[0]);
+ if (m_ArgsInfo.samplingFactor_given) initializer->SetSamplingFactors(m_ArgsInfo.samplingFactor_arg);
+
+ // Initialize
+ initializer->InitializeTransform();
+
+ //-------------------------------------------------------------------------
+ // Initial parameters (passed by reference)
+ //-------------------------------------------------------------------------
+ typedef typename TransformType::ParametersType ParametersType;
+ const unsigned int numberOfParameters = transform->GetNumberOfParameters();
+ ParametersType parameters(numberOfParameters);
+ parameters.Fill( 0.0 );
+ transform->SetParameters( parameters );
+ if (m_ArgsInfo.initCoeff_given) initializer->SetInitialParameters(m_ArgsInfo.initCoeff_arg, parameters);
+
+
+ //=======================================================
+ // Interpolator
+ //=======================================================
+ typedef clitk::GenericInterpolator<args_info_clitkBLUTDIR, FixedImageType, TCoordRep > GenericInterpolatorType;
+ typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
+ genericInterpolator->SetArgsInfo(m_ArgsInfo);
+ typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep > InterpolatorType;
+ typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
+
+
+ //=======================================================
+ // Metric
+ //=======================================================
+ typedef clitk::GenericMetric<args_info_clitkBLUTDIR, FixedImageType,MovingImageType > GenericMetricType;
+ typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
+ genericMetric->SetArgsInfo(m_ArgsInfo);
+ genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
+ if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
+ genericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(0));
+ typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType;
+ typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
+ if (movingMask) metric->SetMovingImageMask(movingMask);
+
+#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
+ if (threadsGiven) {
+ metric->SetNumberOfThreads( threads );
+ if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl;
+ }
+#else
+ if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
+#endif
+
+
+ //=======================================================
+ // Optimizer
+ //=======================================================
+ typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> GenericOptimizerType;
+ GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New();
+ genericOptimizer->SetArgsInfo(m_ArgsInfo);
+ genericOptimizer->SetMaximize(genericMetric->GetMaximize());
+ genericOptimizer->SetNumberOfParameters(transform->GetNumberOfParameters());
+ typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
+ OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer();
+
+
+ //=======================================================
+ // Registration
+ //=======================================================
+ typedef itk::MultiResolutionImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType;
+ typename RegistrationType::Pointer registration = RegistrationType::New();
+ registration->SetMetric( metric );
+ registration->SetOptimizer( optimizer );
+ registration->SetInterpolator( interpolator );
+ registration->SetTransform (transform);
+ if(threadsGiven) {
+ registration->SetNumberOfThreads(threads);
+ if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl;
+ }
+ registration->SetFixedImage( croppedFixedImage );
+ registration->SetMovingImage( movingImage );
+ registration->SetFixedImageRegion( metricRegion );
+ registration->SetFixedImagePyramid( fixedImagePyramid );
+ registration->SetMovingImagePyramid( movingImagePyramid );
+ registration->SetInitialTransformParameters( transform->GetParameters() );
+ registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
+ if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
+
+
+ //================================================================================================
+ // Observers
+ //================================================================================================
+ 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,args_info_clitkBLUTDIR> CommandType;
+ typename CommandType::Pointer command = CommandType::New();
+ command->SetInitializer(initializer);
+ command->SetArgsInfo(m_ArgsInfo);
+ command->SetCommandIterationUpdate(observer);
+ command->SetMaximize(genericMetric->GetMaximize());
+ command->SetMetricRegion(metricRegion);
+ registration->AddObserver( itk::IterationEvent(), command );
+ }
+
+
+ //=======================================================
+ // Let's go
+ //=======================================================
+ if (m_Verbose) std::cout << std::endl << "Starting Registration" << std::endl;
+
+ try
+ {
+ registration->StartRegistration();
+ }
+ catch( itk::ExceptionObject & err )
+ {
+ std::cerr << "ExceptionObject caught while registering!" << std::endl;
+ std::cerr << err << std::endl;
+ return;
+ }
+
+
+ //=======================================================
+ // Get the result
+ //=======================================================
+ OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
+ transform->SetParameters( finalParameters );
+ if (m_Verbose)
+ {
+ std::cout<<"Stop condition description: "
+ <<registration->GetOptimizer()->GetStopConditionDescription()<<std::endl;
+ }
+
+
+ //=======================================================
+ // Get the BSpline coefficient images and write them
+ //=======================================================
+ if (m_ArgsInfo.coeff_given)
+ {
+ typedef typename TransformType::CoefficientImageType CoefficientImageType;
+ typename CoefficientImageType::Pointer coefficientImage =transform->GetCoefficientImage();
+ typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
+ typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
+ coeffWriter->SetInput(coefficientImage);
+ coeffWriter->SetFileName(m_ArgsInfo.coeff_arg);
+ coeffWriter->Update();
+ }
+
+
+
+ //=======================================================
+ // Compute the DVF (only deformable transform)
+ //=======================================================
+ typedef itk::Vector< float, SpaceDimension > DisplacementType;
+ typedef itk::Image< DisplacementType, InputImageType::ImageDimension > DeformationFieldType;
+ typedef itk::TransformToDeformationFieldSource<DeformationFieldType, double> ConvertorType;
+ typename ConvertorType::Pointer filter= ConvertorType::New();
+ filter->SetNumberOfThreads(1);
+ transform->SetBulkTransform(NULL);
+ filter->SetTransform(transform);
+ filter->SetOutputParametersFromImage(fixedImage);
+ filter->Update();
+ typename DeformationFieldType::Pointer field = filter->GetOutput();
+
+
+ //=======================================================
+ // Write the DVF
+ //=======================================================
+ typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;
+ typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
+ fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
+ fieldWriter->SetInput( field );
+ try
+ {
+ fieldWriter->Update();
+ }
+ catch( itk::ExceptionObject & excp )
+ {
+ std::cerr << "Exception thrown writing the DVF" << std::endl;
+ std::cerr << excp << std::endl;