[ITK Community] [Insight-users] Bspline artefact multi resolution - origin offset - Bug?

Emma Saunders emmasaunders123 at gmail.com
Wed Dec 11 17:26:16 EST 2013


Hi everyone

I appear to be getting a step artefact at the boundaries of my registered
image, using a multi-resolution Bspline approach.  The distance of the
artefact co-indices with the origin of the original image.  I can't seem to
find why this occurs.  I set the Bspline grid wrt the original fixed image
so unsure why this is happening.

I would really appreciate any advice:

Below is the code if anyone is interested:

 Many Thanks

Emma






#include "itkMedianImageFilter.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkImageRegistrationMethod.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkSpatialObjectToImageFilter.h"
#include "itkEllipseSpatialObject.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkAffineTransform.h"
#include "itkBSplineDeformableTransform.h" //version 3
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkMemoryProbesCollectorBase.h"
#include <itkHistogramMatchingImageFilter.h>
#include "itkBSplineResampleImageFunction.h"
#include "itkLBFGSBOptimizer.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkLBFGSOptimizer.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"
#include "itkBSplineResampleImageFunction.h"
#include "itkIdentityTransform.h"
#include "itkBSplineDecompositionImageFilter.h"
#include "itkCommand.h"
#include <fstream>
#include "itkCommand.h"

unsigned int Counter = 0;
unsigned int Countold = 0;


//This is for the affine registration
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::RegularStepGradientDescentOptimizer 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( ! itk::IterationEvent().CheckEvent( &event ) )
      {
      return;
      }
     std::cout << optimizer->GetValue() << "   " <<std::endl;

    }
};


