[Insight-users] How I iniatilize the best way
MeanSquareImageToImageMetric ???
javier silva bravo
javier_silva_bravo at hotmail.com
Mon May 16 10:53:33 EDT 2005
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20050516/9c2dfe46/attachment-0001.htm
-------------- next part --------------
#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;
}
More information about the Insight-users
mailing list