ITK  6.0.0
Insight Toolkit
Examples/Statistics/ImageMutualInformation1.cxx
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
// Software Guide : BeginLatex
//
// This example illustrates how to compute the Mutual Information between two
// images using classes from the Statistics framework. Note that you could
// also use for this purpose the ImageMetrics designed for the image
// registration framework.
//
// For example, you could use:
//
// \begin{itemize}
// \item \doxygen{MutualInformationImageToImageMetric}
// \item \doxygen{MattesMutualInformationImageToImageMetric}
// \item \doxygen{MutualInformationHistogramImageToImageMetric}
// \item \doxygen{MutualInformationImageToImageMetric}
// \item \doxygen{NormalizedMutualInformationHistogramImageToImageMetric}
// \item \doxygen{KullbackLeiblerCompareHistogramImageToImageMetric}
// \end{itemize}
//
// Mutual Information as computed in this example, and as commonly used in the
// context of image registration provides a measure of how much uncertainty on
// the value of a pixel in one image is reduced by measuring the homologous
// pixel in the other image. Note that Mutual Information as used here does
// not measure the amount of information that one image provides on the other
// image; this would require us to take into account the spatial
// structures in the images as well as the semantics of the image context in
// terms of an observer.
//
// This implies that there is still an enormous unexploited potential on the
// use of the Mutual Information concept in the domain of medical images,
// among the most interesting of which is the semantic description of
// image in terms of anatomical structures.
//
// \index{Mutual Information!Statistics}
// \index{Statistics!Mutual Information}
// \index{Joint Entropy!Statistics}
// \index{Statistics!Joint Entropy}
// \index{Joint Histogram!Statistics}
// \index{Statistics!Joint Histogram}
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// In this particular example we make use of classes from the Statistics
// framework in order to compute the measure of Mutual Information between two
// images. We assume that both images have the same number of pixels along
// every dimension and that they have the same origin and spacing. Therefore
// the pixels from one image are perfectly aligned with those of the other
// image.
//
// We must start by including the header files of the image, histogram
// filter, reader and Join image filter. We will read both images and use
// the Join image filter in order to compose an image of two components using
// the information of each one of the input images in one component. This is
// the natural way of using the Statistics framework in ITK given that the
// fundamental statistical classes are expecting to receive multi-valued
// measures.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkImage.h"
// Software Guide : EndCodeSnippet
int
main(int argc, char * argv[])
{
if (argc < 3)
{
std::cerr << "Missing command line arguments" << std::endl;
std::cerr << "Usage : ImageMutualInformation1 inputImage1 inputImage2 "
<< std::endl;
return EXIT_FAILURE;
}
// Software Guide : BeginLatex
//
// We define the pixel type and dimension of the images to be read.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using PixelComponentType = unsigned char;
constexpr unsigned int Dimension = 2;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Using the image type we proceed to instantiate the readers for both input
// images. Then, we take their filenames from the command line arguments.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using ReaderType = itk::ImageFileReader<ImageType>;
auto reader1 = ReaderType::New();
auto reader2 = ReaderType::New();
reader1->SetFileName(argv[1]);
reader2->SetFileName(argv[2]);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Using the \doxygen{JoinImageFilter} we use the two input images and put
// them together in an image of two components.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
auto joinFilter = JoinFilterType::New();
joinFilter->SetInput1(reader1->GetOutput());
joinFilter->SetInput2(reader2->GetOutput());
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// At this point we trigger the execution of the pipeline by invoking the
// \code{Update()} method on the Join filter. We must put the call inside a
// try/catch block because the Update() call may potentially result in
// exceptions being thrown.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
try
{
joinFilter->Update();
}
catch (const itk::ExceptionObject & excp)
{
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We now prepare the types to be used for the computation of the joint
// histogram. For this purpose, we take the type of the image resulting from
// the JoinImageFilter and use it as template argument of the
// \doxygen{ImageToHistogramFilter}. We then construct one by invoking the
// \code{New()} method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using VectorImageType = JoinFilterType::OutputImageType;
using HistogramFilterType =
auto histogramFilter = HistogramFilterType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We pass the multiple-component image as input to the histogram filter,
// and setup the marginal scale value that will define the precision to be
// used for classifying values into the histogram bins.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
histogramFilter->SetInput(joinFilter->GetOutput());
histogramFilter->SetMarginalScale(10.0);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We must now define the number of bins to use for each one of the
// components in the joint image. For this purpose we take the
// \code{HistogramSizeType} from the traits of the histogram filter type.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using HistogramSizeType = HistogramFilterType::HistogramSizeType;
HistogramSizeType size(2);
size[0] = 255; // number of bins for the first channel
size[1] = 255; // number of bins for the second channel
histogramFilter->SetHistogramSize(size);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Finally, we must specify the upper and lower bounds for the histogram
// using the \code{SetHistogramBinMinimum()} and
// \code{SetHistogramBinMaximum()} methods. The \code{Update()} method is
// then called in order to trigger the computation of the histogram.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using HistogramMeasurementVectorType =
HistogramFilterType::HistogramMeasurementVectorType;
HistogramMeasurementVectorType binMinimum(3);
HistogramMeasurementVectorType binMaximum(3);
binMinimum[0] = -0.5;
binMinimum[1] = -0.5;
binMinimum[2] = -0.5;
binMaximum[0] = 255.5;
binMaximum[1] = 255.5;
binMaximum[2] = 255.5;
histogramFilter->SetHistogramBinMinimum(binMinimum);
histogramFilter->SetHistogramBinMaximum(binMaximum);
histogramFilter->Update();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The histogram can be recovered from the filter by creating a variable
// with the histogram type taken from the filter traits.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using HistogramType = HistogramFilterType::HistogramType;
const HistogramType * histogram = histogramFilter->GetOutput();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We now walk over all the bins of the joint histogram and compute their
// contribution to the value of the joint entropy. For this purpose we use
// histogram iterators, and the \code{Begin()} and \code{End()} methods.
// Since the values returned from the histogram are measuring frequency we
// must convert them to an estimation of probability by dividing them over
// the total sum of frequencies returned by the \code{GetTotalFrequency()}
// method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
HistogramType::ConstIterator itr = histogram->Begin();
HistogramType::ConstIterator end = histogram->End();
const double Sum = histogram->GetTotalFrequency();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We initialize to zero the variable to use for accumulating the value of
// the joint entropy, and then use the iterator for visiting all the bins of
// the joint histogram. For every bin we compute their contribution to the
// reduction of uncertainty. Note that in order to avoid logarithmic
// operations on zero values, we skip over those bins that have less than
// one count. The entropy contribution must be computed using logarithms in
// base two in order to express entropy in \textbf{bits}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
double JointEntropy = 0.0;
while (itr != end)
{
const double count = itr.GetFrequency();
if (count > 0.0)
{
const double probability = count / Sum;
JointEntropy += -probability * std::log(probability) / std::log(2.0);
}
++itr;
}
// Software Guide : EndCodeSnippet
std::cout << "Joint Entropy = " << JointEntropy << " bits "
<< std::endl;
// Software Guide : BeginLatex
//
// Now that we have the value of the joint entropy we can proceed to
// estimate the values of the entropies for each image independently. This
// can be done by simply changing the number of bins and then recomputing
// the histogram.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
size[0] = 255; // number of bins for the first channel
size[1] = 1; // number of bins for the second channel
histogramFilter->SetHistogramSize(size);
histogramFilter->Update();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We initialize to zero another variable in order to start accumulating the
// entropy contributions from every bin.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
itr = histogram->Begin();
end = histogram->End();
double Entropy1 = 0.0;
while (itr != end)
{
const double count = itr.GetFrequency();
if (count > 0.0)
{
const double probability = count / Sum;
Entropy1 += -probability * std::log(probability) / std::log(2.0);
}
++itr;
}
// Software Guide : EndCodeSnippet
std::cout << "Image1 Entropy = " << Entropy1 << " bits " << std::endl;
// Software Guide : BeginLatex
//
// The same process is used for computing the entropy of the other
// component, simply by swapping the number of bins in the histogram.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
size[0] = 1; // number of bins for the first channel
size[1] = 255; // number of bins for the second channel
histogramFilter->SetHistogramSize(size);
histogramFilter->Update();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The entropy is computed in a similar manner, just by visiting all the
// bins on the histogram and accumulating their entropy contributions.
//
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
itr = histogram->Begin();
end = histogram->End();
double Entropy2 = 0.0;
while (itr != end)
{
const double count = itr.GetFrequency();
if (count > 0.0)
{
const double probability = count / Sum;
Entropy2 += -probability * std::log(probability) / std::log(2.0);
}
++itr;
}
// Software Guide : EndCodeSnippet
std::cout << "Image2 Entropy = " << Entropy2 << " bits " << std::endl;
// Software Guide : BeginLatex
//
// At this point we can compute any of the popular measures of Mutual
// Information. For example
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
double MutualInformation = Entropy1 + Entropy2 - JointEntropy;
// Software Guide : EndCodeSnippet
std::cout << "Mutual Information = " << MutualInformation << " bits "
<< std::endl;
// Software Guide : BeginLatex
//
// or Normalized Mutual Information, where the value of Mutual Information
// is divided by the mean entropy of the input images.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
double NormalizedMutualInformation1 =
2.0 * MutualInformation / (Entropy1 + Entropy2);
// Software Guide : EndCodeSnippet
std::cout << "Normalized Mutual Information 1 = "
<< NormalizedMutualInformation1 << std::endl;
// Software Guide : BeginLatex
//
// A second form of Normalized Mutual Information has been defined as the
// mean entropy of the two images divided by their joint entropy.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
double NormalizedMutualInformation2 = (Entropy1 + Entropy2) / JointEntropy;
// Software Guide : EndCodeSnippet
std::cout << "Normalized Mutual Information 2 = "
<< NormalizedMutualInformation2 << std::endl;
// Software Guide : BeginLatex
//
// You probably will find very interesting how the value of Mutual
// Information is strongly dependent on the number of bins over which the
// histogram is defined.
//
// Software Guide : EndLatex
return EXIT_SUCCESS;
}
itkImageToHistogramFilter.h
itk::JoinImageFilter
Join two images, resulting in an image where each pixel has the components of the first image followe...
Definition: itkJoinImageFilter.h:206
itkImageFileReader.h
itkImage.h
itkJoinImageFilter.h
itk::ImageFileReader
Data source that reads image data from a single file.
Definition: itkImageFileReader.h:75
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:88
New
static Pointer New()
itk::Statistics::ImageToHistogramFilter
This class generates a histogram from an image.
Definition: itkImageToHistogramFilter.h:50
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44