//This is for B spline
   class BCommandIterationUpdate : public itk::Command
  {
  public:
    typedef  BCommandIterationUpdate                     Self;
    typedef  itk::Command                               Superclass;
    typedef  itk::SmartPointer<BCommandIterationUpdate>  Pointer;
    itkNewMacro( BCommandIterationUpdate );
  protected:
    BCommandIterationUpdate() {};

  public:

  typedef itk::LBFGSBOptimizer    OptimizerType;
  typedef   const OptimizerType *                  OptimizerPointer;


  public:

    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)
      {

const OptimizerPointer  optimizer =
                      dynamic_cast< OptimizerPointer >( object );


        if( !(itk::IterationEvent().CheckEvent( &event )) )
          {
          return;
          }

        Countold=Counter;
        Counter = optimizer->GetCurrentIteration();

        if (Counter!=Countold)
        {

        std::cout  <<  optimizer->GetValue() << " " <<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 " << std::endl;
    std::cerr << "   outputImagefile" << std::endl;
     std::cerr << "   outputDeformationField " << std::endl;

    return EXIT_FAILURE;

    }

 // READ AND CREATE THE IMAGES AND DISPLACEMENT FIELD

  const    unsigned int    Dimension = 3;
  typedef float        PixelType;
  typedef itk::Vector< double, 3 >           VectorPixelType;

  typedef itk::Image< PixelType, Dimension >  FixedImageType;
  typedef itk::Image< PixelType, Dimension >  MovingImageType;
  typedef itk::Image< PixelType, Dimension >  OutputImageType;
  typedef itk::Image<  VectorPixelType, 3 > DisplacementFieldType;

  typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
  typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
  typedef itk::ImageFileWriter< OutputImageType >  ImageWriterType;
  typedef itk::ImageFileWriter< DisplacementFieldType >  FieldWriterType;

  typedef itk::ImageFileWriter< DisplacementFieldType >  FieldWriterType;
  FieldWriterType::Pointer FieldWriter = FieldWriterType::New();

  FixedImageReaderType::Pointer  FixedImageReader  =
FixedImageReaderType::New();
  MovingImageReaderType::Pointer MovingImageReader =
MovingImageReaderType::New();
  ImageWriterType::Pointer ImageWriter = ImageWriterType::New();
  FixedImageReader->SetFileName(  argv[1] );
  MovingImageReader->SetFileName( argv[2] );
  ImageWriter->SetFileName( argv[3] );
  FieldWriter->SetFileName( argv[4] );



  FixedImageType::Pointer FixedImage = FixedImageReader->GetOutput();
  MovingImageType::Pointer MovedImage = MovingImageReader->GetOutput();
  OutputImageType::Pointer OutImage = OutputImageType::New();
  DisplacementFieldType::Pointer DeformField = DisplacementFieldType::New();
  std::cout << "First Fixed Image origin is " << FixedImage->GetOrigin() <<
std::endl ;

  FixedImage->Update();
  MovedImage->Update();




  //Median and Gaussian Filter

 typedef itk::MedianImageFilter<FixedImageType, FixedImageType > FilterType;


  // Create and setup a median filter
  FilterType::Pointer medianFilterFixed = FilterType::New();
  FilterType::InputSizeType radius;
  radius.Fill(1);
  medianFilterFixed->SetRadius(radius);
  medianFilterFixed->SetInput( FixedImage );
  medianFilterFixed->Update();

   FilterType::Pointer medianFilterMoved = FilterType::New();
  medianFilterMoved->SetRadius(radius);
  medianFilterMoved->SetInput( MovedImage );
  medianFilterMoved->Update();

  typedef itk::DiscreteGaussianImageFilter<
                                      FixedImageType,
                                      FixedImageType
                                                    > GaussianFilterType;


  GaussianFilterType::Pointer fixedSmoother  = GaussianFilterType::New();
  GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();

  fixedSmoother->SetVariance(1.0 );
  movingSmoother->SetVariance(1.0 );

  fixedSmoother->SetInput(  medianFilterFixed->GetOutput() );
  movingSmoother->SetInput ( medianFilterMoved->GetOutput() );

  fixedSmoother->Update();
  movingSmoother->Update();

   std::cout<< "The fixed origin is " <<
fixedSmoother->GetOutput()->GetOrigin() <<std::endl;
   std::cout << "Fixed Image origin is " << FixedImage->GetOrigin() <<
std::endl ;


//Define the transformations to be used

//PERFROM AN INITIAL AFFINE

  typedef itk::AffineTransform<
                                  double,
                                  Dimension  >     AffineTransformType;


//This will be followed by Bspline


  const unsigned int SpaceDimension = Dimension;
  const unsigned int SplineOrder = 3;
  typedef double CoordinateRepType;

  typedef itk::BSplineDeformableTransform<
                            CoordinateRepType,
                            SpaceDimension,
                            SplineOrder >     DeformableTransformType;
//version 3


 //Now for the Optimizer, Metric and Interpolatror  and Registration

 //which are connected to the registration

 //typedef itk::LBFGSOptimizer       bOptimizerType;
  typedef itk::LBFGSBOptimizer       bOptimizerType;
   bOptimizerType::Pointer      boptimizer     = bOptimizerType::New();

 typedef itk::RegularStepGradientDescentOptimizer       OptimizerType;

  typedef itk::MattesMutualInformationImageToImageMetric<
                                    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();

  // Setup the metric parameters

  metric->ReinitializeSeed( 76926294 );
  metric->SetUseCachingOfBSplineWeights(1);

    const unsigned long numberOfImagePixels =
    FixedImage->GetLargestPossibleRegion().GetNumberOfPixels();

      const unsigned long numberOfSpatialSamplesAffine =
    static_cast<unsigned long>( numberOfImagePixels * 1 );

  metric->SetNumberOfSpatialSamples( numberOfSpatialSamplesAffine );
  metric->SetNumberOfHistogramBins( 40 ); // Reduce the number of bins if a
deformable registration fails. If the number of bins is too large, the
estimated PDFs will be a field of impulses and will inhibit reliable
registration estimation.

  registration->SetMetric(        metric        );
  registration->SetOptimizer(     optimizer     );
  registration->SetInterpolator(  interpolator  );

//create the Affine transform

  AffineTransformType::Pointer  Affinetransform =
AffineTransformType::New();
  registration->SetTransform( Affinetransform );

//input the images

  registration->SetFixedImage(fixedSmoother->GetOutput());
  registration->SetMovingImage(movingSmoother->GetOutput());

//input the region

  registration->SetFixedImageRegion(FixedImage->GetBufferedRegion() );

//set the initial parameters
  typedef RegistrationType::ParametersType AffineParametersType;
  AffineParametersType initialParameters(
Affinetransform->GetNumberOfParameters() );

  Affinetransform->SetIdentity();
  registration->SetInitialTransformParameters(
                             Affinetransform->GetParameters() );


//set scale values for optimizer to search sensibly and set up optimizer

  typedef OptimizerType::ScalesType       OptimizerScalesType;
  OptimizerScalesType optimizerScales(
Affinetransform->GetNumberOfParameters() );

  FixedImageType:: SpacingType fixedspacing = FixedImage->GetSpacing();
  FixedImageType::RegionType region=FixedImage->GetLargestPossibleRegion();
  const unsigned int numberOfPixels = region.GetNumberOfPixels();
  FixedImageType:: SizeType fixedsize =region.GetSize();


  double translationScalea = 1.0 /( 10.0 * fixedspacing[0]  *  fixedsize[0]
 );
  double translationScaleb = 1.0 /( 10.0 * fixedspacing[1]  *  fixedsize[1]
 );
  double translationScalec = 1.0 /( 10.0 * fixedspacing[2]  *  fixedsize[2]
 );


//for 3D use below
  optimizerScales[0] =  1.0;
  optimizerScales[1] =  1.0;
  optimizerScales[2] =  1.0;
  optimizerScales[3] =  1.0;
  optimizerScales[4] =  1.0;
  optimizerScales[5] =  1.0;
  optimizerScales[6] =  1.0;
  optimizerScales[7] =  1.0;
  optimizerScales[8] =  1.0;
  optimizerScales[9]  =  translationScalea;
  optimizerScales[10] =  translationScaleb;
  optimizerScales[11] =  translationScalec;

  //step length was 0.1 larger step length quicker converge but may get
erreaic values in the metric as it decreases
  unsigned int maxNumberOfIterations = 60;
  optimizer->SetMaximumStepLength( 0.05 );
  optimizer->SetMinimumStepLength( 0.001 ); //was 0.00001
  optimizer->SetNumberOfIterations( maxNumberOfIterations );
  optimizer->MaximizeOff();


//Perform Affine Registration put observer on the Affine optimizer

  CommandIterationUpdate::Pointer observera = CommandIterationUpdate::New();
  optimizer->AddObserver( itk::IterationEvent(), observera );

  try
    {
    registration->Update();
    std::cout << " " <<std::endl;
    std::cout << "Optimizer stop condition: "
           << registration->GetOptimizer()->GetStopConditionDescription()
           << std::endl;
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
    }

  std::cout << "Affine Registration completed" << std::endl;
  std::cout << std::endl;
  optimizer->RemoveAllObservers();
  optimizer->RemoveAllObservers();

//Obtain the last Affine transform which will be used as a bulk transform
in the Bspline

  Affinetransform->SetParameters(registration->GetLastTransformParameters()
);


//////////////////////////////////////////////////////////////////////////////////////////////
///NOW FOR THE Bspline Registration
//  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.

//SO must define a transform and a grid for each level

  DeformableTransformType::Pointer  BtransformLow =
DeformableTransformType::New();
  typedef DeformableTransformType::RegionType RegionType;
  RegionType bsplineRegionlow;
  RegionType::SizeType   gridBorderSize;
  RegionType::SizeType   gridSizeOnImagelow;
  RegionType::SizeType   totalGridSizelow;


//.//////////////////////////////////////////////////////////////////////
///GRID SIZE FOR LOW RESOLUTION
/////////////////////////////////////////////////////////////////////////

  gridSizeOnImagelow[0]=12;  //lat coronal, sup-inf //want 5 2 7 10 4 15
  gridSizeOnImagelow[1]=4;    // want
  gridSizeOnImagelow[2]=12;
  gridBorderSize.Fill( 3 );    // Border for spline order = 3 ( 1 lower, 2
upper )
  totalGridSizelow = gridSizeOnImagelow + gridBorderSize;
  bsplineRegionlow.SetSize( totalGridSizelow );
  typedef DeformableTransformType::SpacingType SpacingTypelow;
  SpacingTypelow spacinglow = FixedImage->GetSpacing();
  typedef DeformableTransformType::OriginType OriginTypelow;
  OriginTypelow originlow = FixedImage->GetOrigin();
  FixedImageType::RegionType FixedRegionlow =
