<div dir="ltr"><div>Hi Brain,</div><div><br></div><div>Thanks for your help,</div><div><br></div><div>Below is the code I am using. The images have size [217, 45, 336] and Dimensions [1.48 5.5 1.48 (mm)]. I am using a 3 step multi-resolution framework. Initialized with an affine followed by BsplineSYN. The metric is mattes mutual information, with 40 bins. No pre-processing of the images is performed, I do intend to perform some pre filtering. The metric value appears to take a big jump in the first stage of Bspline at around 100 iterations, using the code below and the result is a mess.</div>
<div><br></div><div style>If you can see anything obvious from my code that would be great.</div><div style><br></div><div style>Thanks for your time I really appreciate it, Kind regards</div><div style><br></div><div style>
Emma</div><div><br></div><div><br></div>
<div><div>#include "itkImageFileReader.h"</div><div>#include "itkImageFileWriter.h"</div><div>#include "itkImageRegistrationMethodv4.h"</div><div>#include "itkAffineTransform.h"</div>
<div>#include "itkMattesMutualInformationImageToImageMetricv4.h"</div><div>#include "itkBSplineSyNImageRegistrationMethod.h"</div><div>#include "itkShrinkImageFilter.h"</div><div>#include "itkBSplineSmoothingOnUpdateDisplacementFieldTransform.h"</div>
<div>#include "itkBSplineSmoothingOnUpdateDisplacementFieldTransformParametersAdaptor.h"</div><div>#include "itkCompositeTransform.h"</div><div>#include "itkVector.h"</div><div>#include <itkRegion.h></div>
<div><br></div><div>template<class TFilter></div><div>class CommandIterationUpdate : public itk::Command</div><div>{</div><div>public:</div><div> typedef CommandIterationUpdate Self;</div><div> typedef itk::Command Superclass;</div>
<div> typedef itk::SmartPointer<Self> Pointer;</div><div> itkNewMacro( Self );</div><div>protected:</div><div> CommandIterationUpdate() {};</div><div><br></div><div>public:</div><div><br></div><div><br></div><div>
void Execute(itk::Object *caller, const itk::EventObject & event)</div><div> {</div><div> Execute( (const itk::Object *) caller, event);</div><div> }</div><div><br></div><div><br></div><div> void Execute(const itk::Object * object, const itk::EventObject & event )</div>
<div> {</div><div> TFilter const * const filter = dynamic_cast<const TFilter *>( object );</div><div><br></div><div> if( typeid( event ) == typeid( itk::InitializeEvent ) )</div><div> {</div><div> // const unsigned int currentLevel = filter->GetCurrentLevel();</div>
<div> }</div><div> else if( typeid( event ) == typeid( itk::IterationEvent ) )</div><div> {</div><div><br></div><div> std::cout << filter->GetCurrentMetricValue() << " " <<std::endl;</div>
<div> </div><div> </div><div> }</div><div> }</div><div>};</div><div><br></div><div><br></div><div>template <unsigned int VImageDimension></div><div>int PerformSimpleImageRegistration( int argc, char *argv[] )</div>
<div>{</div><div> if( argc < 3 )</div><div> {</div><div> std::cout << argv[0] << " fixedImage movingImage outputImage" << std::endl;</div><div> exit( 1 );</div><div> }</div><div>
<br></div><div> typedef float PixelType;</div><div> typedef itk::Image<PixelType, VImageDimension> FixedImageType;</div><div> typedef itk::Image<PixelType, VImageDimension> MovingImageType;</div>
<div><br></div><div> typedef itk::ImageFileReader<FixedImageType> ImageReaderType;</div><div><br></div><div> typename ImageReaderType::Pointer fixedImageReader = ImageReaderType::New();</div><div> fixedImageReader->SetFileName( argv[1] );</div>
<div> fixedImageReader->Update();</div><div> typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();</div><div> fixedImage->Update();</div><div> fixedImage->DisconnectPipeline();</div>
<div><br></div><div> typename ImageReaderType::Pointer movingImageReader = ImageReaderType::New();</div><div> movingImageReader->SetFileName( argv[2] );</div><div> movingImageReader->Update();</div><div> typename MovingImageType::Pointer movingImage = movingImageReader->GetOutput();</div>
<div> movingImage->Update();</div><div> movingImage->DisconnectPipeline();</div><div><br></div><div> typedef itk::AffineTransform<double, VImageDimension> AffineTransformType;</div><div> typedef itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType, AffineTransformType> AffineRegistrationType;</div>
<div> typedef itk::GradientDescentOptimizerv4 GradientDescentOptimizerv4Type;</div><div> typename AffineRegistrationType::Pointer affineSimple = AffineRegistrationType::New();</div><div> affineSimple->SetFixedImage( fixedImage );</div>
<div> affineSimple->SetMovingImage( movingImage );</div><div><br></div><div> // Shrink the virtual domain by specified factors for each level</div><div> //ORIGINAL Image size [217, 45, 336]</div><div> //ORIGINAL Image Dimensions 1.48 5.5 1.48 (mm)</div>
<div> </div><div> typename AffineRegistrationType::ShrinkFactorsPerDimensionContainerType ShrinkFactorsa;</div><div> typename AffineRegistrationType::ShrinkFactorsPerDimensionContainerType ShrinkFactorsb;</div><div>
typename AffineRegistrationType::ShrinkFactorsPerDimensionContainerType ShrinkFactorsc;</div><div><br></div><div> //First level</div><div> ShrinkFactorsa[0] = 3; </div><div> ShrinkFactorsa[1] = 1; </div><div>
ShrinkFactorsa[2] = 3; </div><div><br></div><div> //Second Level </div><div> ShrinkFactorsb[0] = 2; </div><div> ShrinkFactorsb[1] = 1; </div><div> ShrinkFactorsb[2] = 2; </div><div><br></div><div> //Third Level</div>
<div> ShrinkFactorsc[0] = 1;</div><div> ShrinkFactorsc[1] = 1;</div><div> ShrinkFactorsc[2] = 1;</div><div><br></div><div> unsigned int la=0;</div><div> unsigned int lb=1;</div><div> unsigned int lc=2;</div>
<div> </div><div> affineSimple->SetNumberOfLevels(3);</div><div> affineSimple->SetShrinkFactorsPerDimension( la, ShrinkFactorsa );</div><div> affineSimple->SetShrinkFactorsPerDimension( lb, ShrinkFactorsa );</div>
<div> affineSimple->SetShrinkFactorsPerDimension( lc, ShrinkFactorsb );</div><div><br></div><div><br></div><div> //Define the Metric, must use V4 metric with V4 registration</div><div> </div><div> typedef itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,</div>
<div> FixedImageType> MutualInformationMetricType;</div><div><br></div><div> </div><div> typename MutualInformationMetricType::Pointer mutualInformationMetric = MutualInformationMetricType::New();</div>
<div> mutualInformationMetric->SetNumberOfHistogramBins( 40 );</div><div> mutualInformationMetric->SetFixedImage( fixedImage );</div><div> mutualInformationMetric->SetMovingImage( movingImage );</div><div> mutualInformationMetric->SetUseMovingImageGradientFilter( false );</div>
<div> mutualInformationMetric->SetUseFixedImageGradientFilter( false );</div><div> affineSimple->SetMetric( mutualInformationMetric );</div><div><br></div><div><br></div><div> typedef itk::RegistrationParameterScalesFromPhysicalShift<MutualInformationMetricType> AffineScalesEstimatorType;</div>
<div> typename AffineScalesEstimatorType::Pointer scalesEstimator1 = AffineScalesEstimatorType::New();</div><div> scalesEstimator1->SetMetric( mutualInformationMetric );</div><div> scalesEstimator1->SetTransformForward( true );</div>
<div><br></div><div> affineSimple->SmoothingSigmasAreSpecifiedInPhysicalUnitsOn();</div><div> affineSimple->SetSmoothingSigmasAreSpecifiedInPhysicalUnits( true ); //image spacing is used to smooth in physical units</div>
<div> if( affineSimple->GetSmoothingSigmasAreSpecifiedInPhysicalUnits() != true )</div><div> {</div><div> std::cerr << "Returned unexpected value of FALSE." << std::endl;</div><div> return EXIT_FAILURE;</div>
<div> }</div><div><br></div><div> // Smooth by specified gaussian sigmas for each level. These values are specified in</div><div> // physical units. Sigmas of zero cause inconsistency between some platforms.</div><div>
{</div><div> typename AffineRegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;</div><div> smoothingSigmasPerLevel.SetSize( 1 );</div><div> smoothingSigmasPerLevel[0] = 2;</div><div> smoothingSigmasPerLevel[1] = 1;</div>
<div> smoothingSigmasPerLevel[2] = 0.3; //0;</div><div> affineSimple->SetSmoothingSigmasPerLevel( smoothingSigmasPerLevel );</div><div> }</div><div><br></div><div> typedef itk::GradientDescentOptimizerv4 GradientDescentOptimizerv4Type;</div>
<div> typename GradientDescentOptimizerv4Type::Pointer affineOptimizer =</div><div> dynamic_cast<GradientDescentOptimizerv4Type * >( affineSimple->GetOptimizer() );</div><div> if( !affineOptimizer )</div><div>
{</div><div> itkGenericExceptionMacro( "Error dynamic_cast failed" );</div><div> }</div><div> affineOptimizer->SetNumberOfIterations( 75); //need 75</div><div> //When a ScalesEstimator is assigned, the optimizer is enabled by default to estimate learning rate only once, during the first iteration</div>
<div> affineOptimizer->SetDoEstimateLearningRateOnce( false ); //true by default</div><div> affineOptimizer->SetDoEstimateLearningRateAtEachIteration( true );</div><div> affineOptimizer->SetScalesEstimator( scalesEstimator1 );</div>
<div><br></div><div> {</div><div> typedef itk::ImageToImageMetricv4<FixedImageType, MovingImageType> ImageMetricType;</div><div> typename ImageMetricType::Pointer imageMetric = dynamic_cast<ImageMetricType*>( affineSimple->GetMetric() );</div>
<div> //imageMetric->SetUseFloatingPointCorrection(true);</div><div> imageMetric->SetFloatingPointCorrectionResolution(1e4);</div><div> }</div><div><br></div><div> typedef CommandIterationUpdate< GradientDescentOptimizerv4Type > AffineCommandType;</div>
<div> AffineCommandType::Pointer affineobserver = AffineCommandType::New();</div><div> affineOptimizer->AddObserver( itk::IterationEvent(), affineobserver );</div><div><br></div><div> try</div><div> {</div><div> std::cout << "Performing Affine:" << std::endl;</div>
<div> affineSimple->Update();</div><div> </div><div> std::cout << "Last affine iteration is " << affineOptimizer->GetCurrentIteration();</div><div> </div><div> }</div><div> catch( itk::ExceptionObject &e )</div>
<div> {</div><div> std::cerr << "Exception caught: " << e << std::endl;</div><div> return EXIT_FAILURE;</div><div> }</div><div><br></div><div> {</div><div> typedef itk::ImageToImageMetricv4<FixedImageType, MovingImageType> ImageMetricType;</div>
<div> typename ImageMetricType::Pointer imageMetric = dynamic_cast<ImageMetricType*>( affineOptimizer->GetMetric() );</div><div> std::cout << "Affine parameters after registration: " << std::endl</div>
<div> << affineOptimizer->GetCurrentPosition() << std::endl</div><div> << "Last LearningRate: " << affineOptimizer->GetLearningRate() << std::endl</div><div>
<< "Use FltPtCorrex: " << imageMetric->GetUseFloatingPointCorrection() << std::endl</div><div> << "FltPtCorrexRes: " << imageMetric->GetFloatingPointCorrectionResolution() << std::endl</div>
<div> << "Number of threads used: metric: " << imageMetric->GetNumberOfThreadsUsed()</div><div> << std::endl << " optimizer: " << affineOptimizer->GetNumberOfThreads() << std::endl;</div>
<div> }</div><div> </div><div> </div><div> //Now create A composite transform of the input</div><div> </div><div> typedef typename AffineRegistrationType::RealType RealType;</div><div><br></div><div> typedef itk::CompositeTransform<RealType, VImageDimension> CompositeTransformType;</div>
<div> typename CompositeTransformType::Pointer compositeTransform = CompositeTransformType::New();</div><div> compositeTransform->AddTransform( const_cast<typename AffineRegistrationType::OutputTransformType *>( affineSimple->GetOutput()->Get() ) );</div>
<div><br></div><div>////////////////</div><div>//Now for the BsplineSYN method</div><div><br></div><div> //DEFINE THE DISPLACEMENT FIELD</div><div> typedef itk::Vector<RealType, VImageDimension> VectorType;</div><div>
VectorType zeroVector( 0.0 );</div><div> typedef itk::Image<VectorType, VImageDimension> DisplacementFieldType;</div><div> typename DisplacementFieldType::Pointer displacementField = DisplacementFieldType::New();</div>
<div> displacementField->CopyInformation( fixedImage );</div><div> displacementField->SetRegions( fixedImage->GetBufferedRegion() );</div><div> displacementField->Allocate();</div><div> displacementField->FillBuffer( zeroVector );</div>
<div><br></div><div> //DEFINE THE INVERSE DISPLACEMENT FIELD</div><div> typename DisplacementFieldType::Pointer inverseDisplacementField = DisplacementFieldType::New();</div><div> inverseDisplacementField->CopyInformation( fixedImage );</div>
<div> inverseDisplacementField->SetRegions( fixedImage->GetBufferedRegion() );</div><div> inverseDisplacementField->Allocate();</div><div> inverseDisplacementField->FillBuffer( zeroVector );</div><div><br></div>
<div> //DEFINE THE SMOOTHING ON UPDATE FIELD WHICH PARAMETERIZEDS THE DISPLACEMENT FIELD</div><div> typedef itk::BSplineSmoothingOnUpdateDisplacementFieldTransform<RealType, VImageDimension></div><div> BSplineDisplacementFieldTransformType;</div>
<div><br></div><div><br></div><div> //DEFINE THE SYN IMAGE REGISTRATION METHOD</div><div> typedef itk::BSplineSyNImageRegistrationMethod<FixedImageType, MovingImageType, BSplineDisplacementFieldTransformType> DisplacementFieldRegistrationType;</div>
<div> typename DisplacementFieldRegistrationType::Pointer displacementFieldRegistration = DisplacementFieldRegistrationType::New();</div><div><br></div><div> //OPTIMIZER VALUES</div><div><br></div><div> //DEFINE THE OUTPUT TRANSFORM - parameterized by a Bspline///</div>
<div> typename BSplineDisplacementFieldTransformType::Pointer outputDisplacementFieldTransform =</div><div> const_cast<BSplineDisplacementFieldTransformType *>( displacementFieldRegistration->GetOutput()->Get() );</div>
<div><br></div><div> //CREATE THE TRANSFORM ADAPTOR</div><div> typedef itk::BSplineSmoothingOnUpdateDisplacementFieldTransformParametersAdaptor< BSplineDisplacementFieldTransformType > DisplacementFieldTransformAdaptorType;</div>
<div> typename DisplacementFieldRegistrationType::TransformParametersAdaptorsContainerType adaptors;</div><div><br></div><div><br></div><div> //SET THE NUMBER OF LEVELS AND SHRINK FACTORS AND SMOOTHING VALUES FOR THE VIRTUAL DOMAIN</div>
<div> unsigned int numberOfLevels = 3;</div><div><br></div><div> typename DisplacementFieldRegistrationType::NumberOfIterationsArrayType numberOfIterationsPerLevel;</div><div> numberOfIterationsPerLevel.SetSize( 3 );</div>
<div> numberOfIterationsPerLevel[0] = 500;</div><div> numberOfIterationsPerLevel[1] = 300;</div><div> numberOfIterationsPerLevel[2] = 300;</div><div><br></div><div> typename DisplacementFieldRegistrationType::ShrinkFactorsPerDimensionContainerType ShrinkFactorsaa;</div>
<div> typename DisplacementFieldRegistrationType::ShrinkFactorsPerDimensionContainerType ShrinkFactorsbb;</div><div> typename DisplacementFieldRegistrationType::ShrinkFactorsPerDimensionContainerType ShrinkFactorscc;</div>
<div> typename DisplacementFieldRegistrationType::ShrinkFactorsPerDimensionContainerType ShrinkFactorstemp;</div><div><br></div><div> //First level</div><div> ShrinkFactorsaa[0] = 3; </div><div> ShrinkFactorsaa[1] = 1; </div>
<div> ShrinkFactorsaa[2] = 3; </div><div><br></div><div> //Second Level </div><div> ShrinkFactorsbb[0] = 2; </div><div> ShrinkFactorsbb[1] = 1; </div><div> ShrinkFactorsbb[2] = 2; </div><div><br></div><div>
//Third Level</div><div> ShrinkFactorscc[0] = 1;</div><div> ShrinkFactorscc[1] = 1;</div><div> ShrinkFactorscc[2] = 1;</div><div><br></div><div><br></div><div> displacementFieldRegistration->SetNumberOfLevels(3);</div>
<div> displacementFieldRegistration->SetShrinkFactorsPerDimension( la, ShrinkFactorsaa );</div><div> displacementFieldRegistration->SetShrinkFactorsPerDimension( lb, ShrinkFactorsbb );</div><div> displacementFieldRegistration->SetShrinkFactorsPerDimension( lc, ShrinkFactorscc );</div>
<div><br></div><div> typename DisplacementFieldRegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;</div><div> smoothingSigmasPerLevel.SetSize( 3 );</div><div> smoothingSigmasPerLevel[0] = 2;</div><div> smoothingSigmasPerLevel[1] = 1;</div>
<div> smoothingSigmasPerLevel[2] = 0.3;</div><div><br></div><div><br></div><div> for( unsigned int level = 0; level < numberOfLevels; level++ )</div><div> {</div><div> // We use the shrink image filter to calculate the fixed parameters of the virtual</div>
<div> // domain at each level.</div><div><br></div><div> if (level==0)</div><div> </div><div><span class="" style="white-space:pre"> </span>{</div><div> ShrinkFactorstemp=ShrinkFactorsaa;</div><div>
<span class="" style="white-space:pre"> </span> }</div><div><br></div><div> if (level==1)</div><div> </div><div><span class="" style="white-space:pre"> </span>{</div><div> ShrinkFactorstemp=ShrinkFactorsbb;</div>
<div><span class="" style="white-space:pre"> </span> }</div><div><br></div><div> if (level==2)</div><div> </div><div><span class="" style="white-space:pre"> </span>{</div><div> ShrinkFactorstemp=ShrinkFactorscc;</div>
<div><span class="" style="white-space:pre"> </span> }</div><div> </div><div> typedef itk::ShrinkImageFilter<DisplacementFieldType, DisplacementFieldType> ShrinkFilterType;</div><div> typename ShrinkFilterType::Pointer shrinkFilter = ShrinkFilterType::New();</div>
<div> </div><div> std::cout << "using shrink Factor" << ShrinkFactorstemp <<std::endl;</div><div> shrinkFilter->SetShrinkFactors( ShrinkFactorstemp);</div><div> shrinkFilter->SetInput( displacementField );</div>
<div> shrinkFilter->Update();</div><div> std::cout << "New size " << shrinkFilter->GetOutput()->GetBufferedRegion().GetSize() <<std::endl;</div><div> </div><div> //int VirtSize[3];</div>
<div> typename FixedImageType::SizeType VirtSize =shrinkFilter->GetOutput()->GetBufferedRegion().GetSize();</div><div><br></div><div><br></div><div> typename DisplacementFieldTransformAdaptorType::Pointer fieldTransformAdaptor = DisplacementFieldTransformAdaptorType::New();</div>
<div> fieldTransformAdaptor->SetRequiredSpacing( shrinkFilter->GetOutput()->GetSpacing() );</div><div> fieldTransformAdaptor->SetRequiredSize( shrinkFilter->GetOutput()->GetBufferedRegion().GetSize() );</div>
<div> fieldTransformAdaptor->SetRequiredDirection( shrinkFilter->GetOutput()->GetDirection() );</div><div> fieldTransformAdaptor->SetRequiredOrigin( shrinkFilter->GetOutput()->GetOrigin() );</div><div>
fieldTransformAdaptor->SetTransform( outputDisplacementFieldTransform );</div><div><br></div><div> typename BSplineDisplacementFieldTransformType::ArrayType newUpdateMeshSize = newUpdateMeshSize;</div><div> typename BSplineDisplacementFieldTransformType::ArrayType newTotalMeshSize = newTotalMeshSize;</div>
<div> </div><div> if (level==0)</div><div> </div><div> {</div><div><span class="" style="white-space:pre"> </span>//First dimension</div><div><span class="" style="white-space:pre"> </span>newUpdateMeshSize[0] = 0; //Lat</div>
<div><span class="" style="white-space:pre"> </span>newTotalMeshSize[0] = VirtSize[0]/12;</div><div><span class="" style="white-space:pre"> </span>//Second dimension</div><div><span class="" style="white-space:pre"> </span>newUpdateMeshSize[1] = 0; //Ant-post</div>
<div><span class="" style="white-space:pre"> </span>newTotalMeshSize[1] = VirtSize[0]/12; //works better</div><div><span class="" style="white-space:pre"> </span>//Third dimension</div><div><span class="" style="white-space:pre"> </span>newUpdateMeshSize[2] = 0; //Sup-Inf</div>
<div><span class="" style="white-space:pre"> </span>newTotalMeshSize[2] = VirtSize[2]/12;</div><div> </div><div> } </div><div> </div><div> </div><div> //DEFINE THE MESH SIZE FOR THE UPDATE FIELD AND TOTAL FIELD</div>
<div> </div><div> if (level==1)</div><div> </div><div> {</div><div> </div><div><span class="" style="white-space:pre"> </span>//First dimension</div><div><span class="" style="white-space:pre"> </span>newUpdateMeshSize[0] = 0; //Lat</div>
<div><span class="" style="white-space:pre"> </span>newTotalMeshSize[0] = VirtSize[0]/12;</div><div><span class="" style="white-space:pre"> </span>//Second dimension</div><div><span class="" style="white-space:pre"> </span>newUpdateMeshSize[1] = 0; //Ant-post</div>
<div><span class="" style="white-space:pre"> </span>newTotalMeshSize[1] = VirtSize[0]/12;</div><div><span class="" style="white-space:pre"> </span>//Third dimension</div><div><span class="" style="white-space:pre"> </span>newUpdateMeshSize[2] = 0; //Sup-Inf</div>
<div><span class="" style="white-space:pre"> </span>newTotalMeshSize[2] = VirtSize[2]/12;</div><div> </div><div> } </div><div> </div><div> if (level==2)</div><div> </div><div> {</div><div><span class="" style="white-space:pre"> </span>//First dimension</div>
<div><span class="" style="white-space:pre"> </span>newUpdateMeshSize[0] = 0; //Lat</div><div><span class="" style="white-space:pre"> </span>newTotalMeshSize[0] = VirtSize[0]/12;</div><div><span class="" style="white-space:pre"> </span>//Second dimension</div>
<div><span class="" style="white-space:pre"> </span>newUpdateMeshSize[1] = 0; //Ant-post</div><div><span class="" style="white-space:pre"> </span>newTotalMeshSize[1] = VirtSize[0]/12;</div><div><span class="" style="white-space:pre"> </span>//Third dimension</div>
<div><span class="" style="white-space:pre"> </span>newUpdateMeshSize[2] = 0; //Sup-Inf</div><div><span class="" style="white-space:pre"> </span>newTotalMeshSize[2] = VirtSize[2]/12;</div><div> } </div><div> </div>
<div> </div><div> std::cout << "level is " << level <<std::endl;</div><div> std::cout << "newUpdateMeshSize" << newUpdateMeshSize << std::endl;</div><div> std::cout << "newTotalMeshSize" << newTotalMeshSize << std::endl;</div>
<div> </div><div> </div><div> fieldTransformAdaptor->SetMeshSizeForTheUpdateField( newUpdateMeshSize );</div><div> fieldTransformAdaptor->SetMeshSizeForTheTotalField( newTotalMeshSize );</div><div><br>
</div><div> adaptors.push_back( fieldTransformAdaptor.GetPointer() );</div><div> }</div><div><br></div><div> </div><div> typename MutualInformationMetricType::Pointer mutualInformationMetricBspline = MutualInformationMetricType::New();</div>
<div> mutualInformationMetricBspline->SetNumberOfHistogramBins( 40 );</div><div> mutualInformationMetricBspline->SetFixedImage( fixedImage );</div><div> mutualInformationMetricBspline->SetMovingImage( movingImage );</div>
<div> mutualInformationMetricBspline->SetUseMovingImageGradientFilter( false );</div><div> mutualInformationMetricBspline->SetUseFixedImageGradientFilter( false );</div><div><br></div><div> displacementFieldRegistration->SetFixedImage( fixedImage );</div>
<div> displacementFieldRegistration->SetMovingImage( movingImage );</div><div> displacementFieldRegistration->SetNumberOfLevels( 3 );</div><div> displacementFieldRegistration->SetMovingInitialTransform( compositeTransform );</div>
<div> displacementFieldRegistration->SmoothingSigmasAreSpecifiedInPhysicalUnitsOn();</div><div> displacementFieldRegistration->SetSmoothingSigmasPerLevel( smoothingSigmasPerLevel );</div><div> displacementFieldRegistration->SetMetric( mutualInformationMetricBspline );</div>
<div> displacementFieldRegistration->SetLearningRate( 0.25 );</div><div> //in SYN methods learning rate is set initially and modified at each iteration</div><div> //no need for scale as displacement in x,y,z has equal scale</div>
<div> </div><div> displacementFieldRegistration->SetTransformParametersAdaptorsPerLevel( adaptors );</div><div> displacementFieldRegistration->SetNumberOfIterationsPerLevel( numberOfIterationsPerLevel );</div><div>
displacementFieldRegistration->SetTransformParametersAdaptorsPerLevel( adaptors );</div><div> outputDisplacementFieldTransform->SetDisplacementField( displacementField );</div><div> outputDisplacementFieldTransform->SetInverseDisplacementField( inverseDisplacementField );</div>
<div> outputDisplacementFieldTransform->SetSplineOrder( 3 );</div><div><br></div><div> typedef CommandIterationUpdate <DisplacementFieldRegistrationType> DisplacementFieldCommandType;</div><div> typename DisplacementFieldCommandType::Pointer DisplacementFieldObserver = DisplacementFieldCommandType::New();</div>
<div> displacementFieldRegistration->AddObserver( itk::IterationEvent(), DisplacementFieldObserver );</div><div><br></div><div><br></div><div> try</div><div> {</div><div> std::cout << " " << std::endl;</div>
<div> std::cout << "BSpline SyN registration" << std::endl;</div><div> displacementFieldRegistration->Update();</div><div> }</div><div> catch( itk::ExceptionObject &e )</div><div> {</div>
<div> std::cerr << "Exception caught: " << e << std::endl;</div><div> return EXIT_FAILURE;</div><div> }</div><div><br></div><div> compositeTransform->AddTransform( outputDisplacementFieldTransform );</div>
<div><br></div><div><br></div><div> typedef itk::ResampleImageFilter<MovingImageType, FixedImageType> ResampleFilterType;</div><div> typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();</div>
<div> resampler->SetTransform( compositeTransform );</div><div> resampler->SetInput( movingImage );</div><div> resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );</div><div> resampler->SetOutputOrigin( fixedImage->GetOrigin() );</div>
<div> resampler->SetOutputSpacing( fixedImage->GetSpacing() );</div><div> resampler->SetOutputDirection( fixedImage->GetDirection() );</div><div> resampler->SetDefaultPixelValue( 0 );</div><div> resampler->Update();</div>
<div><br></div><div> typedef itk::ImageFileWriter<FixedImageType> WriterType;</div><div> typename WriterType::Pointer writer = WriterType::New();</div><div> writer->SetFileName( argv[3] );</div><div> writer->SetInput( resampler->GetOutput() );</div>
<div> writer->Update();</div><div><br></div><div> return EXIT_SUCCESS;</div><div>}</div><div><br></div><div>int main( int argc, char *argv[] )</div><div>{</div><div> if( argc < 3 )</div><div> {</div><div> std::cout << argv[0] << " fixedImage movingImage outputImage numberOfAffineIterations numberOfDeformableIterations" << std::endl;</div>
<div> exit( 1 );</div><div> }</div><div><br></div><div><br></div><div> const int imgDim=3;</div><div><br></div><div> PerformSimpleImageRegistration<imgDim>( argc, argv );</div><div> </div><div> </div>
<div> return EXIT_SUCCESS;</div><div>}</div><div><br></div></div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Jan 22, 2014 at 5:01 PM, brian avants <span dir="ltr"><<a href="mailto:stnava@gmail.com" target="_blank">stnava@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">multi-resolution set up is usually the key to such issues.<div><br></div><div>but greater detail in your question / code is needed. </div>
<div><br></div><div>syn & variants are applicable for a broad range of images</div>
<div><br></div><div>as you can see in the literature & documentation.</div></div><div class="gmail_extra"><span class="HOEnZb"><font color="#888888"><br clear="all"><div><div><br></div>brian<br><div><br></div><div><br>
</div></div></font></span><div><div class="h5">
<br><br><div class="gmail_quote">On Wed, Jan 22, 2014 at 11:57 AM, Emma Saunders <span dir="ltr"><<a href="mailto:emmasaunders123@gmail.com" target="_blank">emmasaunders123@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr">Thanks Nick and Brian for your input.<div><br></div><div>I was attempting to change the learning rate as I am having poor results using both syn and itkBsplinesyn when registering two MR images of the entire thorax during respiration. Are these methods specifically designed for certain types of images (e.g brain). I initialize my framework with an Affine transform, using mattes MI, the metric value first falls off and then increases. The result is thus only a small amount of movement, around 1/3 of what I would hope for.</div>
<div><br></div><div>I can successfully register the images with a regular Bspline method using ITK, but was hoping to get a SYN method working in order to have a valid inverse.</div><div><br></div><div>I suppose the question is, if I can get regular Bspline to give a descent result should the same hold for the SYN methods or am I missing something?</div>
<div><br></div><div>Thanks for your help</div><div><br></div><div><br></div></div><div><div><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Jan 22, 2014 at 3:30 PM, Emma Saunders <span dir="ltr"><<a href="mailto:emmasaunders123@gmail.com" target="_blank">emmasaunders123@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi all<div><br></div><div>I am trying to use the Syn methods mentioned above.</div><div><br></div><div>
I would simply like to set my framework up to estimate the learning rate at each iteration. So far this has proved problematic as the optimization process is contained within the methods itself and whenever I try to find a solution need to look through ANTS header files and I still have no luck.</div>
<div><br></div><div>Is there a simple way to set the learning rate given a scaleEstimator for these methods (without using the ans registration helper ) or must one simply choose one at the beginning and hope for the best??</div>
<div><br></div><div>I am also finding it difficult to determine which optimizer these methods use as in the ANTS headers there are two options a GradientDescent and a conjugate gradient descent, any ideas??</div>
<div><br></div><div>I initially used the ITK libraries as they provided a method that I could begin using state of the art image processing algorithms without the steep learning curve, is there any pending documentation or software guide for the new ITKV4 framework.</div>
<div><br></div><div>Thanks</div><span><font color="#888888"><div><br></div><div>Emma</div></font></span></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div></div></div>
</blockquote></div><br></div>