[Insight-users] How I iniatilize the best
way MeanSquareImageToImageMetric ???
Luis Ibanez
luis.ibanez at kitware.com
Wed May 18 22:23:21 EDT 2005
Hi Javier,
You shouldn't trust everything that Luis tells you...
but, anyways,
Please post to the list the output of the Command observer
that you connected to the optimizer.
This should include the values of the metric and the
transform parameters at every iteraion of the optimzation
process.
From these values it is easier to figure out what you are
missing in the preparation of the registration.
Thanks
Luis
---------------------------------------------------------------------
javier silva bravo wrote:
> Hello Everybody.
>
> Luis told me when the outputimage is like movil images the method falled.
> So the reason is (he said) I didnt initialize the good way the metric.
> Then How I initialize the metric?
> this is the code.
>
>
> #include "itkMeanSquaresImageToImageMetric.h"
>
> typedef itk::MeanSquaresImageToImageMetric<
>
> FixedImageType,
>
> MovingImageType > MetricType;
>
> MetricType::Pointer metric = MetricType::New();
>
> registration->SetMetric( metric );
>
> metric->SetTransform( transform );
>
> metric->SetFixedImage( fixedImage );
>
> metric->SetMovingImage( movingImageReader->GetOutput() );
>
> metric->SetFixedImageRegion( fixedRegion );
>
> metric->SetTransformParameters( transform->GetParameters() );
>
> optimizer->SetCostFunction(metric);
>
> Thank you.
>
>
> ------------------------------------------------------------------------
> MSN Premium. Protégete, Comunícate y Diviértete Haz clic aquí
> <http://g.msn.com/8HMBESMX/2746??PS=47575> #include
> "itkImageRegistrationMethod.h"
> #include "itkMeanSquaresImageToImageMetric.h"
> //#include "itkLinearInterpolateImageFunction.h"
> #include "itkBSplineInterpolateImageFunction.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 "itkSquaredDifferenceImageFilter.h"
> #include "itkCommand.h"
> #include <fstream>
>
> class CommandIterationUpdate : public itk::Command
> {
> public:
> typedef CommandIterationUpdate Self;
> typedef itk::Command Superclass;
> typedef itk::SmartPointer<Self> Pointer;
> itkNewMacro( Self );
> protected:
> CommandIterationUpdate() {};
> std::ofstream m_file;
>
> public:
> typedef itk::LBFGSBOptimizer OptimizerType;
> typedef const OptimizerType * OptimizerPointer;
>
> void SetFileName(const char* str)
> {
> try{
> m_file.open(str, std::ios::trunc);
> }
> catch(itk::ExceptionObject & err)
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << "Error al tratar de leer el archivo de los valores" <<
> std::endl;
> std::cerr << err << std::endl;
> };
> }
>
>
> 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;
>
> m_file << optimizer->GetCurrentIteration() << " ";
> m_file << optimizer->GetValue() << " ";
> m_file << 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;
> }*/
> if( argc < 3 )
> {
> std::cerr << "Missing Parameters " << std::endl;
> std::cerr << "Usage: " << argv[0];
> std::cerr << " fixedImageFile movingImageFile outputImagefile ";
> std::cerr << " [ValuesField] ";
> return 1;
> }
>
>
> std::cout << "Inicio del programa " << std::endl;
> const unsigned int ImageDimension = 2;
> typedef float 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::MeanSquaresImageToImageMetric<
> FixedImageType,
> MovingImageType > MetricType;
> /*
> typedef itk:: LinearInterpolateImageFunction<
> MovingImageType,
> double > InterpolatorType;*/
> typedef itk::BSplineInterpolateImageFunction<
> MovingImageType,
> double,
> 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();
>
> interpolator->SetSplineOrder(SplineOrder);
>
> metric->SetInterpolator( interpolator );
>
> registration->SetMetric( metric );
> registration->SetOptimizer( optimizer );
> registration->SetInterpolator( interpolator );
>
>
>
> TransformType::Pointer transform = TransformType::New();
> registration->SetTransform( transform );
> metric->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() );
>
> metric->SetFixedImage( fixedImage );
> metric->SetMovingImage( movingImageReader->GetOutput() );
>
> try{
> fixedImageReader->Update();
> }
> catch(itk::ExceptionObject & err)
> {
> std::cerr<<"Error al leer Archivo"<<std::endl;
> std::cerr<<err<<std::endl;
> return -1;
> }
> FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
>
> registration->SetFixedImageRegion( fixedRegion );
> metric->SetFixedImageRegion( fixedRegion );
>
> typedef TransformType::RegionType RegionType;
> RegionType bsplineRegion;
> RegionType::SizeType gridSizeOnImage;
> RegionType::SizeType gridBorderSize;
> RegionType::SizeType totalGridSize;
>
> gridSizeOnImage.Fill( 15 );
> 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() );
> metric->SetTransformParameters( 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->SetCostFunction(metric);
> optimizer->SetCostFunctionConvergenceFactor( 1e+1 );
> optimizer->SetProjectedGradientTolerance( 0.019999 );//1.0 0.155555
> optimizer->SetMaximumNumberOfIterations( 100 );
> optimizer->SetMaximumNumberOfEvaluations( 500 );
> optimizer->SetMaximumNumberOfCorrections( 12 );
> std::cout << " SetProjectedGradientTolerance 0.019999 " << std::endl;
>
> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
> try{
> observer->SetFileName(argv[4]);
> }
> catch(itk::ExceptionObject & err)
> {
> std::cerr <<" Error al leer Archivo "<<std::endl;
> std::cerr <<err<<std::endl;
> }
> optimizer->AddObserver( itk::IterationEvent(), observer );
>
>
>
> itk::TimeProbesCollectorBase collector;
>
> std::cout << std::endl << "Starting Registration" << std::endl;
> std::cout << "Iteracion Valor InfinityNormOfProjectedGradient " <<
> 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;
> }
> collector.Report();
>
> OptimizerType::ParametersType finalParameters =
> registration->GetLastTransformParameters();
>
> transform->SetParameters( finalParameters );
>
> typedef itk::ResampleImageFilter<
> MovingImageType,
> FixedImageType> ResampleFilterType;
>
> ResampleFilterType::Pointer resample = ResampleFilterType::New();
>
> resample->SetTransform( transform );
> resample->SetInput( movingImageReader->GetOutput() );
> resample->SetInterpolator(interpolator);
> 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() );
>
> std::cout << "Creando Imagen De Salida ..." << std::endl;
> 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.
> if( argc >= 7 )
> {
>
> 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;
> }
> }*/
> std::cout << "El programa termino satisfactoriamente" << std::endl;
> 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