mlTSubImage.h

Go to the documentation of this file.

/*************************************************************************************
**
** Copyright 2007, MeVis Medical Solutions AG
**
** The user may use this file in accordance with the license agreement provided with
** the Software or, alternatively, in accordance with the terms contained in a
** written agreement between the user and MeVis Medical Solutions AG.
**
** For further information use the contact form at https://www.mevislab.de/contact
**
**************************************************************************************/
 
#ifndef ML_TSUB_IMAGE_H
#define ML_TSUB_IMAGE_H

 
// ML-includes
#include "mlInitSystemML.h"
#include "mlSubImage.h"
#include "mlTypeTraits.h"
 
#include "mlTemplateSupport.h"
 
#ifdef _MSC_VER
#pragma warning(push)
// warning C4244: 'Argument': Konvertierung von 'typ1' in 'typ2', m�glicher Datenverlust
#pragma warning(disable: 4244)
#endif
 
ML_START_NAMESPACE
 
 
 

template <typename DATATYPE> class TSubImageCursor;
template <typename DATATYPE> class ConstTSubImageCursor;

#define _ML_TSUBIMG_SUBDOT6(pos, offset, stride) \
  (((pos).x - (offset).x) * (stride).x + \
   ((pos).y - (offset).y) * (stride).y + \
   ((pos).z - (offset).z) * (stride).z + \
   ((pos).c - (offset).c) * (stride).c + \
   ((pos).t - (offset).t) * (stride).t + \
   ((pos).u - (offset).u) * (stride).u)

#define _ML_TSUBIMG_SUBDOT3(x, y, z, offset, stride) \
  (((x) - (offset).x) * (stride).x + \
   ((y) - (offset).y) * (stride).y + \
   ((z) - (offset).z) * (stride).z)

#define _ML_TSUBIMG_SUBDOT2(x, y, offset, stride) \
  (((x) - (offset).x) * (stride).x + \
   ((y) - (offset).y) * (stride).y)
 
 
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
template <typename DATATYPE>
class TSubImage : public SubImage
{
 
public:

  typedef DATATYPE ComponentType;

  typedef TSubImageCursor<DATATYPE>      Cursor;
  typedef ConstTSubImageCursor<DATATYPE> ConstCursor;
 
  //---------------------------------------------------------------------------
  //---------------------------------------------------------------------------
  TSubImage(const SubImageBox& box, MLDataType dataType, void* data):SubImage(box, dataType, data)
  {
    // Check whether template type and real data type match if there is valid data.
    if (data && !TypeTraits<DATATYPE>::matches(dataType) )
    {
      ML_PRINT_FATAL_ERROR("TSubImage::TSubImage()", ML_BAD_DATA_TYPE,
        "Mismatch of data type between argument and template typename.");
    }
  }

  TSubImage()
  {
    // Set correct data type enum from template data type.
    setDataType(TypeTraits<DATATYPE>::dataType);
    if constexpr(TypeTraits<DATATYPE>::dataType == ML_INVALID_DATA_TYPE) {
      ML_PRINT_FATAL_ERROR("TSubImage::TSubImage()", ML_BAD_DATA_TYPE,
        "Data type has no registered TypeTraits.");
    }
  }
 

  TSubImage(const SubImage &subImage)
  {
    // Sets box the data is correlated with.
    setBox(subImage.getBox());
 
    // Sets correct data type enum.
    setDataType(subImage.getDataType());
 
    // Stores data pointer (but \e not the memory block handle).
    setData(subImage.getData());
 
    // Sets the source image extent.
    setSourceImageExtent(subImage.getSourceImageExtent());
 
    // Checks whether template type and real data type match if there is valid data.
    // This check only works if TypeTraits has been specified for the argument type.
    if (subImage.getData() && !TypeTraits<DATATYPE>::matches(subImage.getDataType()) )
    {
      ML_PRINT_FATAL_ERROR("TSubImage::TSubImage()", ML_BAD_DATA_TYPE,
        "Mismatch of data type between argument and template typename.");
    }
  }

  inline TSubImage(const TSubImage<DATATYPE>& typedSubImage) :
    SubImage       (typedSubImage)
  {
  }

  inline ~TSubImage() override = default;

  inline TSubImage &operator=(const TSubImage<DATATYPE>& typedSubImage)
  {
    SubImage::operator=(typedSubImage);
    return *this;
  }

 
  //---------------------------------------------------------------------------
  //---------------------------------------------------------------------------

  inline const DATATYPE* getData()    const { return static_cast<DATATYPE*>(_data); }
  inline       DATATYPE* getData()          { return static_cast<DATATYPE*>(_data); }

 
  //---------------------------------------------------------------------------
  //---------------------------------------------------------------------------
  // 64-bit arithmetic is not needed here, because the pointer offset goes from the valid
  // address space to another position inside the valid address space. We do not
  // fall outside 32-bit ranges.
  inline const DATATYPE* getSubImagePointer(const ImageVector& position) const { return static_cast<DATATYPE*>(_data) + ImageVector::dot(position,_stride); }
  inline       DATATYPE* getSubImagePointer(const ImageVector& position)       { return static_cast<DATATYPE*>(_data) + ImageVector::dot(position,_stride); }