FixedImage->GetBufferedRegion();
  FixedImageType::SpacingType   fixedSpacinglow    =
FixedImage->GetSpacing();
  FixedImageType::SizeType FixedImageSizelow = FixedRegionlow.GetSize();

  for(unsigned int r=0; r<Dimension; r++)
    {

    spacinglow[r] = floor( fixedSpacinglow[r] * (FixedImageSizelow[r] - 1)
/ (gridSizeOnImagelow[r] - 1) ); //new way

    }

  FixedImageType::DirectionType gridDirection = FixedImage->GetDirection();
  SpacingTypelow gridOriginOffsetlow = gridDirection * spacinglow;
  OriginTypelow gridOriginlow = originlow - gridOriginOffsetlow;
  BtransformLow->SetGridSpacing( spacinglow );
  BtransformLow->SetGridOrigin( gridOriginlow );
  BtransformLow->SetGridRegion( bsplineRegionlow ); //i.e region of gridsize
  BtransformLow->SetGridDirection( gridDirection );
  std::cout << "grid region for low Bspline is " <<
BtransformLow->GetGridRegion() <<std::endl;


//////////////////////////////////////////////////////////////////////////////
///GRID FOR MEDIUM
////////////////////////////////////////////////////////////////////////////////

  DeformableTransformType::Pointer  BtransformMed=
