Hello all,<br><br>As a further test of the modified VersorRigid3DTransform has given confusing results...<br><br>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".<br>
<br>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.<br>
<br>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.<br><br>I've pasted the modified program code before (sorry for the length) - any endeavors to tease out solutions are most welcome!<br>
<br>Thanks,<br><br>Tim Bhatnagar<br><br>#if defined(_MSC_VER)<br>#pragma warning ( disable : 4786 )<br>#endif<br><br>// Software Guide : BeginCommandLineArgs<br>// INPUTS: {brainweb1e1a10f20.mha}<br>// INPUTS: {brainweb1e1a10f20Rot10Tx15.mha}<br>
// ImageRegistration8Output.mhd<br>// ImageRegistration8DifferenceBefore.mhd<br>// ImageRegistration8DifferenceAfter.mhd<br>// OUTPUTS: {ImageRegistration8Output.png}<br>// OUTPUTS: {ImageRegistration8DifferenceBefore.png}<br>
// OUTPUTS: {ImageRegistration8DifferenceAfter.png}<br>// Software Guide : EndCommandLineArgs<br><br>// Software Guide : BeginLatex<br>//<br>// This example illustrates the use of the \doxygen{VersorRigid3DTransform}<br>
// class for performing registration of two $3D$ images. The example code is<br>// for the most part identical to the code presented in<br>// Section~\ref{sec:RigidRegistrationIn2D}. The major difference is that this<br>
// example is done in $3D$. The class \doxygen{CenteredTransformInitializer} is<br>// used to initialize the center and translation of the transform. The case of<br>// rigid registration of 3D images is probably one of the most commonly found<br>
// cases of image registration.<br>//<br>//<br>// \index{itk::VersorRigid3DTransform}<br>// \index{itk::CenteredTransformInitializer!In 3D}<br>//<br>//<br>// Software Guide : EndLatex <br><br>#include "itkImageRegistrationMethod.h"<br>
<br>#include "itkMeanSquaresImageToImageMetric.h"<br>#include "itkMattesMutualInformationImageToImageMetric.h"<br>#include "itkGradientDifferenceImageToImageMetric.h"<br><br>#include "itkLinearInterpolateImageFunction.h"<br>
#include "itkBSplineInterpolateImageFunction.h"<br><br>#include "itkImage.h"<br><br><br>#include "itkVersorRigid3DTransform.h"<br>#include "itkCenteredTransformInitializer.h"<br><br>
<br>// Software Guide : BeginLatex<br>// <br>// The parameter space of the \code{VersorRigid3DTransform} is not a vector<br>// space, due to the fact that addition is not a closed operation in the space<br>// of versor components. This precludes the use of standard gradient descent<br>
// algorithms for optimizing the parameter space of this transform. A special<br>// optimizer should be used in this registration configuration. The optimizer<br>// designed for this transform is the<br>// \doxygen{VersorRigid3DTransformOptimizer}. This optimizer uses Versor<br>
// composition for updating the first three components of the parameters<br>// array, and Vector addition for updating the last three components of the<br>// parameters array~\cite{Hamilton1866,Joly1905}.<br>//<br>// \index{itk::VersorRigid3DTransformOptimizer!header}<br>
// <br>// Software Guide : EndLatex <br><br>// Software Guide : BeginCodeSnippet<br>#include "itkVersorRigid3DTransformOptimizer.h"<br>// Software Guide : EndCodeSnippet<br><br>#include "itkImageFileReader.h"<br>
#include "itkImageFileWriter.h"<br><br>#include "itkResampleImageFilter.h"<br>#include "itkCastImageFilter.h"<br>#include "itkSubtractImageFilter.h"<br>#include "itkRescaleIntensityImageFilter.h"<br>
#include "itkExtractImageFilter.h"<br><br>#include "itkCommand.h"<br>class CommandIterationUpdate : public itk::Command <br>{<br>public:<br> typedef CommandIterationUpdate Self;<br> typedef itk::Command Superclass;<br>
typedef itk::SmartPointer<Self> Pointer;<br> itkNewMacro( Self );<br>protected:<br> CommandIterationUpdate() {};<br>public:<br> typedef itk::VersorRigid3DTransformOptimizer OptimizerType;<br> typedef const OptimizerType * OptimizerPointer;<br>
<br> void Execute(itk::Object *caller, const itk::EventObject & event)<br> {<br> Execute( (const itk::Object *)caller, event);<br> }<br><br> void Execute(const itk::Object * object, const itk::EventObject & event)<br>
{<br> OptimizerPointer optimizer = <br> dynamic_cast< OptimizerPointer >( object );<br> if( ! itk::IterationEvent().CheckEvent( &event ) )<br> {<br> return;<br> }<br> std::cout << optimizer->GetCurrentIteration() << " ";<br>
std::cout << optimizer->GetValue() << " ";<br> std::cout << optimizer->GetCurrentPosition() << std::endl;<br> }<br>};<br><br><br>int main( int argc, char *argv[] )<br>
{<br> if( argc < 5 )<br> {<br> std::cerr << std::endl;<br> std::cerr << "Missing Parameters " << std::endl;<br> std::cerr << "Usage: " << argv[0];<br> std::cerr << " fixedImageFile movingImageFile ";<br>
std::cerr << " outputImagefile translationInTextFileFormat ";<br> std::cerr << " [differenceBeforeRegistration] [differenceAfterRegistration] ";<br> std::cerr << " [step length] "<< std::endl;<br>
std::cerr << std::endl;<br> return 1;<br> }<br> <br> const unsigned int Dimension = 3;<br> typedef float PixelType;<br><br> typedef itk::Image< PixelType, Dimension > FixedImageType;<br>
typedef itk::Image< PixelType, Dimension > MovingImageType;<br><br><br> // Software Guide : BeginCodeSnippet<br> typedef itk::VersorRigid3DTransform< double > TransformType;<br> // Software Guide : EndCodeSnippet<br>
<br><br><br> typedef itk::VersorRigid3DTransformOptimizer OptimizerType;<br><br> /*<br> MeanSquaresImageToImageMetric<br> MattesMutualInformationImageToImageMetric<br> GradientDifferenceImageToImageMetric<br>
*/<br> typedef itk::MattesMutualInformationImageToImageMetric< <br> FixedImageType, <br> MovingImageType > MetricType;<br><br> /*<br>
LinearInterpolateImageFunction<br> BSplineInterpolateImageFunction<br> */<br> typedef itk::LinearInterpolateImageFunction< <br> MovingImageType,<br> double > InterpolatorType;<br>
<br><br> typedef itk::ImageRegistrationMethod< <br> FixedImageType, <br> MovingImageType > RegistrationType;<br><br> MetricType::Pointer metric = MetricType::New();<br>
OptimizerType::Pointer optimizer = OptimizerType::New();<br> InterpolatorType::Pointer interpolator = InterpolatorType::New();<br> RegistrationType::Pointer registration = RegistrationType::New();<br> <br>
<br> registration->SetMetric( metric );<br> registration->SetOptimizer( optimizer );<br> registration->SetInterpolator( interpolator );<br><br><br> // Software Guide : BeginCodeSnippet<br>
TransformType::Pointer transform = TransformType::New();<br> registration->SetTransform( transform );<br> // Software Guide : EndCodeSnippet<br><br> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;<br>
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;<br><br> FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();<br> MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();<br>
<br> fixedImageReader->SetFileName( argv[1] );<br> movingImageReader->SetFileName( argv[2] );<br><br> // ************************************************************************************************************<br>
<br> // MattesMutualInformationImageToImageMetric<br> metric->UseAllPixelsOn();<br> //metric->SetNumberOfSpatialSamples( 100000 ); // default 100000<br> // This is the number of image samples used to calculate the joint probability distribution. The number<br>
// of spatial samples is clamped to be a minimum of 1. Default value is 50.<br> metric->ComputeGradientOn(); // If true, gradient value is computed. Default value is true.<br> // Select whether the metric will be computed using all the pixels on the fixed image region, or only<br>
// using a set of randomly selected pixels.<br> metric->SetNumberOfHistogramBins( 50 ); // Number of bins to used in the histogram. Typical value is 50.<br><br> // GradientDifferenceImageToImageMetric<br> //metric->ComputeGradientOn(); // If true, gradient value is computed. Default value is true.<br>
<br> //interpolator->SetSplineOrder( 3 ); // default (if not specified) is 3<br> // ************************************************************************************************************<br><br> registration->SetFixedImage( fixedImageReader->GetOutput() );<br>
registration->SetMovingImage( movingImageReader->GetOutput() );<br> fixedImageReader->Update();<br><br> registration->SetFixedImageRegion( <br> fixedImageReader->GetOutput()->GetBufferedRegion() );<br>
<br><br> // Software Guide : BeginCodeSnippet<br> typedef itk::CenteredTransformInitializer< TransformType, <br> FixedImageType, <br> MovingImageType <br>
> TransformInitializerType;<br><br> TransformInitializerType::Pointer initializer = <br> TransformInitializerType::New();<br> // Software Guide : EndCodeSnippet<br>
<br><br><br> // Software Guide : BeginCodeSnippet<br> initializer->SetTransform( transform );<br> initializer->SetFixedImage( fixedImageReader->GetOutput() );<br> initializer->SetMovingImage( movingImageReader->GetOutput() );<br>
// Software Guide : EndCodeSnippet<br><br><br> // Software Guide : BeginCodeSnippet<br> initializer->MomentsOn();<br> // Software Guide : EndCodeSnippet<br><br><br> // Software Guide : BeginCodeSnippet<br> initializer->InitializeTransform();<br>
// Software Guide : EndCodeSnippet<br><br> // Software Guide : BeginCodeSnippet<br> typedef TransformType::VersorType VersorType;<br> typedef VersorType::VectorType VectorType;<br><br> VersorType rotation;<br>
VectorType axis;<br> <br> axis[0] = 0.0;<br> axis[1] = 0.0;<br> axis[2] = 0.1;<br><br> const double angle = 0;<br><br> rotation.Set( axis, angle );<br><br> transform->SetRotation( rotation );<br><br> // Software Guide : EndCodeSnippet<br>
<br><br> // Software Guide : BeginLatex<br> // <br> // We now pass the parameters of the current transform as the initial<br> // parameters to be used when the registration process starts. <br> //<br> // Software Guide : EndLatex <br>
<br> // Software Guide : BeginCodeSnippet<br> registration->SetInitialTransformParameters( transform->GetParameters() );<br> // Software Guide : EndCodeSnippet<br><br> typedef OptimizerType::ScalesType OptimizerScalesType;<br>
OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );<br> const double translationScale = 1.0 / 100000.0;<br><br> optimizerScales[0] = 505;<br> optimizerScales[1] = 505;<br> optimizerScales[2] = 505;<br>
optimizerScales[3] = translationScale;<br> optimizerScales[4] = translationScale;<br> optimizerScales[5] = translationScale;<br><br> optimizer->SetScales( optimizerScales );<br><br> double steplength = 0.1; // 0.35 -> 0.1 ------------------------------------------------------------------<br>
<br> if( argc > 7 )<br> {<br> steplength = atof( argv[7] );<br> }<br><br> optimizer->SetMaximumStepLength( steplength ); <br> optimizer->SetMinimumStepLength( 0.0001 );<br><br> double iteration_number= 300; // ------------------------------------------------------------------------<br>
<br> if( argc > 8 )<br> {<br> iteration_number= atof( argv[8] );<br> }<br><br><br><br> optimizer->SetNumberOfIterations( iteration_number);<br><br><br> // Create the Command observer and register it with the optimizer.<br>
//<br> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();<br> optimizer->AddObserver( itk::IterationEvent(), observer );<br><br><br> try <br> {<br> registration->StartRegistration();<br>
} <br> catch( itk::ExceptionObject & err ) <br> { <br> std::cerr << "ExceptionObject caught !" << std::endl; <br> std::cerr << err << std::endl; <br> return -1;<br> } <br>
<br> OptimizerType::ParametersType finalParameters = <br> registration->GetLastTransformParameters();<br><br><br> const double versorX = finalParameters[0];<br> const double versorY = finalParameters[1];<br>
const double versorZ = finalParameters[2];<br> const double finalTranslationX = finalParameters[3];<br> const double finalTranslationY = finalParameters[4];<br> const double finalTranslationZ = finalParameters[5];<br>
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();<br> const double bestValue = optimizer->GetValue();<br><br> // Print out results<br> //<br> std::cout << std::endl << std::endl;<br>
std::cout << "Result = " << std::endl;<br> std::cout << " versor X = " << versorX << std::endl;<br> std::cout << " versor Y = " << versorY << std::endl;<br>
std::cout << " versor Z = " << versorZ << std::endl;<br> std::cout << " Translation X = " << finalTranslationX << std::endl;<br> std::cout << " Translation Y = " << finalTranslationY << std::endl;<br>
std::cout << " Translation Z = " << finalTranslationZ << std::endl;<br> std::cout << " Iterations = " << numberOfIterations << std::endl;<br> std::cout << " Metric value = " << bestValue << std::endl;<br>
<br> // ------------------------------------------------------------------------------------------------<br> std::ofstream translationFile (argv[4]); // set translation text file<br> translationFile << finalTranslationX << " " << finalTranslationY << " " << finalTranslationZ << "\n";<br>
translationFile.close(); // close translation text file<br> // ------------------------------------------------------------------------------------------------<br><br> // Software Guide : BeginCodeSnippet<br> transform->SetParameters( finalParameters );<br>
<br> TransformType::MatrixType matrix = transform->GetRotationMatrix();<br> TransformType::OffsetType offset = transform->GetOffset();<br><br> std::cout << "Matrix = " << std::endl << matrix << std::endl;<br>
std::cout << "Offset = " << std::endl << offset << std::endl;<br> // Software Guide : EndCodeSnippet<br><br><br> typedef itk::ResampleImageFilter< <br> MovingImageType, <br>
FixedImageType > ResampleFilterType;<br><br> TransformType::Pointer finalTransform = TransformType::New();<br><br> finalTransform->SetCenter( transform->GetCenter() );<br><br> finalTransform->SetParameters( finalParameters );<br>
<br> ResampleFilterType::Pointer resampler = ResampleFilterType::New();<br><br> resampler->SetTransform( finalTransform );<br> resampler->SetInput( movingImageReader->GetOutput() );<br><br> FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();<br>
<br> resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );<br> resampler->SetOutputOrigin( fixedImage->GetOrigin() );<br> resampler->SetOutputSpacing( fixedImage->GetSpacing() );<br>
resampler->SetDefaultPixelValue( 0 );<br> <br> typedef float OutputPixelType;<br><br> typedef itk::Image< OutputPixelType, Dimension > OutputImageType;<br> <br> typedef itk::CastImageFilter< <br> FixedImageType,<br>
OutputImageType > CastFilterType;<br> <br> typedef itk::ImageFileWriter< OutputImageType > WriterType;<br><br><br> WriterType::Pointer writer = WriterType::New();<br>
CastFilterType::Pointer caster = CastFilterType::New();<br><br><br> writer->SetFileName( argv[3] );<br><br><br> caster->SetInput( resampler->GetOutput() );<br> writer->SetInput( caster->GetOutput() );<br>
writer->Update();<br><br><br> typedef itk::SubtractImageFilter< <br> FixedImageType, <br> FixedImageType, <br> FixedImageType > DifferenceFilterType;<br>
<br> DifferenceFilterType::Pointer difference = DifferenceFilterType::New();<br><br><br> typedef itk::RescaleIntensityImageFilter< <br> FixedImageType, <br> OutputImageType > RescalerType;<br>
<br> RescalerType::Pointer intensityRescaler = RescalerType::New();<br> <br> intensityRescaler->SetInput( difference->GetOutput() );<br> intensityRescaler->SetOutputMinimum( 0 );<br> intensityRescaler->SetOutputMaximum( 255 );<br>
<br> difference->SetInput1( fixedImageReader->GetOutput() );<br> difference->SetInput2( resampler->GetOutput() );<br><br> resampler->SetDefaultPixelValue( 0 );<br><br> WriterType::Pointer writer2 = WriterType::New();<br>
writer2->SetInput( intensityRescaler->GetOutput() ); <br> <br><br> // Compute the difference image between the <br> // fixed and resampled moving image.<br> if( argc > 6 )<br> {<br> writer2->SetFileName( argv[6] );<br>
writer2->Update();<br> }<br><br><br> typedef itk::IdentityTransform< double, Dimension > IdentityTransformType;<br> IdentityTransformType::Pointer identity = IdentityTransformType::New();<br><br> // Compute the difference image between the <br>
// fixed and moving image before registration.<br> if( argc > 5 )<br> {<br> resampler->SetTransform( identity );<br> writer2->SetFileName( argv[5] );<br> writer2->Update();<br> }<br><br> return 0;<br>
}<br><br><div class="gmail_quote">On Wed, Mar 23, 2011 at 8:10 PM, Tim Bhatnagar <span dir="ltr"><<a href="mailto:tim.bhatnagar@gmail.com">tim.bhatnagar@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
Hello all,<br><br>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.<br>
<br>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.<br>
<br>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?<br>
<br>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.<br>
<br>I can supply code if needed as well as any other pertinent info deemed necessary. I am using the 3.20 ITK release.<br><br>Thanks so much in advance,<br clear="all"><br>-- <br><font color="#888888">Tim Bhatnagar<br>
</font></blockquote></div><br><br clear="all"><br>-- <br>Tim Bhatnagar<br>