  // 64-bit arithmetic is not needed here, because the pointer offset goes from the valid
  // address space to another position inside the valid address space. We do not
  // fall outside 32-bit ranges.
  inline const DATATYPE* getSubImagePointer(MLint x, MLint y, MLint z) const { return static_cast<DATATYPE*>(_data) + (x*_stride.x + y*_stride.y + _stride.z*z); }
  inline       DATATYPE* getSubImagePointer(MLint x, MLint y, MLint z)       { return static_cast<DATATYPE*>(_data) + (x*_stride.x + y*_stride.y + _stride.z*z); }

  inline const DATATYPE* getImagePointer(const ImageVector& position) const { return static_cast<DATATYPE*>(_data) + _ML_TSUBIMG_SUBDOT6(position, _box.v1, _stride); }
  inline       DATATYPE* getImagePointer(const ImageVector& position)       { return static_cast<DATATYPE*>(_data) + _ML_TSUBIMG_SUBDOT6(position, _box.v1, _stride); }

  inline const DATATYPE* getImagePointer(MLint x, MLint y, MLint z) const { return static_cast<DATATYPE*>(_data) + _ML_TSUBIMG_SUBDOT3(x,y,z, _box.v1, _stride); }
  inline       DATATYPE* getImagePointer(MLint x, MLint y, MLint z)       { return static_cast<DATATYPE*>(_data) + _ML_TSUBIMG_SUBDOT3(x,y,z, _box.v1, _stride); }
 
 
  //-------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------
  inline ImageVector convertPointerToSubImagePosition(DATATYPE* pointer) const
  { return getStride().getVectorPosition(pointer - static_cast<DATATYPE*>(_data)); }

  inline void convertPointerToSubImagePosition(DATATYPE* pointer, MLint* x, MLint* y, MLint* z) const
  {
    const MLint offset = pointer - static_cast<DATATYPE*>(_data);
    const ImageVector stride=getStride();
    const ImageVector imgExt=getImageExtent();
 
    if (x){
      ML_CHECK_THROW(stride.x);
      ML_CHECK_THROW(imgExt.x);
      *x = (offset/stride.x) % imgExt.x;
    }
 
    if (y){
      ML_CHECK_THROW(imgExt.y);
      ML_CHECK_THROW(stride.y);
      *y = (offset/stride.y) % imgExt.y;
    }
 
    if (z){
      ML_CHECK_THROW(imgExt.z);
      ML_CHECK_THROW(stride.z);
      *z = (offset/stride.z) % imgExt.z;
    }
  }
 

  inline ImageVector convertPointerToImagePosition(DATATYPE* pointer) const
  { return convertPointerToSubImagePosition(pointer) + _box.v1; }
 

  inline void convertPointerToImagePosition(DATATYPE* pointer, MLint* x, MLint* y, MLint* z) const
  {
    convertPointerToSubImagePosition(pointer, x, y, z);
    if (x){ *x += _box.v1.x; }
    if (y){ *y += _box.v1.y; }
    if (z){ *z += _box.v1.z; }
  }

 
 
 
  //---------------------------------------------------------------------------
  //---------------------------------------------------------------------------

  inline const DATATYPE& operator[](const ImageVector& position) const    { return getImagePointer(position)[0]; }

  inline DATATYPE& operator[](const ImageVector& position)    { return getImagePointer(position)[0]; }

  inline void setSubImageValue(const ImageVector& position, DATATYPE value) { static_cast<DATATYPE*>(_data)[ ImageVector::dot(position,_stride) ]=value; }

  inline DATATYPE getSubImageValue(const ImageVector& position) const { return static_cast<DATATYPE*>(_data)[ ImageVector::dot(position,_stride) ]; }

  inline void setSubImageValue(MLint x, MLint y, DATATYPE value)      { static_cast<DATATYPE*>(_data)[ x*_stride.x + y*_stride.y ]=value; }
  inline DATATYPE getSubImageValue(MLint x, MLint y)          const { return static_cast<DATATYPE*>(_data)[ x*_stride.x + y*_stride.y ]; }

  inline void setSubImageValue(MLint x, MLint y, MLint z, DATATYPE value){ static_cast<DATATYPE*>(_data)[ x*_stride.x + y*_stride.y+_stride.z*z ]=value; }
  inline DATATYPE getSubImageValue(MLint x, MLint y, MLint z) const { return static_cast<DATATYPE*>(_data)[ x*_stride.x + y*_stride.y+_stride.z*z ]; }

  inline void setImageValue(const ImageVector& position, DATATYPE value) { getImagePointer(position)[0] = value; }
  inline DATATYPE getImageValue(const ImageVector& position) const { return getImagePointer(position)[0]; }

  inline void setImageValue(MLint x, MLint y, DATATYPE value) { static_cast<DATATYPE*>(_data)[_ML_TSUBIMG_SUBDOT2(x,y,_box.v1, _stride)]=value; }
  inline DATATYPE getImageValue(MLint x, MLint y) const { return static_cast<DATATYPE*>(_data)[_ML_TSUBIMG_SUBDOT2(x,y,_box.v1, _stride)]; }

