[Insight-users] ModelToImageRegistration problem
Stefano Diciotti
diciotti at asp.det.unifi.it
Wed Jun 14 02:50:56 EDT 2006
Hi all,
I am working with the ModelToImageRegistration example to
register a model to an image.
I modified that example to register an image with a model constituted by
an ellipse, using translation, rotation and scale transformations.
In particular I modified:
1. The Fixed image is read from file.
2. The Moving image is generated by a single ellipse.
3. The trasform is the CenteredSimilarity2DTransform that
substituted the Euler2DTransform in order to take into account
the scale transformation
4. The OnePlusOneEvolutionaryOptimizer was substituted with the
RegularStepGradientDescentOptimizer as in the ImageRegistration7 Example.
We have tried the program with several fixed images.
For instance:
* an image generated with ImageJ with a white ellipse
with center (58, 68), minor axis = 34 and major axis = 92 in a black
background with noise.
We tried several initial parameters to register the ellipse in the
image. For example
* the model ellipse center in (0,0), minor axis 17 and major axis 46.
Initial Offset (100,100), initial rotation angle 0.5 radians, initial
factor scale 1
For any choise of the parameters, the output always releaved only 1
iteration of the optimizer.
For example this output was displayed
Number of points in the metric = 617
Result =
Scale = 1
Angle (radians) 0.5
Angle (degrees) 28.6479
Center X = 80.7985
Center Y = 161.905
Translation X = 19.2015
Translation Y = -61.9053
Iterations = 1
Metric value = 11783.5
What is going wrong? Why the optimizer do only 1 iteration,
even if the metric value is high??
Thanks,
Stefano
--
Ing. Stefano Diciotti, PhD
Dip. di Elettronica e Telecomunicazioni
Università degli Studi di Firenze
Via S. Marta, 3
50139 Firenze
Tel. 055 4796 443
website: bioingegneria.det.unifi.it/~diciotti
email: stefano.diciotti at unifi.it
diciotti at asp.det.unifi.it
-------------- next part --------------
/*=========================================================================
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
// This example illustrates the use of the \doxygen{SpatialObject} as a
// component of the registration framework in order to perform model based
// registration. The current example creates a geometrical model composed of
// one ellipse. A metric is defined to evaluate the
// fitness between the geometric model and the image.
// Let's look first at the classes required to support
// SpatialObject. In this example we use the
// \doxygen{EllipseSpatialObject} as the basic shape components and we use the
// \doxygen{GroupSpatialObject} to group them together as a representation of
// a more complex shape. Their respective headers are included below.
#include "itkEllipseSpatialObject.h"
#include "itkGroupSpatialObject.h"
// In order to generate the initial synthetic image of the ellipses, we use
// the \doxygen{SpatialObjectToImageFilter} that tests---for every pixel in
// the image---whether the pixel (and hence the spatial object) is
// \emph{inside} or \emph{outside} the geometric model.
#include <itkSpatialObjectToImageFilter.h>
#include "itkImageFileWriter.h"
#include <itkImageToSpatialObjectRegistrationMethod.h>
// A metric is defined to evaluate the fitness between the
// SpatialObject and the Image. The base class for this
// type of metric is the \doxygen{ImageToSpatialObjectMetric}, whose header is
// included below.
#include <itkImageToSpatialObjectMetric.h>
// As in previous registration problems, we have to evaluate the image
// intensity in non-grid positions. The
// \doxygen{LinearInterpolateImageFunction} is used here for this purpose.
#include "itkLinearInterpolateImageFunction.h"
// The SpatialObject is mapped from its own space into the image
// space by using a \doxygen{Transform}. In this
// example, we use the \doxygen{CenteredSimilarity2DTransform}.
#include "itkCenteredTransformInitializer.h"
#include "itkCenteredSimilarity2DTransform.h"
// Registration is fundamentally an optimization problem. Here we include
// the optimizer used to search the parameter space and identify the best
// transformation that will map the shape model on top of the image. The
// optimizer used in this example is the
// \doxygen{RegularStepGradientDescentOptimizer}.
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
// As in previous registration examples, it is important to
// track the evolution of the optimizer as it progresses through the parameter
// space. This is done by using the Command/Observer paradigm. The
// following lines of code implement the \doxygen{Command} observer that
// monitors the progress of the registration. The code is quite
// similar to what we have used in previous registration examples.
#include "itkCommand.h"
class CommandIterationUpdate : public itk::Command
{
public:
typedef CommandIterationUpdate Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
CommandIterationUpdate() {};
// Type defining the optimizer.
public:
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef const OptimizerType * OptimizerPointer;
// Execute method will print data at each iteration
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->GetCurrentPosition() << std::endl;
}
};
// Consider now the most critical component of this new registration
// approach: the metric. This component evaluates the match between the
// SpatialObject and the Image. The
// smoothness and regularity of the metric determine the difficulty of the
// task assigned to the optimizer. In this case, we use a very robust
// optimizer that should be able to find its way even in the most
// discontinuous cost functions. The metric to be implemented should derive
// from the ImageToSpatialObjectMetric class.
//
// The following code implements a simple metric that computes the sum of
// the pixels that are inside the spatial object. In fact, the metric
// maximum is obtained when the model and the image are aligned. The metric
// is templated over the type of the SpatialObject and the type of
// the Image.
template <typename TFixedImage, typename TMovingSpatialObject>
class SimpleImageToSpatialObjectMetric :
public itk::ImageToSpatialObjectMetric<TFixedImage,TMovingSpatialObject>
{
public:
// Standard class typedefs.
typedef SimpleImageToSpatialObjectMetric Self;
typedef itk::ImageToSpatialObjectMetric<TFixedImage,TMovingSpatialObject>
Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef itk::Point<double,2> PointType;
typedef std::list<PointType> PointListType;
typedef TMovingSpatialObject MovingSpatialObjectType;
typedef typename Superclass::ParametersType ParametersType;
typedef typename Superclass::DerivativeType DerivativeType;
typedef typename Superclass::MeasureType MeasureType;
// Method for creation through the object factory.
itkNewMacro(Self);
// Run-time type information (and related methods).
itkTypeMacro(SimpleImageToSpatialObjectMetric, ImageToSpatialObjectMetric);
// Specify the moving spatial object.
void SetMovingSpatialObject( const MovingSpatialObjectType * object)
{
if(!this->m_FixedImage)
{
std::cout << "Please set the image before the moving spatial object" << std::endl;
return;
}
this->m_MovingSpatialObject = object;
m_PointList.clear();
typedef itk::ImageRegionConstIteratorWithIndex<TFixedImage> myIteratorType;
myIteratorType it(this->m_FixedImage,this->m_FixedImage->GetBufferedRegion());
itk::Point<double,2> point;
while(!it.IsAtEnd())
{
for(unsigned int i=0;i<Self::ObjectDimension;i++)
{
point[i]=it.GetIndex()[i];
}
if(this->m_MovingSpatialObject->IsInside(point,99999))
{
m_PointList.push_back(point);
}
++it;
}
std::cout << "Number of points in the metric = " << static_cast<unsigned long>( m_PointList.size() ) << std::endl;
}
unsigned int GetNumberOfParameters(void) const {return 6;};
// Get the Derivatives of the Match Measure
void GetDerivative( const ParametersType &, DerivativeType & ) const
{
return;
}
// The fundamental operation of the metric is its \code{GetValue()} method.
// It is in this method that the fitness value is computed. In our current
// example, the fitness is computed over the points of the
// SpatialObject. For each point, its coordinates are mapped
// through the transform into image space. The resulting point is used
// to evaluate the image and the resulting value is accumulated in a sum.
// Since we are not allowing scale changes, the optimal value of the sum
// will result when all the SpatialObject points are mapped on
// the white regions of the image. Note that the argument for the
// \code{GetValue()} method is the array of parameters of the transform.
//
// \index{ImageToSpatialObjectMetric!GetValue()}
// Get the value for SingleValue optimizers.
MeasureType GetValue( const ParametersType & parameters ) const
{
double value;
this->m_Transform->SetParameters( parameters );
PointListType::const_iterator it = m_PointList.begin();
typename TFixedImage::SizeType size =
this->m_FixedImage->GetBufferedRegion().GetSize();
itk::Index<2> index;
itk::Index<2> start = this->m_FixedImage->GetBufferedRegion().GetIndex();
value = 0;
while(it != m_PointList.end())
{
PointType transformedPoint = this->m_Transform->TransformPoint(*it);
this->m_FixedImage->TransformPhysicalPointToIndex(transformedPoint,index);
if( index[0]> start[0]
&& index[1]> start[1]
&& index[0]< static_cast< signed long >( size[0] )
&& index[1]< static_cast< signed long >( size[1] ) )
{
value += this->m_FixedImage->GetPixel(index);
}
it++;
}
return value;
}
// Get Value and Derivatives for MultipleValuedOptimizers
void GetValueAndDerivative( const ParametersType & parameters,
MeasureType & Value, DerivativeType & Derivative ) const
{
Value = this->GetValue(parameters);
this->GetDerivative(parameters,Derivative);
}
private:
PointListType m_PointList;
};
// Having defined all the registration components we are ready to put the
// pieces together and implement the registration process.
int main( int argc, char *argv[] )
{
if( argc != 2 )
{
std::cerr << "Too many parameters " << std::endl;
std::cerr << "Usage: " << argv[0] << " inputfilename" << std::endl;
return 1;
}
// First we instantiate the GroupSpatialObject and
// EllipseSpatialObject. These two objects are parameterized by
// the dimension of the space. In our current example a $2D$ instantiation
// is created.
typedef itk::GroupSpatialObject< 2 > GroupType;
typedef itk::EllipseSpatialObject< 2 > EllipseType;
// The image is instantiated in the following lines using the pixel
// type and the space dimension. This image uses a \code{float} pixel
// type since we plan to blur it in order to increase the capture radius of
// the optimizer. Images of real pixel type behave better under blurring
// than those of integer pixel type.
typedef itk::Image< float, 2 > ImageType;
typedef float InputPixelType;
typedef itk::Image< InputPixelType, 2 > InputImageType;
typedef itk::ImageFileReader< InputImageType > ReaderType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName( argv[1] );
// Here is where the fun begins! In the following lines we create the
// EllipseSpatialObjects using their \code{New()} methods, and
// assigning the results to SmartPointers. These lines will create
// three ellipses.
EllipseType::Pointer ellipse1 = EllipseType::New();
// Every class deriving from SpatialObject has particular
// parameters enabling the user to tailor its shape. In the case of the
// EllipseSpatialObject, \code{SetRadius()} is used to
// define the ellipse size. An additional \code{SetRadius(Array)} method
// allows the user to define the ellipse axes independently.
EllipseType::ArrayType radius;
radius[0] = 8.5;
radius[1] = 23;
ellipse1->SetRadius( radius );
// The ellipse is created centered in space by default.
// The spatial transform intrinsically associated with the object is
// accessed by the \code{GetTransform()} method. This transform can define
// a translation in space with the \code{SetOffset()} method. We take
// advantage of this feature to place the ellipses at particular
// points in space.
EllipseType::TransformType::OffsetType offset;
offset[ 0 ] = 100.0;
offset[ 1 ] = 100.0;
ellipse1->GetObjectToParentTransform()->SetOffset(offset);
ellipse1->ComputeObjectToWorldTransform();
// Note that after a change has been made in the transform, the
// SpatialObject invokes the method
// \code{ComputeGlobalTransform()} in order to update its global
// transform. The reason for doing this is that SpatialObjects
// can be arranged in hierarchies. It is then possible to change the
// position of a set of spatial objects by moving the parent of the group.
// Now we add the three EllipseSpatialObjects to a
// GroupSpatialObject that will be subsequently passed on to the
// registration method. The GroupSpatialObject facilitates the
// management of the three ellipses as a higher level structure
// representing a complex shape. Groups can be nested any number of levels
// in order to represent shapes with higher detail.
GroupType::Pointer group = GroupType::New();
group->AddSpatialObject( ellipse1 );
// Having the geometric model ready, we proceed to generate the binary
// image representing the imprint of the space occupied by the ellipses.
// The SpatialObjectToImageFilter is used to that end. Note that
// this filter is instantiated over the spatial object used and the image
// type to be generated.
typedef itk::SpatialObjectToImageFilter< GroupType, ImageType >
SpatialObjectToImageFilterType;
// With the defined type, we construct a filter using the \code{New()}
// method. The newly created filter is assigned to a SmartPointer.
SpatialObjectToImageFilterType::Pointer imageFilter =
SpatialObjectToImageFilterType::New();
// The GroupSpatialObject is passed as input to the filter.
imageFilter->SetInput( group );
// The \doxygen{SpatialObjectToImageFilter} acts as a resampling filter.
// Therefore it requires the user to define the size of the desired output
// image. This is specified with the \code{SetSize()} method.
ImageType::SizeType size;
size[ 0 ] = 200;
size[ 1 ] = 200;
imageFilter->SetSize( size );
// Finally we trigger the execution of the filter by calling the
// \code{Update()} method.
try
{
imageFilter->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
// In order to obtain a smoother metric, we blur the image using a
// \doxygen{DiscreteGaussianImageFilter}. This extends the capture radius
// of the metric and produce a more continuous cost function to
// optimize. The following lines instantiate the Gaussian filter and
// create one object of this type using the \code{New()} method.
typedef itk::DiscreteGaussianImageFilter< ImageType, ImageType >
GaussianFilterType;
GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New();
// The output of the SpatialObjectToImageFilter is connected as
// input to the DiscreteGaussianImageFilter.
typedef itk::ImageFileWriter< ImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( "imagefilter.hdr" );
// rescaler->SetInput( filter->GetOutput() );
// writer->SetInput( rescaler->GetOutput() );
writer->SetInput( imageFilter->GetOutput() );
try
{
writer->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
gaussianFilter->SetInput( reader->GetOutput() );
// The variance of the filter is defined as a large value in order to
// increase the capture radius. Finally the execution of the filter is
// triggered using the \code{Update()} method.
const double variance = 20;
gaussianFilter->SetVariance(variance);
try
{
gaussianFilter->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
// The following lines instantiate the type of the
// \doxygen{ImageToSpatialObjectRegistrationMethod} method and instantiate a
// registration object with the \code{New()} method. Note that the
// registration type is templated over the Image and the
// SpatialObject types. The spatial object in this case is the
// group of spatial objects.
typedef itk::ImageToSpatialObjectRegistrationMethod< ImageType, GroupType >
RegistrationType;
RegistrationType::Pointer registration = RegistrationType::New();
// Now we instantiate the metric that is templated over
// the image type and the spatial object type. As usual, the \code{New()}
// method is used to create an object.
typedef SimpleImageToSpatialObjectMetric< ImageType, GroupType > MetricType;
MetricType::Pointer metric = MetricType::New();
// An interpolator will be needed to evaluate the image at non-grid
// positions. Here we instantiate a linear interpolator type.
typedef itk::LinearInterpolateImageFunction< ImageType, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
// The following lines instantiate the optimizer.
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
OptimizerType::Pointer optimizer = OptimizerType::New();
// Next, we instantiate the transform class. In this case we use the
// CenteredSimilarity2DTransform that implements a rigid transform in $2D$
// space.
typedef itk::CenteredSimilarity2DTransform< double > TransformType;
TransformType::Pointer transform = TransformType::New();
// All the components are plugged into the
// ImageToSpatialObjectRegistrationMethod object. The typical
// \code{Set()} methods are used here. Note the use of the
// \code{SetMovingSpatialObject()} method for connecting the spatial object.
// We provide the blurred version of the original synthetic binary
// image as the input image.
double steplength = 1.0; // 1.0
optimizer->SetMaximumStepLength( steplength );
optimizer->SetMinimumStepLength( 0.0001 );
optimizer->SetNumberOfIterations( 500 );
typedef OptimizerType::ScalesType OptimizerScalesType;
OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
const double translationScale = 1.0 / 100.0;
optimizerScales[0] = 10.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = translationScale;
optimizerScales[3] = translationScale;
optimizerScales[4] = translationScale;
optimizerScales[5] = translationScale;
optimizer->SetScales( optimizerScales );
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
registration->SetFixedImage( gaussianFilter->GetOutput() );
registration->SetMovingSpatialObject( group );
registration->SetTransform( transform );
registration->SetInterpolator( interpolator );
registration->SetOptimizer( optimizer );
registration->SetMetric( metric );
typedef itk::CenteredTransformInitializer<
TransformType,
InputImageType,
ImageType > TransformInitializerType;
TransformInitializerType::Pointer initializer = TransformInitializerType::New();
initializer->SetTransform( transform );
initializer->SetFixedImage( reader->GetOutput() );
initializer->SetMovingImage( imageFilter->GetOutput() );
initializer->MomentsOn();
initializer->InitializeTransform();
double initialScale = 1.0;
double initialAngle = 0.5;
transform->SetScale( initialScale );
transform->SetAngle( initialAngle );
// We now pass the parameter of the current transform as the initial
// parameters to be used when the registration process starts.
registration->SetInitialTransformParameters( transform->GetParameters() );
//optimizer->MaximizeOn();
//Finally, we trigger the execution of the registration process with the
// \code{StartRegistration()} method. We place this call in a
// \code{try/catch} block in case any exception is thrown during the
// process.
try
{
registration->StartRegistration();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
OptimizerType::ParametersType finalParameters =
registration->GetLastTransformParameters();
const double finalScale = finalParameters[0];
const double finalAngle = finalParameters[1];
const double finalRotationCenterX = finalParameters[2];
const double finalRotationCenterY = finalParameters[3];
const double finalTranslationX = finalParameters[4];
const double finalTranslationY = finalParameters[5];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
// Print out results
const double finalAngleInDegrees = finalAngle * 45.0 / atan(1.0);
std::cout << std::endl;
std::cout << "Result = " << std::endl;
std::cout << " Scale = " << finalScale << std::endl;
std::cout << " Angle (radians) " << finalAngle << std::endl;
std::cout << " Angle (degrees) " << finalAngleInDegrees << std::endl;
std::cout << " Center X = " << finalRotationCenterX << std::endl;
std::cout << " Center Y = " << finalRotationCenterY << std::endl;
std::cout << " Translation X = " << finalTranslationX << std::endl;
std::cout << " Translation Y = " << finalTranslationY << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
typedef itk::ResampleImageFilter<ImageType,
InputImageType > ResampleFilterType;
TransformType::Pointer finalTransform = TransformType::New();
finalTransform->SetParameters( finalParameters );
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetTransform( finalTransform );
resampler->SetInput( imageFilter->GetOutput() );
ImageType::Pointer fixedImage = reader->GetOutput();
resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resampler->SetOutputOrigin( fixedImage->GetOrigin() );
resampler->SetOutputSpacing( fixedImage->GetSpacing() );
resampler->SetDefaultPixelValue( 0 );
typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, 2> OutputImageType;
typedef itk::CastImageFilter< InputImageType, OutputImageType >
CastFilterType;
typedef itk::ImageFileWriter< OutputImageType > writerType;
writerType::Pointer writer2 = writerType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer2->SetFileName( "reg.hdr" );
caster->SetInput( resampler->GetOutput() );
writer2->SetInput( caster->GetOutput() );
try
{
writer2->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
return 0;
}
-------------- next part --------------
No virus found in this outgoing message.
Checked by AVG Free Edition.
Version: 7.1.394 / Virus Database: 268.8.3/362 - Release Date: 12/06/2006
More information about the Insight-users
mailing list