<div dir="ltr">Hi Nick,<div><br></div><div style>Thanks for your very detailed and helpful response. I really appreciate it. I am making a comparison of the BsplineSYN with SYN. Thus far I have quite a bit of residual motion with BsplineSYN</div>
<div style><br></div><div style>For the run you made, are you able to recall what values of: updateFieldMeshSizeAtBaseLevel and totalFieldMeshSizeAtBaseLevel you chose and the tradeoff's these have?</div><div style><br>
</div><div style>Below is the run I am making in ants for reference</div><div style><br></div><div style><br></div><div style><div>antsRegistration --dimensionality 3 --float 1 --output [BsplineSYN,BsplineSYN.nii, BsplineSYNinv.nii] --interpolation Linear --initial-moving-transform [resp_pos_5.nii,resp_pos_8.nii,1] --transform Rigid[0.1] --metric MI[resp_pos_5.nii,resp_pos_8.nii,1,32,Regular,0.25] --convergence [1000x500x25x10,1e-6,10] --shrink-factors 12x8x4x2 --smoothing-sigmas 4x3x2x1mm --transform Affine[0.1] --metric MI[resp_pos_5.nii,resp_pos_8.nii,1,32,Regular,0.25] --convergence [1000x500x25x10,1e-6,10] --shrink-factors 12x8x4x2 --smoothing-sigmas 4x3x2x1mm --transform BsplineSyN[0.1,10,0,3] --metric CC[resp_pos_5.nii,resp_pos_8.nii,1,4,Regular,0.25] --convergence [100x100x70x50x10,0,10] --shrink-factors 10x6x4x2x1 --smoothing-sigmas 6x4x2x1x0mm --transform BsplineSyN[0.1,10,0,3] --metric Demons[resp_pos_5.nii,resp_pos_8.nii,1,4,Regular,0.25] --convergence [200,0,10] --shrink-factors 1 --smoothing-sigmas 0mm </div>
<div><br></div></div><div style><div style>Many thanks</div><div style><br></div><div style>Emma</div></div><div style><br></div><div style><br></div><div style><br></div><div style><br></div></div><div class="gmail_extra">
<br><br><div class="gmail_quote">On Fri, Jan 24, 2014 at 4:10 PM, Nicholas Tustison <span dir="ltr"><<a href="mailto:ntustison@gmail.com" target="_blank">ntustison@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 style="word-wrap:break-word">Hi Emma,<div><br></div><div>This is a good example of a difficult registration problem and </div><div>illustrates why we often refer people to ANTs instead of the </div><div>simpler registration tests within ITK which demonstrate usage. </div>
<div><br></div><div>There are some preprocessing issues you might want to </div><div>consider such as bias correction and histogram matching. </div><div>More directly related to the image registration part, though, </div>
<div>is the relatively extreme anisotropy and the deformation</div><div>required to register lungs, particularly at the base, during </div><div>different inspiration phases. </div><div><br></div><div>Obviously, you should check that the affine part completed</div>
<div>in an expected way. If the affine registration step failed,</div><div>you can’t expect the deformable part to recover. Assuming</div><div>that the affine registration looks reasonable, I would use the</div><div>neighborhood cross-correlation which is our go-to metric</div>
<div>for deformable registration.</div><div><br></div><div><a href="http://www.itk.org/Doxygen/html/classitk_1_1ANTSNeighborhoodCorrelationImageToImageMetricv4.html" target="_blank">http://www.itk.org/Doxygen/html/classitk_1_1ANTSNeighborhoodCorrelationImageToImageMetricv4.html</a></div>
<div><br></div><div>When I ran your example, I used a multi-resolution </div><div>deformable approach with 4 levels. I also used the CC metric </div><div>with a radius of 4. The radius of 4 is necessary to capture the </div>
<div>basal pulmonary motion. However, the relatively large anisotropy </div><div>doesn’t capture the finer motion at the higher resolution so I ran </div><div>another B-spline SyN image registration stage using the Demons </div>
<div>metric for 200 iterations at the highest level which gave a pretty </div><div>good result. </div><div><br></div><div>So all put together, I needed the following registration stages:</div><div><br></div><div>1. Initialize with a translation transform by aligning the centers </div>
<div> of intensity mass in both images.</div><div>2. Follow-up with a rigid registration using the MI metric, sampled </div><div> at a rate of 25%, for 4 multi-resolution levels (1000,500,250,100</div><div> iterations), learning rate = 0.1.</div>
<div>3. Compose the previous transforms with an affine registration</div><div><div> using the MI metric, sampled at a rate of 25%, for 4 multi-resolution </div><div> levels (1000,500,250,100 iterations) and a learning rate = 0.1.</div>
</div><div>4. Compose the previous transforms with a B-spline SyN registration </div><div> using the CC metric of a radius of 4, with 4 multi-resolution </div><div> levels, with 100,70,50,20 iterations and a learning rate of 0.1</div>
<div>5. Finally, compose the previous transforms with a B-spline SyN</div><div> registration using the Demons metric at the finest resolution level</div><div> with 200 iterations and a learning rate of 0.25.</div>
<div><br></div><div>The first 3 transforms are all linear transforms so I collapsed them to </div><div>a single affine transform file. The last two transforms are displacement</div><div>field transforms with inverses which I also composed to a single</div>
<div>invertible displacement transform.</div><div><br></div><div>There might be other ways of getting a similar (or even better solution)</div><div>but this is something that seems to work and it’s something that you</div>
<div>can put together using classes from the new ITKv4 framework which </div><div>has the flexibility to permit complex registration solutions such as this.</div><div><br></div><div>Putting together an application-level program which takes advantage</div>
<div>of this flexibility can be quite involved which is why we wrote </div><div>antsRegistration.cxx as part of ANTs and often point people to it on</div><div>the ITK list. antsRegistration takes full advantage of the flexibility </div>
<div>we built into the new ITKv4 registration framework which you just would</div><div>not put into test code or simple registration examples.</div><div><br></div><div>Nick</div><div><div class="h5"><div><br></div><div><br>
</div><div><br></div><div><br></div><div><div><div>On Jan 23, 2014, at 8:41 AM, Emma Saunders <<a href="mailto:emmasaunders123@gmail.com" target="_blank">emmasaunders123@gmail.com</a>> wrote:</div><br><blockquote type="cite">
<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>If you can see anything obvious from my code that would be great.</div><div><br></div><div>Thanks for your time I really appreciate it, Kind regards</div><div><br></div><div>
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 style="white-space:pre-wrap"> </span>{</div><div> ShrinkFactorstemp=ShrinkFactorsaa;</div><div>
<span style="white-space:pre-wrap"> </span> }</div><div><br></div><div> if (level==1)</div><div> </div><div><span style="white-space:pre-wrap"> </span>{</div><div> ShrinkFactorstemp=ShrinkFactorsbb;</div>
<div><span style="white-space:pre-wrap"> </span> }</div><div><br></div><div> if (level==2)</div><div> </div><div><span style="white-space:pre-wrap"> </span>{</div><div> ShrinkFactorstemp=ShrinkFactorscc;</div>
<div><span style="white-space:pre-wrap"> </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 style="white-space:pre-wrap"> </span>//First dimension</div><div><span style="white-space:pre-wrap"> </span>newUpdateMeshSize[0] = 0; //Lat</div>
<div><span style="white-space:pre-wrap"> </span>newTotalMeshSize[0] = VirtSize[0]/12;</div><div><span style="white-space:pre-wrap"> </span>//Second dimension</div><div><span style="white-space:pre-wrap"> </span>newUpdateMeshSize[1] = 0; //Ant-post</div>
<div><span style="white-space:pre-wrap"> </span>newTotalMeshSize[1] = VirtSize[0]/12; //works better</div><div><span style="white-space:pre-wrap"> </span>//Third dimension</div><div><span style="white-space:pre-wrap"> </span>newUpdateMeshSize[2] = 0; //Sup-Inf</div>
<div><span style="white-space:pre-wrap"> </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 style="white-space:pre-wrap"> </span>//First dimension</div><div><span style="white-space:pre-wrap"> </span>newUpdateMeshSize[0] = 0; //Lat</div>
<div><span style="white-space:pre-wrap"> </span>newTotalMeshSize[0] = VirtSize[0]/12;</div><div><span style="white-space:pre-wrap"> </span>//Second dimension</div><div><span style="white-space:pre-wrap"> </span>newUpdateMeshSize[1] = 0; //Ant-post</div>
<div><span style="white-space:pre-wrap"> </span>newTotalMeshSize[1] = VirtSize[0]/12;</div><div><span style="white-space:pre-wrap"> </span>//Third dimension</div><div><span style="white-space:pre-wrap"> </span>newUpdateMeshSize[2] = 0; //Sup-Inf</div>
<div><span style="white-space:pre-wrap"> </span>newTotalMeshSize[2] = VirtSize[2]/12;</div><div> </div><div> } </div><div> </div><div> if (level==2)</div><div> </div><div> {</div><div><span style="white-space:pre-wrap"> </span>//First dimension</div>
<div><span style="white-space:pre-wrap"> </span>newUpdateMeshSize[0] = 0; //Lat</div><div><span style="white-space:pre-wrap"> </span>newTotalMeshSize[0] = VirtSize[0]/12;</div><div><span style="white-space:pre-wrap"> </span>//Second dimension</div>
<div><span style="white-space:pre-wrap"> </span>newUpdateMeshSize[1] = 0; //Ant-post</div><div><span style="white-space:pre-wrap"> </span>newTotalMeshSize[1] = VirtSize[0]/12;</div><div><span style="white-space:pre-wrap"> </span>//Third dimension</div>
<div><span style="white-space:pre-wrap"> </span>newUpdateMeshSize[2] = 0; //Sup-Inf</div><div><span style="white-space:pre-wrap"> </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><font color="#888888"><br clear="all"><div><div><br></div>brian<br><div><br></div><div><br>
</div></div></font></span><div><div>
<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 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></blockquote></div><br></div></div></div>
</blockquote></div><br></div>
</blockquote></div><br></div></div></div></div></blockquote></div><br></div>