Hi, <div><br></div><div>I successfully ran DeformableRegistration 3 (symmetric demon force). However, these codes only registration 2D images. Building upon the 2D example, how can I promote the code to account for 3D registration? </div>
<div><br></div><div>Below is the C code that I copy directly from the example. Thanks Stephen.</div><div><br></div><div><div>/*=========================================================================</div><div><br></div>
<div> Program: Insight Segmentation & Registration Toolkit</div><div> Module: $RCSfile: DeformableRegistration3.cxx,v $</div><div> Language: C++</div><div> Date: $Date: 2007-09-07 14:17:42 $</div><div>
Version: $Revision: 1.18 $</div>
<div><br></div><div> Copyright (c) Insight Software Consortium. All rights reserved.</div><div> See ITKCopyright.txt or <a href="http://www.itk.org/HTML/Copyright.htm">http://www.itk.org/HTML/Copyright.htm</a> for details.</div>
<div><br></div><div> This software is distributed WITHOUT ANY WARRANTY; without even </div><div> the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR </div><div> PURPOSE. See the above copyright notices for more information.</div>
<div><br></div><div>=========================================================================*/</div><div>#if defined(_MSC_VER)</div><div>#pragma warning ( disable : 4786 )</div><div>#endif</div><div><br></div><div>#include "itkImageFileReader.h" </div>
<div>#include "itkImageFileWriter.h" </div><div><br></div><div><br></div><div>// Software Guide : BeginLatex</div><div>//</div><div>// This example demonstrates how to use a variant of the ``demons'' algorithm to</div>
<div>// deformably register two images. This variant uses a different formulation</div><div>// for computing the forces to be applied to the image in order to compute the</div><div>// deformation fields. The variant uses both the gradient of the fixed image</div>
<div>// and the gradient of the deformed moving image in order to compute the</div><div>// forces. This mechanism for computing the forces introduces a symmetry with</div><div>// respect to the choice of the fixed and moving images. This symmetry only</div>
<div>// holds during the computation of one iteration of the PDE updates. It is</div><div>// unlikely that total symmetry may be achieved by this mechanism for the</div><div>// entire registration process.</div><div>//</div>
<div>// The first step for using this filter is to include the following header files.</div><div>//</div><div>// Software Guide : EndLatex</div><div><br></div><div>// Software Guide : BeginCodeSnippet</div><div>#include "itkSymmetricForcesDemonsRegistrationFilter.h"</div>
<div>#include "itkHistogramMatchingImageFilter.h"</div><div>#include "itkCastImageFilter.h"</div><div>#include "itkWarpImageFilter.h"</div><div>#include "itkLinearInterpolateImageFunction.h"</div>
<div>// Software Guide : EndCodeSnippet</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div>// The following section of code implements a Command observer</div><div>// that will monitor the evolution of the registration process.</div>
<div>//</div><div> class CommandIterationUpdate : public itk::Command </div><div> {</div><div> public:</div><div> typedef CommandIterationUpdate Self;</div><div> typedef itk::Command Superclass;</div>
<div> typedef itk::SmartPointer<CommandIterationUpdate> Pointer;</div><div> itkNewMacro( CommandIterationUpdate );</div><div> protected:</div><div> CommandIterationUpdate() {};</div><div><br></div><div>
typedef itk::Image< float, 2 > InternalImageType;</div>
<div> typedef itk::Vector< float, 2 > VectorPixelType;</div><div> typedef itk::Image< VectorPixelType, 2 > DeformationFieldType;</div><div><br></div><div> typedef itk::SymmetricForcesDemonsRegistrationFilter<</div>
<div> InternalImageType,</div><div> InternalImageType,</div><div> DeformationFieldType> RegistrationFilterType;</div><div>
<br>
</div><div> public:</div><div><br></div><div> void Execute(itk::Object *caller, const itk::EventObject & event)</div><div> {</div><div> Execute( (const itk::Object *)caller, event);</div><div> }</div>
<div><br></div><div> void Execute(const itk::Object * object, const itk::EventObject & event)</div><div> {</div><div> const RegistrationFilterType * filter = </div><div> dynamic_cast< const RegistrationFilterType * >( object );</div>
<div> if( !(itk::IterationEvent().CheckEvent( &event )) )</div><div> {</div><div> return;</div><div> }</div><div> std::cout << filter->GetMetric() << std::endl;</div>
<div> }</div><div> };</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div>int main( int argc, char *argv[] )</div><div>{</div><div> if( argc < 4 )</div><div> {</div><div> std::cerr << "Missing Parameters " << std::endl;</div>
<div> std::cerr << "Usage: " << argv[0];</div><div> std::cerr << " fixedImageFile movingImageFile ";</div><div> std::cerr << " outputImageFile " << std::endl;</div>
<div> return EXIT_FAILURE;</div><div> }</div><div><br></div><div> // Software Guide : BeginLatex</div><div> //</div><div> // Second, we declare the types of the images.</div><div> //</div><div> // Software Guide : EndLatex</div>
<div><br></div><div> // Software Guide : BeginCodeSnippet</div><div> const unsigned int Dimension = 2;</div><div> typedef unsigned short PixelType;</div><div><br></div><div> typedef itk::Image< PixelType, Dimension > FixedImageType;</div>
<div> typedef itk::Image< PixelType, Dimension > MovingImageType;</div><div> // Software Guide : EndCodeSnippet</div><div><br></div><div> // Set up the file readers</div><div> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;</div>
<div> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;</div><div><br></div><div> FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();</div><div> MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();</div>
<div><br></div><div> fixedImageReader->SetFileName( argv[1] );</div><div> movingImageReader->SetFileName( argv[2] );</div><div><br></div><div><br></div><div> // Software Guide : BeginLatex</div><div> //</div><div>
// Image file readers are set up in a similar fashion to previous examples.</div><div> // To support the re-mapping of the moving image intensity, we declare an</div><div> // internal image type with a floating point pixel type and cast the input</div>
<div> // images to the internal image type.</div><div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div> // Software Guide : BeginCodeSnippet</div><div> typedef float InternalPixelType;</div><div>
typedef itk::Image< InternalPixelType, Dimension > InternalImageType;</div>
<div> typedef itk::CastImageFilter< FixedImageType, </div><div> InternalImageType > FixedImageCasterType;</div><div> typedef itk::CastImageFilter< MovingImageType, </div><div> InternalImageType > MovingImageCasterType;</div>
<div><br></div><div> FixedImageCasterType::Pointer fixedImageCaster = FixedImageCasterType::New();</div><div> MovingImageCasterType::Pointer movingImageCaster = MovingImageCasterType::New();</div><div><br></div><div>
fixedImageCaster->SetInput( fixedImageReader->GetOutput() );</div>
<div> movingImageCaster->SetInput( movingImageReader->GetOutput() ); </div><div> // Software Guide : EndCodeSnippet</div><div><br></div><div> // Software Guide : BeginLatex</div><div> //</div><div> // The demons algorithm relies on the assumption that pixels representing the</div>
<div> // same homologous point on an object have the same intensity on both the</div><div> // fixed and moving images to be registered. In this example, we will</div><div> // preprocess the moving image to match the intensity between the images</div>
<div> // using the \doxygen{HistogramMatchingImageFilter}. </div><div> //</div><div> // \index{itk::Histogram\-Matching\-Image\-Filter}</div><div> //</div><div> // The basic idea is to match the histograms of the two images at a user-specified number of quantile values. For robustness, the histograms are</div>
<div> // matched so that the background pixels are excluded from both histograms.</div><div> // For MR images, a simple procedure is to exclude all gray values that are</div><div> // smaller than the mean gray value of the image.</div>
<div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div> // Software Guide : BeginCodeSnippet</div><div> typedef itk::HistogramMatchingImageFilter<</div><div> InternalImageType,</div>
<div> InternalImageType > MatchingFilterType;</div><div> MatchingFilterType::Pointer matcher = MatchingFilterType::New();</div><div> // Software Guide : EndCodeSnippet</div><div>
<br>
</div><div><br></div><div> // Software Guide : BeginLatex</div><div> // </div><div> // For this example, we set the moving image as the source or input image and</div><div> // the fixed image as the reference image.</div>
<div> //</div><div> // \index{itk::Histogram\-Matching\-Image\-Filter!SetInput()}</div><div> // \index{itk::Histogram\-Matching\-Image\-Filter!SetSourceImage()}</div><div> // \index{itk::Histogram\-Matching\-Image\-Filter!SetReferenceImage()}</div>
<div> //</div><div> // Software Guide : EndLatex </div><div><br></div><div> // Software Guide : BeginCodeSnippet </div><div> matcher->SetInput( movingImageCaster->GetOutput() );</div><div> matcher->SetReferenceImage( fixedImageCaster->GetOutput() );</div>
<div> // Software Guide : EndCodeSnippet</div><div><br></div><div><br></div><div> // Software Guide : BeginLatex</div><div> //</div><div> // We then select the number of bins to represent the histograms and the</div>
<div>
// number of points or quantile values where the histogram is to be</div><div> // matched.</div><div> //</div><div> // \index{itk::Histogram\-Matching\-Image\-Filter!Set\-Number\-Of\-Histogram\-Levels()}</div><div> // \index{itk::Histogram\-Matching\-Image\-Filter!Set\-Number\-Of\-Match\-Points()}</div>
<div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div> // Software Guide : BeginCodeSnippet</div><div> matcher->SetNumberOfHistogramLevels( 1024 );</div><div> matcher->SetNumberOfMatchPoints( 7 );</div>
<div> // Software Guide : EndCodeSnippet</div><div><br></div><div><br></div><div> // Software Guide : BeginLatex</div><div> //</div><div> // Simple background extraction is done by thresholding at the mean</div><div>
// intensity.</div>
<div> //</div><div> // \index{itk::Histogram\-Matching\-Image\-Filter!Threshold\-At\-Mean\-Intensity\-On()}</div><div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div> // Software Guide : BeginCodeSnippet</div>
<div> matcher->ThresholdAtMeanIntensityOn();</div><div> // Software Guide : EndCodeSnippet</div><div><br></div><div><br></div><div> // Software Guide : BeginLatex</div><div> //</div><div> // In the \doxygen{SymmetricForcesDemonsRegistrationFilter}, the deformation field is</div>
<div> // represented as an image whose pixels are floating point vectors.</div><div> //</div><div> // \index{itk::Symmetric\-Forces\-Demons\-Registration\-Filter}</div><div> //</div><div> // Software Guide : EndLatex</div>
<div><br></div><div> // Software Guide : BeginCodeSnippet</div><div> typedef itk::Vector< float, Dimension > VectorPixelType;</div><div> typedef itk::Image< VectorPixelType, Dimension > DeformationFieldType;</div>
<div> typedef itk::SymmetricForcesDemonsRegistrationFilter<</div><div> InternalImageType,</div><div> InternalImageType,</div><div> DeformationFieldType> RegistrationFilterType;</div>
<div> RegistrationFilterType::Pointer filter = RegistrationFilterType::New();</div><div> // Software Guide : EndCodeSnippet</div><div><br></div><div><br></div><div><br></div><div> // Create the Command observer and register it with the registration filter.</div>
<div> //</div><div> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();</div><div> filter->AddObserver( itk::IterationEvent(), observer );</div><div><br></div><div><br></div><div><br></div><div>
<br></div><div> // Software Guide : BeginLatex</div><div> //</div><div> // The input fixed image is simply the output of the fixed image casting</div><div> // filter. The input moving image is the output of the histogram matching</div>
<div> // filter.</div><div> //</div><div> // \index{itk::Symmetric\-Forces\-Demons\-Registration\-Filter!SetFixedImage()}</div><div> // \index{itk::Symmetric\-Forces\-Demons\-Registration\-Filter!SetMovingImage()}</div>
<div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div> // Software Guide : BeginCodeSnippet</div><div> filter->SetFixedImage( fixedImageCaster->GetOutput() );</div><div> filter->SetMovingImage( matcher->GetOutput() );</div>
<div> // Software Guide : EndCodeSnippet</div><div><br></div><div><br></div><div> // Software Guide : BeginLatex</div><div> //</div><div> // The demons registration filter has two parameters: the number of</div><div>
// iterations to be performed and the standard deviation of the Gaussian</div>
<div> // smoothing kernel to be applied to the deformation field after each</div><div> // iteration.</div><div> // \index{itk::Symmetric\-Forces\-Demons\-Registration\-Filter!SetNumberOfIterations()}</div><div> // \index{itk::Symmetric\-Forces\-Demons\-Registration\-Filter!SetStandardDeviations()}</div>
<div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div> // Software Guide : BeginCodeSnippet</div><div> filter->SetNumberOfIterations( 50 );</div><div> filter->SetStandardDeviations( 1.0 );</div>
<div> // Software Guide : EndCodeSnippet</div><div><br></div><div><br></div><div> // Software Guide : BeginLatex</div><div> // </div><div> // The registration algorithm is triggered by updating the filter. The</div><div>
// filter output is the computed deformation field.</div><div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div> // Software Guide : BeginCodeSnippet</div><div> filter->Update();</div><div> // Software Guide : EndCodeSnippet</div>
<div><br></div><div><br></div><div> // Software Guide : BeginLatex</div><div> //</div><div> // The \doxygen{WarpImageFilter} can be used to warp the moving image with</div><div> // the output deformation field. Like the \doxygen{ResampleImageFilter},</div>
<div> // the WarpImageFilter requires the specification of the input image to be</div><div> // resampled, an input image interpolator, and the output image spacing and</div><div> // origin.</div><div> //</div><div> // \index{itk::WarpImageFilter}</div>
<div> // \index{itk::WarpImageFilter!SetInput()}</div><div> // \index{itk::WarpImageFilter!SetInterpolator()}</div><div> // \index{itk::WarpImageFilter!SetOutputSpacing()}</div><div> // \index{itk::WarpImageFilter!SetOutputOrigin()}</div>
<div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div> // Software Guide : BeginCodeSnippet</div><div> typedef itk::WarpImageFilter<</div><div> MovingImageType, </div><div>
MovingImageType,</div><div> DeformationFieldType > WarperType;</div><div> typedef itk::LinearInterpolateImageFunction<</div><div> MovingImageType,</div>
<div> double > InterpolatorType;</div><div> WarperType::Pointer warper = WarperType::New();</div><div> InterpolatorType::Pointer interpolator = InterpolatorType::New();</div>
<div> FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();</div><div><br></div><div> warper->SetInput( movingImageReader->GetOutput() );</div><div> warper->SetInterpolator( interpolator );</div>
<div> warper->SetOutputSpacing( fixedImage->GetSpacing() );</div><div> warper->SetOutputOrigin( fixedImage->GetOrigin() );</div><div> // Software Guide : EndCodeSnippet</div><div><br></div><div><br></div><div>
// Software Guide : BeginLatex</div><div> //</div><div> // Unlike the ResampleImageFilter, the WarpImageFilter</div><div> // warps or transform the input image with respect to the deformation field</div><div> // represented by an image of vectors. The resulting warped or resampled</div>
<div> // image is written to file as per previous examples.</div><div> //</div><div> // \index{itk::WarpImageFilter!SetDeformationField()}</div><div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div>
// Software Guide : BeginCodeSnippet</div><div> warper->SetDeformationField( filter->GetOutput() );</div><div> // Software Guide : EndCodeSnippet</div><div><br></div><div><br></div><div> // Write warped image out to file</div>
<div> typedef unsigned char OutputPixelType;</div><div> typedef itk::Image< OutputPixelType, Dimension > OutputImageType;</div><div> typedef itk::CastImageFilter< </div><div> MovingImageType,</div>
<div> OutputImageType > CastFilterType;</div><div> typedef itk::ImageFileWriter< OutputImageType > WriterType;</div><div><br></div><div> WriterType::Pointer writer = WriterType::New();</div>
<div> CastFilterType::Pointer caster = CastFilterType::New();</div><div><br></div><div> writer->SetFileName( argv[3] );</div><div> </div><div> caster->SetInput( warper->GetOutput() );</div><div> writer->SetInput( caster->GetOutput() );</div>
<div> writer->Update();</div><div><br></div><div><br></div><div> // Software Guide : BeginLatex</div><div> //</div><div> // Let's execute this example using the rat lung data from the previous example.</div><div>
// The associated data files can be found in \code{Examples/Data}:</div><div> //</div><div> // \begin{itemize}</div><div> // \item \code{RatLungSlice1.mha}</div><div> // \item \code{RatLungSlice2.mha}</div><div> // \end{itemize}</div>
<div> //</div><div> // \begin{figure} \center</div><div> // \includegraphics[width=0.44\textwidth]{DeformableRegistration2CheckerboardBefore.eps}</div><div> // \includegraphics[width=0.44\textwidth]{DeformableRegistration2CheckerboardAfter.eps}</div>
<div> // \itkcaption[Demon's deformable registration output]{Checkerboard comparisons</div><div> // before and after demons-based deformable registration.}</div><div> // \label{fig:DeformableRegistration3Output}</div>
<div> // \end{figure}</div><div> // </div><div> // The result of the demons-based deformable registration is presented in</div><div> // Figure \ref{fig:DeformableRegistration3Output}. The checkerboard</div><div> // comparison shows that the algorithm was able to recover the misalignment</div>
<div> // due to expiration.</div><div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div><br></div><div> // Software Guide : BeginLatex</div><div> //</div><div> // It may be also desirable to write the deformation field as an image of</div>
<div> // vectors. This can be done with the following code.</div><div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div> if( argc > 4 ) // if a fourth line argument has been provided...</div><div>
{</div><div><br></div><div> // Software Guide : BeginCodeSnippet</div><div> typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;</div><div><br></div><div> FieldWriterType::Pointer fieldWriter = FieldWriterType::New();</div>
<div> fieldWriter->SetFileName( argv[4] );</div><div> fieldWriter->SetInput( filter->GetOutput() );</div><div><br></div><div> fieldWriter->Update();</div><div> // Software Guide : EndCodeSnippet</div><div>
<br></div><div> // Software Guide : BeginLatex</div><div> //</div><div> // Note that the file format used for writing the deformation field must be</div><div> // capable of representing multiple components per pixel. This is the case</div>
<div> // for the MetaImage and VTK file formats for example.</div><div> //</div><div> // Software Guide : EndLatex</div><div><br></div><div> }</div><div><br></div><div> return EXIT_SUCCESS;</div><div>}</div><div><br>
</div></div>