DeformableTransformType::New();
  RegionType bsplineRegionMed;
  RegionType::SizeType   gridSizeOnImageMed;
  RegionType::SizeType   totalGridSizeMed;

//.//////////////////////////////////////////////////////////////////////
///GRID SIZE FOR MEDIUM RESOLUTION
/////////////////////////////////////////////////////////////////////////
  gridSizeOnImageMed[0]=24;  //lat coronal, sup-inf
  gridSizeOnImageMed[1]=8;    // want
  gridSizeOnImageMed[2]=24;
  totalGridSizeMed = gridSizeOnImageMed + gridBorderSize;
  bsplineRegionMed.SetSize( totalGridSizeMed );
 typedef DeformableTransformType::SpacingType SpacingTypeMed;
  SpacingTypeMed spacingMed = FixedImage->GetSpacing();
  typedef DeformableTransformType::OriginType OriginTypeMed;
  OriginTypeMed originMed = FixedImage->GetOrigin();
  FixedImageType::RegionType FixedRegionMed =
FixedImage->GetBufferedRegion();
  FixedImageType::SpacingType   fixedSpacingMed    =
FixedImage->GetSpacing();
  FixedImageType::SizeType FixedImageSizeMed = FixedRegionMed.GetSize();

      for(unsigned int r=0; r<Dimension; r++)
    {

    spacingMed[r] = floor( fixedSpacingMed[r] * (FixedImageSizeMed[r] - 1)
/ (gridSizeOnImageMed[r] - 1) ); //new way

    }

   SpacingTypeMed gridOriginOffsetMed = gridDirection * spacingMed;
   OriginTypeMed gridOriginMed = originMed - gridOriginOffsetMed;

  BtransformMed->SetGridSpacing( spacingMed );
  BtransformMed->SetGridOrigin( gridOriginMed );
  BtransformMed->SetGridRegion( bsplineRegionMed ); //i.e region of gridsize
  BtransformMed->SetGridDirection( gridDirection );
  std::cout << "grid region for Med Bspline is " <<
BtransformMed->GetGridRegion() <<std::endl;

//////////////////////////////////////////////////////////////////////////////
///GRID FOR HIGH
////////////////////////////////////////////////////////////////////////////////

//Now the Higher Transform

  DeformableTransformType::Pointer  BtransformHigh =
DeformableTransformType::New();
  RegionType bsplineRegionhigh;
  RegionType::SizeType   gridSizeOnImagehigh;
  RegionType::SizeType   totalGridSizehigh;

//.//////////////////////////////////////////////////////////////////////
///GRID SIZE FOR HIGH RESOLUTION
/////////////////////////////////////////////////////////////////////////
  gridSizeOnImagehigh[0]=50;  //lat coronal, sup-inf
  gridSizeOnImagehigh[1]=16;    // want
  gridSizeOnImagehigh[2]=50;
  totalGridSizehigh = gridSizeOnImagehigh + gridBorderSize;
  bsplineRegionhigh.SetSize( totalGridSizehigh );
 typedef DeformableTransformType::SpacingType SpacingTypehigh;
  SpacingTypehigh spacinghigh = FixedImage->GetSpacing();
  typedef DeformableTransformType::OriginType OriginTypehigh;
  OriginTypehigh originhigh = FixedImage->GetOrigin();
  FixedImageType::RegionType FixedRegionhigh =
FixedImage->GetBufferedRegion();
  FixedImageType::SpacingType   fixedSpacinghigh    =
