[Insight-developers] Mattes mutual information metric additions
Karthik Krishnan
Karthik.Krishnan at kitware.com
Tue May 30 14:04:19 EDT 2006
Hi,
We've had PET-CT datasets, (the PET ones in particular) with intensity
distributions where certain regions have very high pixel intensities
that tend to throw the registration off course. I added some code to the
mattes that optionally allows the joint PDF computation of the metric to
be restricted to a fixed range of intensities and it seems to work well.
The API is simple. There are 4 optional methods :
SetFixedImageLowerBound( ), SetMovingImageLowerBound( ),
SetFixedImageUpperBound( ), SetMovingImageUpperBound( ).
I've attached the patch and I would like to commit it, but I want to
confirm with the author or others in case you think its just code
clutter and I should tailor a derived class for my needs.
Thanks
-karthik
PS: Just to clarify, setting the bounds is different from using a
rescale intensity
filter on the fixed and moving images prior to computing the metric.
Rescale intensity filters clamp the intensities at the bounds. Therefore,
you will have an intensity peak at the bounds that also contributes to
estimation of the joint intensity PDF's and distorts the histogram.
-------------- next part --------------
Index: itkMattesMutualInformationImageToImageMetric.h
===================================================================
RCS file: /cvsroot/Insight/Insight/Code/Algorithms/itkMattesMutualInformationImageToImageMetric.h,v
retrieving revision 1.16
diff -u -r1.16 itkMattesMutualInformationImageToImageMetric.h
--- itkMattesMutualInformationImageToImageMetric.h 31 Oct 2005 22:14:35 -0000 1.16
+++ itkMattesMutualInformationImageToImageMetric.h 30 May 2006 17:47:17 -0000
@@ -37,10 +37,6 @@
*
* MattesMutualInformationImageToImageMetric 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
@@ -49,18 +45,20 @@
* 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.
+ * \par Template parameters
+ * This class is templated over the FixedImage type and the MovingImage
+ * type.
*
* \warning This metric assumes that the moving image has already been
* connected to the interpolator outside of this class.
*
+ * \par
* 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.
- *
+ *
+ * \par
* The calculations are based on the method of Mattes et al [1,2]
* where the probability density distribution are estimated using
* Parzen histograms. Since the fixed image PDF does not contribute
@@ -70,42 +68,44 @@
* 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.
- *
- * The option UseAllPixelOn() disables the random sampling and uses
- * all the pixels of the FixedImageRegion in order to estimate the
- * joint intensity PDF.
+ * \par Paramters and Options
+ * \li SetNumberOfSpatialSamples()
+ * \li UseAllPixelsOn()
+ * \li SetNumberOfHistogramBins()
+ * \li SetFixedImageLowerBound()
+ * \li SetFixedImageUpperBound()
+ * \li SetMovingImageLowerBound()
+ * \li SetMovingImageUpperBound()
*
- * During each call of GetValue(), GetDerivatives(),
- * GetValueAndDerivatives(), marginal and joint intensity PDF's
+ * \par
+ * During each call of GetValue(), GetDerivative(),
+ * GetValueAndDerivative(), 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
+ * Once 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.
+ * \par Caveats
+ * \li If a BSplineInterpolationFunction is used, this class obtain
+ * image derivatives from the BSpline interpolator. Otherwise,
+ * image derivatives are computed using central differencing.
+ * \li This class returns the negative mutual information value.
+ * \li 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"
+ * \li "Nonrigid multimodality image registration"
* D. Mattes, 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"
+ * \li "PET-CT Image Registration in the Chest Using Free-form Deformations"
* D. Mattes, 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"
+ * \li "Optimization of Mutual Information for MultiResolution Image
+ * Registration"
* P. Thevenaz and M. Unser
* IEEE Transactions in Image Processing, 9(12) December 2000.
*
@@ -177,7 +177,13 @@
void GetValueAndDerivative( const ParametersType& parameters,
MeasureType& Value, DerivativeType& Derivative ) const;
- /** Number of spatial samples to used to compute metric */
+ /** Number of spatial samples to used to compute metric. On Initialize(),
+ * the FixedImage is uniformly sampled within the FixedImageRegion. The
+ * number of samples used is set via SetNumberOfSpatialSamples().
+ * Typically, the number of spatial samples used should increase with the
+ * image size. A reasonable default is to use 5% of the number of pixels in
+ * the fixed image to estimate the joint intensity PDF. If the option
+ * UseAllPixelsOn() is used, this option has no effect. */
itkSetClampMacro( NumberOfSpatialSamples, unsigned long,
1, NumericTraits<unsigned long>::max() );
itkGetConstReferenceMacro( NumberOfSpatialSamples, unsigned long);
@@ -192,11 +198,21 @@
void ReinitializeSeed(int);
/** 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. */
+ * fixed image region, or only using a set of randomly selected pixels. Its
+ * off by default. */
itkSetMacro(UseAllPixels,bool);
itkGetConstReferenceMacro(UseAllPixels,bool);
itkBooleanMacro(UseAllPixels);
+ /** You might get fixed and moving images with spurious intensity peaks.
+ * Optionally, you may set the lower bounds of the intensities to be
+ * considered for computing the histogram. This option allows to focus the
+ * computation of the Metric in a particular range of intensities that
+ * correspond to features of interest. */
+ void SetFixedImageLowerBound( typename FixedImageType::PixelType );
+ void SetFixedImageUpperBound( typename FixedImageType::PixelType );
+ void SetMovingImageLowerBound( typename MovingImageType::PixelType );
+ void SetMovingImageUpperBound( typename MovingImageType::PixelType );
protected:
@@ -409,7 +425,18 @@
bool m_ReseedIterator;
int m_RandomSeed;
-
+
+ // See documentation of SetFixedImageLowerBound()
+ bool IsWithinFixedImageIntensityBounds( typename FixedImageType::PixelType ) const;
+ bool IsWithinMovingImageIntensityBounds( double ) const;
+ bool m_FixedImageLowerBoundSpecified;
+ bool m_FixedImageUpperBoundSpecified;
+ bool m_MovingImageLowerBoundSpecified;
+ bool m_MovingImageUpperBoundSpecified;
+ typename FixedImageType::PixelType m_FixedImageLowerBound;
+ typename MovingImageType::PixelType m_MovingImageLowerBound;
+ typename FixedImageType::PixelType m_FixedImageUpperBound;
+ typename MovingImageType::PixelType m_MovingImageUpperBound;
};
} // end namespace itk
Index: itkMattesMutualInformationImageToImageMetric.txx
===================================================================
RCS file: /cvsroot/Insight/Insight/Code/Algorithms/itkMattesMutualInformationImageToImageMetric.txx,v
retrieving revision 1.36
diff -u -r1.36 itkMattesMutualInformationImageToImageMetric.txx
--- itkMattesMutualInformationImageToImageMetric.txx 28 Mar 2006 23:47:09 -0000 1.36
+++ itkMattesMutualInformationImageToImageMetric.txx 30 May 2006 17:47:17 -0000
@@ -74,6 +74,12 @@
m_UseAllPixels = false;
m_ReseedIterator = false;
m_RandomSeed = -1;
+
+ // Optional option to specify bounds.. turn it off
+ this->m_FixedImageLowerBoundSpecified = false;
+ this->m_FixedImageUpperBoundSpecified = false;
+ this->m_MovingImageUpperBoundSpecified = false;
+ this->m_MovingImageUpperBoundSpecified = false;
}
@@ -461,7 +467,16 @@
}
// Get sampled fixed image value
- (*iter).FixedImageValue = randIter.Get();
+ typename FixedImageType::PixelType p = randIter.Get();
+
+ // Check if its within any intensity bounds, if specified
+ if (!this->IsWithinFixedImageIntensityBounds( p ))
+ {
+ ++randIter; // jump to another random position
+ continue;
+ }
+ (*iter).FixedImageValue = p;
+
// Translate index to point
(*iter).FixedImagePointValue = inputPoint;
@@ -476,8 +491,18 @@
{
// Get sampled index
FixedImageIndexType index = randIter.GetIndex();
+
// Get sampled fixed image value
- (*iter).FixedImageValue = randIter.Get();
+ typename FixedImageType::PixelType p = randIter.Get();
+
+ // Check if its within any intensity bounds, if specified
+ if (!this->IsWithinFixedImageIntensityBounds( p ))
+ {
+ ++randIter; // jump to another random position
+ continue;
+ }
+ (*iter).FixedImageValue = p;
+
// Translate index to point
this->m_FixedImage->TransformIndexToPhysicalPoint( index,
(*iter).FixedImagePointValue );
@@ -530,7 +555,15 @@
}
// Get sampled fixed image value
- (*iter).FixedImageValue = regionIter.Get();
+ typename FixedImageType::PixelType p = regionIter.Get();
+ if (!this->IsWithinFixedImageIntensityBounds( p ))
+ {
+ ++regionIter;
+ continue;
+ }
+
+ (*iter).FixedImageValue = p;
+
// Translate index to point
(*iter).FixedImagePointValue = inputPoint;
@@ -545,8 +578,17 @@
{
// Get sampled index
FixedImageIndexType index = regionIter.GetIndex();
+
// Get sampled fixed image value
- (*iter).FixedImageValue = regionIter.Get();
+ typename FixedImageType::PixelType p = regionIter.Get();
+ if (!this->IsWithinFixedImageIntensityBounds( p ))
+ {
+ ++regionIter;
+ continue;
+ }
+
+ (*iter).FixedImageValue = p;
+
// Translate index to point
this->m_FixedImage->TransformIndexToPhysicalPoint( index,
(*iter).FixedImagePointValue );
@@ -1266,6 +1308,8 @@
sampleOk = false;
}
}
+
+ sampleOk &= this->IsWithinMovingImageIntensityBounds( movingImageValue );
}
@@ -1405,12 +1449,93 @@
m_WithinSupportRegionArray[counter] = valid;
}
+}
+// Set the lower bound. This is optional. It performs the same functionality
+// as in itk::HistogramImageToImageMetric. The boounds are inclusive. All
+// this does is to prevent samples that lie outside the intensity bounds from
+// contributing to the joint intensity PDF's.
+//
+// Note that setting the bounds is different from using a rescale intensity
+// filter on the fixed and moving images prior to computing the metric.
+// Rescale intensity filters clamp the intensities at the bounds. Therefore,
+// you will have an intensity peak at the bounds that also contributes to
+// estimation of the joint intensity PDF's and distorts the histogram.
+//
+template < class TFixedImage, class TMovingImage >
+void MattesMutualInformationImageToImageMetric<TFixedImage,TMovingImage>
+::SetFixedImageLowerBound( typename FixedImageType::PixelType bound )
+{
+ this->m_FixedImageLowerBoundSpecified = true;
+ this->m_FixedImageLowerBound = bound;
+ this->Modified();
}
+template < class TFixedImage, class TMovingImage >
+void MattesMutualInformationImageToImageMetric<TFixedImage,TMovingImage>
+::SetMovingImageLowerBound( typename MovingImageType::PixelType bound )
+{
+ this->m_MovingImageLowerBoundSpecified = true;
+ this->m_MovingImageLowerBound = bound;
+ this->Modified();
+}
-} // end namespace itk
+template < class TFixedImage, class TMovingImage >
+void MattesMutualInformationImageToImageMetric<TFixedImage,TMovingImage>
+::SetFixedImageUpperBound( typename FixedImageType::PixelType bound )
+{
+ this->m_FixedImageUpperBoundSpecified = true;
+ this->m_FixedImageUpperBound = bound;
+ this->Modified();
+}
+
+template < class TFixedImage, class TMovingImage >
+void MattesMutualInformationImageToImageMetric<TFixedImage,TMovingImage>
+::SetMovingImageUpperBound( typename MovingImageType::PixelType bound )
+{
+ this->m_MovingImageUpperBoundSpecified = true;
+ this->m_MovingImageUpperBound = bound;
+ this->Modified();
+}
+
+template < class TFixedImage, class TMovingImage >
+bool MattesMutualInformationImageToImageMetric<TFixedImage,TMovingImage>
+::IsWithinFixedImageIntensityBounds(
+ typename FixedImageType::PixelType p ) const
+{
+ bool b = true;
+ if (this->m_FixedImageUpperBoundSpecified
+ && (p > this->m_FixedImageUpperBound))
+ {
+ b = false;
+ }
+ if (this->m_FixedImageLowerBoundSpecified
+ && (p < this->m_FixedImageLowerBound))
+ {
+ b = false;
+ }
+ return b;
+}
+
+template < class TFixedImage, class TMovingImage >
+bool MattesMutualInformationImageToImageMetric<TFixedImage,TMovingImage>
+::IsWithinMovingImageIntensityBounds( double p ) const
+{
+ bool b = true;
+ if (this->m_MovingImageUpperBoundSpecified
+ && (p > this->m_MovingImageUpperBound))
+ {
+ b = false;
+ }
+ if (this->m_MovingImageLowerBoundSpecified
+ && (p < this->m_MovingImageLowerBound))
+ {
+ b = false;
+ }
+ return b;
+}
+} // end namespace itk
#endif
More information about the Insight-developers
mailing list