[Insight-users] itkMutualInformationImageToImageMetric

[Maarten Beek]

Hi all,

Is it possible to use the itkMutualInformationImageToImageMetric class by it self, to calculate mutual information values between two data sets?

I currently implemented something like below, but get positive as well as negative values.

Should I set up a full registration pipeline instead (like Examples/Registration/ImageRegistration2.cxx) and get the MI value after one iteration (optimizer->SetNumberOfIterations(1))?

Thanks Maarten



#include <vector>
#include <fstream>

#include "itkImage.h"
#include "itkMutualInformationImageToImageMetric.h"
#include "itkTranslationTransform.h"
#include "itkLinearInterpolateImageFunction.h"

const unsigned int Dimension = 1;
typedef double PixelType;
typedef itk::Image<PixelType, Dimension> ImageType;

ImageType::Pointer readData(const std::string& filename);

int main(int argc, const char* argv[])
{
  typedef itk::MutualInformationImageToImageMetric<ImageType, ImageType>    MetricType;
  typedef itk::TranslationTransform<double, Dimension>                        TransformType;
  typedef itk::LinearInterpolateImageFunction<ImageType, double>            InterpolatorType;

  TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();
  InterpolatorType::Pointer interpolator = InterpolatorType::New();

    std::string file1 = "C:\\Users\\Maarten\\Documents\\Visual Studio 2010\\Projects\\itkMI\\files\\X2.txt";
    std::string file2 = "C:\\Users\\Maarten\\Documents\\Visual Studio 2010\\Projects\\itkMI\\files\\Y2.txt";
    ImageType::Pointer data1 = readData(file1);
    ImageType::Pointer data2 = readData(file2);

    interpolator->SetInputImage(data2);

    MetricType::Pointer metric = MetricType::New();
    metric->SetFixedImage(data1);
    metric->SetFixedImageRegion(data1->GetLargestPossibleRegion());
    metric->SetMovingImage(data2);
    metric->SetTransform(transform);
    metric->SetInterpolator(interpolator);
    //metric->ReinitializeSeed(12345); // with constant seed value -> deterministic
  //metric->SetFixedImageStandardDeviation(  0.4 );
  //metric->SetMovingImageStandardDeviation( 0.4 );
    const unsigned int numberOfPixels = data1->GetBufferedRegion().GetNumberOfPixels();
    //metric->SetNumberOfSpatialSamples( numberOfPixels*0.1 );

    MetricType::MeasureType measure = metric->GetValue(transform->GetParameters());

    std::cout << "MI value: " << measure << std::endl;

    return 1;
}

ImageType::Pointer readData(const std::string& filename)
{
    ImageType::Pointer data = ImageType::New();
    data->Initialize();

    typedef itk::ImageRegionIteratorWithIndex<ImageType> IteratorType;

    std::vector<double> values;
    std::ifstream file(filename);
    if( file.is_open() )
    {
        double value;
        while( !file.eof() )
        {
            file >> value;
            values.push_back(value);
        }
    }
    unsigned int numValues = values.size();
    if( 0 < numValues )
    {
        // maybe use itk::ImportImageFilter?
        ImageType::IndexType start;
        start[0] = 0;
        ImageType::SizeType size;
        size[0] = numValues;
        ImageType::RegionType region;
        region.SetSize(size);
        region.SetIndex(start);
        data->SetRegions(region);
        data->Allocate();

        ImageType::PointType origin;
        origin[0] = 0.0;
        data->SetOrigin(origin);
        ImageType::SpacingType spacing;
        spacing[0] = 1.0;
        data->SetSpacing(spacing);

        IteratorType it(data, data->GetLargestPossibleRegion());
        int i = 0;
        for ( it.GoToBegin(); !it.IsAtEnd(); ++it, ++i)
        {
            it.Set(values[i]);
        }
    }
    return data;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20130626/abfe7c65/attachment.htm>