[Insight-users] using VersorRigid3DTransform, failure to create output image? registration failure?
Tim Bhatnagar
tim.bhatnagar at gmail.com
Thu Mar 24 17:35:40 EDT 2011
Hello all,
As a further test of the modified VersorRigid3DTransform has given confusing
results...
I modified my pre-deformation image-stack by rotating (around one axis) and
translating it by defined parameters (nothing very large, and completely
rigid), and used the transformed image as the "MovingImage".
The program was able to compute the translation fairly accurately, but the
rotation matrix still came out as a singularity with the versor indicating
no rotation. AND the output image was filled with pixelValue = 0.
The images I am using are in the analyze format (.hdr and image file),
16-bit unsigned pixel type, 48 images/stack, 256x256 image resolution.
I've pasted the modified program code before (sorry for the length) - any
endeavors to tease out solutions are most welcome!
Thanks,
Tim Bhatnagar
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
// Software Guide : BeginCommandLineArgs
// INPUTS: {brainweb1e1a10f20.mha}
// INPUTS: {brainweb1e1a10f20Rot10Tx15.mha}
// ImageRegistration8Output.mhd
// ImageRegistration8DifferenceBefore.mhd
// ImageRegistration8DifferenceAfter.mhd
// OUTPUTS: {ImageRegistration8Output.png}
// OUTPUTS: {ImageRegistration8DifferenceBefore.png}
// OUTPUTS: {ImageRegistration8DifferenceAfter.png}
// Software Guide : EndCommandLineArgs
// Software Guide : BeginLatex
//
// This example illustrates the use of the \doxygen{VersorRigid3DTransform}
// class for performing registration of two $3D$ images. The example code is
// for the most part identical to the code presented in
// Section~\ref{sec:RigidRegistrationIn2D}. The major difference is that
this
// example is done in $3D$. The class \doxygen{CenteredTransformInitializer}
is
// used to initialize the center and translation of the transform. The case
of
// rigid registration of 3D images is probably one of the most commonly
found
// cases of image registration.
//
//
// \index{itk::VersorRigid3DTransform}
// \index{itk::CenteredTransformInitializer!In 3D}
//
//
// Software Guide : EndLatex
#include "itkImageRegistrationMethod.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkGradientDifferenceImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkBSplineInterpolateImageFunction.h"
#include "itkImage.h"
#include "itkVersorRigid3DTransform.h"
#include "itkCenteredTransformInitializer.h"
// Software Guide : BeginLatex
//
// The parameter space of the \code{VersorRigid3DTransform} is not a vector
// space, due to the fact that addition is not a closed operation in the
space
// of versor components. This precludes the use of standard gradient
descent
// algorithms for optimizing the parameter space of this transform. A
special
// optimizer should be used in this registration configuration. The
optimizer
// designed for this transform is the
// \doxygen{VersorRigid3DTransformOptimizer}. This optimizer uses Versor
// composition for updating the first three components of the parameters
// array, and Vector addition for updating the last three components of the
// parameters array~\cite{Hamilton1866,Joly1905}.
//
// \index{itk::VersorRigid3DTransformOptimizer!header}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkVersorRigid3DTransformOptimizer.h"
// Software Guide : EndCodeSnippet
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkExtractImageFilter.h"
#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() {};
public:
typedef itk::VersorRigid3DTransformOptimizer 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( ! itk::IterationEvent().CheckEvent( &event ) )
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetCurrentPosition() << std::endl;
}
};
int main( int argc, char *argv[] )
{
if( argc < 5 )
{
std::cerr << std::endl;
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " outputImagefile translationInTextFileFormat ";
std::cerr << " [differenceBeforeRegistration]
[differenceAfterRegistration] ";
std::cerr << " [step length] "<< std::endl;
std::cerr << std::endl;
return 1;
}
const unsigned int Dimension = 3;
typedef float PixelType;
typedef itk::Image< PixelType, Dimension > FixedImageType;
typedef itk::Image< PixelType, Dimension > MovingImageType;
// Software Guide : BeginCodeSnippet
typedef itk::VersorRigid3DTransform< double > TransformType;
// Software Guide : EndCodeSnippet
typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
/*
MeanSquaresImageToImageMetric
MattesMutualInformationImageToImageMetric
GradientDifferenceImageToImageMetric
*/
typedef itk::MattesMutualInformationImageToImageMetric<
FixedImageType,
MovingImageType > MetricType;
/*
LinearInterpolateImageFunction
BSplineInterpolateImageFunction
*/
typedef itk::LinearInterpolateImageFunction<
MovingImageType,
double > InterpolatorType;
typedef itk::ImageRegistrationMethod<
FixedImageType,
MovingImageType > RegistrationType;
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetInterpolator( interpolator );
// Software Guide : BeginCodeSnippet
TransformType::Pointer transform = TransformType::New();
registration->SetTransform( transform );
// Software Guide : EndCodeSnippet
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] );
//
************************************************************************************************************
// MattesMutualInformationImageToImageMetric
metric->UseAllPixelsOn();
//metric->SetNumberOfSpatialSamples( 100000 ); // default 100000
// This is the number of image samples used to calculate the joint
probability distribution. The number
// of spatial samples is clamped to be a minimum of 1. Default value is
50.
metric->ComputeGradientOn(); // If true, gradient value is computed.
Default value is true.
// Select whether the metric will be computed using all the pixels on the
fixed image region, or only
// using a set of randomly selected pixels.
metric->SetNumberOfHistogramBins( 50 ); // Number of bins to used in the
histogram. Typical value is 50.
// GradientDifferenceImageToImageMetric
//metric->ComputeGradientOn(); // If true, gradient value is computed.
Default value is true.
//interpolator->SetSplineOrder( 3 ); // default (if not specified) is 3
//
************************************************************************************************************
registration->SetFixedImage( fixedImageReader->GetOutput() );
registration->SetMovingImage( movingImageReader->GetOutput() );
fixedImageReader->Update();
registration->SetFixedImageRegion(
fixedImageReader->GetOutput()->GetBufferedRegion() );
// Software Guide : BeginCodeSnippet
typedef itk::CenteredTransformInitializer< TransformType,
FixedImageType,
MovingImageType
>
TransformInitializerType;
TransformInitializerType::Pointer initializer =
TransformInitializerType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginCodeSnippet
initializer->SetTransform( transform );
initializer->SetFixedImage( fixedImageReader->GetOutput() );
initializer->SetMovingImage( movingImageReader->GetOutput() );
// Software Guide : EndCodeSnippet
// Software Guide : BeginCodeSnippet
initializer->MomentsOn();
// Software Guide : EndCodeSnippet
// Software Guide : BeginCodeSnippet
initializer->InitializeTransform();
// Software Guide : EndCodeSnippet
// Software Guide : BeginCodeSnippet
typedef TransformType::VersorType VersorType;
typedef VersorType::VectorType VectorType;
VersorType rotation;
VectorType axis;
axis[0] = 0.0;
axis[1] = 0.0;
axis[2] = 0.1;
const double angle = 0;
rotation.Set( axis, angle );
transform->SetRotation( rotation );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We now pass the parameters of the current transform as the initial
// parameters to be used when the registration process starts.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
registration->SetInitialTransformParameters( transform->GetParameters() );
// Software Guide : EndCodeSnippet
typedef OptimizerType::ScalesType OptimizerScalesType;
OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
const double translationScale = 1.0 / 100000.0;
optimizerScales[0] = 505;
optimizerScales[1] = 505;
optimizerScales[2] = 505;
optimizerScales[3] = translationScale;
optimizerScales[4] = translationScale;
optimizerScales[5] = translationScale;
optimizer->SetScales( optimizerScales );
double steplength = 0.1; // 0.35 -> 0.1
------------------------------------------------------------------
if( argc > 7 )
{
steplength = atof( argv[7] );
}
optimizer->SetMaximumStepLength( steplength );
optimizer->SetMinimumStepLength( 0.0001 );
double iteration_number= 300; //
------------------------------------------------------------------------
if( argc > 8 )
{
iteration_number= atof( argv[8] );
}
optimizer->SetNumberOfIterations( iteration_number);
// 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::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
OptimizerType::ParametersType finalParameters =
registration->GetLastTransformParameters();
const double versorX = finalParameters[0];
const double versorY = finalParameters[1];
const double versorZ = finalParameters[2];
const double finalTranslationX = finalParameters[3];
const double finalTranslationY = finalParameters[4];
const double finalTranslationZ = finalParameters[5];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
// Print out results
//
std::cout << std::endl << std::endl;
std::cout << "Result = " << std::endl;
std::cout << " versor X = " << versorX << std::endl;
std::cout << " versor Y = " << versorY << std::endl;
std::cout << " versor Z = " << versorZ << std::endl;
std::cout << " Translation X = " << finalTranslationX << std::endl;
std::cout << " Translation Y = " << finalTranslationY << std::endl;
std::cout << " Translation Z = " << finalTranslationZ << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
//
------------------------------------------------------------------------------------------------
std::ofstream translationFile (argv[4]); // set translation text file
translationFile << finalTranslationX << " " << finalTranslationY << " " <<
finalTranslationZ << "\n";
translationFile.close(); // close translation text file
//
------------------------------------------------------------------------------------------------
// Software Guide : BeginCodeSnippet
transform->SetParameters( finalParameters );
TransformType::MatrixType matrix = transform->GetRotationMatrix();
TransformType::OffsetType offset = transform->GetOffset();
std::cout << "Matrix = " << std::endl << matrix << std::endl;
std::cout << "Offset = " << std::endl << offset << std::endl;
// Software Guide : EndCodeSnippet
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
TransformType::Pointer finalTransform = TransformType::New();
finalTransform->SetCenter( transform->GetCenter() );
finalTransform->SetParameters( finalParameters );
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetTransform( finalTransform );
resampler->SetInput( movingImageReader->GetOutput() );
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resampler->SetOutputOrigin( fixedImage->GetOrigin() );
resampler->SetOutputSpacing( fixedImage->GetSpacing() );
resampler->SetDefaultPixelValue( 0 );
typedef float 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( resampler->GetOutput() );
writer->SetInput( caster->GetOutput() );
writer->Update();
typedef itk::SubtractImageFilter<
FixedImageType,
FixedImageType,
FixedImageType > DifferenceFilterType;
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
typedef itk::RescaleIntensityImageFilter<
FixedImageType,
OutputImageType > RescalerType;
RescalerType::Pointer intensityRescaler = RescalerType::New();
intensityRescaler->SetInput( difference->GetOutput() );
intensityRescaler->SetOutputMinimum( 0 );
intensityRescaler->SetOutputMaximum( 255 );
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( resampler->GetOutput() );
resampler->SetDefaultPixelValue( 0 );
WriterType::Pointer writer2 = WriterType::New();
writer2->SetInput( intensityRescaler->GetOutput() );
// Compute the difference image between the
// fixed and resampled moving image.
if( argc > 6 )
{
writer2->SetFileName( argv[6] );
writer2->Update();
}
typedef itk::IdentityTransform< double, Dimension > IdentityTransformType;
IdentityTransformType::Pointer identity = IdentityTransformType::New();
// Compute the difference image between the
// fixed and moving image before registration.
if( argc > 5 )
{
resampler->SetTransform( identity );
writer2->SetFileName( argv[5] );
writer2->Update();
}
return 0;
}
On Wed, Mar 23, 2011 at 8:10 PM, Tim Bhatnagar <tim.bhatnagar at gmail.com>wrote:
> Hello all,
>
> I am using a somewhat modified (and I'm not entirely certain of the
> modifications as I have inherited this program from a colleague) version of
> VersorRigid3DTransform as a first step in a rigid-then-non-rigid
> registration algorithm (the non-rigid portion is accomplished with 3D
> B-spline deformable registration). The two image sets I have are MR slices
> (48 slices per condition) of a live tissue before and after a relatively
> large deformation. I've segmented only the tissue of interest from the
> images.
>
> After my first attempts with the 'Rigid' aspect of this algorithm, I am
> receiving an output image with all pixel values of '0'. I do receive real
> values for the translation vector, but my rotation matrix is essentially an
> identity matrix (which could be feasible, as I'm not certain how much
> rotation is going on during the deformation). Regardless, the aforementioned
> colleague indicated that a 'zeroed' output image is signifying a
> registration failure.
>
> Can anyone give me some sort of indication of the limits of this rigid
> registration program? If my deformation is too large or affine, should I be
> using a different registration program before the 'non-rigid' part of the
> algorithm?
>
> My downfall is that I am really not well-versed in c++ or much of the
> fundamental workings of ITK, but it seems to be a very powerful tool and
> very much appropriate for my study. I am in the process of 'beefing-up' on
> my understanding of the ITK environment, but in the meantime I would greatly
> appreciate being able to move past this hurdle.
>
> I can supply code if needed as well as any other pertinent info deemed
> necessary. I am using the 3.20 ITK release.
>
> Thanks so much in advance,
>
> --
> Tim Bhatnagar
>
--
Tim Bhatnagar
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110324/e7936196/attachment.htm>
More information about the Insight-users
mailing list