FMEwork/ITK/Sources/ITK/MLITK/ITKSupport/mlITKSupportToolFunctions.h

Go to the documentation of this file.

// **InsertLicense** code
//----------------------------------------------------------------------------------

//----------------------------------------------------------------------------------
#ifndef __mlITKSupportToolFunctions_H
#define __mlITKSupportToolFunctions_H

// Include dll-specific settings.
#include "mlInitSystemITKSupport.h"

#ifndef __mlModuleIncludes_H
#include "mlModuleIncludes.h"
#endif

#ifndef __mlPointList_H
#include "mlPointList.h"
#endif
#ifndef __mlVectorList_H
#include "mlVectorList.h"
#endif
#ifndef __mlXMarkerList_H
#include "mlXMarkerList.h"
#endif


#include <itkImage.h>
#include <itkImportImageFilter.h>
#include <itkMatrix.h>
#include <itkPointSet.h>

ML_START_NAMESPACE

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
extern MLITK_SUPPORT_EXPORT void postITKException(const itk::ExceptionObject &e,
                                                  const Module               *module,
                                                  MLMessageType              messageType,
                                                  const std::string          &handling="");



//---------------------------------------------------------------------------
//------------------------------------------------------------------------------------
extern MLITK_SUPPORT_EXPORT void copyITKDataBufferToMLSubImg(void              *itkData,
                                                             size_t            itkDataTypeSize,
                                                             const SubImageBox &itkBox,
                                                             SubImage          &outSubImg);