  inline void setImageValue(MLint x, MLint y, MLint z, DATATYPE value){ getImagePointer(x,y,z)[0] = value; }
  inline DATATYPE getImageValue(MLint x, MLint y, MLint z)    const { return getImagePointer(x,y,z)[0]; }

 
 
  //---------------------------------------------------------------------------
  //---------------------------------------------------------------------------
  MLint calculateMinMax(DATATYPE& minValue, DATATYPE& maxValue, const SubImageBox * const validBox=nullptr) const
  {
    // If passed box pointer is valid, then compute intersection with
    // this subimage box, otherwise use non intersected subimage box of this.
    const SubImageBox& thisBox = getBox();
    const SubImageBox box     = validBox ? SubImageBox::intersect(thisBox, *validBox) : thisBox;
 
    // Return number of voxels in subimage. Return 0/0 if subimage is empty.
    const MLint numVoxels = box.getNumVoxels();
    if (numVoxels==0){ minValue=0; maxValue=0; return 0; }
 
    // Set min/max to first voxel value.
    minValue = maxValue = getImageValue(box.v1);
 
    // Simply scan all voxels. If box is subimage box, it is faster.
    if (box == thisBox){
      // Not empty. Get pointer to voxel directly after last voxel.
 
      ML_CHECK_THROW(_data);
      const DATATYPE* dataEnd = static_cast<DATATYPE*>(_data) + numVoxels;
 
      // Scan all voxels starting with the second one.
      for (DATATYPE* i=static_cast<DATATYPE*>(_data)+1; i<dataEnd; ++i)
      {
        if (*i < minValue){ minValue = *i; }
        if (*i > maxValue){ maxValue = *i; }
      }
    }
    else{
      // Loop over all rows of the valid subimage region.
      ImageVector p;
      for (MLint u=box.v1.u; u <= box.v2.u; ++u){
        for (MLint t=box.v1.t; t <= box.v2.t; ++t){
          for (MLint c=box.v1.c; c <= box.v2.c; ++c){
            for (MLint z=box.v1.z; z <= box.v2.z; ++z){
              for (MLint y=box.v1.y; y <= box.v2.y; ++y){
                // For inner loop, move pointer.
                p.set(box.v1.x, y, z, c, t, u);
                const DATATYPE *dPtr = getImagePointer(p);
                const MLint pEnd = box.v2.x;
                for (MLint x=box.v1.x; x <= pEnd; ++x, ++dPtr){
                  if (*dPtr < minValue){ minValue = *dPtr; }
                  if (*dPtr > maxValue){ maxValue = *dPtr; }
                }
              }
            }
          }
        }
      }
    }
 
    // Return how many voxel have been scanned.
    return numVoxels;
  };

 
  //---------------------------------------------------------------------------
  //---------------------------------------------------------------------------

#define _ML_CHECK_SUBIMAGE_DATA_POINTERS(FROM_PTR, TO_PTR)                            \
  /* Check for valid data pointers to avoid memory access errors. */                  \
  {                                                                                   \
    if (!FROM_PTR || !TO_PTR){                                                        \
      ML_PRINT_ERROR("TSubImage::copySubImageTyped( )",                                   \
                     ML_BAD_POINTER_OR_0,                                             \
                     "Valid data pointers in source and target subimage required for" \
                     "subimage copying, thus call is ignored.");                      \
      return;                                                                         \
    }                                                                                 \
  }