FixedImage->GetSpacing();
  FixedImageType::SizeType FixedImageSizehigh = FixedRegionhigh.GetSize();

      for(unsigned int r=0; r<Dimension; r++)
    {

    spacinghigh[r] = floor( fixedSpacinghigh[r] * (FixedImageSizehigh[r] -
1) / (gridSizeOnImagehigh[r] - 1) ); //new way

    }

   SpacingTypehigh gridOriginOffsethigh = gridDirection * spacinghigh;
   OriginTypehigh gridOriginhigh = originhigh - gridOriginOffsethigh;

  BtransformHigh->SetGridSpacing( spacinghigh );
  BtransformHigh->SetGridOrigin( gridOriginhigh );
  BtransformHigh->SetGridRegion( bsplineRegionhigh ); //i.e region of
gridsize
  BtransformHigh->SetGridDirection( gridDirection );

  std::cout << "grid region for high Bspline is " <<
BtransformHigh->GetGridRegion() <<std::endl;

///////////////////////////////////////////////////////////////////////
///METRIC NUMBER OF SAMPLES
//////////////////////////////////////////////////////////////////////


//Low
  const unsigned int numberOfBSplineParametersLow =
               BtransformLow->GetNumberOfParameters();

  typedef DeformableTransformType::ParametersType     ParametersTypeLow;
  ParametersTypeLow parametersLow( numberOfBSplineParametersLow );

//Medium

  const unsigned int numberOfBSplineParametersMed =
               BtransformMed->GetNumberOfParameters();

  typedef DeformableTransformType::ParametersType     ParametersTypeMed;
  ParametersTypeMed parametersMed( numberOfBSplineParametersMed );

  const unsigned int numberOfBSplineParametersHigh =
               BtransformHigh->GetNumberOfParameters();

  typedef DeformableTransformType::ParametersType     ParametersTypeHigh;
  ParametersTypeHigh parametersHigh( numberOfBSplineParametersHigh );


     const unsigned long numberOfSamplesBsplineHigh =
       static_cast<unsigned long>(
        vcl_sqrt( static_cast<double>( numberOfBSplineParametersHigh ) *
                static_cast<double>( numberOfPixels ) ) );



///////////////////////////////////////////////////////////////////////
///LOW PARAMETERS
//////////////////////////////////////////////////////////////////////

  typedef DeformableTransformType::ParametersType     ParametersTypelow;

  const unsigned int numberOfBSplineParameterslow =
               BtransformLow->GetNumberOfParameters();

  ParametersTypelow parameterslow( numberOfBSplineParameterslow );
  parameterslow.Fill( 0.0 );
  BtransformLow->SetBulkTransform( Affinetransform );

 //////////////////////////////////////////////////////////////////////////////
///OPTIMIZER FOR LOW RESOLUTION
////////////////////////////////////////////////////////////////////////////////

 // if get IFLAG= -1. LINE SEARCH FAILED
 // increase step lentgth but keep other parameters the same

/////////////////////////////////////////////////////////////////////
// SEE HERE FOR OPTIMIZER ADVICE
//http://www.itk.org/pipermail/insight-users/2011-March/040286.html

///FOR the LBFGSBOptimizer got LOW

  bOptimizerType::BoundSelectionType boundSelectlow(
BtransformLow->GetNumberOfParameters() );
  bOptimizerType::BoundValueType upperBoundlow(
BtransformLow->GetNumberOfParameters() );
  bOptimizerType::BoundValueType lowerBoundlow(
BtransformLow->GetNumberOfParameters() );
  boundSelectlow.Fill( 0 );
  upperBoundlow.Fill( 0.0 );
  lowerBoundlow.Fill( 0.0 );
  boptimizer->SetBoundSelection( boundSelectlow );
  boptimizer->SetUpperBound( upperBoundlow );
  boptimizer->SetLowerBound( lowerBoundlow );
  boptimizer->SetCostFunctionConvergenceFactor( 1.e1 );
  boptimizer->SetProjectedGradientTolerance( 1e-6 );
  boptimizer->SetMaximumNumberOfIterations( 200 ); //was 200
  boptimizer->SetMaximumNumberOfEvaluations( 300 ); //
  boptimizer->SetMaximumNumberOfCorrections( 5 );
