[Insight-users] Setting B-spline coefficients (re-sent and help
needed...)
Karthik Krishnan
Karthik.Krishnan at kitware.com
Tue Jul 11 11:16:39 EDT 2006
The offending lines are
TransformType::Pointer identityTransform = TransformType::New();
identityTransform->SetIdentity();
resample->SetTransform( identityTransform );
You cannot use the transform without setting its parameters first, or at least initializing the size of the parameter array. (This is clearly mentioned in the documentation of the BSPlineDeformableTransform class) I've committed a fix to throw an exception at the SetIdentity() method if you try to do so.
/cvsroot/Insight/Insight/Code/Common/itkBSplineDeformableTransform.txx,v
<-- itkBSplineDeformableTransform.txx
new revision: 1.29; previous revision: 1.28
In any case, you are better off using an itk::IdentityTransform instead
of an Identity BSplinetransform just to resample the moving image to the
fixed image grid. Its just a waste of computation.
typedef itk::IdentityTransform< double, 2 > IdentityTransformType;
IdentitytTransformType::Pointer identityTransform = IdentityTransformType::New();
identityTransform->SetIdentity();
resample->SetTransform( identityTransform );
HTH
-karthik
Xiang, Hong wrote:
> 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;
>}
>
>_______________________________________________
>Insight-users mailing list
>Insight-users at itk.org
>http://www.itk.org/mailman/listinfo/insight-users
>
>
>
More information about the Insight-users
mailing list