  template <typename FROM_DATATYPE>
  void copySubImageReorderColorPlanesToInterleaved(const TSubImage<FROM_DATATYPE>& typedFromImage, const ScaleShiftData& scaleShiftData)
  {
    // Copy in local variables to reduce calculations in getBox().
    const SubImageBox& box =getBox();
    const SubImageBox& tBox=typedFromImage.getBox();
 
    if (box.isEmpty() || tBox.isEmpty()){ return; }
 
    // Calculate overlap of both subimages and copy only if overlap is not empty.
    const SubImageBox intersection = SubImageBox::intersect(typedFromImage.getBox(),getBox());
    if (!intersection.isEmpty()){
 
      // Check data pointers and print errors if necessary.
      _ML_CHECK_SUBIMAGE_DATA_POINTERS(typedFromImage.getData(), getData());
 
      ImageVector p;
      const FROM_DATATYPE* fromPt       = nullptr;
      DATATYPE*            toPt         = nullptr;
 
      const ImageVector toExtent = getExtent();
      const ImageVector toStrides = ImageVector(toExtent.c, toExtent.x, toExtent.y, toExtent.z, toExtent.t, toExtent.u).getStrides();
      const ImageVector offsetReordered(_box.v1.c, _box.v1.x, _box.v1.y, _box.v1.z, _box.v1.t, _box.v1.u);
 
      const MLdouble  scaleDbl  = static_cast<MLdouble>(scaleShiftData.getScale());
      const MLdouble  shiftDbl  = static_cast<MLdouble>(scaleShiftData.getShift());
 
      // Check for identical scaling.
      const bool identicalScaling = MLValuesAreEqual(scaleDbl, static_cast<MLdouble>(1), static_cast<MLdouble>(1)) &&
                                    MLValuesAreEqual(shiftDbl, static_cast<MLdouble>(0), static_cast<MLdouble>(1));
 
      // Traverse all voxels.
      if (identicalScaling) {
        for (p.u=intersection.v1.u; p.u<=intersection.v2.u; ++p.u) {
          for (p.t=intersection.v1.t; p.t<=intersection.v2.t; ++p.t) {
            for (p.c=intersection.v1.c; p.c<=intersection.v2.c; ++p.c) {
              for (p.z=intersection.v1.z; p.z<=intersection.v2.z; ++p.z) {
                for (p.y=intersection.v1.y; p.y<=intersection.v2.y; ++p.y) {
 
                  const ImageVector pos(intersection.v1.x,p.y,p.z,p.c,p.t,p.u);
                  const ImageVector posReordered(p.c, intersection.v1.x,p.y,p.z,p.t,p.u);
                  fromPt= typedFromImage.getImagePointer(pos);
                  toPt  = static_cast<DATATYPE*>(_data) + _ML_TSUBIMG_SUBDOT6(posReordered, offsetReordered, toStrides);
 
                  for (p.x = intersection.v1.x; p.x <= intersection.v2.x; ++p.x, ++fromPt) {
                    *toPt = static_cast<DATATYPE>(*fromPt);
                    toPt += toStrides.y;
                  }
                }
              }
            }
          }
        }
      } else {
        for (p.u=intersection.v1.u; p.u<=intersection.v2.u; ++p.u) {
          for (p.t=intersection.v1.t; p.t<=intersection.v2.t; ++p.t) {
            for (p.c=intersection.v1.c; p.c<=intersection.v2.c; ++p.c) {
              for (p.z=intersection.v1.z; p.z<=intersection.v2.z; ++p.z) {
                for (p.y=intersection.v1.y; p.y<=intersection.v2.y; ++p.y) {
 
                  const ImageVector pos(intersection.v1.x,p.y,p.z,p.c,p.t,p.u);
                  const ImageVector posReordered(p.c, intersection.v1.x,p.y,p.z,p.t,p.u);
                  fromPt= typedFromImage.getImagePointer(pos);
                  toPt  = static_cast<DATATYPE*>(_data) + _ML_TSUBIMG_SUBDOT6(posReordered, offsetReordered, toStrides);
 
                  for (p.x = intersection.v1.x; p.x <= intersection.v2.x; ++p.x, ++fromPt) {
                    *toPt = static_cast<DATATYPE>(static_cast<MLdouble>(*fromPt)*scaleDbl + shiftDbl);
                    toPt += toStrides.y;
                  }
                }
              }
            }
          }
        }
      }
    } // if (!intersection.isEmpty())
  }