///////////////////////////////////////////////////////////////

  metric->SetNumberOfSpatialSamples( numberOfPixels*1 ); //
  RegistrationType::Pointer   bregistration  = RegistrationType::New();
  bregistration->SetNumberOfThreads(1);
  bregistration->SetInitialTransformParameters(
BtransformLow->GetParameters() );
  bregistration->SetTransform(BtransformLow);
  bregistration->SetFixedImageRegion(FixedImage->GetBufferedRegion() );
  bregistration->SetMetric(        metric        );
  bregistration->SetOptimizer(     boptimizer     );
  bregistration->SetInterpolator(  interpolator  );
  bregistration->SetFixedImage(fixedSmoother->GetOutput());
  bregistration->SetMovingImage(movingSmoother->GetOutput());
  bregistration->SetNumberOfThreads(1);


//////////////////////////////////////////////////////////////////////
///REGISTRATION FOR LOW
///////////////////////////////////////////////////////////////////////

  BCommandIterationUpdate::Pointer observerl =
BCommandIterationUpdate::New();
  boptimizer->AddObserver( itk::IterationEvent(), observerl );
  try
    {

    bregistration->Update();
    std::cout << " " <<std::endl;
    std::cout << "Optimizer stop condition: "
           << bregistration->GetOptimizer()->GetStopConditionDescription()
           << std::endl;

    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
    }

 boptimizer->RemoveAllObservers();



/////////////////////////////////////////////////////////////////////
 ///MEDIUM PARAMETERS
////////////////////////////////////////////////////////////////////

  typedef DeformableTransformType::ParametersType     ParametersTypeMed;
  parametersMed.Fill( 0.0 );

//  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.

  unsigned int counterMed = 0;

  for ( unsigned int k = 0; k < SpaceDimension; k++ )
    {
    typedef DeformableTransformType::ImageType ParametersImageTypeMed;
    typedef
itk::ResampleImageFilter<ParametersImageTypeMed,ParametersImageTypeMed>
ResamplerTypeMed;
    ResamplerTypeMed::Pointer upsamplerMed = ResamplerTypeMed::New();

    typedef
itk::BSplineResampleImageFunction<ParametersImageTypeMed,double>
FunctionTypeMed;
    FunctionTypeMed::Pointer functionMed = FunctionTypeMed::New();

    typedef itk::IdentityTransform<double,SpaceDimension>
IdentityTransformTypeMed;
    IdentityTransformTypeMed::Pointer identityMed =
IdentityTransformTypeMed::New();

    upsamplerMed->SetInput( BtransformLow->GetCoefficientImages()[k] );
    upsamplerMed->SetInterpolator( functionMed );
    upsamplerMed->SetTransform( identityMed );
    upsamplerMed->SetSize( BtransformMed->GetGridRegion().GetSize() );
    upsamplerMed->SetOutputSpacing( BtransformMed->GetGridSpacing() );
    upsamplerMed->SetOutputOrigin( BtransformMed->GetGridOrigin() );
    upsamplerMed->SetOutputDirection( FixedImage->GetDirection() );
    upsamplerMed->Update();

    typedef
itk::BSplineDecompositionImageFilter<ParametersImageTypeMed,ParametersImageTypeMed>
      DecompositionTypeMed;
    DecompositionTypeMed::Pointer decompositionMed =
DecompositionTypeMed::New();

    decompositionMed->SetSplineOrder( SplineOrder );
    decompositionMed->SetInput( upsamplerMed->GetOutput() );
    decompositionMed->Update();

    ParametersImageTypeMed::Pointer newCoefficientsMed =
decompositionMed->GetOutput();

    // copy the coefficients into the parameter array
    typedef itk::ImageRegionIterator<ParametersImageTypeMed> IteratorMed;
    IteratorMed it( newCoefficientsMed, BtransformMed->GetGridRegion() );
    while ( !it.IsAtEnd() )
      {
      parametersMed[ counterMed++ ] = it.Get();
      ++it;
      }

    }

  BtransformMed->SetParameters( parametersMed );

////////////////////////////////////////////////////////////////

