[Insight-users] DeformableRegistration6 error
Luis Ibanez
luis.ibanez at kitware.com
Mon Nov 8 20:10:30 EST 2004
Hi Kingaza,
As we explained in the previous email, the error condition
that you are encountering is the result of a premature
termination of the optimizer when it reaches its convergence
criteria.
You want to work in fine tunning the values of the method
SetProjectedGradientTolerance()
The LBFGSBOptimizer will stop when the actual gradient is
smaller than the value that you set. Therefore, if the
optimizer is stopping prematurely that means that your
gradient criteria was too relaxed, e.g. the value was too
large. You probably want to reduce this settinng to values
in the range 1e-6 to 1e-10.
Note that this is closely related to the Image Metric that
you are using, since after all, the derivatives of the metric
are proportional to the magnitude of the metric.
You may also want to fine tune the value of the method
SetCostFunctionConvergenceFactor()
because the optimizer will stop whenever the reduction
in the cost function is less than (factor * machine precision).
The values recommended in the documentation are
1e+12 for low accuracy (may result in early terminatino)
1e+7 for moderate accuracy
1e+1 for high accuracy
Note again that this is proportional to the magnitudes
of the Image Metric that you are using.
It is quite different if you are using Mutual Information
or MeanSquares or NormalizedCrossCorrelation.
Something that may help you at this point is to work with
a very small crooped version of your image. That will
help you run multiple experiments rapidly until you
figure out the parameters that may be appropriate for your
type of image and type of ImageMetric.
Regards,
Luis
------------------------
kingaza at gmail.com wrote:
> I have read the example DeformableRegistration7.cxx, but still have no
> idea about some parameters of LBFGSB optimizer.
> when i run deformableRegistration6 using LBFGSBOptimizer, it seems
> work well in low translation, but in high translation, it do nothing
> but throw a exception. I think this is because of some worng
> configurations to the optimizer, but i am not sure.
> do you think I should read a paper for this optimizer? or is there any
> tips to it?
> thx again!
>
> the attach is my modified code
>
> these are the infomation in command window:
> Starting Registration with low resolution transform
> 0 404.8 7.64756
> 1 241.448 5.85265
> ......
> Starting Registration with high resolution transform
> ExceptionObject caught !
>
> itk::ExceptionObject (0109F7C4)
> Location: "Unknown"
> File: .\itkLBFGSBOptimizer.cxx
> Line: 212
> Descritption: itk::ERROR: LBFGSBOPtimizer(003B67F0): LowerBound array
> does not have sufficient number element
>
>
>
> On Tue, 02 Nov 2004 13:35:10 -0500, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>
>>Hi Kingaza,
>>
>>Please post to the list the error message that you
>>get from the Exception.
>>
>>Did you follow the guidelines of the example
>>
>> DeformableRegistration7.cxx ??
>>
>>Please let us know,
>>
>> Thanks
>>
>>
>>
>>
>> Luis
>>
>>------------------------
>>kingaza at gmail.com wrote:
>>
>>
>>>thx for your help
>>>and now I try to use LBGSBOptimizer inplace of LBGSOptimizer
>>>but it seems that the two optimizers are different in use
>>>for examle,
>>>LBGSBOptimzer does not hava the member function :SetDefaultStepLength()
>>>
>>>and it seems some problems if I do in this way:
>>>
>>>//for transformLow
>>>/* //for LBFGSOptimizer
>>> optimizer->SetGradientConvergenceTolerance( 0.05 );
>>> optimizer->SetLineSearchAccuracy( 0.9 );
>>> optimizer->SetDefaultStepLength( 1.5 );
>>> optimizer->TraceOn();
>>> optimizer->SetMaximumNumberOfFunctionEvaluations( 1000 );
>>>*/
>>> //for LBFGSOptimizer
>>> OptimizerType::BoundSelectionType boundSelect(
>>>transformLow->GetNumberOfParameters() );
>>> OptimizerType::BoundValueType upperBound(
>>>transformLow->GetNumberOfParameters() );
>>> OptimizerType::BoundValueType lowerBound(
>>>transformLow->GetNumberOfParameters() );
>>>
>>> boundSelect.Fill( 0 );
>>> upperBound.Fill( 0.0 );
>>> lowerBound.Fill( 0.0 );
>>>
>>> optimizer->SetBoundSelection( boundSelect );
>>> optimizer->SetUpperBound( upperBound );
>>> optimizer->SetLowerBound( lowerBound );
>>>
>>> optimizer->SetCostFunctionConvergenceFactor( 1e+12 );
>>> optimizer->SetProjectedGradientTolerance( 1.0 );
>>> optimizer->SetMaximumNumberOfIterations( 500 );
>>> optimizer->SetMaximumNumberOfEvaluations( 500 );
>>> optimizer->SetMaximumNumberOfCorrections( 5 );
>>>
>>>//for transformHigh
>>>/* //for LBFGSOptimizer
>>> optimizer->SetGradientConvergenceTolerance( 0.01 );
>>> optimizer->SetDefaultStepLength( 0.25 );
>>>*/
>>> //for LBFGSOptimizer
>>> optimizer->SetProjectedGradientTolerance( 0.01 );
>>>
>>>but it can work only in the first registraion, and throw a exception
>>>in the second registration.
>>>
>>>could you tell me how to use LBGSBOptimizer in a multi-resolution registration?
>>>
>>>On Sat, 30 Oct 2004 11:52:43 -0400, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>>>
>>>
>>>>Hi Kingaza,
>>>>
>>>>This is a known problem of the LBGSOptimizer that tends to happen
>>>>when at the first iteration the optimizer doesn't find a better point
>>>>to move in the parametric space.
>>>>
>>>>It is unrelated to the fact that you are reading DICOM files and it is
>>>>unrelated to the image size. The real cause is that you are inadvertedly
>>>>starting the optimizer in a local optimium.
>>>>
>>>>Please try the new example
>>>>
>>>> DeformableRegistration7.cxx
>>>>
>>>>That uses the modified optimizer
>>>>
>>>> LBFGSBOptimizer
>>>>
>>>>This example was committed recently in CVS, so you will
>>>>have to update your CVS checkout in order to get the file.
>>>>
>>>>http://www.itk.org/cgi-bin/viewcvs.cgi/Examples/Registration/DeformableRegistration7.cxx?rev=1.1&root=Insight&sortby=date&view=log
>>>>
>>>> Regards,
>>>>
>>>> Luis
>>>>
>>>>------------------------------------
>>>>
>>>>
>>>>kingaza at gmail.com wrote:
>>>>
>>>>
>>>>
>>>>>hi all,
>>>>>
>>>>>I try DeformableRegistration6 using images which size is 256*256, but
>>>>>an exception occures:
>>>>>
>>>>>itk::ExceptionObject (0108F840)
>>>>>Location: "Unknown"
>>>>>File: .\itkLBFGSOptimizer.cxx
>>>>>LIne: 252
>>>>>Description: itk::ERROR: LBFGSOpimeizer(00376680): Error occured in optimization
>>>>>
>>>>>and then i debug deeply, and find in itkLBFGSOptimizer.cxx
>>>>>parameters.size = 0 while initialPosition.size() =128
>>>>>this is where the exception is caught.
>>>>>I had thought it is because of the dicom file format , so I try it
>>>>>again using png file, and the same appears.
>>>>>But if I use images which size is 64*64, anything is OK, and the
>>>>>output is perfect.
>>>>>please help me!
>>>>>
>>>>>
>>>>>
>>>>>itkLBFGSOptimizer.cxx
>>>>>/////////////////////////////////////////////////////////////////////////////////////////////
>>>>>......
>>>>>m_VnlOptimizer->minimize( parameters );
>>>>>
>>>>>if ( parameters.size() != initialPosition.size() )
>>>>> {
>>>>> // set current position to initial position and throw an exception
>>>>> this->SetCurrentPosition( initialPosition );
>>>>> itkExceptionMacro( << "Error occured in optimization" );
>>>>> }
>>>>>......
>>>>>/////////////////////////////////////////////////////////////////////////////////////////////
>>>>>_______________________________________________
>>>>>Insight-users mailing list
>>>>>Insight-users at itk.org
>>>>>http://www.itk.org/mailman/listinfo/insight-users
>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>
>>
>
>
> ------------------------------------------------------------------------
>
> /*=========================================================================
>
> Program: Insight Segmentation & Registration Toolkit
> Module: $RCSfile: DeformableRegistration6.cxx,v $
> Language: C++
> Date: $Date: 2004/08/05 15:04:02 $
> Version: $Revision: 1.3 $
>
> Copyright (c) Insight Software Consortium. All rights reserved.
> See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
>
> This software is distributed WITHOUT ANY WARRANTY; without even
> the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
> PURPOSE. See the above copyright notices for more information.
>
> =========================================================================*/
>
> // Software Guide : BeginLatex
> //
> // This example illustrates the use of the \doxygen{BSplineDeformableTransform}
> // class in a manually controlled multi-resolution scheme. Here we define two
> // transforms at two different resolution levels. A first registration is
> // perfomed with the spline grid of low resolution, and the results are then
> // used for initializing a higher resolution grid. Since this example is quite
> // similar to the previous example on the use of the
> // \code{BSplineDeformableTransform} we omit here most of the details already
> // discussed and will focus on the aspects related to the multi-resolution
> // approach.
> //
> // \index{itk::BSplineDeformableTransform}
> // \index{itk::BSplineDeformableTransform!DeformableRegistration}
> // \index{itk::LBFGSOptimizer}
> //
> //
> // Software Guide : EndLatex
>
> #include "itkImageRegistrationMethod.h"
> #include "itkMeanSquaresImageToImageMetric.h"
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkImage.h"
>
> // Software Guide : BeginLatex
> //
> // We include the header files for the transform and the optimizer.
> //
> // \index{itk::BSplineDeformableTransform!header}
> // \index{itk::LBFGSOptimizer!header}
> //
> // Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> #include "itkBSplineDeformableTransform.h"
> #include "itkLBFGSOptimizer.h"
> #include "itkLBFGSBOptimizer.h"
> // Software Guide : EndCodeSnippet
>
>
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
>
> #include "itkResampleImageFilter.h"
> #include "itkCastImageFilter.h"
> #include "itkSquaredDifferenceImageFilter.h"
>
> #include "itkBSplineResampleImageFunction.h"
> #include "itkIdentityTransform.h"
> #include "itkBSplineDecompositionImageFilter.h"
>
> // NOTE: the LBFGSOptimizer does not invoke events
>
> // The following section of code implements a Command observer
> // used to monitor the evolution of the registration process.
> //
> #include "itkCommand.h"
> class CommandIterationUpdate : public itk::Command
> {
> public:
> typedef CommandIterationUpdate Self;
> typedef itk::Command Superclass;
> typedef itk::SmartPointer<Self> Pointer;
> itkNewMacro( Self );
> protected:
> CommandIterationUpdate() {};
> public:
> typedef itk::LBFGSBOptimizer OptimizerType;
> typedef const OptimizerType * OptimizerPointer;
>
> void Execute(itk::Object *caller, const itk::EventObject & event)
> {
> Execute( (const itk::Object *)caller, event);
> }
>
> void Execute(const itk::Object * object, const itk::EventObject & event)
> {
> OptimizerPointer optimizer =
> dynamic_cast< OptimizerPointer >( object );
> if( typeid( event ) != typeid( itk::IterationEvent ) )
> {
> return;
> }
> std::cout << optimizer->GetCurrentIteration() << " ";
> std::cout << optimizer->GetValue() << " ";
> std::cout << optimizer->GetInfinityNormOfProjectedGradient() << std::endl;
> }
> };
>
>
> int main( int argc, char *argv[] )
> {
> if( argc < 4 )
> {
> std::cerr << "Missing Parameters " << std::endl;
> std::cerr << "Usage: " << argv[0];
> std::cerr << " fixedImageFile movingImageFile outputImagefile ";
> std::cerr << " [differenceOutputfile] [differenceBeforeRegistration] ";
> std::cerr << " [deformationField] ";
> return 1;
> }
>
> const unsigned int ImageDimension = 2;
> typedef float PixelType;
>
> typedef itk::Image< PixelType, ImageDimension > FixedImageType;
> typedef itk::Image< PixelType, ImageDimension > MovingImageType;
>
>
> // Software Guide : BeginLatex
> //
> // We instantiate now the type of the \code{BSplineDeformableTransform} using
> // as template parameters the type for coordinates representation, the
> // dimension of the space, and the order of the BSpline.
> //
> // \index{BSplineDeformableTransform|New}
> // \index{BSplineDeformableTransform|Instantiation}
> //
> // Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> const unsigned int SpaceDimension = ImageDimension;
> const unsigned int SplineOrder = 3;
> typedef double CoordinateRepType;
>
> typedef itk::BSplineDeformableTransform<
> CoordinateRepType,
> SpaceDimension,
> SplineOrder > TransformType;
> // Software Guide : EndCodeSnippet
>
>
> // typedef itk::LBFGSOptimizer OptimizerType;
> typedef itk::LBFGSBOptimizer OptimizerType;
>
> typedef itk::MeanSquaresImageToImageMetric<
> FixedImageType,
> MovingImageType > MetricType;
>
> typedef itk:: LinearInterpolateImageFunction<
> MovingImageType,
> double > InterpolatorType;
>
> typedef itk::ImageRegistrationMethod<
> FixedImageType,
> MovingImageType > RegistrationType;
>
> MetricType::Pointer metric = MetricType::New();
> OptimizerType::Pointer optimizer = OptimizerType::New();
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> RegistrationType::Pointer registration = RegistrationType::New();
>
>
> registration->SetMetric( metric );
> registration->SetOptimizer( optimizer );
> registration->SetInterpolator( interpolator );
>
>
> // Software Guide : BeginLatex
> //
> // We construct two transform objects, each one will be configured for a resolution level.
> // Notice than in this multi-resolution scheme we are not modifying the
> // resolution of the image, but rather the flexibility of the deformable
> // transform itself.
> //
> // \index{itk::RegistrationMethod!SetTransform()}
> //
> // Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> TransformType::Pointer transformLow = TransformType::New();
> registration->SetTransform( transformLow );
> // Software Guide : EndCodeSnippet
>
> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
>
> FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
> MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
>
> fixedImageReader->SetFileName( argv[1] );
> movingImageReader->SetFileName( argv[2] );
>
> FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
>
> registration->SetFixedImage( fixedImage );
> registration->SetMovingImage( movingImageReader->GetOutput() );
>
> fixedImageReader->Update();
>
> FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
>
> registration->SetFixedImageRegion( fixedRegion );
>
> // Software Guide : BeginLatex
> //
> // Here we define the parameters of the BSplineDeformableTransform grid. We
> // arbitrarily decide to use a grid with $5 \times 5$ nodes within the image.
> // The reader should note that the BSpline computation requires a
> // finite support region ( 1 grid node at the lower borders and 2
> // grid nodes at upper borders). Therefore in this example, we set
> // the grid size to be $8 \times 8$ and place the grid origin such that
> // grid node (1,1) coinicides with the first pixel in the fixed image.
> //
> // \index{BSplineDeformableTransform}
> //
> // Software Guide : EndLatex
>
>
> typedef TransformType::RegionType RegionType;
> RegionType bsplineRegionLow;
> RegionType::SizeType gridBorderSize;
> RegionType::SizeType totalGridSize;
>
> gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower, 2 upper )
>
> // Software Guide : BeginLatex
> //
> // Here we define the parameters of the BSpline transform at low resolution
> //
> // \index{BSplineDeformableTransform}
> //
> // Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> RegionType::SizeType gridLowSizeOnImage;
> gridLowSizeOnImage.Fill( 5 );
> totalGridSize = gridLowSizeOnImage + gridBorderSize;
>
> RegionType bsplineRegion;
> bsplineRegion.SetSize( totalGridSize );
>
> typedef TransformType::SpacingType SpacingType;
> SpacingType spacingLow = fixedImage->GetSpacing();
>
> typedef TransformType::OriginType OriginType;
> OriginType originLow = fixedImage->GetOrigin();;
>
> FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();
>
> for(unsigned int r=0; r<ImageDimension; r++)
> {
> spacingLow[r] *= floor( static_cast<double>(fixedImageSize[r] - 1) /
> static_cast<double>(gridLowSizeOnImage[r] - 1) );
> originLow[r] -= spacingLow[r];
> }
>
> transformLow->SetGridSpacing( spacingLow );
> transformLow->SetGridOrigin( originLow );
> transformLow->SetGridRegion( bsplineRegion );
>
> typedef TransformType::ParametersType ParametersType;
>
> const unsigned int numberOfParameters =
> transformLow->GetNumberOfParameters();
>
> ParametersType parametersLow( numberOfParameters );
>
> parametersLow.Fill( 0.0 );
>
> transformLow->SetParameters( parametersLow );
> // Software Guide : EndCodeSnippet
>
>
>
> // Software Guide : BeginLatex
> //
> // We now pass the parameters of the current transform as the initial
> // parameters to be used when the registration process starts.
> //
> // Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> registration->SetInitialTransformParameters( transformLow->GetParameters() );
>
> /*
> optimizer->SetGradientConvergenceTolerance( 0.05 );
> optimizer->SetLineSearchAccuracy( 0.9 );
> optimizer->SetDefaultStepLength( 1.5 );
> optimizer->TraceOn();
> optimizer->SetMaximumNumberOfFunctionEvaluations( 1000 );
> */
>
> OptimizerType::BoundSelectionType boundSelect( transformLow->GetNumberOfParameters() );
> OptimizerType::BoundValueType upperBound( transformLow->GetNumberOfParameters() );
> OptimizerType::BoundValueType lowerBound( transformLow->GetNumberOfParameters() );
>
> boundSelect.Fill( 0 );
> upperBound.Fill( 0.0 );
> lowerBound.Fill( 0.0 );
>
> optimizer->SetBoundSelection( boundSelect );
> optimizer->SetUpperBound( upperBound );
> optimizer->SetLowerBound( lowerBound );
>
> optimizer->SetCostFunctionConvergenceFactor( 1e+12 );
> optimizer->SetProjectedGradientTolerance( 1.0 );
> optimizer->SetMaximumNumberOfIterations( 500 );
> optimizer->SetMaximumNumberOfEvaluations( 500 );
> optimizer->SetMaximumNumberOfCorrections( 5 );
>
> // Create the Command observer and register it with the optimizer.
> //
> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
> optimizer->AddObserver( itk::IterationEvent(), observer );
>
> std::cout << "Starting Registration with low resolution transform" << std::endl;
>
> try
> {
> registration->StartRegistration();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return -1;
> }
> // Software Guide : EndCodeSnippet
>
>
> // Software Guide : BeginLatex
> //
> // Once the registration has finished with the low resolution grid, we
> // proceed to instantiate a higher resolution
> // \code{BSplineDeformableTransform}.
> //
> // Software Guide : EndLatex
>
> TransformType::Pointer transformHigh = TransformType::New();
>
> RegionType::SizeType gridHighSizeOnImage;
> gridHighSizeOnImage.Fill( 10 );
> totalGridSize = gridHighSizeOnImage + gridBorderSize;
>
> bsplineRegion.SetSize( totalGridSize );
>
> SpacingType spacingHigh = fixedImage->GetSpacing();
> OriginType originHigh = fixedImage->GetOrigin();;
>
> for(unsigned int rh=0; rh<ImageDimension; rh++)
> {
> spacingHigh[rh] *= floor( static_cast<double>(fixedImageSize[rh] - 1) /
> static_cast<double>(gridHighSizeOnImage[rh] - 1) );
> originHigh[rh] -= spacingHigh[rh];
> }
>
> transformHigh->SetGridSpacing( spacingHigh );
> transformHigh->SetGridOrigin( originHigh );
> transformHigh->SetGridRegion( bsplineRegion );
>
> ParametersType parametersHigh( transformHigh->GetNumberOfParameters() );
> parametersHigh.Fill( 0.0 );
>
> // Software Guide : BeginLatex
> //
> // Now we need to initialize the BSpline coefficients of the higher resolution
> // transform. This is done by first computing the actual deformation field
> // at the higher resolution from the lower resolution BSpline coefficients.
> // Then a BSpline decomposition is done to obtain the BSpline coefficient of
> // the higher resolution transform.
> //
> // Software Guide : EndLatex
>
> unsigned int counter = 0;
>
> for ( unsigned int k = 0; k < SpaceDimension; k++ )
> {
> typedef TransformType::ImageType ParametersImageType;
> typedef itk::ResampleImageFilter<ParametersImageType,ParametersImageType> ResamplerType;
> ResamplerType::Pointer upsampler = ResamplerType::New();
>
> typedef itk::BSplineResampleImageFunction<ParametersImageType,double> FunctionType;
> FunctionType::Pointer function = FunctionType::New();
>
> typedef itk::IdentityTransform<double,SpaceDimension> IdentityTransformType;
> IdentityTransformType::Pointer identity = IdentityTransformType::New();
>
> upsampler->SetInput( transformLow->GetCoefficientImage()[k] );
> upsampler->SetInterpolator( function );
> upsampler->SetTransform( identity );
> upsampler->SetSize( transformHigh->GetGridRegion().GetSize() );
> upsampler->SetOutputSpacing( transformHigh->GetGridSpacing() );
> upsampler->SetOutputOrigin( transformHigh->GetGridOrigin() );
>
> typedef itk::BSplineDecompositionImageFilter<ParametersImageType,ParametersImageType>
> DecompositionType;
> DecompositionType::Pointer decomposition = DecompositionType::New();
>
> decomposition->SetSplineOrder( SplineOrder );
> decomposition->SetInput( upsampler->GetOutput() );
> decomposition->Update();
>
> ParametersImageType::Pointer newCoefficients = decomposition->GetOutput();
>
> // copy the coefficients into the parameter array
> typedef itk::ImageRegionIterator<ParametersImageType> Iterator;
> Iterator it( newCoefficients, transformHigh->GetGridRegion() );
> while ( !it.IsAtEnd() )
> {
> parametersHigh[ counter++ ] = it.Get();
> ++it;
> }
>
> }
>
> transformHigh->SetParameters( parametersHigh );
>
> // Software Guide : BeginLatex
> //
> // We now pass the parameters of the high resolution transform as the initial
> // parameters to be used in a second stage of the registration process.
> //
> // Software Guide : EndLatex
>
> std::cout << "Starting Registration with high resolution transform" << std::endl;
>
> // Software Guide : BeginCodeSnippet
> registration->SetInitialTransformParameters( transformHigh->GetParameters() );
> registration->SetTransform( transformHigh );
>
> // Software Guide : BeginLatex
> //
> // Typically, we will also want to tighten the optimizer parameters
> // when we move from lower to higher resolution grid.
> //
> // Software Guide : EndLatex
>
> /* optimizer->SetGradientConvergenceTolerance( 0.01 );
> optimizer->SetDefaultStepLength( 0.25 );
> */
> optimizer->SetProjectedGradientTolerance( 0.01 );
>
> try
> {
> registration->StartRegistration();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return -1;
> }
> // Software Guide : EndCodeSnippet
>
>
>
> // Finally we use the last transform parameters in order to resample the image.
> //
> transformHigh->SetParameters( registration->GetLastTransformParameters() );
>
> typedef itk::ResampleImageFilter<
> MovingImageType,
> FixedImageType > ResampleFilterType;
>
> ResampleFilterType::Pointer resample = ResampleFilterType::New();
>
> resample->SetTransform( transformHigh );
> resample->SetInput( movingImageReader->GetOutput() );
>
> resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
> resample->SetOutputOrigin( fixedImage->GetOrigin() );
> resample->SetOutputSpacing( fixedImage->GetSpacing() );
> resample->SetDefaultPixelValue( 100 );
>
> typedef unsigned char OutputPixelType;
>
> typedef itk::Image< OutputPixelType, ImageDimension > OutputImageType;
>
> typedef itk::CastImageFilter<
> FixedImageType,
> OutputImageType > CastFilterType;
>
> typedef itk::ImageFileWriter< OutputImageType > WriterType;
>
>
> WriterType::Pointer writer = WriterType::New();
> CastFilterType::Pointer caster = CastFilterType::New();
>
>
> writer->SetFileName( argv[3] );
>
>
> caster->SetInput( resample->GetOutput() );
> writer->SetInput( caster->GetOutput() );
>
>
> try
> {
> writer->Update();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return -1;
> }
>
>
>
> typedef itk::SquaredDifferenceImageFilter<
> FixedImageType,
> FixedImageType,
> OutputImageType > DifferenceFilterType;
>
> DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
>
> WriterType::Pointer writer2 = WriterType::New();
> writer2->SetInput( difference->GetOutput() );
>
>
> // Compute the difference image between the
> // fixed and resampled moving image.
> if( argc >= 5 )
> {
> difference->SetInput1( fixedImageReader->GetOutput() );
> difference->SetInput2( resample->GetOutput() );
> writer2->SetFileName( argv[4] );
> try
> {
> writer2->Update();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return -1;
> }
> }
>
>
> // Compute the difference image between the
> // fixed and moving image before registration.
> if( argc >= 6 )
> {
> writer2->SetFileName( argv[5] );
> difference->SetInput1( fixedImageReader->GetOutput() );
> difference->SetInput2( movingImageReader->GetOutput() );
> try
> {
> writer2->Update();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return -1;
> }
> }
>
>
>
> // Generate the explicit deformation field resulting from
> // the registration.
>
> typedef itk::Vector< float, ImageDimension > VectorType;
> typedef itk::Image< VectorType, ImageDimension > DeformationFieldType;
>
> DeformationFieldType::Pointer field = DeformationFieldType::New();
> field->SetRegions( fixedRegion );
> field->SetOrigin( fixedImage->GetOrigin() );
> field->SetSpacing( fixedImage->GetSpacing() );
> field->Allocate();
>
> typedef itk::ImageRegionIterator< DeformationFieldType > FieldIterator;
> FieldIterator fi( field, fixedRegion );
>
> fi.GoToBegin();
>
> TransformType::InputPointType fixedPoint;
> TransformType::OutputPointType movingPoint;
> DeformationFieldType::IndexType index;
>
> VectorType displacement;
>
> while( ! fi.IsAtEnd() )
> {
> index = fi.GetIndex();
> field->TransformIndexToPhysicalPoint( index, fixedPoint );
> movingPoint = transformHigh->TransformPoint( fixedPoint );
> displacement[0] = movingPoint[0] - fixedPoint[0];
> displacement[1] = movingPoint[1] - fixedPoint[1];
> fi.Set( displacement );
> ++fi;
> }
>
>
>
> typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;
> FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
>
> fieldWriter->SetInput( field );
>
> if( argc >= 7 )
> {
> fieldWriter->SetFileName( argv[6] );
> try
> {
> fieldWriter->Update();
> }
> catch( itk::ExceptionObject & excp )
> {
> std::cerr << "Exception thrown " << std::endl;
> std::cerr << excp << std::endl;
> return EXIT_FAILURE;
> }
> }
>
> return 0;
> }
>
More information about the Insight-users
mailing list