//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename SIZE_TYPE>
SIZE_TYPE ITKSizeFromMLVector(const ImageVector &vec)
{
  ML_TRACE_IN("template ITKSizeFromMLVector()");

  SIZE_TYPE sizeObj;
  if (sizeObj.GetSizeDimension() > static_cast<unsigned int>(vec.dim())){
    ML_PRINT_FATAL_ERROR("ITKSizeFromMLVector", ML_BAD_DIMENSION, "Too high ITK dimension. The ML can handle only 6 dimensions.");
  }
  else{
    for (unsigned int i=0; i<sizeObj.GetSizeDimension(); i++){ sizeObj[i] = vec[i]; }
  }
  return sizeObj;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename INDEX_TYPE>
INDEX_TYPE ITKIndexFromMLVector(const ImageVector &vec)
{
  ML_TRACE_IN("template ITKIndexFromMLVector()");

  INDEX_TYPE sizeObj;
  if (sizeObj.GetIndexDimension() > static_cast<unsigned int>(vec.dim())){
    ML_PRINT_FATAL_ERROR("ITKIndexFromMLVector", ML_BAD_DIMENSION, "Too high ITK dimension. The ML can handle only 6 dimensions.");
  }
  else{
    for (unsigned int i=0; i<sizeObj.GetIndexDimension(); i++){ sizeObj[i] = vec[i]; }
  }
  return sizeObj;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename INDEX_TYPE>
ImageVector MLVectorFromITKSize(const INDEX_TYPE &sizeObj, MLint defaultVal)
{
  ML_TRACE_IN("template MLVectorFromITKSize()");

  ImageVector vec(0);
  unsigned int sizeDim = sizeObj.GetSizeDimension();
  unsigned int vecDim  = static_cast<unsigned int>(vec.dim());
  if (sizeDim > vecDim){
    sizeDim = vecDim;
    ML_PRINT_WARNING("MLVectorFromITKSize", ML_BAD_DIMENSION, "Too high ITK dimension. Highest dimensions ignored.");
  }

  for (unsigned int i=0; i < vecDim; i++){ vec[i] = (i < sizeDim) ? sizeObj[i] : defaultVal; }
  return vec;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename INDEX_TYPE>
ImageVector MLVectorFromITKIndex(const INDEX_TYPE &indexObj, MLint defaultVal)
{
  ML_TRACE_IN("template MLVectorFromITKIndex()");

  ImageVector vec(0);
  unsigned int idxDim = indexObj.GetIndexDimension();
  unsigned int vecDim = static_cast<unsigned int>(vec.dim());
  if (idxDim > vecDim){
    idxDim = vecDim;
    ML_PRINT_WARNING("MLVectorFromITKIndex", ML_BAD_DIMENSION, "Too high ITK dimension. Highest dimensions ignored.");
  }

  for (unsigned int i=0; i < vecDim; i++){ vec[i] = (i < idxDim) ? indexObj[i] : defaultVal; }
  return vec;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename VECTOR_TYPE>
VECTOR_TYPE ITKVectorFromMLVec3(const Vector3 &vec, MLdouble defaultComp)
{
  ML_TRACE_IN("template ITKVectorFromMLVec3()");

  VECTOR_TYPE sizeObj;
  for (unsigned int i=0; i<sizeObj.GetVectorDimension(); i++){
    sizeObj[i] = i < 3 ? vec[i] : defaultComp;
  }
  return sizeObj;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename VECTOR_TYPE>
VECTOR_TYPE ITKVectorFromMLVec4(const Vector4 &vec, MLdouble defaultComp)
{
  ML_TRACE_IN("template ITKVectorFromMLVec4()");

  VECTOR_TYPE sizeObj;
  for (unsigned int i=0; i<sizeObj.GetVectorDimension(); i++){
    sizeObj[i] = i < 4 ? vec[i] : defaultComp;
  }
  return sizeObj;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename POINT_TYPE>
POINT_TYPE ITKPointFromMLVec3(const Vector3 &vec, MLdouble defaultComp)
{
  ML_TRACE_IN("template ITKPointFromMLVec3()");

  POINT_TYPE sizeObj;
  for (unsigned int i=0; i<sizeObj.GetPointDimension(); i++){
    sizeObj[i] = i < 3 ? vec[i] : defaultComp;
  }
  return sizeObj;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename POINT_TYPE>
POINT_TYPE ITKPointFromMLVec4(const Vector4 &vec, MLdouble defaultComp)
{
  ML_TRACE_IN("template ITKPointFromMLVec4()");

  POINT_TYPE sizeObj;
  for (unsigned int i=0; i<sizeObj.GetPointDimension(); i++){
    sizeObj[i] = i < 4 ? vec[i] : defaultComp;
  }
  return sizeObj;
}


//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename REGION_PARENT_TYPE>
typename REGION_PARENT_TYPE::RegionType ITKRegionFromMLSubImgBox(const SubImageBox &subImgBox)
{
  ML_TRACE_IN("template ITKRegionFromMLSubImgBox()");

  typename REGION_PARENT_TYPE::RegionType region;

  // Set ITK region origin to the ML SubImageBox origin.
  region.SetIndex(ITKIndexFromMLVector<ITKML_TYPENAME REGION_PARENT_TYPE::IndexType>(subImgBox.v1));

  // Set ITK region extent to the ML SubImageBox extent.
  region.SetSize(ITKSizeFromMLVector<ITKML_TYPENAME REGION_PARENT_TYPE::SizeType>(subImgBox.getExtent()));
  return region;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename REGION_PARENT_TYPE>
SubImageBox MLSubImgBoxFromITKRegion(const typename REGION_PARENT_TYPE::RegionType &region)
{
  ML_TRACE_IN("template MLSubImgBoxFromITKRegion()");

  // Get ITK region origin into the ML SubImageBox origin.
  // For dimensions not specified by ITK use coordinate 0.
  SubImageBox box;
  box.v1 = MLVectorFromITKIndex(region.GetIndex(), 0);

  // Get ITK region extent into the ML SubImageBox extent.
  // For dimensions not specified by ITK use extent 1.
  ImageVector ext = MLVectorFromITKSize(region.GetSize(), 1);

  // Return the constructed SubImageBox.
  box.v2 = box.v1 + (ext - ImageVector(1));
  return box;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename ARRAY_TYPE, typename STL_VECTOR>
ARRAY_TYPE ITKArrayFromSTLVector(const STL_VECTOR &stlVec)
{
  ML_TRACE_IN("template ITKArrayFromSTLVector()");

  const size_t mFieldSize = stlVec.size();
  ARRAY_TYPE retArray(mFieldSize);
  for (size_t i=0; i<mFieldSize; ++i){ retArray[i] = stlVec[i]; }
  return retArray;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename STL_CONTAINER, typename ITK_ARRAY_TYPE>
const STL_CONTAINER STLVectorFromITKArray(const ITK_ARRAY_TYPE &itkArray)
{
  ML_TRACE_IN("template STLVectorFromITKArray()");

  const size_t itkFieldSize = itkArray.size();
  STL_CONTAINER stlVector;
  for (size_t i=0; i<itkFieldSize; ++i){ stlVector.push_back(itkArray[i]); }
  return stlVector;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename INDEX_TYPE>
INDEX_TYPE ITKIndexFromMLVec6(const Vector6 &vec)
{
  ML_TRACE_IN("template ITKIndexFromMLVector()");

  INDEX_TYPE sizeObj;
  if (sizeObj.GetIndexDimension() > 6){
    ML_PRINT_FATAL_ERROR("ITKIndexFromMLVec6", ML_BAD_DIMENSION, "Too high ITK dimension. The ML can handle only 6 dimensions.");
  }
  else{
    for (unsigned int i=0; i<sizeObj.GetIndexDimension(); i++){ sizeObj[i] = vec[i]; }
  }
  return sizeObj;
}


//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<class FilterType>
typename FilterType::NodeContainer::Pointer ITKNodeContainerFromBasePointer(Base *baseVal)
{
  ML_TRACE_IN("ITKNodeContainerFromBasePointer()");

  XMarkerListContainer *xmlc = NULL;
  XMarkerList          *xml  = NULL;
  PointList            *pl   = NULL;
  VectorList           *vl   = NULL;

  // Get base field.
  if (baseVal){
    // Check for different base object types.
    if (ML_BASE_IS_A(baseVal, XMarkerListContainer)){
      xmlc = static_cast<XMarkerListContainer*>(baseVal);
    }
    else if (ML_BASE_IS_A(baseVal, XMarkerList)){
      xml = static_cast<XMarkerList*>(baseVal);
    }
    else if (ML_BASE_IS_A(baseVal, PointList)){
      pl = static_cast<PointList*>(baseVal);
    }
    else if (ML_BASE_IS_A(baseVal, VectorList)){
      vl = static_cast<VectorList*>(baseVal);
    }
    else{
      // No valid type in base field. That can happen.
    }
  }

  // Get number of points from list.
  MLssize_t numVals = (xmlc ? static_cast<MLssize_t>(xmlc->getList()->size()) :
                      (xml  ? static_cast<MLssize_t>(xml->size())             :
                      (pl   ? pl->getNum() :
                      (vl   ? vl->getNum() : 0))));


  if (numVals > 0){
    // Create a NodeContainer and initialize it.
    typename FilterType::NodeContainer::Pointer nc = FilterType::NodeContainer::New();
    nc->Initialize();

    // Insert all list elements into the NodeContainer.
    for (MLssize_t c=0; c < numVals; ++c){
      // Get point list value if we have such a list.
      Vector6 plVal(0);
      if (pl){
        float px=0, py=0, pz=0;
        pl->getValue(c, px, py, pz);
        plVal[0] = px;
        plVal[1] = py;
        plVal[2] = pz;
      }

      // Get point list value if we have such a list.
      Vector6 vlVal(0);
      int vecType=0;
      if (vl){
        float px=0, py=0, pz=0;
        vl->getPoint(c, vecType, px, py, pz);
        vlVal[0]=px;
        vlVal[1]=py;
        vlVal[2]=pz;
      }

      // get position from any of the lists.
      Vector6 mlVec((xmlc ? (static_cast<XMarker*>(xmlc->getList()->getItemAt(c)))->pos :
                    (xml  ? (static_cast<XMarker*>(xml->getItemAt(c)))->pos             :
                    (pl   ? plVal                                            :
                    (vl   ? vlVal                                            : Vector6(0))))));

      typename FilterType::NodeType pnt;
      pnt.SetValue(0);
      pnt.SetIndex(ITKIndexFromMLVec6<ITKML_TYPENAME FilterType::NodeType::IndexType>(mlVec));
      nc->InsertElement(c, pnt);
    }

    return nc;
  }
  else{
    return NULL;
  }
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename DTYPE, unsigned int ROW_DIM, unsigned int COL_DIM>
typename itk::Matrix<DTYPE, ROW_DIM, COL_DIM> ITKMatrixFromMLMatrix(const Matrix4 &mat)
{
  ML_TRACE_IN("ITKMatrixFromMLMatrix()");

  // Create default matrix.
  itk::Matrix<DTYPE, ROW_DIM, COL_DIM> retMat;
  retMat.SetIdentity();

  // Do not allow more than four dimensions.
  int maxRowDim = ROW_DIM;
  if (maxRowDim > 4){
    maxRowDim = 4;
    ML_PRINT_WARNING("ITKMatrixFromMLMatrix", ML_BAD_DIMENSION, "Too high row dimension of ITK matrix. Only 4 dimensions will be converted.");
  }
  int maxColDim = COL_DIM;
  if (maxColDim > 4){
    maxColDim = 4;
    ML_PRINT_WARNING("ITKMatrixFromMLMatrix", ML_BAD_DIMENSION, "Too high col dimension of ITK matrix. Only 4 dimensions will be converted.");
  }

  // Copy elements.
  for (int row=0; row < maxRowDim; ++row){
    for (int col=0; col < maxColDim; ++col){
      retMat[row][col] = mat[row][col];
    }
  }

  return retMat;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<typename DTYPE, unsigned int ROW_DIM, unsigned int COL_DIM>
Matrix4 MLMatrixFromITKMatrix(const ITKML_TYPENAME itk::Matrix<DTYPE, ROW_DIM, COL_DIM> &mat, bool fillWithID=false)
{
  ML_TRACE_IN("MLMatrixFromITKMatrix()");

  // Do not allow more than four dimensions.
  int maxRowDim = ROW_DIM;
  if (maxRowDim > 4){
    maxRowDim = 4;
    ML_PRINT_WARNING("MLMatrixFromITKMatrix", ML_BAD_DIMENSION, "Too high row dimension of ITK matrix. Only 4 dimensions will be converted.");
  }
  int maxColDim = COL_DIM;
  if (maxColDim > 4){
    maxColDim = 4;
    ML_PRINT_WARNING("MLMatrixFromITKMatrix", ML_BAD_DIMENSION, "Too high col dimension of ITK matrix. Only 4 dimensions will be converted.");
  }

  // Copy elements.
  Matrix4 retMat;
  if (fillWithID){ retMat = Matrix4::getIdentity(); } else { retMat.set(0); }
  for (int row=0; row < maxRowDim; ++row){
    for (int col=0; col < maxColDim; ++col){
      retMat[row][col] = mat[row][col];
    }
  }

  return retMat;
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template<class POINTSETTYPE>
typename POINTSETTYPE::Pointer ITKPointSetFromBasePointer(Base *baseVal)
{
  ML_TRACE_IN("ITKPointSetFromBasePointer()");

  XMarkerListContainer *xmlc = NULL;
  XMarkerList          *xml  = NULL;
  PointList            *pl   = NULL;
  VectorList           *vl   = NULL;

  // Get base field.
  if (baseVal){
    // Check for different base object types.
    if (ML_BASE_IS_A(baseVal, XMarkerListContainer)){ xmlc = static_cast<XMarkerListContainer*>(baseVal); }
    else if (ML_BASE_IS_A(baseVal, XMarkerList))    { xml  = static_cast<XMarkerList*>(baseVal);          }
    else if (ML_BASE_IS_A(baseVal, PointList))      { pl   = static_cast<PointList*>(baseVal);            }
    else if (ML_BASE_IS_A(baseVal, VectorList))     { vl   = static_cast<VectorList*>(baseVal);           }
    else{ /* No valid type in base field. That can happen. */ }
  }

  // Get number of points from list.
  MLssize_t numVals = (xmlc ? static_cast<MLssize_t>(xmlc->getList()->size()) :
                      (xml  ? static_cast<MLssize_t>(xml->size())             :
                      (pl   ? pl->getNum()                                    :
                      (vl   ? vl->getNum()                                    : 0))));


  // Determine maximum point dimension, clamp it to 6.
  int maxDim = POINTSETTYPE::PointDimension;
  if (maxDim > 6){
    maxDim = 6;
    ML_PRINT_WARNING("ITKPointSetFromBasePointer", 
                     ML_BAD_DIMENSION, 
                     "Too high dimension of ITK PointSet. Only 6 dimensions will be converted.");
  }

  if (numVals > 0){

    // Create point set object.
    typename POINTSETTYPE::Pointer outputPointSet = POINTSETTYPE::New();

    typedef typename POINTSETTYPE::PointDataContainer DataContainer;
    outputPointSet->SetPointData( DataContainer::New());
    outputPointSet->GetPoints()->Reserve( numVals );
    outputPointSet->GetPointData()->Reserve( numVals );

    typename POINTSETTYPE::PointIdentifier  pointId = 0;
    typename POINTSETTYPE::PointType  point;

    // Insert all list elements into the NodeContainer.
    for (MLssize_t c=0; c < numVals; ++c){
      // Get point list value if we have such a list.
      Vector6 plVal(0);
      if (pl){
        float px=0, py=0, pz=0;
        pl->getValue(c, px, py, pz);
        plVal[0] = px;
        plVal[1] = py;
        plVal[2] = pz;
      }

      // Get point list value if we have such a list.
      Vector6 vlVal(0);
      int vecType=0;
      if (vl){
        float px=0, py=0, pz=0;
        vl->getPoint(c, vecType, px, py, pz);
        vlVal[0]=px;
        vlVal[1]=py;
        vlVal[2]=pz;
      }

      // get position from any of the lists.
      Vector6 mlVec((xmlc ? (static_cast<XMarker*>(xmlc->getList()->getItemAt(c)))->pos :
                    (xml  ? (static_cast<XMarker*>(xml->getItemAt(c)))->pos             :
                    (pl   ? plVal                                                       :
                    (vl   ? vlVal                                                       : Vector6(0))))));

      // Copy point components.
      for (int j=0; j < maxDim ;j++){ point[j] = mlVec[j]; }
      outputPointSet->SetPoint(pointId++, point );
    }

    return outputPointSet;
  }
  else{
    // Create and return default (empty) point set object.
    typename POINTSETTYPE::Pointer outputPointSet = POINTSETTYPE::New();
    return outputPointSet;
  }
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template <typename ITK_IMPORT_IMAGE_FILTER_TYPE>
void ITKSetOriginFromVec3(itk::ImportImageFilter<typename ITK_IMPORT_IMAGE_FILTER_TYPE::OutputImagePixelType, 
                                                 ITK_IMPORT_IMAGE_FILTER_TYPE::OutputImageType::ImageDimension>* importImageFilter, 
                          const Vector3 &orig)
{
  importImageFilter->SetOrigin(ITKPointFromMLVec3<ITKML_TYPENAME ITK_IMPORT_IMAGE_FILTER_TYPE::OriginType>(orig,0));
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template <typename ITK_IMAGE_TYPE>
void ITKSetOriginFromVec3(itk::Image<typename ITK_IMAGE_TYPE::PixelType, 
                                     ITK_IMAGE_TYPE::ImageDimension>* image, 
                          const Vector3 &orig)
{
  image->SetOrigin(ITKPointFromMLVec3<ITKML_TYPENAME ITK_IMAGE_TYPE::PointType>(orig,0));
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template <typename ITK_CLASS_TYPE>
void setITKWorldFromMedicalImageProperty(MedicalImageProperties props,
                                         ITK_CLASS_TYPE *image,
                                         bool correctSVS=true)
{
  ML_TRACE_IN("template setITKWorldFromMedicalImageProperty()");

  if (!image){ return; }

  // Get ML voxel size/scaling.
  const Vector3 scales(props.getVoxelSize());

  // Set voxel spacing, origin and direction cosines of output image and use 1 for unknown higher dimensional components.
  image->SetSpacing(ITKVectorFromMLVec3<ITKML_TYPENAME ITK_CLASS_TYPE::SpacingType>(scales, 1));

  // Correct the sub voxel position by translating it to the center.
  // In the ML the voxel coordinate is in the corner, in itk and vtk
  // it is in the center.
  if (correctSVS){ props.translateVoxelToWorldMatrix(Vector3( 0.5, 0.5, 0.5)); }

  Matrix4 mat = props.getVoxelToWorldMatrix();

  // Determine image origin for itk image:
  // Extract (sub voxel) corrected translation vector.
  const Vector3 orig(mat[0][3], mat[1][3], mat[2][3]);
  ITKSetOriginFromVec3<ITK_CLASS_TYPE>(image, orig);
  
  // Image direction/orientation:
  // Multiply mat with inverse voxel scale factor to get normalized direction cosines/orientation matrix.
  const Vector3 invScales(MLValueIs0WOM(scales[0]) ? 1 : 1/scales[0],
                          MLValueIs0WOM(scales[1]) ? 1 : 1/scales[1],
                          MLValueIs0WOM(scales[2]) ? 1 : 1/scales[2]);

  mat[0][0] = mat[0][0]*invScales[0];
  mat[1][0] = mat[1][0]*invScales[0];
  mat[2][0] = mat[2][0]*invScales[0];
  mat[3][0] = 0;

  mat[0][1] = mat[0][1]*invScales[1];
  mat[1][1] = mat[1][1]*invScales[1];
  mat[2][1] = mat[2][1]*invScales[1];
  mat[3][1] = 0;

  mat[0][2] = mat[0][2]*invScales[2];
  mat[1][2] = mat[1][2]*invScales[2];
  mat[2][2] = mat[2][2]*invScales[2];
  mat[3][2] = 0;

  mat[0][3] = 0;
  mat[1][3] = 0;
  mat[2][3] = 0;
  mat[3][3] = 1;

  typename ITK_CLASS_TYPE::DirectionType dirCosines = ITKMatrixFromMLMatrix<MLdouble,
                                                                            ITK_CLASS_TYPE::DirectionType::RowDimensions,
                                                                            ITK_CLASS_TYPE::DirectionType::ColumnDimensions>(mat);
  image->SetDirection(dirCosines);
}


//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template <typename ITK_CLASS_TYPE>
void setMLWorldFromITKScaleOriginAndOrientation(const ITK_CLASS_TYPE *image,
                                                MedicalImageProperties &props,
                                                bool correctSVS=true)
{
  ML_TRACE_IN("template setMLWorldFromITKScaleOriginAndOrientation()");

  if (!image){ return; }
  if (ITK_CLASS_TYPE::ImageDimension < 2){
    ML_PRINT_ERROR("setMLWorldFromITKScaleOriginAndOrientation", 
                   ML_BAD_PARAMETER,
                   "1 D itk world matrix conversion is not supported; "
                   "using identity as return value.");
    props.setVoxelToWorldMatrix(Matrix4::getIdentity());
  }

  // Get location, voxel scaling and orientation of itk image or importer.
  const typename ITK_CLASS_TYPE::SpacingType   spacing   = image->GetSpacing();
  const typename ITK_CLASS_TYPE::PointType     origin    = image->GetOrigin();
  const typename ITK_CLASS_TYPE::DirectionType direction = image->GetDirection();

  // Copy 3x3 itk matrix components to corresponding ones in id ML matrix.
  const Matrix4 dirCosines = MLMatrixFromITKMatrix<double,
                                                   ITK_CLASS_TYPE::ImageDimension,
                                                   ITK_CLASS_TYPE::ImageDimension>(direction, true);

  // Create the ML world matrix as id and set up scaling, orientation and origin provided by the itk image or importer.
  Matrix4 mat = Matrix4::getIdentity();

  // Compose scale, directionCosines and translation to the new ML world matrix.
  mat[0][0] = dirCosines[0][0]*spacing[0];
  mat[1][0] = dirCosines[1][0]*spacing[0];
  mat[2][0] = dirCosines[2][0]*spacing[0];
  mat[3][0] = 0;

  mat[0][1] = dirCosines[0][1]*spacing[1];
  mat[1][1] = dirCosines[1][1]*spacing[1];
  mat[2][1] = dirCosines[2][1]*spacing[1];
  mat[3][1] = 0;

  // Avoid invalid accesses to spacing if we operate on 2D images.
  const double spacing2 = (ITK_CLASS_TYPE::ImageDimension > 2) ? spacing[2] : 1;
  mat[0][2] = dirCosines[0][2]*spacing2;
  mat[1][2] = dirCosines[1][2]*spacing2;
  mat[2][2] = dirCosines[2][2]*spacing2;
  mat[3][2] = 0;

  mat[0][3] = origin[0];
  mat[1][3] = origin[1];
  mat[2][3] = (ITK_CLASS_TYPE::ImageDimension > 2) ? origin[2] : 0;
  mat[3][3] = 1;

  props.setVoxelToWorldMatrix(mat);

  // Subtract half voxel shift - in MedicalImageProperties the matrix considers
  // voxel positions at the corner and not in the center.
  if (correctSVS){ props.translateVoxelToWorldMatrix(Vector3(-0.5, -0.5, -0.5)); }
}


//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template <typename ITK_INDATATYPE, unsigned int DIM>
typename  itk::ImportImageFilter<ITK_INDATATYPE, DIM>::Pointer getITKImportImageFromSubImg(const SubImage               &inSubImg,
                                                                                           const MedicalImageProperties &props,
                                                                                           bool                         correctSVS=true)
{
  ML_TRACE_IN("template getITKImageFromSubImg()");

  // Create a new ImportImageFilter.
  typedef itk::ImportImageFilter<ITK_INDATATYPE, DIM> ImportFilterType;
  typedef typename ImportFilterType::Pointer ImportFilterPointerType;

  ImportFilterPointerType importer = ImportFilterType::New();

  // Get box from the input subimage, convert it to a ITK region and
  // set it as output image region of the image import filter.
  SubImageBox inSubImgBox(inSubImg.getBox());
  importer->SetRegion(ITKRegionFromMLSubImgBox<ImportFilterType>(inSubImgBox));

  // Set voxel spacing, origin and orientation of itk image.
  setITKWorldFromMedicalImageProperty<ImportFilterType>(props, &(*importer), correctSVS);

  // Set pointer to the data to be imported by the image import filter.
  importer->SetImportPointer(static_cast<ITK_INDATATYPE*>(inSubImg.getData()),
                             inSubImgBox.getNumVoxels(),
                             false); // we don't want to let ITK managing memory
  importer->Update();

  // Return the importer object.
  return importer;
}


  //----------------------------------------------------------------------------------
  //----------------------------------------------------------------------------------
  template <typename RETURN_TYPE_PTR, typename FILTER_TYPE, typename VOXEL_TYPE>
  RETURN_TYPE_PTR getInputAsItkImportImageAndSubImg(Module                &op,
                                                    int                   inIdx,
                                                    TSubImage<VOXEL_TYPE> &dataSubImg,
                                                    bool                  correctSVS=true)
  {
    // Get and load connected image only if it's valid and not a a redirected one
    // (which may occur in cases of optional mask image inputs).
    PagedImage* pInImg = op.getUpdatedInputImage(inIdx);
    if (pInImg){
      // Get 3D box input image.
      const ImageVector _inImgExt = pInImg->getImageExtent();
      const SubImageBox inImgBox(ImageVector(_inImgExt.x, _inImgExt.y, _inImgExt.z,1,1,1));

      // Check whether TSubImage contains enough memory, if not we need to free and reallocate it.
      if (dataSubImg.getData() && (dataSubImg.getBox().getExtent() != inImgBox.getExtent())){
        dataSubImg.free();
      }
      // Set up data box and allocate data if still not there.
      dataSubImg.setBox(inImgBox);
      if (!dataSubImg.getData()){
        dataSubImg.allocateAsMemoryBlockHandle(ML_FATAL_MEMORY_ERROR);
      }
      if (!dataSubImg.getData()){
        return NULL;
      }

      // Get image data from input.
      MLErrorCode err = pInImg->getTile(dataSubImg);
      if (ML_RESULT_OK == err){
        // Create a new ImportImageFilter.
        typename FILTER_TYPE::Pointer importer = FILTER_TYPE::New();

        // Get box from the input subimage, convert it to a ITK region and
        // set it as output image region of the image import filter.
        importer->SetRegion(ITKRegionFromMLSubImgBox<FILTER_TYPE>(inImgBox));

        // Pass world transformation of ML image to world transformation of itk image or importer.
        setITKWorldFromMedicalImageProperty<FILTER_TYPE>(*pInImg, &(*importer), correctSVS);

        // Set pointer to the data to be imported by the image import filter.
        // Pass true to make ITK manage the memory, i.e. if the last import object
        // disappears the data also disappears.
        importer->SetImportPointer(static_cast<VOXEL_TYPE*>(dataSubImg.getData()), inImgBox.getNumVoxels(), false);
        importer->Update();
        return importer->GetOutput();
      }
      else{
        ML_PRINT_ERROR("mlITKSupportToolFunctions::getInputAsItkImportImageAndSubImg()",
                       err,
                       "Failed to request input image data, probably subsequent operations will fail");
      }
    }

    // Error, no image could be loaded.
    return NULL;
  }


ML_END_NAMESPACE

#endif