使用互信息进行图像配准

熵

熵是描述一个事件的不确定性的，在英文中的描述是disorder，也就是无序，杂乱。写一个事件的熵可以表示为：

$$E(A)=-\sum _A{p}_i{\log }_2{p}_i$$

当$p_i=0$的时候，$p_i\log p_i=0$。也就是事件所有可能情况的概率与其对应对数的乘积就和的相反数。假设一个事件$A(X)$存在 { x = 1 , x = 2 , x = 3 , x = 4 } \{x=1,x=2,x=3,x=4\} {x=1,x=2,x=3,x=4}这几种情况，并且经过多次实验，发现这几种情况出现的频数分别为：

x	1	2	3	4
频数	8	3	7	2

那么可以计算出，$p_1=\frac{8}{20}=0.4,p_2=\frac{3}{20}=0.15,p_3=\frac{7}{20}=0.35,p_4=\frac{2}{20}=0.1$。则事件A的熵为：

$$E(A)=-(p_1{\log }_2{p}_1+p_2{\log }_2{p}_2+p_3{\log }_2{p}_3+p_4{\log }_2{p}_4)=1.8016$$

在熵的理论中，认为在自然条件下，所有的事情都会慢慢变得杂乱无序，而熵是描述杂乱无序的程度，所以熵都是会变大的。而熵的最大和最小值分别是多少？

还是事件$A$，重复做了20次实验，发现所有的结果都是$x=1$的情况，这个时候事件是非常确定的，是最不无序的状态，这个时候熵是最小的。

x	1	2	3	4
频数	20	0	0	0

$$E_{min}(A)=-p_1{\log }_2{p}_1=0$$

而当做了20次实验，发现所有情况出现的概率是均等的，也就是各出现了5次的话，这个时候是最杂乱无序的，熵是最大的。

x	1	2	3	4
频数	5	5	5	5

$$E_{max}(A)=-4\times \frac{1}{4}{\log }_2(\frac{1}{4})=2$$

