1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576
|
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: ImageRegistration2o.cxx,v $
Language: C++
Date: $Date: 2009-06-24 12:09:08 $
Version: $Revision: 1.8 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
// Software Guide : BeginLatex
//
// Some of the most challenging cases of image registration arise when images
// of different modalities are involved. In such cases, metrics based on
// direct comparison of gray levels are not applicable. It has been
// extensively shown that metrics based on the evaluation of mutual
// information is the best way to overcome the difficulties of multi-modality
// registration.
//
// \index{itk::Image\-Registration\-Method!Multi-Modality}
//
// The following simple example illustrates how multiple imaging modalities
// can be registered using the ITK registration framework. The first
// difference is the use of the \doxygen{MutualInformationImageToImageMetric} as
// the cost-function to be optimized and the second difference is the use of
// the \doxygen{GradientDescentOptimizer}. Due to the stochastic nature of the
// metric computation, the values are too noisy to work
// successfully with the \doxygen{RegularStepGradientDescentOptimizer}.
// Therefore, we will use the simpler GradientDescentOptimizer with
// a user defined learning rate. The following headers declare the basic
// components of this registration method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkImageRegistrationMethod.h"
#include "itkTranslationTransform.h"
#include "itkMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkGradientDescentOptimizer.h"
#include "itkOrientedImage.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// One way to simplify the computation of the mutual information is
// to normalize the statistical distribution of the two input images. The
// \doxygen{NormalizeImageFilter} is the perfect tool for this task.
// It rescales the intensities of the input images in order to produce an
// output image with zero mean and unit variance. This filter has been
// discussed in Section \ref{sec:CastingImageFilters}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkNormalizeImageFilter.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Additionally, low-pass filtering of the images to be registered will also
// increase robustness against noise. In this example, we will use the
// \doxygen{DiscreteGaussianImageFilter} for the purpose. The
// characteristics of this filter have been discussed in Section
// \ref{sec:BlurringFilters}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkDiscreteGaussianImageFilter.h"
// Software Guide : EndCodeSnippet
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
// The following section of code implements a Command observer
// that will monitor the evolution of the registration process.
//
#include "itkCommand.h"
class CommandIterationUpdate : public itk::Command
{
public:
typedef CommandIterationUpdate Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
CommandIterationUpdate() {};
public:
typedef itk::GradientDescentOptimizer 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 < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << "outputImagefile [differenceImage]" << std::endl;
return EXIT_FAILURE;
}
// Software Guide : BeginLatex
//
// The moving and fixed images should be instantiated first.
//
// Software Guide : EndLatex
//
// Software Guide : BeginCodeSnippet
const unsigned int Dimension = 2;
typedef unsigned short PixelType;
typedef itk::OrientedImage< PixelType, Dimension > FixedImageType;
typedef itk::OrientedImage< PixelType, Dimension > MovingImageType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// It is convenient to work with an internal image type because mutual
// information will perform better on images with a normalized statistical
// distribution. The fixed and moving images will be normalized and
// converted to this internal type.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef float InternalPixelType;
typedef itk::OrientedImage< InternalPixelType, Dimension > InternalImageType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The rest of the image registration components are instantiated as
// illustrated in Section \ref{sec:IntroductionImageRegistration} with
// the use of the \code{InternalImageType}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::TranslationTransform< double, Dimension > TransformType;
typedef itk::GradientDescentOptimizer OptimizerType;
typedef itk::LinearInterpolateImageFunction<
InternalImageType,
double > InterpolatorType;
typedef itk::ImageRegistrationMethod<
InternalImageType,
InternalImageType > RegistrationType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The mutual information metric type is instantiated using the image
// types.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::MutualInformationImageToImageMetric<
InternalImageType,
InternalImageType > MetricType;
// Software Guide : EndCodeSnippet
TransformType::Pointer transform = TransformType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetOptimizer( optimizer );
registration->SetTransform( transform );
registration->SetInterpolator( interpolator );
// Software Guide : BeginLatex
//
// The metric is created using the \code{New()} method and then
// connected to the registration object.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
MetricType::Pointer metric = MetricType::New();
registration->SetMetric( metric );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The metric requires a number of parameters to be selected, including
// the standard deviation of the Gaussian kernel for the fixed image
// density estimate, the standard deviation of the kernel for the moving
// image density and the number of samples use to compute the densities
// and entropy values. Details on the concepts behind the computation of
// the metric can be found in Section
// \ref{sec:MutualInformationMetric}. Experience has
// shown that a kernel standard deviation of $0.4$ works well for images
// which have been normalized to a mean of zero and unit variance. We
// will follow this empirical rule in this example.
//
// \index{itk::Mutual\-Information\-Image\-To\-Image\-Metric!SetFixedImageStandardDeviation()}
// \index{itk::Mutual\-Information\-Image\-To\-Image\-Metric!SetMovingImageStandardDeviation()}
// \index{itk::Mutual\-Information\-Image\-To\-Image\-Metric!SetNumberOfSpatialSamples()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
metric->SetFixedImageStandardDeviation( 0.4 );
metric->SetMovingImageStandardDeviation( 0.4 );
// 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] );
// Software Guide : BeginLatex
//
// The normalization filters are instantiated using the fixed and moving
// image types as input and the internal image type as output.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::NormalizeImageFilter<
FixedImageType,
InternalImageType
> FixedNormalizeFilterType;
typedef itk::NormalizeImageFilter<
MovingImageType,
InternalImageType
> MovingNormalizeFilterType;
FixedNormalizeFilterType::Pointer fixedNormalizer =
FixedNormalizeFilterType::New();
MovingNormalizeFilterType::Pointer movingNormalizer =
MovingNormalizeFilterType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The blurring filters are declared using the internal image type as both
// the input and output types. In this example, we will set the variance
// for both blurring filters to 2.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::DiscreteGaussianImageFilter<
InternalImageType,
InternalImageType
> GaussianFilterType;
GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
fixedSmoother->SetVariance( 2.0 );
movingSmoother->SetVariance( 2.0 );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The output of the readers becomes the input to the normalization
// filters. The outputs of the normalization filters is connected as
// input to the blurring filters. The input to the registration method
// is taken from the blurring filters.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
fixedNormalizer->SetInput( fixedImageReader->GetOutput() );
movingNormalizer->SetInput( movingImageReader->GetOutput() );
fixedSmoother->SetInput( fixedNormalizer->GetOutput() );
movingSmoother->SetInput( movingNormalizer->GetOutput() );
registration->SetFixedImage( fixedSmoother->GetOutput() );
registration->SetMovingImage( movingSmoother->GetOutput() );
// Software Guide : EndCodeSnippet
fixedNormalizer->Update();
FixedImageType::RegionType fixedImageRegion =
fixedNormalizer->GetOutput()->GetBufferedRegion();
registration->SetFixedImageRegion( fixedImageRegion );
typedef RegistrationType::ParametersType ParametersType;
ParametersType initialParameters( transform->GetNumberOfParameters() );
initialParameters[0] = 0.0; // Initial offset in mm along X
initialParameters[1] = 0.0; // Initial offset in mm along Y
registration->SetInitialTransformParameters( initialParameters );
// Software Guide : BeginCodeSnippet
const unsigned int numberOfPixels = fixedImageRegion.GetNumberOfPixels();
const unsigned int numberOfSamples =
static_cast< unsigned int >( numberOfPixels * 0.01 );
metric->SetNumberOfSpatialSamples( numberOfSamples );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Since larger values of mutual information indicate better matches than
// smaller values, we need to maximize the cost function in this example.
// By default the GradientDescentOptimizer class is set to minimize the
// value of the cost-function. It is therefore necessary to modify its
// default behavior by invoking the \code{MaximizeOn()} method.
// Additionally, we need to define the optimizer's step size using the
// \code{SetLearningRate()} method.
//
// \index{itk::Gradient\-Descent\-Optimizer!MaximizeOn()}
// \index{itk::Image\-Registration\-Method!Maximize vs Minimize}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
optimizer->SetLearningRate( 15.0 );
optimizer->SetNumberOfIterations( 200 );
optimizer->MaximizeOn();
// Software Guide : EndCodeSnippet
// Create the Command observer and register it with the optimizer.
//
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
try
{
registration->StartRegistration();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch( itk::ExceptionObject & err )
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
ParametersType finalParameters = registration->GetLastTransformParameters();
double TranslationAlongX = finalParameters[0];
double TranslationAlongY = finalParameters[1];
unsigned int numberOfIterations = optimizer->GetCurrentIteration();
double bestValue = optimizer->GetValue();
// Print out results
//
std::cout << "Result = " << std::endl;
std::cout << " Translation X = " << TranslationAlongX << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
// Software Guide : BeginLatex
//
// Let's execute this example over two of the images provided in
// \code{Examples/Data}:
//
// \begin{itemize}
// \item \code{BrainT1SliceBorder20.png}
// \item \code{BrainProtonDensitySliceShifted13x17y.png}
// \end{itemize}
//
// \begin{figure}
// \center
// \includegraphics[width=0.44\textwidth]{BrainT1SliceBorder20.eps}
// \includegraphics[width=0.44\textwidth]{BrainProtonDensitySliceShifted13x17y.eps}
// \itkcaption[Multi-Modality Registration Inputs]{A T1 MRI (fixed image) and a proton
// density MRI (moving image) are provided as input to the registration method.}
// \label{fig:FixedMovingImageRegistration2}
// \end{figure}
//
// The second image is the result of intentionally translating the image
// \code{BrainProtonDensitySliceBorder20.png} by $(13,17)$ millimeters. Both
// images have unit-spacing and are shown in Figure
// \ref{fig:FixedMovingImageRegistration2}. The registration is
// stopped at 200 iterations and produces as result the
// parameters:
//
// \begin{verbatim}
// Translation X = 12.8804
// Translation Y = 16.7718
// \end{verbatim}
// These values are approximately within half a pixel of
// the true misaligment introduced in the moving image.
//
// Software Guide : EndLatex
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
TransformType::Pointer finalTransform = TransformType::New();
finalTransform->SetParameters( finalParameters );
finalTransform->SetFixedParameters( transform->GetFixedParameters() );
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform( finalTransform );
resample->SetInput( movingImageReader->GetOutput() );
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resample->SetOutputOrigin( fixedImage->GetOrigin() );
resample->SetOutputSpacing( fixedImage->GetSpacing() );
resample->SetOutputDirection( fixedImage->GetDirection() );
resample->SetDefaultPixelValue( 100 );
typedef unsigned char OutputPixelType;
typedef itk::OrientedImage< OutputPixelType, Dimension > OutputImageType;
typedef itk::CastImageFilter<
FixedImageType,
OutputImageType > CastFilterType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName( argv[3] );
caster->SetInput( resample->GetOutput() );
writer->SetInput( caster->GetOutput() );
writer->Update();
// Software Guide : BeginLatex
//
// \begin{figure}
// \center
// \includegraphics[width=0.32\textwidth]{ImageRegistration2Output.eps}
// \includegraphics[width=0.32\textwidth]{ImageRegistration2CheckerboardBefore.eps}
// \includegraphics[width=0.32\textwidth]{ImageRegistration2CheckerboardAfter.eps}
// \itkcaption[Multi-Modality Registration outputs]{Mapped moving image (left)
// and composition of fixed and moving images before (center) and after
// (right) registration.}
// \label{fig:ImageRegistration2Output}
// \end{figure}
//
// The moving image after resampling is presented on the left
// side of Figure \ref{fig:ImageRegistration2Output}. The center and right
// figures present a checkerboard composite of the fixed and
// moving images before and after registration.
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// \begin{figure}
// \center
// \includegraphics[width=0.44\textwidth]{ImageRegistration2TraceTranslations.eps}
// \includegraphics[width=0.44\textwidth]{ImageRegistration2TraceTranslations2.eps}
// \itkcaption[Multi-Modality Registration plot of translations]{Sequence of
// translations during the registration process. On the left are iterations 0 to
// 200. On the right are iterations 150 to 200.}
// \label{fig:ImageRegistration2TraceTranslations}
// \end{figure}
//
// Figure \ref{fig:ImageRegistration2TraceTranslations} shows the sequence
// of translations followed by the optimizer as it searched the parameter
// space. The left plot shows iterations $0$ to $200$ while the right
// figure zooms into iterations $150$ to $200$. The area covered by the
// right figure has been highlighted by a rectangle in the left image. It
// can be seen that after a certain number of iterations the optimizer
// oscillates within one or two pixels of the true solution. At this
// point it is clear that more iterations will not help. Instead it is
// time to modify some of the parameters of the registration process, for
// example, reducing the learning rate of the optimizer and continuing the
// registration so that smaller steps are taken.
//
// \begin{figure}
// \center
// \includegraphics[width=0.44\textwidth]{ImageRegistration2TraceMetric.eps}
// \includegraphics[width=0.44\textwidth]{ImageRegistration2TraceMetric2.eps}
// \itkcaption[Multi-Modality Registration plot of metrics]{The sequence of metric
// values produced during the registration process. On the left are iterations 0 to 200.
// On the right are iterations 100 to 200.}
// \label{fig:ImageRegistration2TraceMetric}
// \end{figure}
//
// Figure \ref{fig:ImageRegistration2TraceMetric} shows the sequence of
// metric values computed as the optimizer searched the parameter space.
// The left plot shows values when iterations are extended from $0$ to
// $200$ while the right figure zooms into iterations $100$ to $200$. The
// fluctuations in the metric value are due to the stochastic nature in
// which the measure is computed. At each call of \code{GetValue()}, two
// new sets of intensity samples are randomly taken from the image to
// compute the density and entropy estimates. Even with the fluctuations,
// the measure initially increases overall with the number of iterations.
// After about 150 iterations, the metric value merely oscillates without further
// noticeable convergence. The trace plots in Figure
// \ref{fig:ImageRegistration2TraceMetric} highlight one of the
// difficulties associated with this particular metric: the stochastic
// oscillations make it difficult to determine convergence and limit the
// use of more sophisticated optimization methods. As explained above,
// the reduction of the learning rate as the registration progresses is
// very important in order to get precise results.
//
// This example shows the importance of tracking the evolution of the
// registration method in order to obtain insight into the characteristics
// of the particular problem at hand and the components being used. The
// behavior revealed by these plots usually helps to identify possible
// improvements in the setup of the registration parameters.
//
// Software Guide : EndLatex
return EXIT_SUCCESS;
}
|