[Insight-users] Setting B-spline coefficients (re-sent and help
needed...)
Xiang, Hong
HXIANG at LROC.HARVARD.EDU
Tue Jul 11 09:47:07 EDT 2006
7. Setting B-spline coefficients (as prompted in the
OutputWindow) (Xiang, Hong)
------------------------------
Message: 7
Date: Fri, 7 Jul 2006 18:49:57 -0400
From: "Xiang, Hong" <HXIANG at LROC.HARVARD.EDU>
Subject: Setting B-spline coefficients (as prompted in the OutputWindow)
To: <insight-users at itk.org>
Cc: Karthik.Krishnan at kitware.com, luis.ibanez at kitware.com, "Xiang,
Hong" <HXIANG at LROC.HARVARD.EDU>
Hello ITK users,
While running a modified code from DeformableRegistration8.cxx, I ran
into a warning OutputWindow with a rolling message of
"...\code\common\itkBSplineDeformableTransform.txx, line 5... B-spline
coefficients have not been set".
It appears that it did not afftect the results of the registration (as
confirmed from the output files -- checkerboard displays and the
deformation field), but it caused the job to run overtime of more than
10 minutes after seeing the iterations were finished (as prompted in
the command line).
I pasted the modified code below and hope to get a fix to the problem,
any clue would be much appreciated.
Regards,
Hong
P.S. the code,
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkImageRegistrationMethod.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkImage.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkBSplineDeformableTransform.h"
#include "itkLBFGSBOptimizer.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCheckerBoardImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
//
// 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; // line 24
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( !(itk::IterationEvent().CheckEvent( &event )) )
{
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 << " [CheckerBoardBefore] [CheckerBoardAfter] ";
std::cerr << " [DeformationField] ";
return 1;
}
const unsigned int ImageDimension = 2;
typedef signed short PixelType;
typedef itk::Image< PixelType, ImageDimension > FixedImageType;
typedef itk::Image< PixelType, ImageDimension > MovingImageType;
const unsigned int SpaceDimension = ImageDimension;
const unsigned int SplineOrder = 3;
typedef double CoordinateRepType;
typedef itk::BSplineDeformableTransform<
CoordinateRepType,
SpaceDimension,
SplineOrder > TransformType;
typedef itk::LBFGSBOptimizer 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();
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetInterpolator( interpolator );
TransformType::Pointer transform = TransformType::New();
registration->SetTransform( transform );
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 );
typedef TransformType::RegionType RegionType;
RegionType bsplineRegion;
RegionType::SizeType gridSizeOnImage;
RegionType::SizeType gridBorderSize;
RegionType::SizeType totalGridSize;
gridSizeOnImage.Fill( 12 );
gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower, 2 upper
)
totalGridSize = gridSizeOnImage + gridBorderSize;
bsplineRegion.SetSize( totalGridSize );
typedef TransformType::SpacingType SpacingType;
SpacingType spacing = fixedImage->GetSpacing();
typedef TransformType::OriginType OriginType;
OriginType origin = fixedImage->GetOrigin();;
FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();
for(unsigned int r=0; r<ImageDimension; r++)
{
spacing[r] *= floor( static_cast<double>(fixedImageSize[r] - 1) /
static_cast<double>(gridSizeOnImage[r] - 1) );
origin[r] -= spacing[r];
}
transform->SetGridSpacing( spacing );
transform->SetGridOrigin( origin );
transform->SetGridRegion( bsplineRegion );
typedef TransformType::ParametersType ParametersType;
const unsigned int numberOfParameters =
transform->GetNumberOfParameters();
ParametersType parameters( numberOfParameters );
parameters.Fill( 0.0 );
transform->SetParameters( parameters );
registration->SetInitialTransformParameters( transform->GetParameters() );
OptimizerType::BoundSelectionType boundSelect(
transform->GetNumberOfParameters() );
OptimizerType::BoundValueType upperBound( transform->GetNumberOfParameters()
);
OptimizerType::BoundValueType lowerBound( transform->GetNumberOfParameters()
);
boundSelect.Fill( 0 );
upperBound.Fill( 0.0 );
lowerBound.Fill( 0.0 );
optimizer->SetBoundSelection( boundSelect );
optimizer->SetUpperBound( upperBound );
optimizer->SetLowerBound( lowerBound );
optimizer->SetCostFunctionConvergenceFactor( 1e+7 );
optimizer->SetProjectedGradientTolerance( 1e-4 );
optimizer->SetMaximumNumberOfIterations( 500 );
optimizer->SetMaximumNumberOfEvaluations( 500 );
optimizer->SetMaximumNumberOfCorrections( 12 );
//
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
metric->SetNumberOfHistogramBins( 50 );
const unsigned int numberOfSamples = fixedRegion.GetNumberOfPixels() / 10;
metric->SetNumberOfSpatialSamples( numberOfSamples );
metric->ReinitializeSeed( 76926294 );
// Add a time probe
itk::TimeProbesCollectorBase collector;
std::cout << std::endl << "Starting Registration" << std::endl;
try
{
collector.Start( "Registration" );
registration->StartRegistration();
collector.Stop( "Registration" );
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
OptimizerType::ParametersType finalParameters =
registration->GetLastTransformParameters();
// std::cout << "Last Transform Parameters" << std::endl;
// std::cout << finalParameters << std::endl;
// Report the time taken by the registration
collector.Report();
transform->SetParameters( finalParameters );
//
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform( transform );
resample->SetInput( movingImageReader->GetOutput() );
resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resample->SetOutputOrigin( fixedImage->GetOrigin() );
resample->SetOutputSpacing( fixedImage->GetSpacing() );
resample->SetDefaultPixelValue( 100 );
typedef signed short 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() );
writer->Update();
//
// Generate checkerboards before and after registration
//
typedef itk::CheckerBoardImageFilter< FixedImageType > CheckerBoardFilterType;
CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
checker->SetInput1( fixedImage );
checker->SetInput2( resample->GetOutput() );
caster->SetInput( checker->GetOutput() );
writer->SetInput( caster->GetOutput() );
resample->SetDefaultPixelValue( 0 );
////
// Before registration
TransformType::Pointer identityTransform = TransformType::New();
identityTransform->SetIdentity();
resample->SetTransform( identityTransform );
if( argc > 4 )
{
writer->SetFileName( argv[4] );
writer->Update();
}
// After registration
resample->SetTransform( transform );
if( argc > 5 )
{
writer->SetFileName( argv[5] );
writer->Update();
}
////
// Generate the explicit deformation field resulting from
// the registration.
if( argc > 6 )
{
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 = transform->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 );
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