更一般地来讲，当一个事件有n个可能出现的情况，那么最大和最小的熵分别是：
$$\{\begin{array}{ll}{E}_{max}(A)=-{\log }_2(\frac{1}{n})={\log }_2(n)& (当所有情况等概率出现)\\ E_{min}(A)=0& (当只有一种情况出现)\end{array}$$

互信息

互信息(Mutual Information)是熵理论中的概念，描述了当一个事件的引入，对原来事件不确定性的减少量。计算互信息需要使用到，熵，联合熵，条件熵等概率，方便理解，可以使用到韦恩图的方式来理解。

• 熵：事件的不确定度，韦恩图中蓝色和红色圆圈分别表示事件$X,Y$的熵。
• 联合熵：两个事件共同的不确定度，当两个事件完全独立的时候，联合熵等于两个事件熵的和，不完全独立的时候，等于两个事件韦恩图的并集。并且有$H(X,Y)\le H(X)+H(Y)$，$H(X,Y)=-{\sum }_x{\sum }_{y}p(x,y){\log }_{2}p(x,y)$。
• 条件熵：事件X在事件Y引入的情况下，剩余的不确定度，在图中表示为两个圆去除中间共同的部分。
• 互信息：事件X在引入事件Y后不确定度的减少量，即为中间公共的部分,计算公式为：$I(X,Y)={\sum }_x{\sum }_{y}p(x,y)\log \frac{p(x,y)}{p(x)p(y)}$，但是一般地我们更喜欢使用下面的公式来进行互信息的计算。

$$I(X,Y)=H(X)+H(Y)-H(X,Y)\text{\hspace{1em}(互信息计算公式)}$$

显然，当互信息越大，表示两个事件的"重合"越高，"相关性"越大。

使用互信息的图像配准

这部分代码来自 ITK Examples/RegistrationITKv4/ImageRegistration17.cxx。数据来自 ITK Testing Data。由于ITK给的代码读取部分有问题，添加了部分修改。

• 在include部分添加了 #include <itkPNGImageIO.h>的包含头文件。

• 在读取文件之前添加ImageIO。

  ......
  fixedImageReader->SetFileName(argv[1]);
  fixedImageReader->SetImageIO(itk::PNGImageIO::New());
  movingImageReader->SetFileName(argv[2]);
  movingImageReader->SetImageIO(itk::PNGImageIO::New());
  ......
  

• 并且，不需要保存信息，将Print out results以后的结果都注释掉了。

输入参数分别为：

C:/depolytools/InsightToolkit-5.2.0/source/Examples/Data/BrainProtonDensitySliceBorder20.png C:/depolytools/InsightToolkit-5.2.0/source/Examples/Data/BrainProtonDensitySliceShifted13x17y.png 0 0

需要替换以下图片的路径即可。

#include "itkImageRegistrationMethod.h"

#include "itkTranslationTransform.h"
#include "itkMutualInformationHistogramImageToImageMetric.h"
#include "itkAmoebaOptimizer.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include <itkPNGImageIO.h>
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCommand.h"
class CommandIterationUpdate : public itk::Command
{
public:
    using Self = CommandIterationUpdate;
    using Superclass = itk::Command;
    using Pointer = itk::SmartPointer<Self>;
    itkNewMacro(Self);

protected:
    CommandIterationUpdate() { m_IterationNumber = 0; }

public:
    using OptimizerType = itk::AmoebaOptimizer;
    using OptimizerPointer = const OptimizerType*;

    void
        Execute(itk::Object* caller, const itk::EventObject& event) override
    {
        Execute((const itk::Object*)caller, event);
    }

    void
        Execute(const itk::Object* object, const itk::EventObject& event) override
    {
        auto optimizer = static_cast<OptimizerPointer>(object);
        if (!itk::IterationEvent().CheckEvent(&event))
        {
            return;
        }
        std::cout << m_IterationNumber++ << "   ";
        std::cout << optimizer->GetCachedValue() << "   ";
        std::cout << optimizer->GetCachedCurrentPosition() << std::endl;
    }

private:
    unsigned long m_IterationNumber;
};


// C:/depolytools/InsightToolkit-5.2.0/source/Examples/Data/BrainProtonDensitySliceBorder20.png C:/depolytools/InsightToolkit-5.2.0/source/Examples/Data/BrainProtonDensitySliceShifted13x17y.png 0 0
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 ";
        std::cerr << " [initialTx] [initialTy]" << std::endl;
        return EXIT_FAILURE;
    }

    constexpr unsigned int Dimension = 2;
    using PixelType = unsigned char;

    using FixedImageType = itk::Image<PixelType, Dimension>;
    using MovingImageType = itk::Image<PixelType, Dimension>;

    using TransformType = itk::TranslationTransform<double, Dimension>;

    using OptimizerType = itk::AmoebaOptimizer;
    using InterpolatorType =
        itk::LinearInterpolateImageFunction<MovingImageType, double>;
    using RegistrationType =
        itk::ImageRegistrationMethod<FixedImageType, MovingImageType>;


    using MetricType =
        itk::MutualInformationHistogramImageToImageMetric<FixedImageType,
        MovingImageType>;


    auto transform = TransformType::New();
    auto optimizer = OptimizerType::New();
    auto interpolator = InterpolatorType::New();
    auto registration = RegistrationType::New();

    registration->SetOptimizer(optimizer);
    registration->SetTransform(transform);
    registration->SetInterpolator(interpolator);


    auto metric = MetricType::New();
    registration->SetMetric(metric);

    using HistogramSizeType = MetricType::HistogramSizeType;

    HistogramSizeType histogramSize;

    histogramSize.SetSize(2);

    histogramSize[0] = 256;
    histogramSize[1] = 256;

    metric->SetHistogramSize(histogramSize);
    metric->ComputeGradientOff();


    const unsigned int numberOfParameters = transform->GetNumberOfParameters();


    using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
    using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;

    auto fixedImageReader = FixedImageReaderType::New();
    auto movingImageReader = MovingImageReaderType::New();

    fixedImageReader->SetFileName(argv[1]);
    fixedImageReader->SetImageIO(itk::PNGImageIO::New());
    movingImageReader->SetFileName(argv[2]);
    movingImageReader->SetImageIO(itk::PNGImageIO::New());

    registration->SetFixedImage(fixedImageReader->GetOutput());
    registration->SetMovingImage(movingImageReader->GetOutput());

    fixedImageReader->Update();
    movingImageReader->Update();

    FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();

    registration->SetFixedImageRegion(fixedImage->GetBufferedRegion());


    transform->SetIdentity();

    using ParametersType = RegistrationType::ParametersType;

    ParametersType initialParameters = transform->GetParameters();

    initialParameters[0] = 0.0;
    initialParameters[1] = 0.0;

    if (argc > 5)
    {
        initialParameters[0] = std::stod(argv[4]);
        initialParameters[1] = std::stod(argv[5]);
    }

    registration->SetInitialTransformParameters(initialParameters);

    std::cout << "Initial transform parameters = ";
    std::cout << initialParameters << std::endl;

    OptimizerType::ParametersType simplexDelta(numberOfParameters);
    simplexDelta.Fill(5.0);

    optimizer->AutomaticInitialSimplexOff();
    optimizer->SetInitialSimplexDelta(simplexDelta);
    optimizer->MaximizeOn();
    optimizer->SetParametersConvergenceTolerance(0.1); // 1/10th pixel
    optimizer->SetFunctionConvergenceTolerance(0.001); // 0.001 bits
    optimizer->SetMaximumNumberOfIterations(200);

    // Create the Command observer and register it with the optimizer.
    //
    auto observer = CommandIterationUpdate::New();
    optimizer->AddObserver(itk::IterationEvent(), observer);


    try
    {
        registration->Update();
        std::cout << "Optimizer stop condition: "
            << registration->GetOptimizer()->GetStopConditionDescription()
            << std::endl;
    }
    catch (const itk::ExceptionObject& err)
    {
        std::cout << "ExceptionObject caught !" << std::endl;
        std::cout << err << std::endl;
        return EXIT_FAILURE;
    }

    ParametersType finalParameters = registration->GetLastTransformParameters();

    const double finalTranslationX = finalParameters[0];
    const double finalTranslationY = finalParameters[1];

    double bestValue = optimizer->GetValue();
    
    // Print out results
    //
    std::cout << "Result = " << std::endl;
    std::cout << " Translation X = " << finalTranslationX << std::endl;
    std::cout << " Translation Y = " << finalTranslationY << std::endl;
    std::cout << " Metric value  = " << bestValue << std::endl;

    return EXIT_SUCCESS;
}

这里使用到得优化器是Amoeba，不需要梯度信息，所以在metric里可以设置metric->ComputeGradientOff();，而配准结果越好，互信息越大，所以需要在优化器设置optimizer->MaximizeOn();。直接运行代码，可以很快得到比较好的输出。

Result =
 Translation X = 13.011
 Translation Y = 17.0249
 Metric value  = -2.49226

理论上来讲，这里的Metric value的值应该是一个正数，并且不断变大，但是实际上输出却是一个不断变小的负数，这是由于优化器本身是做最小化优化的，为了直接用于最大化优化，在Metric的值上添加一个负号就变成了最小化优化的结果。