[Insight-users] Re: Problems with Insight's registration application
Terry Tang
ttang at homeworking.org
Sun Oct 22 22:49:41 EDT 2006
Hello Luis,
Thanks for your reply. I have tried the IntensityBased2D3DRegistration application with the same commands too, but there was only 1 iteration for each level; the final parameters were tx=ty=tz=0. The problem is that when regGradDescOptimizer calculates the transformedGradient, it is zero, so the optimization ended without further iterations.
Did you get zero's for the final parameters too? If not, could it be that we are using different versions? I noticed that there were some bugs with the gradientDifferenceMetric, is that true still? (I'm using revision 12.1.05 of this metric, with the same IntensityBased2D3DRegistration.cxx as provided by InsightApplications-2.8).
I also tried with a different metric (Mattes MI), but I would have no (0) samples and the exception is thrown (see 2D3DRegistration.cxx).
My other questions are:
1. The size and resolution of the projectedImage in GenerateProjection.cxx is fixed to 256 and 1.0 respectively, what is the rationale behind fixing these parameters?
2. How can we estimate the initial parameters for different transforms?
a) How can we estimate the rotation scales/ translation scales for VersorRigid3DTransform?
b) I have tried to use VersorRigid3DTransform in place of EulerCenteredTransform; if I were to use the same set of brain images, would I need to translate the 3d image to [-100,-100,-100] as well? Particularly, if the images I need to register have the following res/ sizes: 2D: [256,256], [2.21,2.21]/ 3D: [128,128,90], [4.8,4.8,4.8], could I assume that the center of the 3D is an ideal start position? (see attached images)
Thanks so much for your help.
Ying
--- luis.ibanez at kitware.com wrote:
From: Luis Ibanez <luis.ibanez at kitware.com>
To: ttang at homeworking.org
CC: Insight Users <insight-users at itk.org>
Subject: Re: Problems with Insight's registration application
Date: Sun, 22 Oct 2006 19:15:52 -0400
Hi Ying,
I'm assuming that you are talking about the application:
InsightApplications/IntensityBased2D3DRegistration/
Is that right ?
[In the future, you may want to be more specific when you
post your questions to the mailing list]
---
I just retried the commands that are described in the README.txt
file and they work for me. I'm not getting any error related
to number of parameters.
Did you changed the Transform of this code ?
If so, could you please post your code to the list ?
Thanks
Luis
PS: Here are the commands I tried, just as described in
the README.txt file:
==========================================================
A) GenerateProjection
A.1) Translate the image by (-100,-100,-100) then project.
./GenerateProjection -t -100 -100 -100
/data/BrainWeb/brainweb165a10f17.mha projection.mhd
A.2) Translate the image by (-100,-100,-100),
rotate 10 degrees around x, then project.
./GenerateProjection -t -100 -100 -100 -rx 10
/data/BrainWeb/brainweb165a10f17.mha projection.mhd
B) IntensityBased2D3DRegistration
B.1) Register the image that was generated.
./IntensityBased2D3DRegistration -t -100 -100 -100 projection.mhd
/data/BrainWeb/brainweb165a10f17.mha
B.2) Register the image that was generated with rotation.
./IntensityBased2D3DRegistration -t -100 -100 -100 -rx 10
projection.mhd /data/BrainWeb/brainweb165a10f17.mha
===================================================
-----------------------------
ttang at homeworking.org wrote:
> Hello Luis,
>
>
>
> I would like to perform 2d-3d registration using mutual information and possibly other metrics. I have tried running the Insight example but experience the same problem as Rafi did (see below):
>
>
>
>
> by rafi Feb 15, 2006; 03:26am :: Rate this Message: (use ratings to moderate[?])
>
> Reply | Reply to Author | Show Only this Message | Link to this Message
>
>
> Hi Karthik,
>
> Thank you very much; it is what I was looking for.
>
> However, it seems that the code is not working.
>
> I did what it says in the README file and it fails at run time with an exception claiming the number of parameters defined was 6 while the transformation expected 9...
>
> After fixing this, it still does not work and finish after the first iteration claiming the gradient difference is zero (I followed it in the debugger) while the 2 images are not the same.
>
>>From the debugger, it seems that the gradient difference is always zero no mater how you change the parameters
>
> Can you advice what is missing there?
>
> Thanks,
>
> Rafi
>
>
>
>
> I tried to use Mattes' MI instead of graident difference but then I had problems with calculating the derivative. What is the status of this application? Could you please give us some suggestions/ directions on this?
>
>
>
> Thanks for your help,
> Ying
>
>
_____________________________________________________________
Visit http://www.homeworking.com ~ To help you work at home
-------------- next part --------------
// This program implements an intensity based 2D-3D registration
// algorithm using the \code{RayCastInterpolateImageFunction} class
// and \code{GradientDifferenceImageToImageMetric} similarity measure.
//#include <iostream>
#include "itkImageRegistrationMethod.h"
#include "itkMultiResolutionImageRegistrationMethod.h"
#include "itkRecursiveMultiResolutionPyramidImageFilter.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkGradientDifferenceImageToImageMetric.h"
#include "itkRayCastInterpolateImageFunction.h"
#include "itkVersorRigid3DTransformOptimizer.h"
#include "itkVersorRigid3DTransform.h"
#include "itkCenteredEuler3DTransform.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkVersorTransformOptimizer.h"
#include "itkVersorTransform.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"
#include "itkCommand.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
class CommandIterationUpdate : public itk::Command
{
public:
typedef CommandIterationUpdate Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
CommandIterationUpdate() {};
public:
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef const OptimizerType *OptimizerPointer;
typedef itk::VersorRigid3DTransformOptimizer Versor3DOptimizerType;
typedef const Versor3DOptimizerType *Versor3DOptimizerPointer;
typedef itk::VersorTransformOptimizer VersorOptimizerType;
typedef const VersorOptimizerType *VersorOptimizerPointer;
typedef Versor3DOptimizerType::ParametersType Versor3DParametersType;
typedef VersorOptimizerType::ParametersType VersorParametersType;
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 );
Versor3DOptimizerPointer versor3DOptimizer = dynamic_cast< Versor3DOptimizerPointer >( object );
VersorOptimizerPointer versorOptimizer = dynamic_cast< VersorOptimizerPointer >( object );
if( typeid( event ) != typeid( itk::IterationEvent ) )
{
return;
}
if(optimizerID == 0)
{
std::cout << optimizer->GetCurrentIteration();
std::cout << ") " << optimizer->GetValue();
std::cout << " " << optimizer->GetCurrentPosition() << std::endl;
}
else if(optimizerID==1)
{
Versor3DParametersType p( 6 );
double a,b,c,d,e,f;
p = optimizer->GetCurrentPosition();
a = p[0]; b= p[1]; c = p[2]; d= p[3]; e= p[4]; f= p[5];
std::cout << optimizer->GetCurrentIteration();
std::cout << ") " << optimizer->GetValue();
std::cout << "[" << a << " " << b << " " << c << " " << d << " " << e << " " << f << "]\n";
}
else
{
}
}
unsigned int optimizerID;
};
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// A second \emph{Command/Observer} is used to reduce the minimum
// registration step length on each occasion that the resolution of
// the multi-scale registration is increased (see
// \doxygen{MultiResImageRegistration1} for more info).
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
template <typename TRegistration>
class RegistrationInterfaceCommand : public itk::Command
{
public:
typedef RegistrationInterfaceCommand Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
RegistrationInterfaceCommand() {};
public:
unsigned int optimizerID;
typedef TRegistration MultiRegistrationType;
typedef MultiRegistrationType * MultiRegistrationPointer;
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef OptimizerType * OptimizerPointer;
typedef itk::VersorRigid3DTransformOptimizer Versor3DOptimizerType;
typedef Versor3DOptimizerType * Versor3DOptimizerPointer;
typedef itk::VersorTransformOptimizer VersorOptimizerType;
typedef VersorOptimizerType * VersorOptimizerPointer;
void Execute(itk::Object * object, const itk::EventObject & event)
{
if( typeid( event ) != typeid( itk::IterationEvent ) )
{
return;
}
MultiRegistrationPointer reg = dynamic_cast<MultiRegistrationPointer>( object );
OptimizerPointer optimizer = dynamic_cast< OptimizerPointer >( reg->GetOptimizer() );
Versor3DOptimizerPointer versor3DOptimizer = dynamic_cast< Versor3DOptimizerPointer >( reg->GetOptimizer() );
if ( reg->GetCurrentLevel() != 0 )
{
if(optimizerID==0)
{
optimizer->SetMaximumStepLength( optimizer->GetCurrentStepLength() );
optimizer->SetMinimumStepLength( optimizer->GetMinimumStepLength() / 4.0 );
}
else
{
versor3DOptimizer->SetMaximumStepLength( optimizer->GetCurrentStepLength() );
versor3DOptimizer->SetMinimumStepLength( optimizer->GetMinimumStepLength() /4.0);
}
}
// optimizer->Print( std::cout );
}
void Execute(const itk::Object * , const itk::EventObject & )
{ return; }
};
void usage()
{
std::cerr << "\n";
std::cerr << "Usage: itkMultiResIntensityBasedRegn2D3D <options> Image2D Volume3D\n";
std::cerr << " Registers a 3D volume to a 3D/2D image. \n\n";
std::cerr << " where <options> is one or more of the following:\n\n";
std::cerr << " <-h> Display (this) usage information\n";
std::cerr << " <-v> Verbose output [default: yes]\n";
std::cerr << " <-dbg> Debugging output [default: no]\n";
std::cerr << " <-mode> Registration mode [default: 1] 1) 3D-3D and 2) 3D-2D \n";
std::cerr << " <-n int> The number of scales to apply [default: 2]\n";
std::cerr << " <-maxScale int> The scale factor corresponding to max resolution [default: 1]\n";
std::cerr << " <-step float float> Maximum and minimum step sizes [default: 4 and 0.01]\n";
std::cerr << " <-fl float> Focal length or source to image distance [default: 400mm]\n";
std::cerr << " <-t float float float> Translation parameter of the camera \n";
std::cerr << " <-trans> Transform [0 = versor3D, 1 = CenteredEuler3D] ";
std::cerr << " <-transScales float float float> Translation scales\n";
std::cerr << " <-rx float> Rotation around x,y,z axis in degrees \n";
std::cerr << " <-ry float>\n";
std::cerr << " <-rz float>\n";
std::cerr << " <-normal float float> The 2D projection normal position [default: 0x0mm]\n";
std::cerr << " <-cor float float float> The centre of rotation relative to centre of volume\n";
std::cerr << " <-threshold float> Threshold [default: 0]\n";
std::cerr << " <-metric int> Metric [default: 0 = Mattes MI]\n";
std::cerr << " <-nBins int> Number of bins [default: 20]\n";
std::cerr << " <-nSamples int> Number of samples [default: 10000]\n";
std::cerr << " <-diff file> Difference image filename\n";
std::cerr << " <-o file> Output image filename\n\n";
std::cerr << " by john.hipwell at kcl.ac.uk\n";
std::cerr << " and thomas at hartkens.de\n";
std::cerr << " (Imaging Sciences KCL London)\n\n";
//exit(1);
}
int main( int argc, char *argv[] )
{
char *fileImage2D = NULL;
char *fileVolume3D = NULL;
char *fileOutput = NULL;
// char *fileDifference;
char myFileImage2D[100];
char myFileVolume3D[100];;
char myFileOutput[100];
char input[100];
bool ok;
bool verbose = true;
bool debug = false;
unsigned int nScales = 2;
int maxScale = 1;
double rx = 0.;
double ry = 0.;
double rz = 0.;
double tx = 0.;
double ty = 0.;
double tz = 0.;
double cx = 0.;
double cy = 0.;
double cz = 0.;
double focalLength = 400.;
double maxStepSize = 4.;
double minStepSize = 1.;
double transXScale, transYScale, transZScale ;
transXScale = transYScale = transZScale = 0.000135;
double o2Dx = 0;
double o2Dy = 0;
double threshold=0;
int metricID=0;
int transformID=0;
int nSamples=10000;
int nBins=30;
int isDebug = 0;
int mode=1;
if(argc<=1)
{
isDebug = 1;
std::cout << "\nEnter filename prefix of 2d image:" ;
std::cin >> input;
sprintf(myFileImage2D, "%s", input);
std::cout << "\nEnter filename prefix of 3d image:" ;
std::cin >> input;
sprintf(myFileVolume3D, "%s", input);
std::cout << "\nEnter filename prefix of output image:" ;
std::cin >> input;
sprintf(myFileOutput, "%s", input);
std::cout << "\nEnter metric to be used:" ;
std::cin >> input;
metricID = atoi(input);
std::cout << "\nEnter transform to be used:" ;
std::cin >> input;
transformID = atoi(input);
}
// Parse command line parameters
while (argc > 1)
{
ok = false;
if ((ok == false) && (strcmp(argv[1], "-h") == 0))
{
argc--; argv++;
ok = true;
usage();
}
if ((ok == false) && (strcmp(argv[1], "-v") == 0))
{
argc--; argv++;
ok = true;
verbose = true;
}
if ((ok == false) && (strcmp(argv[1], "-dbg") == 0))
{
argc--; argv++;
ok = true;
debug = true;
}
if ((ok == false) && (strcmp(argv[1], "-n") == 0))
{
argc--; argv++;
ok = true;
nScales=atoi(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-maxScale") == 0))
{
argc--; argv++;
ok = true;
maxScale=atoi(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-fl") == 0))
{
argc--; argv++;
ok = true;
focalLength=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-step") == 0))
{
argc--; argv++;
ok = true;
maxStepSize=atof(argv[1]);
argc--; argv++;
minStepSize=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-t") == 0))
{
argc--; argv++;
ok = true;
tx=atof(argv[1]);
argc--; argv++;
ty=atof(argv[1]);
argc--; argv++;
tz=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-trans") == 0))
{
argc--; argv++;
ok = true;
transformID=atoi(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-transScales") == 0))
{
argc--; argv++;
ok = true;
transXScale = atof(argv[1]);
argc--; argv++;
transYScale = atof(argv[1]);
argc--; argv++;
transZScale = atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-rx") == 0))
{
argc--; argv++;
ok = true;
rx=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-ry") == 0))
{
argc--; argv++;
ok = true;
ry=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-rz") == 0))
{
argc--; argv++;
ok = true;
rz=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-normal") == 0))
{
argc--; argv++;
ok = true;
o2Dx=atof(argv[1]);
argc--; argv++;
o2Dy=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-cor") == 0))
{
argc--; argv++;
ok = true;
cx=atof(argv[1]);
argc--; argv++;
cy=atof(argv[1]);
argc--; argv++;
cz=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-threshold") == 0))
{
argc--; argv++;
ok = true;
threshold=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-mode") == 0))
{
argc--; argv++;
ok = true;
mode=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-metric") == 0))
{
argc--; argv++;
ok = true;
metricID=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-nBins") == 0))
{
argc--; argv++;
ok = true;
nBins=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-nSamples") == 0))
{
argc--; argv++;
ok = true;
nSamples=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-diff") == 0))
{
argc--; argv++;
ok = true;
// fileDifference = argv[1];
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-o") == 0))
{
argc--; argv++;
ok = true;
fileOutput = argv[1];
argc--; argv++;
}
if (ok == false)
{
if (fileImage2D == NULL)
{
fileImage2D = argv[1];
argc--;
argv++;
}
else if (fileVolume3D == NULL)
{
fileVolume3D = argv[1];
argc--;
argv++;
}
else
{
std::cerr << "ERROR: Can not parse argument " << argv[1] << std::endl;
usage();
}
}
}
if (verbose)
{
if (fileImage2D) std::cout << "Input 2D image: " << fileImage2D << std::endl;
if (fileVolume3D) std::cout << "Input 3D image: " << fileVolume3D << std::endl;
if (fileOutput) std::cout << "Output image: " << fileOutput << std::endl;
}
// Software Guide : BeginLatex
//
// We begin the program proper by defining the 2D and 3D images. The
// \code{MultiResolutionImageRegistrationMethod} requires that both
// images have the same dimension so the 2D image is given
// dimension 3 and the size of the \emph{z} dimension is set to unity.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
const unsigned int Dimension = 3;
typedef short PixelType;
typedef itk::Image< PixelType, Dimension > ImageType2D;
typedef itk::Image< PixelType, Dimension > ImageType3D;
//
// Due to the recursive nature of the process by which the
// downsampled images are computed by the image pyramids of the
// multi-resolution registration, the output images are required to
// have real pixel types. We declare this internal image type to be
// \code{InternalPixelType}:
//
typedef float InternalPixelType;
typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
typedef itk::CenteredEuler3DTransform< double > TransformType;
typedef itk::VersorRigid3DTransform< double > Versor3DTransformType;
typedef itk::VersorRigid3DTransformOptimizer Versor3DOptimizerType;
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef itk::GradientDifferenceImageToImageMetric< InternalImageType, InternalImageType > GradientMetricType;
typedef itk::MattesMutualInformationImageToImageMetric< InternalImageType, InternalImageType > MattesMetricType;
typedef itk::RayCastInterpolateImageFunction< InternalImageType, double > InterpolatorType;
typedef itk::ImageRegistrationMethod< InternalImageType, InternalImageType > RegistrationType;
typedef itk::MultiResolutionImageRegistrationMethod< InternalImageType, InternalImageType > MultiRegistrationType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// In the multi-resolution framework, a
// \doxygen{MultiResolutionPyramidImageFilter} is used to create a
// pyramid of downsampled images. The size of each downsampled image
// is specified by the user in the form of a schedule of shrink
// factors. A description of the filter and the format of the
// schedules are found in Section \ref{sec:ImagePyramids}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::MultiResolutionPyramidImageFilter< InternalImageType, InternalImageType > ImagePyramidType2D;
typedef itk::MultiResolutionPyramidImageFilter< InternalImageType, InternalImageType > ImagePyramidType3D;
typedef itk::RecursiveMultiResolutionPyramidImageFilter< InternalImageType, InternalImageType > FixedImagePyramidType;
typedef itk::RecursiveMultiResolutionPyramidImageFilter< InternalImageType, InternalImageType > MovingImagePyramidType;
ImagePyramidType3D::Pointer imagePyramid3Db = ImagePyramidType3D::New();
ImagePyramidType3D::Pointer imagePyramid3D = ImagePyramidType3D::New();
ImagePyramidType2D::Pointer imagePyramid2D = ImagePyramidType2D::New();
MovingImagePyramidType::Pointer movingPyramid = MovingImagePyramidType::New();
FixedImagePyramidType::Pointer fixedPyramid = FixedImagePyramidType::New();
GradientMetricType::Pointer GradientMetric = GradientMetricType::New();
MattesMetricType::Pointer MattesMetric = MattesMetricType::New();
Versor3DTransformType::Pointer versor3DTransform = Versor3DTransformType::New();
Versor3DOptimizerType::Pointer versor3DOptimizer = Versor3DOptimizerType::New();
TransformType::Pointer transform = TransformType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
MultiRegistrationType::Pointer multiRegistration = MultiRegistrationType::New();
RegistrationType::Pointer registration = RegistrationType::New();
if(metricID==1)
{
registration->SetMetric( GradientMetric );
multiRegistration->SetMetric( GradientMetric );
}
else
{
registration->SetMetric( MattesMetric );
multiRegistration->SetMetric( MattesMetric );
}
MattesMetric->SetNumberOfSpatialSamples( nSamples );
MattesMetric->SetNumberOfHistogramBins( nBins );
multiRegistration->SetInterpolator( interpolator );
registration->SetInterpolator( interpolator );
std::cout << "\nNumber of params used: " << transform->GetNumberOfParameters();
unsigned int shrinkFactors[3], shrinkFactors2[3];
unsigned int xStartFactor, yStartFactor, zStartFactor, xStartFactor2, yStartFactor2, zStartFactor2;
xStartFactor = yStartFactor = zStartFactor= xStartFactor2 = yStartFactor2 = zStartFactor2 = 2;
shrinkFactors[0] = xStartFactor;
shrinkFactors[1] = yStartFactor;
shrinkFactors[2] = zStartFactor;
shrinkFactors2[0] = xStartFactor2;
shrinkFactors2[1] = yStartFactor2;
shrinkFactors2[2] = zStartFactor2;
fixedPyramid->SetStartingShrinkFactors( shrinkFactors );
movingPyramid->SetStartingShrinkFactors( shrinkFactors2 );
fixedPyramid->SetNumberOfLevels( nScales );
movingPyramid->SetNumberOfLevels( nScales );
multiRegistration->SetFixedImagePyramid( fixedPyramid );
multiRegistration->SetMovingImagePyramid( movingPyramid );
if (debug)
{
//metric->DebugOn();
//transform->DebugOn();
if( transformID ==1)
optimizer->DebugOn();
else if( transformID ==0)
versor3DOptimizer->DebugOn();
//interpolator->DebugOn();
//registration->DebugOn();
}
// Software Guide : BeginLatex
//
// The 2- and 3-D images are read from files,
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::ImageFileReader< ImageType2D > ImageReaderType2D;
typedef itk::ImageFileReader< ImageType3D > ImageReaderType3D;
ImageReaderType2D::Pointer imageReader2D = ImageReaderType2D::New();
ImageReaderType3D::Pointer imageReader3D = ImageReaderType3D::New();
if(isDebug)
{
imageReader2D->SetFileName( myFileImage2D );
imageReader3D->SetFileName( myFileVolume3D );
}
else
{
imageReader2D->SetFileName( fileImage2D );
imageReader3D->SetFileName( fileVolume3D );
}
typedef itk::CastImageFilter< ImageType2D, InternalImageType > CastFilterType2D;
typedef itk::CastImageFilter< ImageType3D, InternalImageType > CastFilterType3D;
CastFilterType2D::Pointer caster2D = CastFilterType2D::New();
CastFilterType3D::Pointer caster3D = CastFilterType3D::New();
caster2D->SetInput( imageReader2D->GetOutput() );
caster3D->SetInput( imageReader3D->GetOutput() );
caster2D->Update();
caster3D->Update();
registration->SetFixedImage( caster2D->GetOutput() );
registration->SetMovingImage( caster3D->GetOutput() );
registration->SetFixedImageRegion( caster2D->GetOutput()->GetBufferedRegion() );
multiRegistration->SetFixedImage( caster2D->GetOutput() );
multiRegistration->SetMovingImage( caster3D->GetOutput() );
multiRegistration->SetFixedImageRegion( caster2D->GetOutput()->GetBufferedRegion() );
// Initialise the transform
// ~~~~~~~~~~~~~~~~~~~~~~~~
transform->SetComputeZYX(true);
// Software Guide : BeginLatex
//
// The transform is initialised with the translation [tx,ty,tz] and
// rotation [rx,ry,rz] specified on the command line
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
TransformType::OutputVectorType translation;
translation[0] = tx;
translation[1] = ty;
translation[2] = tz;
transform->SetTranslation(translation);
transform->SetRotation(M_PI/180.0*rx, M_PI/180.0*ry, M_PI/180.0*rz);
const Versor3DTransformType::VersorType versor;
double vx = sin(rx/2.0);
double vy = sin(ry/2.0);
double vz = sin(rz/2.0);
double w = 1.0;
// versor.Set( vx , vy, vz, w);
versor3DTransform->SetTranslation(translation);
versor3DTransform->SetRotation( versor );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The centre of rotation is set by default to the centre of the 3D
// volume but can be offset from this position using a command
// line specified translation [cx,cy,cz]
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
double origin3D[ Dimension ];
const itk::Vector<double, 3> resolution3D = imageReader3D->GetOutput()->GetSpacing();
typedef ImageType3D::RegionType ImageRegionType3D;
typedef ImageRegionType3D::SizeType SizeType3D;
ImageRegionType3D region3D = caster3D->GetOutput()->GetBufferedRegion();
SizeType3D size3D = region3D.GetSize();
origin3D[0] = resolution3D[0]*((double) size3D[0])/2.;
origin3D[1] = resolution3D[1]*((double) size3D[1])/2.;
origin3D[2] = resolution3D[2]*((double) size3D[2])/2.;
TransformType::InputPointType center;
center[0] = cx + origin3D[0];
center[1] = cy + origin3D[1];
center[2] = cz + origin3D[2];
transform->SetCenter(center);
versor3DTransform->SetCenter(center);
if (verbose)
{
std::cout << "3D image size: "
<< size3D[0] << ", " << size3D[1] << ", " << size3D[2] << std::endl
<< " resolution: "
<< resolution3D[0] << ", " << resolution3D[1] << ", " << resolution3D[2] << std::endl
<< "Transform: " << transform << std::endl;
}
// Set the origin of the 2D image
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Software Guide : BeginLatex
//
// For correct (perspective) projection of the 3D volume, the 2D
// image needs to be placed at a certain distance (the focal length
// or source-to-image distance \emph{focalLength}) from the focal
// point, and the normal from the imaging plane to the focal point
// needs to be specified.
//
// By default, ie. if [tx,ty,tz] is [0,0,0], the 3D volume is placed
// half-way between the focal point and the imaging plane. Similarly
// the imaging plane normal is set by default to the center of the
// 2D image but may be modified from this using the command line
// parameters [o2Dx,o2Dy].
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
double origin2D[ Dimension ];
const itk::Vector<double, 3> resolution2D = imageReader2D->GetOutput()->GetSpacing();
typedef ImageType2D::RegionType ImageRegionType2D;
typedef ImageRegionType2D::SizeType SizeType2D;
ImageRegionType2D region2D = caster2D->GetOutput()->GetBufferedRegion();
SizeType2D size2D = region2D.GetSize();
origin2D[0] = origin3D[0] + o2Dx - resolution2D[0]*((double) size2D[0] - 1.)/2.;
origin2D[1] = origin3D[1] + o2Dy - resolution2D[1]*((double) size2D[1] - 1.)/2.;
origin2D[2] = origin3D[2] + focalLength/2.;
imageReader2D->GetOutput()->SetOrigin( origin2D );
if (verbose)
{
std::cout << "2D image size: "
<< size2D[0] << ", " << size2D[1] << ", " << size2D[2] << std::endl
<< " resolution: "
<< resolution2D[0] << ", " << resolution2D[1] << ", " << resolution2D[2] << std::endl
<< " and position: "
<< origin2D[0] << ", " << origin2D[1] << ", " << origin2D[2] << std::endl;
}
// Set the pyramid schedule for the 2D image
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Software Guide : BeginLatex
//
// The multi-resolution hierachy is specified with respect to the 2D
// image. The number of scales or levels defaults to two, with each level
// differing in resolution by a factor of two. The highest resolution (final)
// level corresponds to the resolution of the input 2D image but
// this can be altered by the user via the \emph{maxScale} command
// line parameter. The number of scales \emph{nScales} can also be
// specified by the user.
//
// Software Guide : EndLatex
std::vector<double> sampledRes2D(nScales, 0);
std::vector<ImageRegionType2D> imageRegionPyramid2D(nScales);
imagePyramid2D->SetNumberOfLevels( nScales );
sampledRes2D.reserve( nScales );
imageRegionPyramid2D.reserve( nScales );
typedef ImagePyramidType2D::ScheduleType ScheduleType2D;
ScheduleType2D schedule2D = imagePyramid2D->GetSchedule();
for ( unsigned int dim = 0; dim < ImageType2D::ImageDimension; dim++ )
{
schedule2D[ nScales-1 ][ dim ] = maxScale;
}
for ( int level=nScales-2; level >= 0; level-- )
{
for ( unsigned int dim = 0; dim < ImageType2D::ImageDimension; dim++ )
{
schedule2D[ level ][ dim ] = 2*schedule2D[ level+1 ][ dim ];
}
}
std::cout << "\nSchedule2d:\n"<< schedule2D;
// Compute the 2D ImageRegion corresponding to each level of the
// pyramid.
typedef ImageRegionType2D::IndexType IndexType2D;
IndexType2D inputStart2D = region2D.GetIndex();
for ( unsigned int level=0; level < nScales; level++ )
{
SizeType2D size;
IndexType2D start;
sampledRes2D[ level ] = 0.;
for ( unsigned int dim = 0; dim < ImageType2D::ImageDimension; dim++ )
{
float scaleFactor = static_cast<float>( schedule2D[ level ][ dim ] );
size[ dim ] = static_cast<SizeType2D::SizeValueType>(
floor( static_cast<float>( size2D[ dim ] ) / scaleFactor ) );
if( size[ dim ] < 1 )
{
size[ dim ] = 1;
schedule2D[ level ][ dim ] = 1;
}
std::cout << level << " " << dim << " "
<< size[ dim ] << " " << schedule2D[ level ][ dim ]
<< std::endl;
scaleFactor = static_cast<float>( schedule2D[ level ][ dim ] );
start[ dim ] = static_cast<IndexType2D::IndexValueType>(
ceil( static_cast<float>( inputStart2D[ dim ] ) / scaleFactor ) );
if ( dim < 2 )
{
sampledRes2D[ level ] += scaleFactor*resolution2D[ dim ]
*scaleFactor*resolution2D[ dim ];
}
}
sampledRes2D[ level ] = sqrt( sampledRes2D[ level ] );
imageRegionPyramid2D[ level ].SetSize( size );
imageRegionPyramid2D[ level ].SetIndex( start );
}
imagePyramid2D->SetSchedule( schedule2D );
std::cout << "\nSchedule2d:\n"<< schedule2D;
// Set the pyramid schedule for the 3D image
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Software Guide : BeginLatex
//
// The 3D image pyramid naturally contains the same number of levels
// as the multi-resolution schedule for the 2D image. In addition
// the sub-sampling factors for each level are set such that
// resolution of the 3D image in each case is reduced to no greater
// than the corresponding 2D image resolution at that scale. This
// ensures that the 3D volume is reduced in size as far as possible,
// minimising ray-casting computation time, whilst retaining
// sufficient information to ensure accurate production of the DRR
// at the resolution of the 2D image.
//
// Software Guide : EndLatex
std::vector<ImageRegionType3D> imageRegionPyramid3D(nScales);
imagePyramid3D->SetNumberOfLevels( nScales );
imageRegionPyramid3D.reserve( nScales );
typedef ImagePyramidType3D::ScheduleType ScheduleType3D;
ScheduleType3D schedule3D = imagePyramid3D->GetSchedule();
// Compute the 2D image pyramid schedule such that the 3D volume
// resolution is no greater than the 2D image resolution.
for ( unsigned int level=0; level < nScales; level++ )
{
for ( unsigned int dim = 0; dim < ImageType3D::ImageDimension; dim++)
{
double res = resolution3D[ dim ];
schedule3D[ level ][ dim ] = 1;
while ( res*2. < sampledRes2D[ level ] )
{
schedule3D[ level ][ dim ] *= 2;
res *= 2.;
}
}
}
// Compute the 3D ImageRegion corresponding to each level of the
// pyramid.
typedef ImageRegionType3D::IndexType IndexType3D;
IndexType3D inputStart3D = region3D.GetIndex();
for ( unsigned int level=0; level < nScales; level++ )
{
SizeType3D size;
IndexType3D start;
for ( unsigned int dim = 0; dim < ImageType3D::ImageDimension; dim++)
{
float scaleFactor = static_cast<float>( schedule3D[ level ][ dim ] );
size[ dim ] = static_cast<SizeType3D::SizeValueType>(
floor( static_cast<float>( size3D[ dim ] ) / scaleFactor ) );
if( size[ dim ] < 1 )
{
size[ dim ] = 1;
schedule3D[ level ][ dim ] = 1;
}
scaleFactor = static_cast<float>( schedule3D[ level ][ dim ] );
start[ dim ] = static_cast<IndexType3D::IndexValueType>(
ceil( static_cast<float>( inputStart3D[ dim ] ) / scaleFactor ) );
}
imageRegionPyramid3D[ level ].SetSize( size );
imageRegionPyramid3D[ level ].SetIndex( start );
}
imagePyramid3D->SetSchedule( schedule3D );
std::cout << "\nSchedule3d:\n"<< schedule3D;
// Initialise the ray cast interpolator
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Software Guide : BeginLatex
//
// The ray cast interpolator is used to project the 3D volume. It
// does this by casting rays from the (transformed) focal point to
// each (transformed) pixel coordinate in the 2D image.
//
// In addition a threshold may be specified to ensure that only
// intensities greater than a given value contribute to the
// projected volume. This can be used, for instance, to remove soft
// tissue from projections of CT data and force the registration
// to find a match which aligns bony structures in the images.
//
InterpolatorType::InputPointType focalPoint;
focalPoint[0] = origin3D[0];
focalPoint[1] = origin3D[1];
focalPoint[2] = origin3D[2] - focalLength/2.;
interpolator->SetFocalPoint(focalPoint);
interpolator->SetThreshold(threshold);
interpolator->SetTransform(transform);
// Software Guide : EndCodeSnippet
// NB. Interpolator input image set in itkImageToImageMatric
if (verbose)
{
std::cout << "Focal point: "
<< focalPoint[0] << ", " << focalPoint[1] << ", " << focalPoint[2] << std::endl
<< "Threshold: " << threshold << std::endl;
}
// Set up the transform and start position
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// The maximum and minimum step lengths are specified. The default
// values for these parameters are 4 and 1 (degrees or mm, depending
// on the parameter) respectively, for the coarsest resolution (so
// step sizes of 4, 2 and 1 will be used). For finer scales the
// minimum step size is reduced by a factor of four, so the second
// scale will be performed with step sizes of 1, 0.5 and 0.25 mm or
// degrees.
//
if(metricID==0)
{
optimizer->SetMaximize( true );
versor3DOptimizer->MinimizeOff();
}
else
{
optimizer->SetMinimize( true );
versor3DOptimizer->MinimizeOn();
}
versor3DOptimizer->SetMaximumStepLength( maxStepSize );
versor3DOptimizer->SetMinimumStepLength( minStepSize );
versor3DOptimizer->SetNumberOfIterations( 200 );
itk::VersorRigid3DTransformOptimizer::ScalesType versor3DScales( 6 );
versor3DScales[0] = 1.0;
versor3DScales[1] = 1.0;
versor3DScales[2] = 1.0;
versor3DScales[3] = transXScale;
versor3DScales[4] = transYScale;
versor3DScales[5] = transZScale;
versor3DOptimizer->SetScales( versor3DScales );
//
// The optimizer weightings are set such that one degree equates to
// one millimeter.
optimizer->SetMaximumStepLength( maxStepSize );
optimizer->SetMinimumStepLength( minStepSize );
optimizer->SetNumberOfIterations( 200 );
itk::Optimizer::ScalesType regScales( 9 );
regScales[0] = 180./M_PI;
regScales[1] = 180./M_PI;
regScales[2] = 180./M_PI;
regScales[3] = 1.;
regScales[4] = 1.;
regScales[5] = 1.;
regScales[6] = 1.;
regScales[7] = 1.;
regScales[8] = 1.;
optimizer->SetScales( regScales );
if (verbose)
{
optimizer->Print( std::cout );
}
// Create the observers
// ~~~~~~~~~~~~~~~~~~~~
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
versor3DOptimizer->AddObserver( itk::IterationEvent(), observer );
typedef RegistrationInterfaceCommand<MultiRegistrationType> CommandType;
CommandType::Pointer command = CommandType::New();
multiRegistration->AddObserver( itk::IterationEvent(), command );
multiRegistration->SetNumberOfLevels( nScales );
//imagePyramid3D->Print(std::cout);
//imagePyramid2D->Print(std::cout);
if(transformID==1)
{
registration->SetOptimizer( optimizer );
registration->SetTransform( transform );
std::cout << "\nTransform params: \n" << transform->GetParameters();
registration->SetInitialTransformParameters( transform->GetParameters() );
multiRegistration->SetOptimizer( optimizer );
multiRegistration->SetTransform( transform );
multiRegistration->SetInitialTransformParameters( transform->GetParameters() );
observer->optimizerID = 0;
}
else if( transformID ==0)
{
registration->SetOptimizer( versor3DOptimizer );
registration->SetTransform( versor3DTransform );
registration->SetInitialTransformParameters( versor3DTransform->GetParameters() );
multiRegistration->SetOptimizer( versor3DOptimizer );
multiRegistration->SetTransform( versor3DTransform );
multiRegistration->SetInitialTransformParameters( versor3DTransform->GetParameters() );
observer->optimizerID = 1;
}
// Start the registration
// ~~~~~~~~~~~~~~~~~~~~~~
if (verbose)
{
std::cout << "Starting the registration now..." << std::endl;
}
try
{
if(nScales>1)
multiRegistration->StartRegistration();
else
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;
if(nScales>1)
finalParameters= multiRegistration->GetLastTransformParameters();
else
finalParameters= registration->GetLastTransformParameters();
const double TranslationAlongX = finalParameters[0];
const double TranslationAlongY = finalParameters[1];
const double TranslationAlongZ = finalParameters[2];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
std::cout << "Result = " << std::endl;
std::cout << " Translation X = " << TranslationAlongX << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << std::endl;
std::cout << " Translation Z = " << TranslationAlongZ << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
// Write out the projection image at the registration position
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
typedef itk::ResampleImageFilter< InternalImageType, InternalImageType > FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetInput( caster3D->GetOutput() );
filter->SetDefaultPixelValue( 0 );
if(transformID==0)
filter->SetTransform( versor3DTransform );
else
filter->SetTransform( transform );
filter->SetInterpolator( interpolator );
filter->SetSize( imageReader2D->GetOutput()->GetLargestPossibleRegion().GetSize() );
filter->SetOutputOrigin( imageReader2D->GetOutput()->GetOrigin() );
filter->SetOutputSpacing( imageReader2D->GetOutput()->GetSpacing() );
try {
std::cout << "Writing image: " << fileOutput << std::endl;
filter->Update();
}
catch( itk::ExceptionObject & err ) {
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
std::cout << "\n1st filter updated";
typedef itk::CastImageFilter< InternalImageType, ImageType2D > OutputCastFilterType;
OutputCastFilterType::Pointer outputCaster = OutputCastFilterType::New();
outputCaster->SetInput( filter->GetOutput() );
typedef itk::ImageFileWriter< ImageType2D > WriterType;
WriterType::Pointer writer = WriterType::New();
if(isDebug)
writer->SetFileName( myFileOutput );
else
writer->SetFileName( fileOutput );
writer->SetInput( outputCaster->GetOutput() );
try {
std::cout << "Writing image: " << fileOutput << std::endl;
writer->Update();
}
catch( itk::ExceptionObject & err ) {
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
return 0;
}
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: IntensityBased2D3DRegistration.cxx,v $
Language: C++
Date: $Date: 2005/01/02 20:23:52 $
Version: $Revision: 1.4 $
Copyright (c) 2002 Insight 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 program implements an intensity based 2D-3D registration
// algorithm using the \code{RayCastInterpolateImageFunction} class
// and \code{GradientDifferenceImageToImageMetric} similarity measure.
//
// Software Guide : EndLatex
#include <iostream>
// Software Guide : BeginLatex
//
// The algorithm performs a multi-resolution coarse to fine
// registration using the
// \code{MultiResolutionImageRegistrationMethod} class.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkMultiResolutionImageRegistrationMethod.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The transformation used is a rigid 3D Euler transform with the
// provision of a center of rotation which defaults to the center of
// the 3D volume. In practice the center of the particular feature of
// interest in the 3D volume should be used.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkCenteredEuler3DTransform.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We have chosen to implement the registration using the gradient
// difference similarity measure.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkGradientDifferenceImageToImageMetric.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// This is an intensity based registration algorith so ray casting is
// used to project the 3D volume onto pixels in the target 2D image.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkRayCastInterpolateImageFunction.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Finally a gradient descent optimizer is used to search for the
// registration position.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkRegularStepGradientDescentOptimizer.h"
// Software Guide : EndCodeSnippet
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"
#include "itkCommand.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
// Software Guide : BeginLatex
//
// First we define the command class to allow us to monitor the registration.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
class CommandIterationUpdate : public itk::Command
{
public:
typedef CommandIterationUpdate Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
CommandIterationUpdate() {};
public:
typedef itk::RegularStepGradientDescentOptimizer 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;
}
std::cout << "Iteration: " << optimizer->GetCurrentIteration() << std::endl;
std::cout << "Similarity: " << optimizer->GetValue() << std::endl;
std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl;
}
};
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// A second \emph{Command/Observer} is used to reduce the minimum
// registration step length on each occasion that the resolution of
// the multi-scale registration is increased (see
// \doxygen{MultiResImageRegistration1} for more info).
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
template <typename TRegistration>
class RegistrationInterfaceCommand : public itk::Command
{
public:
typedef RegistrationInterfaceCommand Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
RegistrationInterfaceCommand() {};
public:
typedef TRegistration RegistrationType;
typedef RegistrationType * RegistrationPointer;
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef OptimizerType * OptimizerPointer;
void Execute(itk::Object * object, const itk::EventObject & event)
{
if( typeid( event ) != typeid( itk::IterationEvent ) )
{
return;
}
RegistrationPointer registration =
dynamic_cast<RegistrationPointer>( object );
// If this is the first resolution level we assume that the
// maximum step length (representing the first step size) and the
// minimum step length (representing the convergence criterion)
// have already been set. At each subsequent resolution level, we
// will reduce the minimum step length by a factor of four in order
// to allow the optimizer to focus on progressively smaller
// regions. The maximum step length is set to the current step
// length. In this way, when the optimizer is reinitialized at the
// beginning of the registration process for the next level, the
// step length will simply start with the last value used for the
// previous level. This will guarantee the continuity of the path
// taken by the optimizer through the parameter space.
OptimizerPointer optimizer = dynamic_cast< OptimizerPointer >(
registration->GetOptimizer() );
if ( registration->GetCurrentLevel() != 0 )
{
optimizer->SetMaximumStepLength( optimizer->GetCurrentStepLength() );
optimizer->SetMinimumStepLength( optimizer->GetMinimumStepLength() / 4.0 );
}
optimizer->Print( std::cout );
}
void Execute(const itk::Object * , const itk::EventObject & )
{ return; }
};
// Software Guide : EndCodeSnippet
void usage()
{
std::cerr << "\n";
std::cerr << "Usage: itkMultiResIntensityBasedRegn2D3D <options> Image2D Volume3D\n";
std::cerr << " Registers a 3D volume to a 2D image. \n\n";
std::cerr << " where <options> is one or more of the following:\n\n";
std::cerr << " <-h> Display (this) usage information\n";
std::cerr << " <-v> Verbose output [default: no]\n";
std::cerr << " <-dbg> Debugging output [default: no]\n";
std::cerr << " <-n int> The number of scales to apply [default: 2]\n";
std::cerr << " <-maxScale int> The scale factor corresponding to max resolution [default: 1]\n";
std::cerr << " <-step float float> Maximum and minimum step sizes [default: 4 and 0.01]\n";
std::cerr << " <-fl float> Focal length or source to image distance [default: 400mm]\n";
std::cerr << " <-t float float float> Translation parameter of the camera \n";
std::cerr << " <-rx float> Rotation around x,y,z axis in degrees \n";
std::cerr << " <-ry float>\n";
std::cerr << " <-rz float>\n";
std::cerr << " <-normal float float> The 2D projection normal position [default: 0x0mm]\n";
std::cerr << " <-cor float float float> The centre of rotation relative to centre of volume\n";
std::cerr << " <-threshold float> Threshold [default: 0]\n";
std::cerr << " <-diff file> Difference image filename\n";
std::cerr << " <-o file> Output image filename\n\n";
std::cerr << " by john.hipwell at kcl.ac.uk\n";
std::cerr << " and thomas at hartkens.de\n";
std::cerr << " (Imaging Sciences KCL London)\n\n";
exit(1);
}
int main( int argc, char *argv[] )
{
char *fileImage2D = NULL;
char *fileVolume3D = NULL;
char *fileOutput = NULL;
// char *fileDifference;
bool ok;
bool verbose = false;
bool debug = false;
unsigned int nScales = 2;
int maxScale = 1;
double rx = 0.;
double ry = 0.;
double rz = 0.;
double tx = 0.;
double ty = 0.;
double tz = 0.;
double cx = 0.;
double cy = 0.;
double cz = 0.;
double focalLength = 400.;
double maxStepSize = 4.;
double minStepSize = 1.;
double o2Dx = 0;
double o2Dy = 0;
double threshold=0;
// Parse command line parameters
if (argc <= 1)
usage();
while (argc > 1)
{
ok = false;
if ((ok == false) && (strcmp(argv[1], "-h") == 0))
{
argc--; argv++;
ok = true;
usage();
}
if ((ok == false) && (strcmp(argv[1], "-v") == 0))
{
argc--; argv++;
ok = true;
verbose = true;
}
if ((ok == false) && (strcmp(argv[1], "-dbg") == 0))
{
argc--; argv++;
ok = true;
debug = true;
}
if ((ok == false) && (strcmp(argv[1], "-n") == 0))
{
argc--; argv++;
ok = true;
nScales=atoi(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-maxScale") == 0))
{
argc--; argv++;
ok = true;
maxScale=atoi(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-fl") == 0))
{
argc--; argv++;
ok = true;
focalLength=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-step") == 0))
{
argc--; argv++;
ok = true;
maxStepSize=atof(argv[1]);
argc--; argv++;
minStepSize=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-t") == 0))
{
argc--; argv++;
ok = true;
tx=atof(argv[1]);
argc--; argv++;
ty=atof(argv[1]);
argc--; argv++;
tz=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-rx") == 0))
{
argc--; argv++;
ok = true;
rx=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-ry") == 0))
{
argc--; argv++;
ok = true;
ry=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-rz") == 0))
{
argc--; argv++;
ok = true;
rz=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-normal") == 0))
{
argc--; argv++;
ok = true;
o2Dx=atof(argv[1]);
argc--; argv++;
o2Dy=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-cor") == 0))
{
argc--; argv++;
ok = true;
cx=atof(argv[1]);
argc--; argv++;
cy=atof(argv[1]);
argc--; argv++;
cz=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-threshold") == 0))
{
argc--; argv++;
ok = true;
threshold=atof(argv[1]);
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-diff") == 0))
{
argc--; argv++;
ok = true;
// fileDifference = argv[1];
argc--; argv++;
}
if ((ok == false) && (strcmp(argv[1], "-o") == 0))
{
argc--; argv++;
ok = true;
fileOutput = argv[1];
argc--; argv++;
}
if (ok == false)
{
if (fileImage2D == NULL)
{
fileImage2D = argv[1];
argc--;
argv++;
}
else if (fileVolume3D == NULL)
{
fileVolume3D = argv[1];
argc--;
argv++;
}
else
{
std::cerr << "ERROR: Can not parse argument " << argv[1] << std::endl;
usage();
}
}
}
if (verbose)
{
if (fileImage2D) std::cout << "Input 2D image: " << fileImage2D << std::endl;
if (fileVolume3D) std::cout << "Input 3D image: " << fileVolume3D << std::endl;
if (fileOutput) std::cout << "Output image: " << fileOutput << std::endl;
}
// Software Guide : BeginLatex
//
// We begin the program proper by defining the 2D and 3D images. The
// \code{MultiResolutionImageRegistrationMethod} requires that both
// images have the same dimension so the 2D image is given
// dimension 3 and the size of the \emph{z} dimension is set to unity.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
const unsigned int Dimension = 3;
typedef short PixelType;
typedef itk::Image< PixelType, Dimension > ImageType2D;
typedef itk::Image< PixelType, Dimension > ImageType3D;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The following lines define each of the components used in the
// registration: The transform, optimizer, metric, interpolator and
// the multi-resolution registration method itself.
//
// Due to the recursive nature of the process by which the
// downsampled images are computed by the image pyramids of the
// multi-resolution registration, the output images are required to
// have real pixel types. We declare this internal image type to be
// \code{InternalPixelType}:
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef float InternalPixelType;
typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
typedef itk::CenteredEuler3DTransform< double > TransformType;
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef itk::GradientDifferenceImageToImageMetric<
InternalImageType,
InternalImageType > MetricType;
typedef itk::RayCastInterpolateImageFunction<
InternalImageType,
double > InterpolatorType;
typedef itk::MultiResolutionImageRegistrationMethod<
InternalImageType,
InternalImageType > RegistrationType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// In the multi-resolution framework, a
// \doxygen{MultiResolutionPyramidImageFilter} is used to create a
// pyramid of downsampled images. The size of each downsampled image
// is specified by the user in the form of a schedule of shrink
// factors. A description of the filter and the format of the
// schedules are found in Section \ref{sec:ImagePyramids}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::MultiResolutionPyramidImageFilter<
InternalImageType,
InternalImageType > ImagePyramidType2D;
typedef itk::MultiResolutionPyramidImageFilter<
InternalImageType,
InternalImageType > ImagePyramidType3D;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Each of the registration components are instantiated in the
// usual way...
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
ImagePyramidType2D::Pointer imagePyramid2D =
ImagePyramidType2D::New();
ImagePyramidType3D::Pointer imagePyramid3D =
ImagePyramidType3D::New();
MetricType::Pointer metric = MetricType::New();
TransformType::Pointer transform = TransformType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// and passed to the registration method:
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetTransform( transform );
registration->SetInterpolator( interpolator );
registration->SetFixedImagePyramid( imagePyramid2D );
registration->SetMovingImagePyramid( imagePyramid3D );
// Software Guide : EndCodeSnippet
if (debug)
{
//metric->DebugOn();
//transform->DebugOn();
optimizer->DebugOn();
//interpolator->DebugOn();
//registration->DebugOn();
}
// Software Guide : BeginLatex
//
// The 2- and 3-D images are read from files,
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::ImageFileReader< ImageType2D > ImageReaderType2D;
typedef itk::ImageFileReader< ImageType3D > ImageReaderType3D;
ImageReaderType2D::Pointer imageReader2D = ImageReaderType2D::New();
ImageReaderType3D::Pointer imageReader3D = ImageReaderType3D::New();
imageReader2D->SetFileName( fileImage2D );
imageReader3D->SetFileName( fileVolume3D );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// but before connecting these images to the registration we need
// to cast them to the internal image type using
// \doxygen{CastImageFilters}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::CastImageFilter<
ImageType2D, InternalImageType > CastFilterType2D;
typedef itk::CastImageFilter<
ImageType3D, InternalImageType > CastFilterType3D;
CastFilterType2D::Pointer caster2D = CastFilterType2D::New();
CastFilterType3D::Pointer caster3D = CastFilterType3D::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The output of the readers is connected as input to the cast
// filters. The inputs to the registration method are taken from the
// cast filters.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
caster2D->SetInput( imageReader2D->GetOutput() );
caster3D->SetInput( imageReader3D->GetOutput() );
caster2D->Update();
caster3D->Update();
registration->SetFixedImage( caster2D->GetOutput() );
registration->SetMovingImage( caster3D->GetOutput() );
registration->SetFixedImageRegion(
caster2D->GetOutput()->GetBufferedRegion() );
// Software Guide : EndCodeSnippet
// Initialise the transform
// ~~~~~~~~~~~~~~~~~~~~~~~~
transform->SetComputeZYX(true);
// Software Guide : BeginLatex
//
// The transform is initialised with the translation [tx,ty,tz] and
// rotation [rx,ry,rz] specified on the command line
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
TransformType::OutputVectorType translation;
translation[0] = tx;
translation[1] = ty;
translation[2] = tz;
transform->SetTranslation(translation);
transform->SetRotation(M_PI/180.0*rx, M_PI/180.0*ry, M_PI/180.0*rz);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The centre of rotation is set by default to the centre of the 3D
// volume but can be offset from this position using a command
// line specified translation [cx,cy,cz]
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
double origin3D[ Dimension ];
const itk::Vector<double, 3> resolution3D = imageReader3D->GetOutput()->GetSpacing();
typedef ImageType3D::RegionType ImageRegionType3D;
typedef ImageRegionType3D::SizeType SizeType3D;
ImageRegionType3D region3D = caster3D->GetOutput()->GetBufferedRegion();
SizeType3D size3D = region3D.GetSize();
origin3D[0] = resolution3D[0]*((double) size3D[0])/2.;
origin3D[1] = resolution3D[1]*((double) size3D[1])/2.;
origin3D[2] = resolution3D[2]*((double) size3D[2])/2.;
TransformType::InputPointType center;
center[0] = cx + origin3D[0];
center[1] = cy + origin3D[1];
center[2] = cz + origin3D[2];
transform->SetCenter(center);
// Software Guide : EndCodeSnippet
if (verbose)
{
std::cout << "3D image size: "
<< size3D[0] << ", " << size3D[1] << ", " << size3D[2] << std::endl
<< " resolution: "
<< resolution3D[0] << ", " << resolution3D[1] << ", " << resolution3D[2] << std::endl
<< "Transform: " << transform << std::endl;
}
// Set the origin of the 2D image
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Software Guide : BeginLatex
//
// For correct (perspective) projection of the 3D volume, the 2D
// image needs to be placed at a certain distance (the focal length
// or source-to-image distance \emph{focalLength}) from the focal
// point, and the normal from the imaging plane to the focal point
// needs to be specified.
//
// By default, ie. if [tx,ty,tz] is [0,0,0], the 3D volume is placed
// half-way between the focal point and the imaging plane. Similarly
// the imaging plane normal is set by default to the center of the
// 2D image but may be modified from this using the command line
// parameters [o2Dx,o2Dy].
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
double origin2D[ Dimension ];
const itk::Vector<double, 3> resolution2D = imageReader2D->GetOutput()->GetSpacing();
typedef ImageType2D::RegionType ImageRegionType2D;
typedef ImageRegionType2D::SizeType SizeType2D;
ImageRegionType2D region2D = caster2D->GetOutput()->GetBufferedRegion();
SizeType2D size2D = region2D.GetSize();
origin2D[0] = origin3D[0] + o2Dx - resolution2D[0]*((double) size2D[0] - 1.)/2.;
origin2D[1] = origin3D[1] + o2Dy - resolution2D[1]*((double) size2D[1] - 1.)/2.;
origin2D[2] = origin3D[2] + focalLength/2.;
imageReader2D->GetOutput()->SetOrigin( origin2D );
registration->SetFixedImageRegion( caster2D->GetOutput()->GetBufferedRegion() );
// Software Guide : EndCodeSnippet
if (verbose)
{
std::cout << "2D image size: "
<< size2D[0] << ", " << size2D[1] << ", " << size2D[2] << std::endl
<< " resolution: "
<< resolution2D[0] << ", " << resolution2D[1] << ", " << resolution2D[2] << std::endl
<< " and position: "
<< origin2D[0] << ", " << origin2D[1] << ", " << origin2D[2] << std::endl;
}
// Set the pyramid schedule for the 2D image
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Software Guide : BeginLatex
//
// The multi-resolution hierachy is specified with respect to the 2D
// image. The number of scales or levels defaults to two, with each level
// differing in resolution by a factor of two. The highest resolution (final)
// level corresponds to the resolution of the input 2D image but
// this can be altered by the user via the \emph{maxScale} command
// line parameter. The number of scales \emph{nScales} can also be
// specified by the user.
//
// Software Guide : EndLatex
std::vector<double> sampledRes2D;
std::vector<ImageRegionType2D> imageRegionPyramid2D;
imagePyramid2D->SetNumberOfLevels( nScales );
sampledRes2D.reserve( nScales );
imageRegionPyramid2D.reserve( nScales );
typedef ImagePyramidType2D::ScheduleType ScheduleType2D;
ScheduleType2D schedule2D = imagePyramid2D->GetSchedule();
for ( unsigned int dim = 0; dim < ImageType2D::ImageDimension; dim++ )
{
schedule2D[ nScales-1 ][ dim ] = maxScale;
}
for ( int level=nScales-2; level >= 0; level-- )
{
for ( unsigned int dim = 0; dim < ImageType2D::ImageDimension; dim++ )
{
schedule2D[ level ][ dim ] = 2*schedule2D[ level+1 ][ dim ];
}
}
// Compute the 2D ImageRegion corresponding to each level of the
// pyramid.
typedef ImageRegionType2D::IndexType IndexType2D;
IndexType2D inputStart2D = region2D.GetIndex();
for ( unsigned int level=0; level < nScales; level++ )
{
SizeType2D size;
IndexType2D start;
sampledRes2D[ level ] = 0.;
for ( unsigned int dim = 0; dim < ImageType2D::ImageDimension; dim++ )
{
float scaleFactor = static_cast<float>( schedule2D[ level ][ dim ] );
size[ dim ] = static_cast<SizeType2D::SizeValueType>(
floor( static_cast<float>( size2D[ dim ] ) / scaleFactor ) );
if( size[ dim ] < 1 )
{
size[ dim ] = 1;
schedule2D[ level ][ dim ] = 1;
}
std::cout << level << " " << dim << " "
<< size[ dim ] << " " << schedule2D[ level ][ dim ]
<< std::endl;
scaleFactor = static_cast<float>( schedule2D[ level ][ dim ] );
start[ dim ] = static_cast<IndexType2D::IndexValueType>(
ceil( static_cast<float>( inputStart2D[ dim ] ) / scaleFactor ) );
if ( dim < 2 )
{
sampledRes2D[ level ] += scaleFactor*resolution2D[ dim ]
*scaleFactor*resolution2D[ dim ];
}
}
sampledRes2D[ level ] = sqrt( sampledRes2D[ level ] );
imageRegionPyramid2D[ level ].SetSize( size );
imageRegionPyramid2D[ level ].SetIndex( start );
}
imagePyramid2D->SetSchedule( schedule2D );
// Set the pyramid schedule for the 3D image
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Software Guide : BeginLatex
//
// The 3D image pyramid naturally contains the same number of levels
// as the multi-resolution schedule for the 2D image. In addition
// the sub-sampling factors for each level are set such that
// resolution of the 3D image in each case is reduced to no greater
// than the corresponding 2D image resolution at that scale. This
// ensures that the 3D volume is reduced in size as far as possible,
// minimising ray-casting computation time, whilst retaining
// sufficient information to ensure accurate production of the DRR
// at the resolution of the 2D image.
//
// Software Guide : EndLatex
std::vector<ImageRegionType3D> imageRegionPyramid3D;
imagePyramid3D->SetNumberOfLevels( nScales );
imageRegionPyramid3D.reserve( nScales );
typedef ImagePyramidType3D::ScheduleType ScheduleType3D;
ScheduleType3D schedule3D = imagePyramid3D->GetSchedule();
// Compute the 2D image pyramid schedule such that the 3D volume
// resolution is no greater than the 2D image resolution.
for ( unsigned int level=0; level < nScales; level++ )
{
for ( unsigned int dim = 0; dim < ImageType3D::ImageDimension; dim++)
{
double res = resolution3D[ dim ];
schedule3D[ level ][ dim ] = 1;
while ( res*2. < sampledRes2D[ level ] )
{
schedule3D[ level ][ dim ] *= 2;
res *= 2.;
}
}
}
// Compute the 3D ImageRegion corresponding to each level of the
// pyramid.
typedef ImageRegionType3D::IndexType IndexType3D;
IndexType3D inputStart3D = region3D.GetIndex();
for ( unsigned int level=0; level < nScales; level++ )
{
SizeType3D size;
IndexType3D start;
for ( unsigned int dim = 0; dim < ImageType3D::ImageDimension; dim++)
{
float scaleFactor = static_cast<float>( schedule3D[ level ][ dim ] );
size[ dim ] = static_cast<SizeType3D::SizeValueType>(
floor( static_cast<float>( size3D[ dim ] ) / scaleFactor ) );
if( size[ dim ] < 1 )
{
size[ dim ] = 1;
schedule3D[ level ][ dim ] = 1;
}
scaleFactor = static_cast<float>( schedule3D[ level ][ dim ] );
start[ dim ] = static_cast<IndexType3D::IndexValueType>(
ceil( static_cast<float>( inputStart3D[ dim ] ) / scaleFactor ) );
}
imageRegionPyramid3D[ level ].SetSize( size );
imageRegionPyramid3D[ level ].SetIndex( start );
}
imagePyramid3D->SetSchedule( schedule3D );
// Initialise the ray cast interpolator
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Software Guide : BeginLatex
//
// The ray cast interpolator is used to project the 3D volume. It
// does this by casting rays from the (transformed) focal point to
// each (transformed) pixel coordinate in the 2D image.
//
// In addition a threshold may be specified to ensure that only
// intensities greater than a given value contribute to the
// projected volume. This can be used, for instance, to remove soft
// tissue from projections of CT data and force the registration
// to find a match which aligns bony structures in the images.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
InterpolatorType::InputPointType focalPoint;
focalPoint[0] = origin3D[0];
focalPoint[1] = origin3D[1];
focalPoint[2] = origin3D[2] - focalLength/2.;
interpolator->SetFocalPoint(focalPoint);
interpolator->SetThreshold(threshold);
interpolator->SetTransform(transform);
// Software Guide : EndCodeSnippet
// NB. Interpolator input image set in itkImageToImageMatric
if (verbose)
{
std::cout << "Focal point: "
<< focalPoint[0] << ", " << focalPoint[1] << ", " << focalPoint[2] << std::endl
<< "Threshold: " << threshold << std::endl;
}
// Set up the transform and start position
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Software Guide : BeginLatex
//
// The registration start position is intialised using the
// transformation parameters.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
registration->SetInitialTransformParameters( transform->GetParameters() );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We wish to maximise the gradient difference similarity measure.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
optimizer->SetMaximize( true );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The maximum and minimum step lengths are specified. The default
// values for these parameters are 4 and 1 (degrees or mm, depending
// on the parameter) respectively, for the coarsest resolution (so
// step sizes of 4, 2 and 1 will be used). For finer scales the
// minimum step size is reduced by a factor of four, so the second
// scale will be performed with step sizes of 1, 0.5 and 0.25 mm or
// degrees.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
optimizer->SetMaximumStepLength( maxStepSize );
optimizer->SetMinimumStepLength( minStepSize );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The maximum number of iterations is set to 200.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
optimizer->SetNumberOfIterations( 200 );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The optimizer weightings are set such that one degree equates to
// one millimeter.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
itk::Optimizer::ScalesType weightings( 6 );
weightings[0] = 180./M_PI;
weightings[1] = 180./M_PI;
weightings[2] = 180./M_PI;
weightings[3] = 1.;
weightings[4] = 1.;
weightings[5] = 1.;
optimizer->SetScales( weightings );
// Software Guide : EndCodeSnippet
if (verbose)
{
optimizer->Print( std::cout );
}
// Create the observers
// ~~~~~~~~~~~~~~~~~~~~
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
// Software Guide : BeginLatex
//
// Once all the registration components are in place we can create
// an instance of the interface command and connect it to the
// registration object using the \code{AddObserver()} method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
CommandType::Pointer command = CommandType::New();
registration->AddObserver( itk::IterationEvent(), command );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Finally we set the number of multi-resolution levels:
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
registration->SetNumberOfLevels( nScales );
// Software Guide : EndCodeSnippet
imagePyramid3D->Print(std::cout);
imagePyramid2D->Print(std::cout);
// Start the registration
// ~~~~~~~~~~~~~~~~~~~~~~
if (verbose)
{
std::cout << "Starting the registration now..." << std::endl;
}
try
{
// Software Guide : BeginLatex
//
// and start the registration.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
registration->StartRegistration();
// Software Guide : EndCodeSnippet
}
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();
const double TranslationAlongX = finalParameters[0];
const double TranslationAlongY = finalParameters[1];
const double TranslationAlongZ = finalParameters[2];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
std::cout << "Result = " << std::endl;
std::cout << " Translation X = " << TranslationAlongX << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << std::endl;
std::cout << " Translation Z = " << TranslationAlongZ << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
// Write out the projection image at the registration position
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
typedef itk::ResampleImageFilter< InternalImageType, InternalImageType > FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetInput( caster3D->GetOutput() );
filter->SetDefaultPixelValue( 0 );
filter->SetTransform( transform );
filter->SetInterpolator( interpolator );
filter->SetSize( imageReader2D->GetOutput()->GetLargestPossibleRegion().GetSize() );
filter->SetOutputOrigin( imageReader2D->GetOutput()->GetOrigin() );
filter->SetOutputSpacing( imageReader2D->GetOutput()->GetSpacing() );
typedef itk::CastImageFilter< InternalImageType, ImageType2D > OutputCastFilterType;
OutputCastFilterType::Pointer outputCaster = OutputCastFilterType::New();
outputCaster->SetInput( filter->GetOutput() );
typedef itk::ImageFileWriter< ImageType2D > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( fileOutput );
writer->SetInput( outputCaster->GetOutput() );
try {
std::cout << "Writing image: " << fileOutput << std::endl;
writer->Update();
}
catch( itk::ExceptionObject & err ) {
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
return 0;
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: brainProjection.hdr
Type: application/octet-stream
Size: 348 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20061022/e8627990/brainProjection-0001.obj
-------------- next part --------------
A non-text attachment was scrubbed...
Name: NM_tomo003.hdr
Type: application/octet-stream
Size: 348 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20061022/e8627990/NM_tomo003-0001.obj
-------------- next part --------------
A non-text attachment was scrubbed...
Name: BPPOS001_DS.hdr
Type: application/octet-stream
Size: 348 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20061022/e8627990/BPPOS001_DS-0002.obj
-------------- next part --------------
A non-text attachment was scrubbed...
Name: BPPOS001_DS.img
Type: application/octet-stream
Size: 131072 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20061022/e8627990/BPPOS001_DS-0003.obj
More information about the Insight-users
mailing list