[Insight-users] Mutual Information Error
Jordan Soet
jds at interchange.ubc.ca
Mon Jun 6 14:54:38 EDT 2005
Thanks for your advice, Luis, that's what I've been thinking too,
however, I still can't get it working (although, it doesn't really
matter anymore as I've been able to get it working with the
VersorRigid3DTransform, which should be enough for what I'm doing...
however I'd still like to get the AffineTransform working, so as not to
be restricted by the rigid transform). First I have another question
about the Versor transform though. I'm just wondering about the
optimizer selection for it. Obviously, this transform has the special
optimizer for it, the VersorRigid3DTransformOptimizer. In the software
guide it says that this is made because gradient descent optimizers
don't work with this transform, but what I'm wondering is would the
OnePlusOneEvolutionaryOptimizer work with it? Because it says in the
software guide that this optimizer is very good at working with the
Mutual Information metric, so I kind of wanted to work with it. Which do
you think would be better?
Anyways, back to my initial problem, I thought that was the problem too,
and I'm pretty sure I had originally been initializing it completely
wrong (setting the parameters as all 0's for an identity transform),
however, I've since changed that and it still doesn't work. I originally
tried setting the parameters for an identity transform, [1 0 0 0 1 0 0 0
1], then I tried using the ->SetIdentity() method of the transform, then
I tried using the CenteredTransformInitializer, then I tried using the
CenteredAffineTransform with it and none of it worked (note, I tried all
of these using the same fixed and moving image, and also a fixed and
moving image with the same origin and spacing). Anyways, I attached my
code this time, if you can figure anything out, that'd be great, if not,
thanks for all your help, you've been extremely helpful.
Thanks,
Jordan
Luis Ibanez wrote:
> Hi Jordan,
>
> The most likely cause is that you are not providing a
> good initialization for the Affine Transform.
>
> You probably want to use the CenteredTransformInitializer
> as described in the ITK Software Guide
>
> http://www.itk.org/ItkSoftwareGuide.pdf
>
>
>
> Regards,
>
>
> Luis
>
>
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: ImageRegistration11.cxx,v $
Language: C++
Date: $Date: 2004/12/28 14:42:48 $
Version: $Revision: 1.6 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
// Software Guide : BeginLatex
//
// This example illustrates how to combine the MutualInformation metric with an
// Evolutionary algorithm for optimization. Evolutionary algorithms are
// naturally well-suited for optimizing the Mutual Information metric given its
// random and noisy behavior.
//
// The structure of the example is almost identical o the one illustrated in
// ImageRegistration4. Therefore we will focus here on the setup that is
// specifically required for the evolutionary optimizer.
//
//
// \index{itk::ImageRegistrationMethod!Multi-Modality}
// \index{itk::OnePlusOneEvolutionaryOptimizer!Multi-Modality}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkImageRegistrationMethod.h"
#include "itkCenteredAffineTransform.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkOnePlusOneEvolutionaryOptimizer.h"
#include "itkNormalVariateGenerator.h"
#include "itkCenteredTransformInitializer.h"
#include "itkImage.h"
// Software Guide : EndCodeSnippet
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.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;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
CommandIterationUpdate() { m_LastMetricValue = 0.0; };
public:
typedef itk::OnePlusOneEvolutionaryOptimizer 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( typeid( event ) != typeid( itk::IterationEvent ) )
{
return;
}
double currentValue = optimizer->GetValue();
// Only print out when the Metric value changes
if( fabs( m_LastMetricValue - currentValue ) > 1e-7 )
{
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << currentValue << " ";
std::cout << optimizer->GetCurrentPosition() << std::endl;
m_LastMetricValue = currentValue;
}
}
private:
double m_LastMetricValue;
};
int main( int argc, char *argv[] )
{
if( argc < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << "outputImagefile "<< std::endl;
return 1;
}
const unsigned int Dimension = 3;
typedef unsigned short PixelType;
typedef itk::Image< PixelType, Dimension > FixedImageType;
typedef itk::Image< PixelType, Dimension > MovingImageType;
typedef itk::CenteredAffineTransform< double, Dimension > TransformType;
typedef itk::OnePlusOneEvolutionaryOptimizer OptimizerType;
typedef itk::LinearInterpolateImageFunction<
MovingImageType,
double > InterpolatorType;
typedef itk::ImageRegistrationMethod<
FixedImageType,
MovingImageType > RegistrationType;
typedef itk::MattesMutualInformationImageToImageMetric<
FixedImageType,
MovingImageType > MetricType;
typedef itk::CenteredTransformInitializer< TransformType,
FixedImageType,
MovingImageType
> TransformInitializerType;
TransformInitializerType::Pointer initializer = TransformInitializerType::New();
TransformType::Pointer transform = TransformType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetOptimizer( optimizer );
registration->SetTransform( transform );
registration->SetInterpolator( interpolator );
MetricType::Pointer metric = MetricType::New();
registration->SetMetric( metric );
metric->ReinitializeSeed(12321);
metric->SetNumberOfHistogramBins( 50 );
metric->UseAllPixelsOn();
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] );
registration->SetFixedImage( fixedImageReader->GetOutput() );
registration->SetMovingImage( movingImageReader->GetOutput() );
initializer->SetTransform( transform );
initializer->SetFixedImage( fixedImageReader->GetOutput() );
initializer->SetMovingImage( movingImageReader->GetOutput() );
initializer->MomentsOn();
try{
initializer->InitializeTransform();
}
catch( itk::ExceptionObject & err )
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return -1;
}
registration->SetInitialTransformParameters( transform->GetParameters() );
registration->SetFixedImageRegion(
fixedImageReader->GetOutput()->GetBufferedRegion() );
typedef itk::Statistics::NormalVariateGenerator GeneratorType;
GeneratorType::Pointer generator = GeneratorType::New();
generator->Initialize(12345);
optimizer->MaximizeOff();
optimizer->SetNormalVariateGenerator( generator );
optimizer->Initialize( 10 );
optimizer->SetEpsilon( 1.0 );
optimizer->SetMaximumIteration( 4000 );
// Translation scale should be 1/(image diagonal)
// = 1/sqrt( 159^2 + 163^2 + 101^2)
float translationScale = 1.0 / 250.0;
typedef OptimizerType::ScalesType OptimizerScalesType;
OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
optimizerScales[0] = 1.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = 1.0;
optimizerScales[3] = 1.0;
optimizerScales[4] = 1.0;
optimizerScales[5] = 1.0;
optimizerScales[6] = 1.0;
optimizerScales[7] = 1.0;
optimizerScales[8] = 1.0;
optimizerScales[9] = translationScale;
optimizerScales[10] = translationScale;
optimizerScales[11] = translationScale;
optimizerScales[12] = translationScale;
optimizerScales[13] = translationScale;
optimizerScales[14] = translationScale;
optimizer->SetScales( optimizerScales );
// Software Guide : EndCodeSnippet
// Create the Command observer and register it with the optimizer.
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
try
{
registration->StartRegistration();
}
catch( itk::ExceptionObject & err )
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return -1;
}
typedef RegistrationType::ParametersType ParametersType;
ParametersType finalParameters = registration->GetLastTransformParameters();
unsigned int numberOfIterations = optimizer->GetCurrentIteration();
double bestValue = optimizer->GetValue();
// Print out results
//
std::cout << "Result: " << std::endl;
std::cout << " Parameters = " << finalParameters << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
TransformType::Pointer finalTransform = TransformType::New();
finalTransform->SetParameters( finalParameters );
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform( finalTransform );
resample->SetInput( movingImageReader->GetOutput() );
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resample->SetOutputOrigin( fixedImage->GetOrigin() );
resample->SetOutputSpacing( fixedImage->GetSpacing() );
resample->SetDefaultPixelValue( 0 );
typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, Dimension > 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();
return 0;
}
More information about the Insight-users
mailing list