///FOR the LBFGSBOptimizer MEDIUM

  bOptimizerType::BoundSelectionType boundSelectMed(
BtransformMed->GetNumberOfParameters() );
  bOptimizerType::BoundValueType upperBoundMed(
BtransformMed->GetNumberOfParameters() );
  bOptimizerType::BoundValueType lowerBoundMed(
BtransformMed->GetNumberOfParameters() );
  boundSelectMed.Fill( 0 );
  upperBoundMed.Fill( 0.0 );
  lowerBoundMed.Fill( 0.0 );
  boptimizer->SetBoundSelection( boundSelectMed );
  boptimizer->SetUpperBound( upperBoundMed );
  boptimizer->SetLowerBound( lowerBoundMed );
  boptimizer->SetCostFunctionConvergenceFactor( 1.e12 ); //was 1.e7
 //set/Get the CostFunctionConvergenceFactor. Algorithm terminates
 //  when the reduction in cost function is less than factor * epsmcj
 // where epsmch is the machine precision.
 // Typical values for factor: 1e+12 for low accuracy;
 // 1e+7 for moderate accuracy and 1e+1 for extremely high accuracy.
  boptimizer->SetProjectedGradientTolerance( 1e-6 );
  //Gradient tolerance determines a lower threshold which cuases the
optimization to stop if its
//value reached, recommended values for gradient tolerance are 1e-6 to
1e-10, higher values might result in premature elimination of the
optimization process
  boptimizer->SetMaximumNumberOfIterations( 200); //was 200
  boptimizer->SetMaximumNumberOfEvaluations( 300 );
  boptimizer->SetMaximumNumberOfCorrections( 5 );

  //M Number of correction- number of function/gradient pairs used to build
 model
  //This parameter is passed to the function which is used to create
optimizer object.
  // Increase of number of corrections decreases number of function
evaluations
  //and iterations needed to converge, but increases computational overhead
associated with iteration.
  // On well-conditioned problems can be as small as 3-10.
  //If function changes rapidly in some directions and slowly in other
ones,
  //then you can try increasing M in order to increase convergenc

///////////////////////////////////////////////////////////////////////
///REGISTRATION FOR MEDIUM
/////////////////////////////////////////////////////////////////////
  //For time probe

   metric->SetNumberOfSpatialSamples( numberOfPixels*1 ); //
  bregistration->SetInitialTransformParameters(
BtransformMed->GetParameters() );
  bregistration->SetTransform( BtransformMed );

  itk::TimeProbesCollectorBase chronometerm;
  itk::MemoryProbesCollectorBase memorymeterm;

  BCommandIterationUpdate::Pointer observerm =
BCommandIterationUpdate::New();
  boptimizer->AddObserver( itk::IterationEvent(), observerm );


  try
    {

    bregistration->Update();


    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
    }

  boptimizer->RemoveAllObservers();

 BtransformMed->SetParameters( bregistration->GetLastTransformParameters()
);

//NOW HIGH RESOLUTION
/////////////////////////////////////////////////////////////////////
 ///HIGH PARAMETERS
////////////////////////////////////////////////////////////////////

  parametersHigh.Fill( 0.0 );

  unsigned int counter = 0;

  for ( unsigned int k = 0; k < SpaceDimension; k++ )
    {
    typedef DeformableTransformType::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( BtransformMed->GetCoefficientImages()[k] );
    upsampler->SetInterpolator( function );
    upsampler->SetTransform( identity );
    upsampler->SetSize( BtransformHigh->GetGridRegion().GetSize() );
    upsampler->SetOutputSpacing( BtransformHigh->GetGridSpacing() );
    upsampler->SetOutputOrigin( BtransformHigh->GetGridOrigin() );
    upsampler->SetOutputDirection( FixedImage->GetDirection() );
    upsampler->Update();


    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, BtransformHigh->GetGridRegion() );
    while ( !it.IsAtEnd() )
      {
      parametersHigh[ counter++ ] = it.Get();
      ++it;
      }

    }

  BtransformHigh->SetParameters( parametersHigh );


/////////////////////////////////////////////////////////////////////////////////////////

///FOR the LBFGSBOptimizer HIGH

  bOptimizerType::BoundSelectionType boundSelecthigh(
BtransformHigh->GetNumberOfParameters() );
  bOptimizerType::BoundValueType upperBoundhigh(
BtransformHigh->GetNumberOfParameters() );
  bOptimizerType::BoundValueType lowerBoundhigh(
BtransformHigh->GetNumberOfParameters() );
  boundSelecthigh.Fill( 0 );
  upperBoundhigh.Fill( 0.0 );
  lowerBoundhigh.Fill( 0.0 );
  boptimizer->SetBoundSelection( boundSelecthigh );
  boptimizer->SetUpperBound( upperBoundhigh );
  boptimizer->SetLowerBound( lowerBoundhigh );
  boptimizer->SetCostFunctionConvergenceFactor( 1.e12 );
  boptimizer->SetProjectedGradientTolerance( 1e-6 );//was 1e-6
  boptimizer->SetMaximumNumberOfIterations( 1000 ); //was 200
  boptimizer->SetMaximumNumberOfEvaluations( 300 );
  boptimizer->SetMaximumNumberOfCorrections( 5 );



