Is there any way to convert an Eigen::Matrix back

问题:

I used Eigen library to convert several itk::image images into matrices, and do some dense linear algebra computations on them. Finally, I have the output as a matrix, but I need it in itk::image form. Is there any way to do this?

    const unsigned int numberOfPixels = importSize[0] * importSize[1];
    float* array1 = inverseU.data();
    float* localBuffer = new float[numberOfPixels];
    std::memcpy(localBuffer, array1, numberOfPixels);
    const bool importImageFilterWillOwnTheBuffer = true;
    importFilter->SetImportPointer(localBuffer,numberOfPixels,importImageFilterWillOwnTheBuffer);
    importFilter->Update();

inverseU is the Eigen library matrix (float), importSize is the size of this matrix. When I give importFilter->GetOutput(), and write the result to file, the image I get is like this, which is not correct.

This is the matrix inverseU. https://drive.google.com/file/d/0B3L9EtRhN11QME16SGtfSDJzSWs/view?usp=sharing . It is supposed to give a retinal fundus image in image form, I got the matrix after doing deblurring.

回答1:

Take a look at the ImportImageFilter of itk. In particular, it may be used to build an itk::Image starting from a C-style array (example).

Someone recently asked how to convert a CImg image to ITK image. My answer might be a starting point...

A way to get the array out of a matrix A from Eigen may be found here :

double* array=A.data();  

EDIT : here is a piece of code to turn a matrix of float into a png image saved with ITK. First, the matrix is converted to an itk Image of float. Then, this image is rescaled an cast to a image of unsigned char, using the RescaleIntensityImageFilter as explained here. Finally, the image is saved in png format.

#include <iostream>
#include <itkImage.h>

using namespace itk;
using namespace std;

#include <Eigen/Dense>
using Eigen::MatrixXf;

#include <itkImportImageFilter.h>
#include <itkImageFileWriter.h>
#include "itkRescaleIntensityImageFilter.h"

void eigen_To_ITK (MatrixXf mat)
{

    const unsigned int Dimension = 2;

    typedef itk::Image<unsigned char, Dimension>  UCharImageType;
    typedef itk::Image< float, Dimension > FloatImageType;
    typedef itk::ImportImageFilter< float, Dimension > ImportFilterType;
    ImportFilterType::Pointer importFilter = ImportFilterType::New();

    typedef itk::RescaleIntensityImageFilter< FloatImageType, UCharImageType > RescaleFilterType;
    RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();

    typedef itk::ImageFileWriter<  UCharImageType  > WriterType;
    WriterType::Pointer writer = WriterType::New();


    FloatImageType::SizeType imsize;
    imsize[0] = mat.rows();
    imsize[1] = mat.cols();

    ImportFilterType::IndexType start;
    start.Fill( 0 );
    ImportFilterType::RegionType region;
    region.SetIndex( start );
    region.SetSize( imsize );
    importFilter->SetRegion( region );

    const itk::SpacePrecisionType origin[ Dimension ] = { 0.0, 0.0 };
    importFilter->SetOrigin( origin );

    const itk::SpacePrecisionType spacing[ Dimension ] = { 1.0, 1.0 };
    importFilter->SetSpacing( spacing );

    const unsigned int numberOfPixels = imsize[0] * imsize[1];

    const bool importImageFilterWillOwnTheBuffer = true;

    float * localBuffer = new float[ numberOfPixels ];
    float * it = localBuffer;

    memcpy(it, mat.data(), numberOfPixels*sizeof(float));
    importFilter->SetImportPointer( localBuffer, numberOfPixels,importImageFilterWillOwnTheBuffer );

    rescaleFilter ->SetInput(importFilter->GetOutput());
    rescaleFilter->SetOutputMinimum(0);
    rescaleFilter->SetOutputMaximum(255);


    writer->SetFileName( "output.png" );
    writer->SetInput(rescaleFilter->GetOutput() );
    writer->Update();

}

int main()
{

    const int rows = 42;
    const int cols = 90;
    MatrixXf mat1(rows, cols);
    mat1.topLeftCorner(rows/2, cols/2) = MatrixXf::Zero(rows/2, cols/2);
    mat1.topRightCorner(rows/2, cols/2) = MatrixXf::Identity(rows/2, cols/2);
    mat1.bottomLeftCorner(rows/2, cols/2) = -MatrixXf::Identity(rows/2, cols/2);
    mat1.bottomRightCorner(rows/2, cols/2) = MatrixXf::Zero(rows/2, cols/2);

    mat1+=0.1*MatrixXf::Random(rows,cols);

    eigen_To_ITK (mat1);

    cout<<"running fine"<<endl;
    return 0;
}

The program is build using CMake. Here is the CMakeLists.txt :

cmake_minimum_required(VERSION 2.8 FATAL_ERROR)
project(ItkTest)

find_package(ITK REQUIRED)
include(${ITK_USE_FILE})

# to include eigen. This path may need to be changed
include_directories(/usr/local/include/eigen3)

add_executable(MyTest main.cpp)
target_link_libraries(MyTest ${ITK_LIBRARIES})