  //
  //Note: The following function definition must be implemented in the class definition
  //      See BUG: LNK2001 on Member Function When Use Nested Class Template, Article ID: Q128789
  template <typename FROM_DATATYPE>
  void copySubImageTyped(const TSubImage<FROM_DATATYPE> &typedFromImg, const ScaleShiftData& scaleShiftData)
  {
    // Any box empty? Then no copy is necessary.
    if (getBox().isEmpty() || typedFromImg.getBox().isEmpty()){
      return;
    }
 
    if (scaleShiftData.getReorderMode() == ScaleShiftData::ReorderColorPlanesToInterleaved) {
      copySubImageReorderColorPlanesToInterleaved(typedFromImg, scaleShiftData);
      return;
    }
 
    // Get double and long double versions of scale and shift values.
    const MLdouble  scaleDbl  = scaleShiftData.getScale();
    const MLdouble  shiftDbl  = scaleShiftData.getShift();
 
    // Check for identical scaling and use default copying then.
    if (MLValuesAreEqual(scaleDbl, 1.0, 1.0) && MLValuesAreEqual(shiftDbl, 0.0, 1.0)){
      copySubImage(typedFromImg);
      return;
    }
 
    // Get data type properties. If we calculate with normal datatypes then we calculate in double values;
    // on long double built-in types we use long double, and if we use extended types we remain in extended types.
    const MLint isDTScalarType  = MLIsScalarTypePtr(static_cast<DATATYPE*>(nullptr));
    const MLint isFDTScalarType = MLIsScalarTypePtr(static_cast<FROM_DATATYPE*>(nullptr));
 
    if (getDataType()==typedFromImg.getDataType()){
      if (getBox()==typedFromImg.getBox()){
        // Images have same data type and same position/extent:
        //-----------------------------------------------------
        const FROM_DATATYPE  *fromPt  = typedFromImg.getData();
        DATATYPE       *toPt    = getData();
 
        // Check data pointers and print errors if necessary.
        _ML_CHECK_SUBIMAGE_DATA_POINTERS(fromPt, toPt);
 
        // Get number and the size in bytes of all voxels of this subimage.
        const MLint size     = getNumVoxels();
        const MLint bytesize = size * MLSizeOf(getDataType());
 
        // Copy the memory to the target buffer, it must fit since types are identical.
        memcpy( toPt, fromPt, bytesize );
 
        // When we work with extended types we must avoid to change to double during calculations since
        // that causes information loss. In normal case we remain in double.
        if (!isDTScalarType){
          // NOTE/TODO:
          // Currently extended types must support adding and multiplying with doubles
          // for this to work.
          for (DATATYPE* toPtEnd=toPt+size; toPt<toPtEnd; ++toPt){
            *toPt=static_cast<DATATYPE>(static_cast<DATATYPE>(*toPt)*scaleDbl + shiftDbl);
          }
        }
        else{
          for (DATATYPE* toPtEnd=toPt+size; toPt<toPtEnd; ++toPt){
            *toPt=static_cast<DATATYPE>( static_cast<MLdouble>(*toPt)*scaleDbl + shiftDbl);
          }
        }
      }
      else{
        // Images have same data type:
        //----------------------------
        // Calculate intersection of both subimages. Copy only if intersection is not empty.
        const SubImageBox intersection = SubImageBox::intersect(typedFromImg.getBox(),getBox());
        if (!intersection.isEmpty()){
 
          // Check data pointers and print errors if necessary.
          _ML_CHECK_SUBIMAGE_DATA_POINTERS(typedFromImg.getData(), getData());
 
          ImageVector p;
          const FROM_DATATYPE *fromPt  = nullptr;
          DATATYPE       *toPt         = nullptr,
                         *toPt_        = nullptr,
                         *toPtEnd      = nullptr;
          const ImageVector    fromStr      = typedFromImg.getStride();
          const ImageVector    toStr        = getStride();
          const MLint dist       = intersection.v2.x - intersection.v1.x;
          const MLint size_x     = dist+1;
          const MLint bytesize_x = size_x * MLSizeOf(getDataType());
 
          // Traverse all 6 dimensions.
          for (p.u=intersection.v1.u; p.u<=intersection.v2.u; ++p.u){
            for (p.t=intersection.v1.t; p.t<=intersection.v2.t; ++p.t){
              for (p.c=intersection.v1.c; p.c<=intersection.v2.c; ++p.c){
                for (p.z=intersection.v1.z; p.z<=intersection.v2.z; ++p.z){
 
                  // Get start position to each slice of data.
                  const ImageVector pos(intersection.v1.x, intersection.v1.y, p.z, p.c, p.t, p.u);
                  fromPt = typedFromImg.getImagePointer(pos);
                  toPt   =              getImagePointer(pos);
                  ML_CHECK_THROW(fromPt);
                  ML_CHECK_THROW(toPt);
 
                  for (p.y=intersection.v1.y; p.y<=intersection.v2.y; ++p.y,fromPt+=fromStr.y,toPt+=toStr.y){
 
                    // Copy buffer to target buffer. Below we work only on one buffer for rescaling.
                    memcpy(toPt, fromPt, bytesize_x);
                    toPt_ = toPt;
 
                    if (!isDTScalarType){
                      // See description in upper case "if (getDataType()==typedFromImg.getDataType())".
                      for (toPtEnd = toPt_+size_x;
                           toPt_<toPtEnd;
                           ++toPt_){
                        *toPt_=static_cast<DATATYPE>(static_cast<DATATYPE>(*toPt_)*scaleDbl + shiftDbl);
                      }
                    }
                    else{
                      for (toPtEnd = toPt_+size_x;
                           toPt_<toPtEnd;
                           ++toPt_){
                        *toPt_=static_cast<DATATYPE>(static_cast<MLdouble>(*toPt_)*scaleDbl + shiftDbl);
                      }
                    }
                  }
                }
              }
            }
          }
        } // if (!intersection.isEmpty())
      } // else
    } // if (getDataType()==typedFromImg.getDataType())
 
    // Images have different data type:
    //---------------------------------
    else{
      // Copy overlap of both subimages and copy only if it is not empty.
      const SubImageBox intersection = SubImageBox::intersect(typedFromImg.getBox(),getBox());
      if (!intersection.isEmpty()){
 
        // Check data pointers and print errors if necessary.
        _ML_CHECK_SUBIMAGE_DATA_POINTERS(typedFromImg.getData(), getData());
 
        ImageVector p;
        const FROM_DATATYPE  *fromPt = nullptr;
        DATATYPE       *toPt         = nullptr;
 
        // Traverse all 6 dimensions.
        for (p.u=intersection.v1.u; p.u<=intersection.v2.u; ++p.u){
          for (p.t=intersection.v1.t; p.t<=intersection.v2.t; ++p.t){
            for (p.c=intersection.v1.c; p.c<=intersection.v2.c; ++p.c){
              for (p.z=intersection.v1.z; p.z<=intersection.v2.z; ++p.z){
                for (p.y=intersection.v1.y; p.y<=intersection.v2.y; ++p.y){
 
                  // Get start position to each slice of data.
                  const ImageVector pos(intersection.v1.x,p.y,p.z,p.c,p.t,p.u);
                  fromPt = typedFromImg.getImagePointer(pos);
                  toPt   =              getImagePointer(pos);
                  const MLint startX = intersection.v1.x;
                  const MLint endX   = intersection.v2.x;
 
                  if (!isFDTScalarType){
                    // When we work with extended types we must avoid to change to double during calculations since
                    // that causes information loss. In normal case we remain in double.
                    for (p.x = startX;  p.x <= endX;  ++p.x, ++fromPt, ++toPt){
                      *toPt=static_cast<DATATYPE>(static_cast<FROM_DATATYPE>(*fromPt)*scaleDbl + shiftDbl);
                    }
                  }
                  else{
                    for (p.x = startX;  p.x <= endX;  ++p.x, ++fromPt, ++toPt){
                      *toPt=static_cast<DATATYPE>(static_cast<MLdouble>(*fromPt)*scaleDbl + shiftDbl);
                    }
                  }
                }
              }
            }
          }
        } // for (p.u)
 
      } // if (!intersection.isEmpty())
 
    } // else of if (getDataType()==typedFromImg.getDataType())
  }
 

 
 
  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
 