///////////////////////////////////////////////////////////////////////
///REGISTRATION FOR HIGH
/////////////////////////////////////////////////////////////////////
  //For time probe

   metric->SetNumberOfSpatialSamples( numberOfPixels*1 ); //
  bregistration->SetInitialTransformParameters(
BtransformHigh->GetParameters() );
  bregistration->SetTransform( BtransformHigh );

  itk::TimeProbesCollectorBase chronometerh;
  itk::MemoryProbesCollectorBase memorymeterh;

  BCommandIterationUpdate::Pointer observerh =
BCommandIterationUpdate::New();
  boptimizer->AddObserver( itk::IterationEvent(), observerh );

 try
    {

    bregistration->Update();

   std::cout << " " <<std::endl;
    std::cout << "Optimizer stop condition: "
           << bregistration->GetOptimizer()->GetStopConditionDescription()
           << std::endl;
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
    }

  // Finally we use the last transform parameters in order to resample the
image.
  //

  boptimizer->RemoveAllObservers();

 BtransformHigh->SetParameters( bregistration->GetLastTransformParameters()
);


  typedef itk::ResampleImageFilter<
                            MovingImageType,
                            FixedImageType >    ResampleFilterType;

  ResampleFilterType::Pointer resample = ResampleFilterType::New();


  resample->SetTransform( BtransformHigh );
  resample->SetInput( MovedImage );
  resample->SetSize(    FixedImage->GetLargestPossibleRegion().GetSize() );
  resample->SetOutputOrigin(  FixedImage->GetOrigin() );
  resample->SetOutputSpacing( FixedImage->GetSpacing() );
  resample->SetOutputDirection( FixedImage->GetDirection() );
  resample->SetDefaultPixelValue( 0 );
  resample->Update();

  ImageWriter->SetInput(resample->GetOutput());

    try
      {
      ImageWriter->Update();
      }
    catch( itk::ExceptionObject & excp )
      {
      std::cerr << "Exception thrown " << std::endl;
      std::cerr << excp << std::endl;
      return EXIT_FAILURE;
      }


  FixedImageType::RegionType FixedRegion = FixedImage->GetBufferedRegion();
  DeformField->SetRegions( FixedRegion );
  DeformField->SetOrigin( FixedImage->GetOrigin() );
  DeformField->SetSpacing( FixedImage->GetSpacing() );
  FixedImageType:: DirectionType Fixeddirection =FixedImage->GetDirection();
  DeformField->SetDirection(Fixeddirection);
  DeformField->Allocate();
  DeformField->Update();

//WRITE THE DEFORMATION FIELD

    typedef itk::ImageRegionIterator< DisplacementFieldType > FieldIterator;
    FieldIterator fi( DeformField, FixedRegion );

    fi.GoToBegin();

    DeformableTransformType::InputPointType  fixedPointa;
    DeformableTransformType::OutputPointType movingPointa;
    DisplacementFieldType::IndexType indexa;

    VectorPixelType displacement;

    while( ! fi.IsAtEnd() )
      {
      indexa = fi.GetIndex();
      DeformField->TransformIndexToPhysicalPoint( indexa, fixedPointa );
      movingPointa = BtransformHigh->TransformPoint( fixedPointa );
      displacement = movingPointa - fixedPointa;
      fi.Set( displacement );
      ++fi;
      }

    FieldWriter->SetInput( DeformField );
    //std::cout << "Writing deformation field ...";

    try
      {
      FieldWriter->Update();
      }
    catch( itk::ExceptionObject & excp )
      {
      std::cerr << "Exception thrown " << std::endl;
      std::cerr << excp << std::endl;
      return EXIT_FAILURE;
      }


  return EXIT_SUCCESS;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20131211/d0d69676/attachment.html>
-------------- next part --------------
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://www.itk.org/mailman/listinfo/insight-users


More information about the Community mailing list