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 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703
|
/*=========================================================================
Program: Advanced Normalization Tools
Module: $RCSfile: itkAvantsMutualInformationRegistrationFunction.h,v $
Language: C++
Date: $Date: 2009/01/08 15:14:48 $
Version: $Revision: 1.20 $
Copyright (c) ConsortiumOfANTS. All rights reserved.
See accompanying COPYING.txt or
http://sourceforge.net/projects/advants/files/ANTS/ANTSCopyright.txt 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.
=========================================================================*/
#ifndef __itkAvantsMutualInformationRegistrationFunction_h
#define __itkAvantsMutualInformationRegistrationFunction_h
#include "itkImageFileWriter.h"
#include "itkImageToImageMetric.h"
#include "itkAvantsPDEDeformableRegistrationFunction.h"
#include "itkCovariantVector.h"
#include "itkPoint.h"
#include "itkIndex.h"
#include "itkBSplineKernelFunction.h"
#include "itkBSplineDerivativeKernelFunction.h"
#include "itkCentralDifferenceImageFunction.h"
#include "itkBSplineInterpolateImageFunction.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkBSplineDeformableTransform.h"
#include "itkTranslationTransform.h"
#include "itkArray2D.h"
#include "itkImageBase.h"
#include "itkTransform.h"
#include "itkInterpolateImageFunction.h"
#include "itkSingleValuedCostFunction.h"
#include "itkExceptionObject.h"
#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkSpatialObject.h"
#include "itkConstNeighborhoodIterator.h"
namespace itk
{
/** \class AvantsMutualInformationRegistrationFunction
* \brief Computes the mutual information between two images to be
* registered using the method of Avants et al.
*
* AvantsMutualInformationRegistrationFunction computes the mutual
* information between a fixed and moving image to be registered.
*
* This class is templated over the FixedImage type and the MovingImage
* type.
*
* The fixed and moving images are set via methods SetFixedImage() and
* SetMovingImage(). This metric makes use of user specified Transform and
* Interpolator. The Transform is used to map points from the fixed image to
* the moving image domain. The Interpolator is used to evaluate the image
* intensity at user specified geometric points in the moving image.
* The Transform and Interpolator are set via methods SetTransform() and
* SetInterpolator().
*
* If a BSplineInterpolationFunction is used, this class obtain
* image derivatives from the BSpline interpolator. Otherwise,
* image derivatives are computed using central differencing.
*
* \warning This metric assumes that the moving image has already been
* connected to the interpolator outside of this class.
*
* The method GetValue() computes of the mutual information
* while method GetValueAndDerivative() computes
* both the mutual information and its derivatives with respect to the
* transform parameters.
*
* The calculations are based on the method of Avants et al [1,2]
* where the probability density distribution are estimated using
* Parzen histograms. Since the fixed image PDF does not contribute
* to the derivatives, it does not need to be smooth. Hence,
* a zero order (box car) BSpline kernel is used
* for the fixed image intensity PDF. On the other hand, to ensure
* smoothness a third order BSpline kernel is used for the
* moving image intensity PDF.
*
* On Initialize(), the FixedImage is uniformly sampled within
* the FixedImageRegion. The number of samples used can be set
* via SetNumberOfSpatialSamples(). Typically, the number of
* spatial samples used should increase with the image size.
*
* During each call of GetValue(), GetDerivatives(),
* GetValueAndDerivatives(), marginal and joint intensity PDF's
* values are estimated at discrete position or bins.
* The number of bins used can be set via SetNumberOfHistogramBins().
* To handle data with arbitray magnitude and dynamic range,
* the image intensity is scale such that any contribution to the
* histogram will fall into a valid bin.
*
* One the PDF's have been contructed, the mutual information
* is obtained by doubling summing over the discrete PDF values.
*
*
* Notes:
* 1. This class returns the negative mutual information value.
* 2. This class in not thread safe due the private data structures
* used to the store the sampled points and the marginal and joint pdfs.
*
* References:
* [1] "Nonrigid multimodality image registration"
* D. Avants, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank
* Medical Imaging 2001: Image Processing, 2001, pp. 1609-1620.
* [2] "PET-CT Image Registration in the Chest Using Free-form Deformations"
* D. Avants, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank
* IEEE Transactions in Medical Imaging. Vol.22, No.1,
January 2003. pp.120-128.
* [3] "Optimization of Mutual Information for MultiResolution Image
* Registration"
* P. Thevenaz and M. Unser
* IEEE Transactions in Image Processing, 9(12) December 2000.
*
* \ingroup RegistrationMetrics
* \ingroup ThreadUnSafe
*/
template <class TFixedImage,class TMovingImage , class TDeformationField>
class ITK_EXPORT AvantsMutualInformationRegistrationFunction :
public AvantsPDEDeformableRegistrationFunction< TFixedImage, TMovingImage , TDeformationField>
{
public:
/** Standard class typedefs. */
typedef AvantsMutualInformationRegistrationFunction Self;
typedef AvantsPDEDeformableRegistrationFunction< TFixedImage,
TMovingImage, TDeformationField > Superclass;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro( AvantsMutualInformationRegistrationFunction,
AvantsPDEDeformableRegistrationFunction );
/** MovingImage image type. */
typedef typename Superclass::MovingImageType MovingImageType;
typedef typename Superclass::MovingImagePointer MovingImagePointer;
/** FixedImage image type. */
typedef typename Superclass::FixedImageType FixedImageType;
typedef typename Superclass::FixedImagePointer FixedImagePointer;
typedef typename FixedImageType::IndexType IndexType;
typedef typename FixedImageType::SizeType SizeType;
typedef typename FixedImageType::SpacingType SpacingType;
/** Deformation field type. */
typedef typename Superclass::VectorType VectorType;
typedef typename Superclass::DeformationFieldType DeformationFieldType;
typedef typename Superclass::DeformationFieldTypePointer
DeformationFieldTypePointer;
/** Inherit some enums from the superclass. */
itkStaticConstMacro(ImageDimension, unsigned int,Superclass::ImageDimension);
/** Inherit some enums from the superclass. */
typedef typename Superclass::PixelType PixelType;
typedef typename Superclass::RadiusType RadiusType;
typedef typename Superclass::NeighborhoodType NeighborhoodType;
typedef typename Superclass::FloatOffsetType FloatOffsetType;
typedef typename Superclass::TimeStepType TimeStepType;
/** Interpolator type. */
typedef double CoordRepType;
typedef // // LinearInterpolateImageFunction<MovingImageType,CoordRepType>
BSplineInterpolateImageFunction<MovingImageType,CoordRepType>
InterpolatorType;
typedef typename InterpolatorType::Pointer InterpolatorPointer;
typedef typename InterpolatorType::PointType PointType;
typedef InterpolatorType DefaultInterpolatorType;
// typedef LinearInterpolateImageFunction<MovingImageType,CoordRepType>
//DefaultInterpolatorType;
/** Covariant vector type. */
typedef CovariantVector<double,itkGetStaticConstMacro(ImageDimension)> CovariantVectorType;
/** Gradient calculator type. */
typedef CentralDifferenceImageFunction<FixedImageType> GradientCalculatorType;
typedef typename GradientCalculatorType::Pointer GradientCalculatorPointer;
/** Set the moving image interpolator. */
void SetMovingImageInterpolator( InterpolatorType * ptr )
{ m_MovingImageInterpolator = ptr; }
/** Get the moving image interpolator. */
InterpolatorType * GetMovingImageInterpolator(void)
{ return m_MovingImageInterpolator; }
/** This class uses a constant timestep of 1. */
virtual TimeStepType ComputeGlobalTimeStep(void *itkNotUsed(GlobalData)) const
{ return 1; }
/** Return a pointer to a global data structure that is passed to
* this object from the solver at each calculation. */
virtual void *GetGlobalDataPointer() const
{
GlobalDataStruct *global = new GlobalDataStruct();
// global->m_SumOfSquaredDifference = 0.0;
/// global->m_NumberOfPixelsProcessed = 0L;
// global->m_SumOfSquaredChange = 0;
return global;
}
/** Release memory for global data structure. */
virtual void ReleaseGlobalDataPointer( void *GlobalData ) const
{
delete (GlobalDataStruct *) GlobalData;
}
/** Set the object's state before each iteration. */
virtual void InitializeIteration();
typedef double CoordinateRepresentationType;
/** Types inherited from Superclass. */
typedef TranslationTransform<CoordinateRepresentationType,
// itkGetStaticConstMacro(ImageDimension),
itkGetStaticConstMacro(ImageDimension)> TransformType;
typedef ImageToImageMetric< TFixedImage, TMovingImage > Metricclass;
typedef typename TransformType::Pointer TransformPointer;
typedef typename Metricclass::TransformJacobianType TransformJacobianType;
// typedef typename Metricclass::InterpolatorType InterpolatorType;
typedef typename Metricclass::MeasureType MeasureType;
typedef typename Metricclass::DerivativeType DerivativeType;
typedef typename TransformType::ParametersType ParametersType;
typedef typename Metricclass::FixedImageConstPointer FixedImageConstPointer;
typedef typename Metricclass::MovingImageConstPointer MovingImageCosntPointer;
// typedef typename TransformType::CoordinateRepresentationType CoordinateRepresentationType;
/** Index and Point typedef support. */
typedef typename FixedImageType::IndexType FixedImageIndexType;
typedef typename FixedImageIndexType::IndexValueType FixedImageIndexValueType;
typedef typename MovingImageType::IndexType MovingImageIndexType;
typedef typename TransformType::InputPointType FixedImagePointType;
typedef typename TransformType::OutputPointType MovingImagePointType;
/** The marginal PDFs are stored as std::vector. */
typedef float PDFValueType;
// typedef std::vector<PDFValueType> MarginalPDFType;
typedef Image<PDFValueType,1> MarginalPDFType;
typedef typename MarginalPDFType::IndexType MarginalPDFIndexType;
/** Typedef for the joint PDF and PDF derivatives are stored as ITK Images. */
typedef Image<PDFValueType,2> JointPDFType;
typedef Image<PDFValueType,3> JointPDFDerivativesType;
typedef JointPDFType::IndexType JointPDFIndexType;
typedef JointPDFType::PixelType JointPDFValueType;
typedef JointPDFType::RegionType JointPDFRegionType;
typedef JointPDFType::SizeType JointPDFSizeType;
typedef JointPDFDerivativesType::IndexType JointPDFDerivativesIndexType;
typedef JointPDFDerivativesType::PixelType JointPDFDerivativesValueType;
typedef JointPDFDerivativesType::RegionType JointPDFDerivativesRegionType;
typedef JointPDFDerivativesType::SizeType JointPDFDerivativesSizeType;
/** Get the value and derivatives for single valued optimizers. */
double GetValueAndDerivative( IndexType index,
MeasureType& Value, DerivativeType& Derivative1 , DerivativeType& Derivative2 ) ;
/** Get the value and derivatives for single valued optimizers. */
double GetValueAndDerivativeInv( IndexType index,
MeasureType& Value, DerivativeType& Derivative1 , DerivativeType& Derivative2 ) ;
/** Number of spatial samples to used to compute metric */
// itkSetClampMacro( NumberOfSpatialSamples, unsigned long,
// 1, NumericTraits<unsigned long>::max() );
//itkGetConstReferenceMacro( NumberOfSpatialSamples, unsigned long);
/** Number of bins to used in the histogram. Typical value is 50. */
// itkSetClampMacro( NumberOfHistogramBins, unsigned long,
// 1, NumericTraits<unsigned long>::max() );
// itkGetConstReferenceMacro( NumberOfHistogramBins, unsigned long);
void SetNumberOfHistogramBins(unsigned long nhb) { m_NumberOfHistogramBins=nhb+2*this->m_Padding;}
unsigned long GetNumberOfHistogramBins() {return m_NumberOfHistogramBins;}
void SetTransform(TransformPointer t){m_Transform=t;}
TransformPointer GetTransform(){return m_Transform;}
void SetInterpolator(InterpolatorPointer t){m_Interpolator=t;}
InterpolatorPointer GetInterpolator(){ return m_Interpolator; }
void GetProbabilities();
double ComputeMutualInformation()
{
/** Ray Razlighi changes : May 3, 2010: improves convergence.
1- The padding is lowered to 2.
2- The MI energy is changed to actual MI, please note the value of this MI is in the range of [0 min(H(x),H(y))]. in case you need it to be normalized.
3- The Natural logarithm is changed to log2.
4- In the ComputeMutualInformation() the iterator range has been changed to cover the entire PDF.
5- In the FitIndexInBins() the truncation by floor has been modified to round.
6- The normalization is done based on NomberOfHistogramBins-1 instead of NomberOfHistogramBins. */
float px,py,pxy;
double mival = 0;
double mi;
unsigned long ct = 0;
typename JointPDFType::IndexType index;
for (unsigned int ii=0; ii<m_NumberOfHistogramBins; ii++)
{
MarginalPDFIndexType mind;
mind[0]=ii;
px = m_FixedImageMarginalPDF->GetPixel(mind);
for (unsigned int jj=0; jj<m_NumberOfHistogramBins; jj++)
{
mind[0]=jj;
py = m_MovingImageMarginalPDF->GetPixel(mind);
float denom = px*py;
index[0]=ii;
index[1]=jj;
//pxy=m_JointPDF->GetPixel(index);
JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() +
( ii* m_NumberOfHistogramBins );
// Move the pointer to the first affected bin
int pdfMovingIndex = static_cast<int>( jj );
pdfPtr += pdfMovingIndex;
pxy=*(pdfPtr);
mi=0;
if (fabs(denom) > 0 )
{
if (pxy/denom > 0)
{
//true mi
mi = pxy*log(pxy/denom);
//test mi
//mi = 1.0 + log(pxy/denom);
ct++;
}
}
mival += mi;
}
// std::cout << " II " << ii << " JJ " << ii << " pxy " << pxy << " px " << px << std::endl;
}
this->m_Energy = -1.0*mival/log(2);
return this->m_Energy;
}
virtual VectorType ComputeUpdateInv(const NeighborhoodType &neighborhood,
void *globalData,
const FloatOffsetType &offset = FloatOffsetType(0.0))
{
VectorType update;
update.Fill(0.0);
IndexType oindex = neighborhood.GetIndex();
FixedImageType* img =const_cast<FixedImageType *>(this->Superclass::m_FixedImage.GetPointer());
if (!img) return update;
typename FixedImageType::SpacingType spacing=img->GetSpacing();
typename FixedImageType::SizeType imagesize=img->GetLargestPossibleRegion().GetSize();
for (unsigned int dd=0; dd<ImageDimension; dd++)
{
if ( oindex[dd] < 1 ||
oindex[dd] >= static_cast<typename IndexType::IndexValueType>(imagesize[dd]-1) )
return update;
}
CovariantVectorType fixedGradient;
double loce=0.0;
ParametersType fdvec1(ImageDimension);
ParametersType fdvec2(ImageDimension);
fdvec1.Fill(0);
fdvec2.Fill(0);
fixedGradient = m_FixedImageGradientCalculator->EvaluateAtIndex( oindex );
double nccm1=0;
loce=this->GetValueAndDerivativeInv(oindex,nccm1,fdvec1,fdvec2);
float eps=10;
if ( loce > eps ) loce=eps;
if ( loce < eps*(-1.0) ) loce=eps*(-1.0);
for (int imd=0; imd<ImageDimension; imd++) update[imd]=loce*fixedGradient[imd]*spacing[imd]*(1);
if (this->m_MetricImage) this->m_MetricImage->SetPixel(oindex,loce);
return update;
}
/* Normalizing the image to the range of [0 1] */
double GetMovingParzenTerm( double intensity )
{
double windowTerm = static_cast<double>( intensity ) - this->m_MovingImageTrueMin;
windowTerm = windowTerm / ( this->m_MovingImageTrueMax - this->m_MovingImageTrueMin ) ;
return windowTerm ;
}
double GetFixedParzenTerm( double intensity )
{
double windowTerm = static_cast<double>( intensity ) - this->m_FixedImageTrueMin;
windowTerm = windowTerm / ( this->m_FixedImageTrueMax - this->m_FixedImageTrueMin ) ;
return windowTerm ;
}
/* find the image index in the number of histogram bins */
unsigned int FitIndexInBins( double windowTerm )
{
unsigned int movingImageParzenWindowIndex =
static_cast<unsigned int>( this->m_Padding + round( windowTerm * (float)(this->m_NumberOfHistogramBins - 1 - this->m_Padding))) ;
// Make sure the extreme values are in valid bins
if ( movingImageParzenWindowIndex < this->m_Padding )
{
movingImageParzenWindowIndex = this->m_Padding -1 ;
}
else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - this->m_Padding ) )
{
movingImageParzenWindowIndex = m_NumberOfHistogramBins - this->m_Padding - 1;
}
return movingImageParzenWindowIndex;
}
double FitContIndexInBins( double windowTerm ) {
return ( (double) this->m_Padding + windowTerm*(float)(this->m_NumberOfHistogramBins-this->m_Padding));
}
virtual VectorType ComputeUpdate(const NeighborhoodType &neighborhood,
void *globalData,
const FloatOffsetType &offset = FloatOffsetType(0.0))
{
VectorType update;
update.Fill(0.0);
IndexType oindex = neighborhood.GetIndex();
FixedImageType* img =const_cast<FixedImageType *>(this->Superclass::m_MovingImage.GetPointer());
if (!img) return update;
typename FixedImageType::SpacingType spacing=img->GetSpacing();
typename FixedImageType::SizeType imagesize=img->GetLargestPossibleRegion().GetSize();
for (unsigned int dd=0; dd<ImageDimension; dd++)
{
if ( oindex[dd] < 1 ||
oindex[dd] >= static_cast<typename IndexType::IndexValueType>(imagesize[dd]-1) )
return update;
}
CovariantVectorType movingGradient;
double loce=0.0;
ParametersType fdvec1(ImageDimension);
ParametersType fdvec2(ImageDimension);
fdvec1.Fill(0);
fdvec2.Fill(0);
movingGradient = m_MovingImageGradientCalculator->EvaluateAtIndex( oindex );
double nccm1=0;
loce=this->GetValueAndDerivative(oindex,nccm1,fdvec1,fdvec2);
float eps=10;
if ( loce > eps ) loce=eps;
if ( loce < eps*(-1.0) ) loce=eps*(-1.0);
for (int imd=0; imd<ImageDimension; imd++) update[imd]=loce*movingGradient[imd]*spacing[imd]*(1);
return update;
}
void WriteImages()
{
if (this->m_MetricImage)
{
typedef ImageFileWriter<FixedImageType> writertype;
typename writertype::Pointer w= writertype::New();
w->SetInput( this->m_MetricImage);
w->SetFileName("ZZmetric.nii.gz");
w->Write();
}
}
void SetOpticalFlow(bool b){ m_OpticalFlow = b; }
typename JointPDFType::Pointer GetJointPDF() { return m_JointPDF; }
typename JointPDFType::Pointer GetJointHist() { return m_JointHist; }
void SetFixedImageMask( FixedImageType* img) {m_FixedImageMask=img; }
/** FixedImage image neighborhood iterator type. */
typedef ConstNeighborhoodIterator<FixedImageType> FixedImageNeighborhoodIteratorType;
/** A global data type for this class of equation. Used to store
* iterators for the fixed image. */
struct GlobalDataStruct
{
FixedImageNeighborhoodIteratorType m_FixedImageIterator;
};
/** The fixed image marginal PDF. */
mutable MarginalPDFType::Pointer m_FixedImageMarginalPDF;
/** The moving image marginal PDF. */
mutable MarginalPDFType::Pointer m_MovingImageMarginalPDF;
/** The joint PDF and PDF derivatives. */
typename JointPDFType::Pointer m_JointPDF;
unsigned long m_NumberOfSpatialSamples;
unsigned long m_NumberOfParameters;
/** Variables to define the marginal and joint histograms. */
unsigned long m_NumberOfHistogramBins;
double m_MovingImageNormalizedMin;
double m_FixedImageNormalizedMin;
double m_MovingImageTrueMin;
double m_MovingImageTrueMax;
double m_FixedImageTrueMin;
double m_FixedImageTrueMax;
double m_FixedImageBinSize;
double m_MovingImageBinSize;
protected:
AvantsMutualInformationRegistrationFunction();
virtual ~AvantsMutualInformationRegistrationFunction() {};
void PrintSelf(std::ostream& os, Indent indent) const;
private:
AvantsMutualInformationRegistrationFunction(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
typename JointPDFType::Pointer m_JointHist;
typename JointPDFDerivativesType::Pointer m_JointPDFDerivatives;
/** Typedefs for BSpline kernel and derivative functions. */
typedef BSplineKernelFunction<3> CubicBSplineFunctionType;
typedef BSplineDerivativeKernelFunction<3>
CubicBSplineDerivativeFunctionType;
/** Cubic BSpline kernel for computing Parzen histograms. */
typename CubicBSplineFunctionType::Pointer m_CubicBSplineKernel;
typename CubicBSplineDerivativeFunctionType::Pointer
m_CubicBSplineDerivativeKernel;
/**
* Types and variables related to image derivative calculations.
* If a BSplineInterpolationFunction is used, this class obtain
* image derivatives from the BSpline interpolator. Otherwise,
* image derivatives are computed using central differencing.
*/
typedef CovariantVector< double,
itkGetStaticConstMacro(ImageDimension) > ImageDerivativesType;
/** Boolean to indicate if the interpolator BSpline. */
bool m_InterpolatorIsBSpline;
// boolean to determine if we use mono-modality assumption
bool m_OpticalFlow;
/** Typedefs for using BSpline interpolator. */
typedef
BSplineInterpolateImageFunction<MovingImageType,
CoordinateRepresentationType> BSplineInterpolatorType;
/** Pointer to BSplineInterpolator. */
typename BSplineInterpolatorType::Pointer m_BSplineInterpolator;
/** Typedefs for using central difference calculator. */
typedef CentralDifferenceImageFunction<MovingImageType,
CoordinateRepresentationType> DerivativeFunctionType;
/** Pointer to central difference calculator. */
typename DerivativeFunctionType::Pointer m_DerivativeCalculator;
/**
* Types and variables related to BSpline deformable transforms.
* If the transform is of type third order BSplineDeformableTransform,
* then we can speed up the metric derivative calculation by
* only inspecting the parameters within the support region
* of a mapped point.
*/
/** Boolean to indicate if the transform is BSpline deformable. */
bool m_TransformIsBSpline;
/** The number of BSpline parameters per image dimension. */
long m_NumParametersPerDim;
/**
* The number of BSpline transform weights is the number of
* of parameter in the support region (per dimension ). */
unsigned long m_NumBSplineWeights;
/**
* Enum of the deformabtion field spline order.
*/
enum { DeformationSplineOrder = 3 };
/**
* Typedefs for the BSplineDeformableTransform.
*/
typedef BSplineDeformableTransform<
CoordinateRepresentationType,
::itk::GetImageDimension<FixedImageType>::ImageDimension,
DeformationSplineOrder> BSplineTransformType;
typedef typename BSplineTransformType::WeightsType
BSplineTransformWeightsType;
typedef typename BSplineTransformType::ParameterIndexArrayType
BSplineTransformIndexArrayType;
/**
* Variables used when transform is of type BSpline deformable.
*/
typename BSplineTransformType::Pointer m_BSplineTransform;
/**
* Cache pre-transformed points, weights, indices and
* within support region flag.
*/
typedef typename BSplineTransformWeightsType::ValueType WeightsValueType;
typedef Array2D<WeightsValueType> BSplineTransformWeightsArrayType;
typedef typename BSplineTransformIndexArrayType::ValueType IndexValueType;
typedef Array2D<IndexValueType> BSplineTransformIndicesArrayType;
typedef std::vector<MovingImagePointType> MovingImagePointArrayType;
typedef std::vector<bool> BooleanArrayType;
BSplineTransformWeightsArrayType m_BSplineTransformWeightsArray;
BSplineTransformIndicesArrayType m_BSplineTransformIndicesArray;
MovingImagePointArrayType m_PreTransformPointsArray;
BooleanArrayType m_WithinSupportRegionArray;
typename TFixedImage::SpacingType m_FixedImageSpacing;
typename TFixedImage::PointType m_FixedImageOrigin;
typedef FixedArray<unsigned long,
::itk::GetImageDimension<FixedImageType>::ImageDimension> ParametersOffsetType;
ParametersOffsetType m_ParametersOffset;
mutable TransformPointer m_Transform;
InterpolatorPointer m_Interpolator;
InterpolatorPointer m_FixedImageInterpolator;
InterpolatorPointer m_MovingImageInterpolator;
GradientCalculatorPointer m_FixedImageGradientCalculator;
GradientCalculatorPointer m_MovingImageGradientCalculator;
FixedImagePointer m_FixedImageMask;
MovingImagePointer m_MovingImageMask;
double m_NormalizeMetric;
float m_Normalizer;
typedef BSplineInterpolateImageFunction<JointPDFType,double> pdfintType;
typename pdfintType::Pointer pdfinterpolator;
typedef BSplineInterpolateImageFunction<JointPDFDerivativesType,double> dpdfintType;
typename dpdfintType::Pointer dpdfinterpolator;
typedef BSplineInterpolateImageFunction<MarginalPDFType,double> pdfintType2;
typename pdfintType2::Pointer pdfinterpolator2;
typename pdfintType2::Pointer pdfinterpolator3;
unsigned int m_Padding;
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkAvantsMutualInformationRegistrationFunction.cxx"
#endif
#endif
|