  //---------------------------------------------------------------------------
  //---------------------------------------------------------------------------
  inline void fill(DATATYPE value)
  {
    ML_CHECK_THROW(_data);
 
    const DATATYPE* dataEnd = static_cast<DATATYPE*>(_data) + getNumVoxels();
 
    for (DATATYPE* i = static_cast<DATATYPE*>(_data); i<dataEnd; ++i){ *i=value; }
  };
 
 
  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  inline void fillBordersWithValue(const SubImageBox& box, DATATYPE fillValue)
  {
    SubImage::fillBordersWithTypeData(box, reinterpret_cast<const MLTypeData*>(&fillValue));
  }
 
  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  inline void fillInvalidRegionWithValue(DATATYPE value) {
    fillBordersWithValue(getValidRegion(), value);
  }
 
  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  inline void fillInvalidRegionWithBorderValues() {
    fillBordersWithBorderValues(getValidRegion());
  }

};
 
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template <typename DATATYPE>
class TSubImageCursorBase {
 
protected:
 
  TSubImageCursorBase(const TSubImage<DATATYPE>& subImage) {
    _data = const_cast<DATATYPE*>(subImage.getData());
    _subImageOffset = subImage.getBox().v1;
    _stride = subImage.getStride();
    _cursor = _data;
  }
 
public:
 
  //---------------------------------------------------------------------------
  //---------------------------------------------------------------------------
  inline void      setSubImagePosition(const ImageVector& position) { _cursor = _data+ImageVector::dot(position,_stride); }
  inline void      setSubImagePosition(MLint x, MLint y, MLint z) { _cursor = _data+_stride.x*x+_stride.y*y+_stride.z*z; }
  inline void      setImagePosition(const ImageVector& position) { _cursor = _data + _ML_TSUBIMG_SUBDOT6(position,_subImageOffset,_stride); }
  inline void      setImagePosition(MLint x, MLint y, MLint z)    { _cursor = _data+ + _ML_TSUBIMG_SUBDOT3(x,y,z,_subImageOffset,_stride); }
  inline void      setPosition(const DATATYPE *pointer) { _cursor = const_cast<DATATYPE*>(pointer); }
  inline void      moveByOffset(const ImageVector& offset) { _cursor += ImageVector::dot(offset,_stride); }
  inline void      moveByOffset(MLint x, MLint y, MLint z) { _cursor += _stride.x*x+_stride.y*y+_stride.z*z; }
  inline void      moveX() { _cursor += _stride.x; }
  inline void      moveY() { _cursor += _stride.y; }
  inline void      moveZ() { _cursor += _stride.z; }
  inline void      moveC() { _cursor += _stride.c; }
  inline void      moveT() { _cursor += _stride.t; }
  inline void      moveU() { _cursor += _stride.u; }
  inline void      reverseMoveX() { _cursor -= _stride.x; }
  inline void      reverseMoveY() { _cursor -= _stride.y; }
  inline void      reverseMoveZ() { _cursor -= _stride.z; }
  inline void      reverseMoveC() { _cursor -= _stride.c; }
  inline void      reverseMoveT() { _cursor -= _stride.t; }
  inline void      reverseMoveU() { _cursor -= _stride.u; }

 
  //---------------------------------------------------------------------------
  // Note: All these functions do not need explicit 64-bit arithmetic, as
  //       pointer arithmetic always moves from valid address area to another
  //       valid address area even on 32-bit systems, and calculations do not fall
  //       outside valid ranges.
  //---------------------------------------------------------------------------
  inline DATATYPE getValue() const { return *(_cursor); }

  inline DATATYPE getValueWithOffset(const ImageVector& offset) const { return *(_cursor+ImageVector::dot(offset,_stride)); }

  inline DATATYPE getValueWithOffset(MLint dx, MLint dy, MLint dz) const { return *(_cursor+_stride.x*dx+_stride.y*dy+_stride.z*dz); }

 
protected:
  DATATYPE* _data;
  ImageVector _subImageOffset;
  ImageVector _stride;
  DATATYPE* _cursor;
};
 
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template <typename DATATYPE>
class ConstTSubImageCursor : public TSubImageCursorBase<DATATYPE> {
 
typedef TSubImageCursorBase<DATATYPE> T;
 
public:
  ConstTSubImageCursor(const TSubImage<DATATYPE>& subImage) : TSubImageCursorBase<DATATYPE>(subImage) {
  }
 
  //---------------------------------------------------------------------------
  inline const DATATYPE* getPointer() const { return T::_cursor; }
  inline const DATATYPE* getPointerWithOffset(const ImageVector& offset) const { return T::_cursor+ImageVector::dot(offset,T::_stride); }
  inline const DATATYPE* getPointerWithOffset(MLint dx, MLint dy, MLint dz) const { return (T::_cursor+T::_stride.x*dx+T::_stride.y*dy+T::_stride.z*dz); }
};
 
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template <typename DATATYPE>
class TSubImageCursor : public TSubImageCursorBase<DATATYPE> {
 
typedef TSubImageCursorBase<DATATYPE> T;
 
public:
 
  TSubImageCursor(TSubImage<DATATYPE>& subImage) : TSubImageCursorBase<DATATYPE>(subImage) {
  }
 
  //---------------------------------------------------------------------------
  inline DATATYPE* getPointer() const { return T::_cursor; }
  inline DATATYPE* getPointerWithOffset(const ImageVector& offset) const { return T::_cursor+ImageVector::dot(offset,T::_stride); }
  inline DATATYPE* getPointerWithOffset(MLint dx, MLint dy, MLint dz) const { return (T::_cursor+T::_stride.x*dx+T::_stride.y*dy+T::_stride.z*dz); }
 
  //---------------------------------------------------------------------------
  // Note: All these functions do not need explicit 64-bit arithmetic, because
  //       pointer arithmetic always moves from valid address area to another
  //       valid address area even on 32-bit systems, and calculations do not fall
  //       outside valid ranges.
  //---------------------------------------------------------------------------
  inline void setValue(DATATYPE value) const { *(T::_cursor) = value; }
  inline void setValueWithOffset(const ImageVector& offset, DATATYPE value) const { *(T::_cursor+ImageVector::dot(offset,T::_stride))=value; }
  inline void setValueWithOffset(MLint dx, MLint dy, MLint dz, DATATYPE value) const { *(T::_cursor+T::_stride.x*dx+T::_stride.y*dy+T::_stride.z*dz)=value; }
};
 
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
template <typename DATATYPE>
class TSubImageWithCursor : public TSubImage<DATATYPE>
{
  typedef TSubImage<DATATYPE> T;
 
public:
  TSubImageWithCursor():TSubImage<DATATYPE>(), _cursor(nullptr)
  {
  }

  TSubImageWithCursor(const TSubImage<DATATYPE>& subImage):TSubImage<DATATYPE>(subImage), _cursor(nullptr)
  {
  }

  TSubImageWithCursor(const SubImage& subImage):TSubImage<DATATYPE>(subImage), _cursor(nullptr)
  {
  }

  TSubImageWithCursor(const TSubImageWithCursor& subImage):TSubImage<DATATYPE>(subImage),_cursor(subImage._cursor)
  {
  }

  inline TSubImageWithCursor &operator=(const TSubImageWithCursor<DATATYPE> &tSubImg)
  {
    TSubImage<DATATYPE>::operator=(tSubImg);
    _cursor = tSubImg._cursor;
    return *this;
  }
 
  //---------------------------------------------------------------------------
  //---------------------------------------------------------------------------
  inline void setCursorSubImagePosition(const ImageVector& position) { _cursor = static_cast<DATATYPE*>(T::_data)+ImageVector::dot(position,T::_stride); }
  inline void setCursorSubImagePosition(MLint x, MLint y, MLint z) { _cursor = static_cast<DATATYPE*>(T::_data)+T::_stride.x*x+T::_stride.y*y+T::_stride.z*z; }
  inline void setCursorImagePosition(const ImageVector& position) { _cursor = static_cast<DATATYPE*>(T::_data)+ _ML_TSUBIMG_SUBDOT6(position, T::_box.v1 ,T::_stride); }
  inline void setCursorImagePosition(MLint x, MLint y, MLint z) { _cursor = static_cast<DATATYPE*>(T::_data)+ _ML_TSUBIMG_SUBDOT3(x,y,z, T::_box.v1 ,T::_stride); }
  inline void setCursorPosition(const DATATYPE* pointer) { _cursor = const_cast<DATATYPE*>(pointer); }
  inline void moveCursorByOffset(const ImageVector& offset) { _cursor += ImageVector::dot(offset,T::_stride); }
  inline void moveCursorByOffset(MLint x, MLint y, MLint z) { _cursor += T::_stride.x*x+T::_stride.y*y+T::_stride.z*z; }
  inline void moveCursorX() { _cursor += T::_stride.x; }
  inline void moveCursorY() { _cursor += T::_stride.y; }
  inline void moveCursorZ() { _cursor += T::_stride.z; }
  inline void moveCursorC() { _cursor += T::_stride.c; }
  inline void moveCursorT() { _cursor += T::_stride.t; }
  inline void moveCursorU() { _cursor += T::_stride.u; }
  inline void reverseMoveCursorX() { _cursor -= T::_stride.x; }
  inline void reverseMoveCursorY() { _cursor -= T::_stride.y; }
  inline void reverseMoveCursorZ() { _cursor -= T::_stride.z; }
  inline void reverseMoveCursorC() { _cursor -= T::_stride.c; }
  inline void reverseMoveCursorT() { _cursor -= T::_stride.t; }
  inline void reverseMoveCursorU() { _cursor -= T::_stride.u; }

 
  //---------------------------------------------------------------------------
  // Note: All these functions do not need explicit 64-bit arithmetic, because
  //       pointer arithmetic always goes from valid address area to another
  //       valid address area even on 32-bit systems, and calculations do not fall
  //       outside valid ranges.
  //---------------------------------------------------------------------------
  inline DATATYPE getCursorValue() const { return *(_cursor); }
  inline DATATYPE getCursorValueWithOffset(const ImageVector& offset) const { return *(_cursor+ImageVector::dot(offset,T::_stride)); }
  inline DATATYPE getCursorValueWithOffset(MLint dx, MLint dy, MLint dz) const { return *(_cursor+T::_stride.x*dx+T::_stride.y*dy+T::_stride.z*dz); }

 
  //---------------------------------------------------------------------------
  inline DATATYPE* getCursorPointer() const { return _cursor; }
  inline DATATYPE* getCursorPointerWithOffset(const ImageVector& offset) const { return _cursor+ImageVector::dot(offset,T::_stride); }
  inline DATATYPE* getCursorPointerWithOffset(MLint dx, MLint dy, MLint dz) const { return (_cursor+T::_stride.x*dx+T::_stride.y*dy+T::_stride.z*dz); }
 
  //---------------------------------------------------------------------------
  // Note: All these functions do not need explicit 64-bit arithmetic, because
  //       pointer arithmetic always goes from valid address area to another
  //       valid address area even on 32-bit systems, and calculations do not fall
  //       outside valid ranges.
  //---------------------------------------------------------------------------
  inline void setCursorValue(DATATYPE value) { *(_cursor) = value; }
  inline void setCursorValueWithOffset(const ImageVector& offset, DATATYPE value) { *(_cursor+ImageVector::dot(offset,T::_stride))=value; }
  inline void setCursorValueWithOffset(MLint dx, MLint dy, MLint dz, DATATYPE value) { *(_cursor+T::_stride.x*dx+T::_stride.y*dy+T::_stride.z*dz)=value; }

 
private:
  DATATYPE* _cursor;
};
 
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
template <typename T>
TSubImage<T>& tsubimage_cast(SubImage& subImg)
{
  // Check whether template type and real data type match if there is valid data.
  if (subImg.getData() && !TypeTraits<T>::matches(subImg.getDataType()) )
  {
    ML_PRINT_FATAL_ERROR("tsubimage_cast", ML_BAD_DATA_TYPE,
      "Mismatch of data type between argument and template typename.");
  }
  return *(static_cast<TSubImage<T>* >(&subImg));
}
 
template <typename T>
const TSubImage<T>& tsubimage_cast(const SubImage& subImg)
{
  // Check whether template type and real data type match if there is valid data.
  if (subImg.getData() && !TypeTraits<T>::matches(subImg.getDataType()) )
  {
    ML_PRINT_FATAL_ERROR("tsubimage_cast", ML_BAD_DATA_TYPE,
      "Mismatch of data type between argument and template typename.");
  }
  return *(static_cast<const TSubImage<T>* >(&subImg));
}
 
template <typename T>
TSubImage<T>* tsubimage_cast(SubImage* subImg)
{
  // Check whether template type and real data type match if there is valid data.
  if (subImg->getData() && !TypeTraits<T>::matches(subImg->getDataType()) )
  {
    ML_PRINT_FATAL_ERROR("tsubimage_cast", ML_BAD_DATA_TYPE,
      "Mismatch of data type between argument and template typename.");
  }
  return static_cast<TSubImage<T>* >(subImg);
}
 
template <typename T>
const TSubImage<T>* tsubimage_cast(const SubImage* subImg)
{
  // Check whether template type and real data type match if there is valid data.
  if (subImg->getData() && !TypeTraits<T>::matches(subImg->getDataType()) )
  {
    ML_PRINT_FATAL_ERROR("tsubimage_cast", ML_BAD_DATA_TYPE,
      "Mismatch of data type between argument and template typename.");
  }
  return static_cast<const TSubImage<T>* >(subImg);
}

 
 
ML_END_NAMESPACE
 
namespace std {
 
  //---------------------------------------------------------------------------
  //         Debug printing for TSubImage.
  //---------------------------------------------------------------------------

  template <typename DATATYPE>
  inline ostream &operator<<(ostream &ostr, const ML_NAMESPACE::TSubImage<DATATYPE> &v)
  {
    // Get box of subimage and traverse all of its voxels.
    const ML_NAMESPACE::SubImageBox& box = v.getBox();
    ML_NAMESPACE::ImageVector p;
    for (p.u=box.v1.u; p.u<=box.v2.u; ++p.u){
      for (p.t=box.v1.t; p.t<=box.v2.t; ++p.t){
        for (p.c=box.v1.c; p.c<=box.v2.c; ++p.c){
          for (p.z=box.v1.z; p.z<=box.v2.z; ++p.z){
            for (p.y=box.v1.y; p.y<=box.v2.y; ++p.y){
              // Print row of voxels after a single row start position.
              p.x = box.v1.x;
              ostr << p << ": ";
              for (; p.x<=box.v2.x; p.x++){
                ostr << v.getImageValue(p) << " ";
              }
              ostr << std::endl;
            }
          }
        }
      }
    }
    return ostr;
  }
 
}
 
#ifdef _MSC_VER
#pragma warning(pop)
#endif
 
 
#endif //of __mlTSubImage_H