mlKernel.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_KERNEL_H
#define ML_KERNEL_H
 
 
 
// ML-includes
#include "mlInitSystemKernel.h"
#include <mlMemory.h>
#include <cerrno>
 
ML_START_NAMESPACE
 
  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE> class TKernel {
 
  public:
 
    //--------------------------------------------------------------------------------------------------------
    //--------------------------------------------------------------------------------------------------------
    enum KernModifier {
      KERN_SET       = 0,
      KERN_MULT         ,
      KERN_ADD          ,
      KERN_INVDIV       ,
      KERN_INVSUB       ,
      KERN_SQR          ,
      KERN_SQRT         ,
      KERN_POW          ,
      KERN_LOG          ,
      KERN_NORMALIZE    ,
      KERN_GAUSS        ,
      KERN_SPHERE       ,
      KERN_MIRROR       ,
      KERN_MIRROR_X     ,
      KERN_MIRROR_Y     ,
      KERN_MIRROR_Z     ,
      KERN_MIRROR_C     ,
      KERN_MIRROR_T     ,
      KERN_MIRROR_U     ,
 
      NUM_KERN_MODIFYERS
    };
 
    inline TKernel();
 
    inline TKernel(const TKernel &kern);
 
    inline virtual ~TKernel();
 
    const TKernel& operator=(const TKernel &kern);
 
    void reset();
 
    inline size_t getTabSize() const;
 
    inline const ImageVector* getCoordTab(size_t dim = 0) const;
 
    inline const KDATATYPE* getValueTab(size_t dim = 0) const;
 
    KDATATYPE getValueTabSum() const;
 
    KDATATYPE getMinValue() const;
 
    KDATATYPE getMaxValue() const;
 
    KDATATYPE getNegativeValueSum() const;
 
    KDATATYPE getPositiveValueSum() const;
 
    void manipulateKernelElements(KernModifier mode, KDATATYPE v);
 
    inline const ImageVector &getExtent() const;
 
    inline const ImageVector &getNegativeExtent() const;
 
    inline const ImageVector &getPositiveExtent() const;
 
    static inline MLint coordToIndex(MLint x, MLint y, MLint z, MLint c, MLint t, MLint u, const ImageVector &size);
 
    static inline MLint coordToIndex(const ImageVector &p, const ImageVector &size);
 
    static inline ImageVector indexToCoord(MLint idx, const ImageVector &ext);
 
    static inline ImageVector calculateNegativeExtentFromExtent(const ImageVector &ext);
 
    static inline ImageVector calculatePositiveExtentFromExtent(const ImageVector &ext);
 
    MLint findIndex(const ImageVector &pos) const;
 
    void setPartialKernel(size_t numElems, ImageVector * const coords, KDATATYPE * const values=nullptr);
 
    std::string getKernel(bool asLines=false, MLint fieldWidth=0, MLint precision=10) const;
 
    std::string setKernel(const std::string &kernString);
 
    void setKernel(const ImageVector &ext, const KDATATYPE * const values=nullptr, bool * const mask=nullptr);
 
    template <typename KDATATYPE2>
      void setKernel(const ImageVector &ext,
                     const KDATATYPE2 * const xaxis, const KDATATYPE2 * const yaxis, const KDATATYPE2 * const zaxis,
                     const KDATATYPE2 * const caxis, const KDATATYPE2 * const taxis, const KDATATYPE2 * const uaxis,
                     bool normalize)
    {
      // Clear current separability members.
      _invalidateSeparableInfos();
 
      MLint b=0;
 
      // Check for valid parameters. Return empty kernel otherwise.
      if (!xaxis || !yaxis || !zaxis || !caxis || !taxis || !uaxis){
        _clearTables();
        _calculateRealKernelExt();
        return;
      }
 
      // Pointer to matrix which contains the products of the corresponding axis-values
      KDATATYPE *valuematrix =
       static_cast<KDATATYPE*>(Memory::allocateMemory(sizeof(KDATATYPE) * ext.compMul(), ML_FATAL_MEMORY_ERROR));
 
      for (MLint u=0; u<ext.u; ++u){
        for (MLint t=0; t<ext.t; ++t){
          for (MLint c=0; c<ext.c; ++c){
            for (MLint z=0; z<ext.z; ++z){
              for (MLint y=0; y<ext.y; ++y){
                for (MLint x=0; x<ext.x; ++x){
                  valuematrix[b++] = static_cast<KDATATYPE>(xaxis[x]*yaxis[y]*zaxis[z]*caxis[c]*taxis[t]*uaxis[u]);
                }
      }}}}}
 
      setKernel(ext,valuematrix);
 
      Memory::freeMemory(valuematrix);
      valuematrix=nullptr;
 
      // If desired then normalize all kernel elements so that their sum becomes 1.
      if (normalize){ manipulateKernelElements(KERN_NORMALIZE, 0); }
    }
 
    void setSeparableKernel(const std::vector<typename std::vector<KDATATYPE> > &separableRows);
 
    void resizeKernel(const ImageVector &ext, KDATATYPE newVal=0);
 
    void fillUndefinedElements(KDATATYPE newVal=0);
 
    void makeCircular();
 
    void mirror(int dimension=-1);
 
    static MLldouble binomialcoeff(MLint n, MLint k);
 
    static std::vector<KDATATYPE> get1DGauss(size_t numSamples, bool normalize=true);
 
    void setGauss(const ImageVector &ext);
 
 
    //-------------------------------------------------------------------------------------------
    //
    //
    // If the kernel is interpreted as a separable kernel then the first six rows
    // of the first slice of the kernel matrix are taken.
    // These six rows are interpreted as the six 1-dimensional axis in x,y,z,c,t, and u
    // dimensions which are used as 1-D filter kernels for a 6D separable filtering.
    //
    //
    // For example:
    //
    //  | 1 3 5 6 4 |
    //  | 2 7 2     |
    //  | 4 1 9 5   |
    //  |           |
    //  | 3 8 3     |
    //  |           |
    //
    // describes a separable kernel with the following axes extent:
    // X=5, Y=3, Z=4, C=0, T=3, U=0.
    // An extent of 0 (for example C and U dimension) is used to describe
    // that filtering with a 1-D kernel in that dimension is not applied.
    //
    //-------------------------------------------------------------------------------------------
    inline void         setSeparable(bool isSeparableVal);
 
    inline bool         isSeparable() const;
 
    inline ImageVector  getSeparableDimEntries() const;
 
    inline ImageVector  getSeparableOneDimExtents() const;
 
    size_t              getSeparableDimIndex(size_t dim=0) const;
 
    const std::vector<ImageVector> &getSeparableCoordTab() const;
 
 
  protected:
    //-------------------------------------------------------------------------------------------
    //
    //                           Protected methods.
    //
    //-------------------------------------------------------------------------------------------
 
    void        _init();
 
    void        _clearTables();
 
    void        _addCoordinate(const ImageVector &pos, KDATATYPE value, bool update);
 
    void        _calculateRealKernelExt();
 
    void        _validateSeparabilityInfos() const;
 
    void        _invalidateSeparableInfos();
 
 
 
  private:
    //-------------------------------------------------------------------------------------------
    //
    //                            Member variables.
    //
    //-------------------------------------------------------------------------------------------
 
    ImageVector    _ext;
 
    ImageVector    _negExt;
 
    ImageVector    _posExt;
 
    size_t         _tabSize;
 
    KDATATYPE*     _valueTab;
 
    ImageVector*   _coordTab;
 
 
 
    bool           _isSeparable;
 
 
    mutable bool   _areSeparabilityInfosUpToDate;
 
    mutable ImageVector _separableDimEntries;
 
    mutable ImageVector _separableOneDimExt;
 
    mutable ImageVector _separableExt;
 
    mutable ImageVector _separableNegExt;
 
    mutable ImageVector _separablePosExt;
 
    mutable std::vector<ImageVector> _separableCoordTab;
  };  // class TKernel
 
 
 
 
 
  //------------------------------------------------------------------------
  // IMPLEMENTATION PART: No user specific things below.
  //------------------------------------------------------------------------
 
 
  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    TKernel<KDATATYPE>::TKernel()
  {
    _init();
  }
 
  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    TKernel<KDATATYPE>::TKernel(const TKernel &kern)
  {
    _init();
    operator=(kern);
  }
 
  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    TKernel<KDATATYPE>::~TKernel()
  {
    _clearTables();
  }
 
  //-------------------------------------------------------------------------------------------
  // Assignment operator. Sets current instance to contents of \c kern.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    const TKernel<KDATATYPE>& TKernel<KDATATYPE>::operator=(const TKernel &kern)
  {
    // Assign kernel to this if it's not this.
    if (&kern != this){
      // Clean up current kernel table.
      _clearTables();
 
      // Set new kernel table size, the kernel voxels and the kernel values.
      setPartialKernel(kern._tabSize, kern._coordTab, kern._valueTab);
 
      // Copy mutable members for separable kernel support.
      _isSeparable                  = kern._isSeparable;
      _areSeparabilityInfosUpToDate = kern._areSeparabilityInfosUpToDate;
      _separableDimEntries          = kern._separableDimEntries;
      _separableOneDimExt           = kern._separableOneDimExt;
      _separableExt                 = kern._separableExt;
      _separableNegExt              = kern._separableNegExt;
      _separablePosExt              = kern._separablePosExt;
      _separableCoordTab            = kern._separableCoordTab;
    }
 
    // Return copy.
    return *this;
  }
 
  //-------------------------------------------------------------------------------------------
  // Reset kernel to construction state.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::reset()
  {
    // Remove dynamic structures.
    _clearTables();
 
    // Call constructor initialization. We can do that here, because
    // all internal tables have been cleared.
    _init();
  }
 
  //-------------------------------------------------------------------------------------------
  // Returns the current number of kernel elements.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    size_t TKernel<KDATATYPE>::getTabSize() const
  {
    if (_isSeparable){
      _validateSeparabilityInfos();
      return getSeparableCoordTab().size();
    }
    else{
      return _tabSize;
    }
  }
 
  //-------------------------------------------------------------------------------------------
  // Gets a table of coordinates pointing to all kernel elements
  // which are currently defined.
  // The size of the returned table is given by getTabSize().
  // In the case of separable filtering it returns those which are
  // valid for the current pass of separable kernel filtering.
  // It is a subset of the array given by \c getSeparableCoordTab().
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    const ImageVector* TKernel<KDATATYPE>::getCoordTab(size_t dim) const
  {
    return _isSeparable ? &(getSeparableCoordTab()[getSeparableDimIndex(dim)]) : _coordTab;
  }
 
  //-------------------------------------------------------------------------------------------
  // Gets the table of kernel element values which are currently defined.
  // The size of the returned table is given by getTabSize().
  // In the case of separable filtering the subset of kernel elements
  // for the pass \c dim can be requested.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    const KDATATYPE* TKernel<KDATATYPE>::getValueTab(size_t dim) const
  {
    return _isSeparable ? (_valueTab + getSeparableDimIndex(dim)) : _valueTab;
  }
 
  //-------------------------------------------------------------------------------------------
  // Return the sum of all kernel element values
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    KDATATYPE TKernel<KDATATYPE>::getValueTabSum() const
  {
    // Scan all elements of the value table and sum them up.
    KDATATYPE v=0;
    for (size_t i=0; i<_tabSize; i++){ v+=_valueTab[i]; }
    return v;
  }
 
  //-------------------------------------------------------------------------------------------
  // Return the minimum value of all kernel element values
  // If kernel value table is empty then 0 is returned.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    KDATATYPE TKernel<KDATATYPE>::getMinValue() const
  {
    // Scan all elements of the value table and search the minimum.
    KDATATYPE v = (_tabSize==0) ? static_cast<KDATATYPE>(0) : _valueTab[0];
    for (size_t i=1; i<_tabSize; i++){ if (_valueTab[i] < v){ v = _valueTab[i]; } }
    return v;
  }
 
  //-------------------------------------------------------------------------------------------
  // Return the maximum value of all kernel element values
  // If kernel value table is empty then 1 is returned.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    KDATATYPE TKernel<KDATATYPE>::getMaxValue() const
  {
    // Scan all elements of the value table and search the minimum.
    KDATATYPE v = (_tabSize==0) ? static_cast<KDATATYPE>(1) : _valueTab[0];
    for (size_t i=1; i<_tabSize; i++){ if (_valueTab[i] > v){ v = _valueTab[i]; } }
    return v;
  }
 
 
  //-------------------------------------------------------------------------------------------
  // Return the sum of all negative kernel element values.
  // If kernel value table is empty then 0 is returned.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    KDATATYPE TKernel<KDATATYPE>::getNegativeValueSum() const
  {
    KDATATYPE v = 0;
    for (size_t i=0; i<_tabSize; i++){
      // Make checks for not "> 0" instead of "< 0" to avoid warnings on unsigned KDATATYPEs.
      if (!(_valueTab[i] > 0)){ v += _valueTab[i]; }
    }
    return v;
  }
 
  //-------------------------------------------------------------------------------------------
  // Return the sum of all positive kernel element values.
  // If kernel value table is empty then 0 is returned.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    KDATATYPE TKernel<KDATATYPE>::getPositiveValueSum() const
  {
    KDATATYPE v = 0;
    for (size_t i=0; i<_tabSize; i++){ if (_valueTab[i] > 0){ v += _valueTab[i]; } }
    return v;
  }
 
 
  //-------------------------------------------------------------------------------------------
  // Modify all kernel element values with the value \c v.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::manipulateKernelElements(KernModifier mode, KDATATYPE v)
  {
    // Invalidate current information for separable kernels.
    _invalidateSeparableInfos();
 
    switch (mode){
    case KERN_SET:{
      // Set each kernel element to v.
      for (size_t i=0; i<_tabSize; i++){ _valueTab[i] = v; }
    }
    break;
 
    case KERN_MULT:{
      // Multiply each kernel element with v.
      for (size_t i=0; i<_tabSize; i++){ _valueTab[i] *= v; }
    }
    break;
 
    case KERN_ADD:{
      // Add v to each kernel element.
      for (size_t i=0; i<_tabSize; i++){ _valueTab[i] += v; }
    }
    break;
 
    case KERN_INVDIV:{
      // Set each kernel element to v / kernelElement. Zero kernel elements are left unchanged.
      for (size_t i=0; i<_tabSize; i++){ if (!MLValueIs0WOM(_valueTab[i])) {_valueTab[i] = v/_valueTab[i]; } }
    }
    break;
 
    case KERN_INVSUB:{
      // Set each kernel element by v - kernelElement.
      for (size_t i=0; i<_tabSize; i++){ _valueTab[i] = v-_valueTab[i]; }
    }
    break;
 
    case KERN_SQR:{
      // Compute all squares of _valueTab[i].
      for (size_t i=0; i<_tabSize; i++){ _valueTab[i] *= _valueTab[i]; }
    }
    break;
 
    case KERN_SQRT:{
      // Compute all square roots of _valueTab[i].  Negative kernel elements
      // are left unchanged.
      for (size_t i=0; i<_tabSize; i++){
        // We do not use >= because but two comparisons to avoid warnings
        // in case of unsigned KDATATYPE.
        if ((_valueTab[i]>0) || (MLValueIs0WOM(_valueTab[i]))){
          // As documented in class description we cast in case of integer KDATATYPES.
          _valueTab[i] = static_cast<KDATATYPE>(sqrt(static_cast<double>(_valueTab[i])));
        }
      }
    }
    break;
 
    case KERN_POW:{
      // Compute all _valueTab[i] raised to the power of v.
      for (size_t i=0; i<_tabSize; i++){
        // Linux man page says: Error on "The  argument x is negative and y is not an integral value.
        // This would result in a complex number." This then will set errno (EDOM).
        // So we leave kernel elements unchanged in that case. We do not use >= because but
        // two comparisons to avoid warnings in case of unsigned KDATATYPE.
        if ( (_valueTab[i] > static_cast<KDATATYPE>(0)) || MLValueIs0WOM(_valueTab[i]) ){
          // As documented in class description we cast in case of integer KDATATYPES.
          errno = 0;
          KDATATYPE buf = static_cast<KDATATYPE>(powl(_valueTab[i], v));
          if (0 == errno){
            _valueTab[i] = buf;
          }
        }
      }
    }
    break;
 
    case KERN_LOG:{
      // Compute logarithm of base v of all kernel elements.
      // Do not apply operations on kernel elements which are not permitted,
      // e.g. negative bases or negative values or base == 1. In that case
      // leave kernel elements unchanged.
      // For LOGv(K) we use log(K)/log(v) instead.
      if ((v > 0) && (MLValuesDifferWOM(v, static_cast<KDATATYPE>(1))) ){
        MLdouble invBase=1.0/log(static_cast<double>(v));
        for (size_t i=0; i<_tabSize; i++){
          KDATATYPE val = _valueTab[i];
          // As documented in class description we cast in case of integer KDATATYPES.
          if (val > 0){ _valueTab[i] = static_cast<KDATATYPE>(log(static_cast<double>(val))*invBase); }
        }
      }
    }
    break;
 
    case KERN_NORMALIZE:{
      // Multiply all kernel element values with a value so that their sum is 1.
      // If sum is zero then values are left unchanged.
      KDATATYPE sum=getValueTabSum();
      // As documented in class description we cast in case of integer KDATATYPES.
      if (!MLValueIs0WOM(sum)){ manipulateKernelElements(KERN_MULT, static_cast<KDATATYPE>(1./sum)); }
    }
    break;
 
    case KERN_GAUSS:{
      // Set all kernel elements to binomial values and normalize.
      setGauss(_ext);
    }
    break;
 
    case KERN_SPHERE:{
      // Throw away kernel corners to make it approximately spherical.
      makeCircular();
    }
    break;
 
    case KERN_MIRROR:{
      // Apply mirroring operation.
      mirror();
    }
    break;
 
    case KERN_MIRROR_X:{
      // Apply mirroring operation in x dimension.
      mirror(0);
    }
    break;
 
    case KERN_MIRROR_Y:{
      // Apply mirroring operation in y dimension.
      mirror(1);
    }
    break;
 
    case KERN_MIRROR_Z:{
      // Apply mirroring operation in z dimension.
      mirror(2);
    }
    break;
 
    case KERN_MIRROR_C:{
      // Apply mirroring operation in c dimension.
      mirror(3);
    }
    break;
 
    case KERN_MIRROR_T:{
      // Apply mirroring operation in t dimension.
      mirror(4);
    }
    break;
 
    case KERN_MIRROR_U:{
      // Apply mirroring operation in u dimension.
      mirror(5);
    }
    break;
 
    case NUM_KERN_MODIFYERS:{
      ML_PRINT_ERROR("TKernel<KDATATYPE>::manipulateKernelElements(mlKernModifier mode, KDATATYPE v)",
                     ML_PROGRAMMING_ERROR, "NUM_KERN_MODIFYERS is an invalid parameter. Ignoring kernel manipulation.");
    }
    break;
 
    default:{
      ML_PRINT_ERROR("TKernel<KDATATYPE>::manipulateKernelElements(mlKernModifier mode, KDATATYPE v)",
                     ML_BAD_PARAMETER, "Invalid parameter passed. Ignoring kernel manipulation.");
    }
    break;
 
    }
  }
 
  //-------------------------------------------------------------------------------------------
  // Returns the kernel extents in 6D. It defines the rectangular region in
  // which all coordinates returned by \c getCoordTab() are found.
  // Note that the returned region might be larger than required
  // e.g. after removing elements from kernel.
  // In the case of separable filtering it returns a vector where the
  // components 0,...,5 contain the extent of the region spanned by the
  // 6 separated 1-D kernels. Note that components for those dimensions
  // where no filtering takes place are set to 1.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    const ImageVector &TKernel<KDATATYPE>::getExtent() const
  {
    if (_isSeparable){
      _validateSeparabilityInfos();
      return _separableExt;
    }
    else{
      return _ext;
    }
  }
 
  //-------------------------------------------------------------------------------------------
  // The extent of the kernel to both sides. The sum of both +1 is the extent of the kernel.
  // Using \c getNegativeExtent() as negative extent of an image and \c getPositiveExtent() as positive extent
  // increment for the image then the kernel can be placed correctly on all normal image
  // voxels without having voxel accesses out of range.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    const ImageVector &TKernel<KDATATYPE>::getNegativeExtent() const
  {
    if (_isSeparable){
      _validateSeparabilityInfos();
      return _separableNegExt;
    }
    else{
      return _negExt;
    }
  }
 
  //-------------------------------------------------------------------------------------------
  // See \c getNegativeExtent().
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    const ImageVector &TKernel<KDATATYPE>::getPositiveExtent() const
  {
    if (_isSeparable){
      _validateSeparabilityInfos();
      return _separablePosExt;
    }
    else{
      return _posExt;
    }
  }
 
  //-------------------------------------------------------------------------------------------
  // Converts the coordinate (x,y,z,c,t,u) into the kernel to an index into an array
  // with 6D extents given by \c size.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    MLint TKernel<KDATATYPE>::coordToIndex(MLint x, MLint y, MLint z, MLint c, MLint t, MLint u, const ImageVector &size)
  {
    return SubImage::coordToIndex(x,y,z,c,t,u, size);
  }
 
  //-------------------------------------------------------------------------------------------
  // Converts the coordinate \p into the kernel with extents \c size to an index.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    MLint TKernel<KDATATYPE>::coordToIndex(const ImageVector &p, const ImageVector &size)
  {
    return SubImage::coordToIndex(p, size);
  }
 
  //-------------------------------------------------------------------------------------------
  // Converts an index into an array with extents \c ext to a coordinate.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    ImageVector TKernel<KDATATYPE>::indexToCoord(MLint idx, const ImageVector &ext)
  {
    return SubImage::indexToCoord(idx, ext);
  }
 
  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    ImageVector TKernel<KDATATYPE>::calculateNegativeExtentFromExtent(const ImageVector &ext)
  {
    return ImageVector(ext.x/2, ext.y/2, ext.z/2, ext.c/2, ext.t/2, ext.u/2);
  }
 
  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    ImageVector TKernel<KDATATYPE>::calculatePositiveExtentFromExtent(const ImageVector &ext)
  {
    return ImageVector(ext.x-ext.x/2-1,  ext.y-ext.y/2-1,  ext.z-ext.z/2-1,
                  ext.c-ext.c/2-1,  ext.t-ext.t/2-1,  ext.u-ext.u/2-1);
  }
 
  //-------------------------------------------------------------------------------------------
  // Return index to kernel element if it exists; otherwise return -1.
  // Note that this method needs to search; so it's not very efficient.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    MLint TKernel<KDATATYPE>::findIndex(const ImageVector &pos) const
  {
    // Search coordinate in _coordtab.
    for (size_t c=0; c < _tabSize; c++){
      if (_coordTab[c] == pos){
        return c;
      }
    }
 
    return -1;
  }
 
  //-------------------------------------------------------------------------------------------
  // Defines a set of local kernel coordinates and optionally a set
  // of values which define the kernel elements.
  // \c numElems defines the number of coordinates.
  // If desired the corresponding set of kernel values can be passed in \c values.
  // Otherwise all kernel values are set to (KDATATYPE)(1.0/(KDATATYPE)_tabSize).
  // Note that \c coords and \c values have to contain \c numElems entries if passed.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::setPartialKernel(size_t numElems, ImageVector * const coords, KDATATYPE * const values)
  {
    // Clear all dynamic tables.
    _clearTables();
 
    // Define the new number of kernel voxels.
    _tabSize = numElems;
 
    // Recreate the coordinate and value tables
    if (_tabSize > 0)      {
      _coordTab = new ImageVector[_tabSize];
 
      _valueTab = new KDATATYPE[_tabSize];
 
      // Copy all passed coordinates into the coordinate table.
      for (size_t c=0; c < _tabSize; c++){ _coordTab[c] = coords[c]; }
 
      // If we have a user defined value table then copy it.
      // Otherwise fill table with values 1.0/static_cast<KDATATYPE>(_tabSize).
      if (values){
        memcpy(_valueTab, values, _tabSize*sizeof(KDATATYPE));
      }
      else{
        // As documented in class description we cast in case of integer KDATATYPES.
        KDATATYPE v = static_cast<KDATATYPE>(1.0/static_cast<KDATATYPE>(_tabSize));
        for (size_t d=0; d<_tabSize; d++){ _valueTab[d] = v; }
      }
    }
 
    // Calculate the real kernel size, i.e. the bounding box around all
    // used kernel elements.
    _calculateRealKernelExt();
  }
 
 
  //-------------------------------------------------------------------------------------------
  // Returns the current kernel elements and values as string:
  // The string needs has the following format:
  //
  // \code
  //
  // (x_0, y_0, z_0, c_0, t_0, u_0):v_0
  // ...
  // (x_n, y_n, z_n, c_n, t_n, u_n):v_n
  //
  // \endcode
  //
  // where the coordinates left from ':' specify the 6d coordinate
  // of the kernel element; the value v_x right from ':' specifies the
  // value of the kernel element.
  // If asLines is true then a second string format is used:
  //
  // \code
  //
  // (*, y_0, z_0, c_0, t_0, u_0):x_0 ... x_n
  // ...
  // (*, y_n, z_n, c_n, t_n, u_n):y_n ... y_n
  //
  // \endcode
  //
  // So all kernel elements of a row are saved in a string line; so
  // the x coordinates becomes invalid and is set as an asterisk.
  // To have a minimum field width pass \c fieldWidth and the digits after
  // the period if given by \c precision. Note that these settings are used
  // only if asLines is true.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    std::string TKernel<KDATATYPE>::getKernel(bool asLines, MLint fieldWidth, MLint precision) const
  {
    std::string str="";
    char strLine[256]="";
 
    // Clamp field width to sensible values.
    if (fieldWidth > 128){ fieldWidth=128; }
    if (fieldWidth <   0){ fieldWidth=0;   }
 
    if (asLines){
      // Print kernel as lines.
      for (MLint u=0; u < getExtent().u; u++){
        for (MLint t=0; t < getExtent().t; t++){
          for (MLint c=0; c < getExtent().c; c++){
            for (MLint z=0; z < getExtent().z; z++){
              for (MLint y=0; y < getExtent().y; y++){
                // Print line start into string.
                snprintf(strLine, 255, "(*,%lld,%lld,%lld,%lld,%lld):", y, z, c, t, u);
 
                // To compensate line start differences if any component has more than
                // one digit do append spaces until string has at least 16 chars. So
                // two components may have more than one digit or one can have two more digits.
                while (strlen(strLine) < 16){ strcat(strLine, " "); }
 
                // Add Line start to existing kernel string.
                str += strLine;
 
                // Add all row elements of the kernel to the string.
                for (MLint x=0; x < getExtent().x; x++){
                  // Does the element exist?
                  MLint idx = findIndex(ImageVector(x,y,z,c,t,u));
                  if (idx < 0){
                    // Kernel element does not exist.
                    // Fill field with fieldWidth spaces and add ','. Then terminate string.
                    // For the first field do not add a comma. So we have one more space then
                    // in those cases where we have a comma in it. So we need to add one space
                    // more in cases with comma.
                    strLine[0] = x==0 ? ' ' : ',';
                    MLint w= x==0 ? fieldWidth : fieldWidth+1;
                    for (MLint i=1; i < w; i++){ strLine[i] = ' '; }
                    strLine[w] = 0;
                  }
                  else{
                    // Print comma and element if it's not the last one; otherwise print only the element.
                    // integrate field width for float field of aligned size. Always cast
                    // element values to double because string reading and parsing of extended
                    // types would lead to further problems.
                    // Hence we store only scalar values.
                    char format[64];
                    snprintf(format, 63, "%s%%%d.%dlf", "%s", static_cast<MLint32>(fieldWidth), static_cast<MLint32>(precision));
                    snprintf(strLine, 255, format, 0==x ? "" : ",", static_cast<MLdouble>(_valueTab[idx]));
                  }
 
                  // Append (empty) element to the end of the line.
                  str += strLine;
                }
                // Terminate line with '\n'
                str += '\n';
              }
              // Add an empty line between xy planes of the kernel.
              str += '\n';
            }  // z
          }  // c
        }  // t
      } // u
    }
    else{
      // Print each kernel element into a string and append it to the entire kernel string.
      for (size_t c=0; c < _tabSize; c++){
        ImageVector coord = _coordTab[c];
        KDATATYPE value = _valueTab[c];
        snprintf(strLine, 255, "%s:%lf\n", coord.print("(", ",", ")").c_str(), static_cast<MLdouble>(value));
        str += strLine;
      }
    }
 
    return str;
  }
 
 
 
  //-------------------------------------------------------------------------------------------
  // Defines elements and values of the kernel matrix by a string.
  // The string needs to have the following format:
  //
  // \code
  //
  // (x_0, y_0, z_0, c_0, t_0, u_0):v_0
  // ...
  // (x_n, y_n, z_n, c_n, t_n, u_n):v_n
  //
  // \endcode
  //
  // where the coordinates left from ':' specify the 6d coordinate
  // of the kernel element; the value v_x right from ':' specifies the
  // value of the kernel element. Only one coordinate entry is scanned per line.
  // So big kernels with a few elements can be set easily. As long as lines with
  // kernel elements are found elements are set in the kernel table.
  // Return string is empty on successful scan. On errors it contains the number
  // of error lines.
  //
  // A second way to specify more kernel elements at once is with the
  // following string line format:
  //
  // \code
  //
  // (*, y_0, z_0, c_0, t_0, u_0):x_0, ... ,x_n
  // ...
  // (*, y_n, z_n, c_n, t_n, u_n):y_n, ... ,y_n
  //
  // \endcode
  //
  // So all kernel elements of a row are found in a string line. Note
  // that the x coordinates becomes invalid and must set as an asterisk.
  // Empty elements in the kernel can be left empty before between and
  // after commas.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    std::string TKernel<KDATATYPE>::setKernel(const std::string &kernString)
  {
    // Initialize an error string to return error line numbers.
    std::string errStr = "";
    char *copy = nullptr;
 
    // Make kernel empty.
    reset();
 
    // Scan only if string contains enough characters to specify one kernel element.
    if (kernString.length() < 15){
      // String too short to specify a valid kernel. Return 0 as error line.
      errStr = "0";
      return errStr;
    }
    else{
      // Get a modifiable copy of the string.
      std::string str = kernString;
 
      // Assure that a line feed terminates the string.
      if (str[str.length()-1] != '\n'){ str += '\n'; }
 
      // Scan line by line until no valid line is found.
      MLint lineNumber = 0;
      while (str != ""){
        // Scan line for following format : "(x,y,z,c,t,u):v"
        MLint x=0, y=0, z=0, c=0, t=0, u=0, num=0;
        MLdouble fValue = 0;
        num = sscanf(str.c_str(), "(%lld,%lld,%lld,%lld,%lld,%lld):%lf", &x, &y, &z, &c, &t, &u, &fValue);
        KDATATYPE value = static_cast<KDATATYPE>(fValue);
        if (num == 7){
          // Okay, we found a valid line. Add values to kernel coordinate table.
          _addCoordinate(ImageVector(x,y,z,c,t,u), value, false);
        }
        else{
          // Line is invalid. Test for the line format.
          num = sscanf(str.c_str(), "(*,%lld,%lld,%lld,%lld,%lld):", &y, &z, &c, &t, &u);
          if (5 == num){
            // Yes, the line format is ok. So search ":".
            const size_t cPos = str.find(':');
            if (cPos != std::string::npos){
              // Yes, found. First part of line has already been scanned and so removed it.
              str.erase(0, cPos+1);
 
              // Flag to stop line scan and counter for number of scanned elements.
              bool go = true;
              MLint xCoord = 0;
 
              // Scan until we reach at end of line or an error occurs.
              while (go && (str.length() != 0)){
                const char ch = str[0];
 
                // Kill space characters and continue.
                if (' ' == ch){
                  str.erase(0, 1);
                }
 
                // Line end? Then finish line scan.
                else if ('\n' == ch) {
                  // Terminate loop and leave '\n' because later it's searched.
                  go = false;
                }
 
                // Colon? Then increase element number count and kill colon.
                else if (',' == ch) {
                  // Go forward to next kernel element in row and remove comma.
                  ++xCoord;
                  str.erase(0, 1);
                }
                else{
                  // No space, no comma, so try to scan a number.
 
                  // First char where a number incompatible character is found.
                  char *stopPos=nullptr;
 
                  // Copy string.
                  copy = Memory::duplicateString(str.c_str(), ML_FATAL_MEMORY_ERROR);
                  if (copy){
                    // Scan string for number.
                    const MLdouble val = strtod(copy, &stopPos);
 
                    // More than 0 characters scanned successfully?
                    if (stopPos != copy){
                      // Yes, a number at beginning. Add it at the right position in the kernel.
                      _addCoordinate(ImageVector(xCoord, y,z,c,t,u), static_cast<KDATATYPE>(val), false);
 
                      // Remove number at start of string. Use pointer arithmetics to
                      // determine number of scanned number characters.
                      str.erase(0, stopPos-copy);
                    }
                    else{
                      // No parsing possible, there must be an error. So terminate scan of this line.
                      go = false;
 
                      // Add line number to error string.
                      char err[30]="";
                      snprintf(err, 29, " %lld", lineNumber);
                      errStr += err;
                    }
                  }
 
                  // Remove copy of line.
                  if (copy){
                    Memory::freeMemory(copy);
                    copy = nullptr;
                  }
                } // else
              } // while
            } // if (cPos != std::string::npos)
          } // if (num == 5)
        } // else
 
        // Last line?
        const size_t pos = str.find('\n');
        if (pos != std::string::npos){
          // Not the last line, because scanned linefeed is not at end of string.
          // So remove line from string.
          str.erase(0, pos+1);
        }
        else{
          // Last line. Make string empty to terminate loop.
          str = "";
        }
 
        // Increase number of scanned lines.
        lineNumber++;
      }  // while (str != "")
    } // else
 
    // Update the real kernel size.
    _calculateRealKernelExt();
 
    return errStr;
  }
 
 
 
  //-------------------------------------------------------------------------------------------
  // Defines a complete kernel matrix whose extents are defined by \c ext.
  // If desired the set of kernel values can be passed in \c values.
  // Otherwise all kernel values are set to 1.0/(KDATATYPE)ext.compMul(), i.e.
  // to the number of entries of the matrix.
  // If desired a set of mask values can be specified. If \c mask[n] is true then
  // tab[n] is included in the kernel; otherwise it's not part of the kernel.
  // Internally the specified kernel matrix is handled as a table of kernel elements.
  // Note that \c values and \c mask must have ext.compMul() elements if they are defined.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::setKernel(const ImageVector &ext, const KDATATYPE * const values, bool * const mask)
  {
    // Clear all dynamic tables.
    _clearTables();
 
    // Compute maximum size of value and kernel table.
    size_t maxElems = ext.compMul();
 
    // Do we have a mask ? Then count number of enabled kernel elements.
    _tabSize=0;
    if (mask){
      for (size_t c=0; c<maxElems; c++){ if (mask[c]){ _tabSize++; } }
    }
    else{
      // No mask set. All kernel elements are used.
      _tabSize = maxElems;
    }
 
    // Recreate the coordinate and value tables with the size of used elements.
    if (_tabSize > 0){
      _coordTab = new ImageVector[_tabSize];
      _valueTab = new KDATATYPE[_tabSize];
      if (!_coordTab || !_valueTab){
        _clearTables();
        ML_PRINT_FATAL_ERROR("TKernel<KDATATYPE>::setKernel(...):", ML_NO_MEMORY,
                             "Cannot create kernel filter tables, so tables are reset. This will probably cause other errors.");
      }
 
      if (_coordTab && _valueTab){
        // If we have a user defined value table then copy it.
        // Otherwise fill table with values 1.0/(KDATATYPE)_tabSize.
        if (values){
          memcpy(_valueTab, values, _tabSize*sizeof(KDATATYPE));
        }
        else{
          // As documented in class description we cast in case of integer KDATATYPES.
          KDATATYPE v = static_cast<KDATATYPE>(1.0/static_cast<KDATATYPE>(_tabSize));
          for (size_t d=0; d<_tabSize; d++){ _valueTab[d] = v; }
        }
 
        // Define all kernel coordinates if no mask exists or if the mask exists and the
        // entry for the current kernel voxel is true.
        size_t idx=0;
        for (size_t e=0; e<maxElems; e++){
          if (!mask || (mask && mask[e])){
            _coordTab[idx] = indexToCoord(e, ext);
            idx++;
          }
        }
 
        if (idx !=_tabSize){
          ML_PRINT_FATAL_ERROR("TKernel<KDATATYPE>::setKernel(...):", ML_CALCULATION_ERROR,
            "Internal error. Cannot create kernel filter tables. This will probably cause other errors.");
          _clearTables();
        }
      } // if (_coordTab && _valueTab)
    }
 
    // Calculate the real kernel size, i.e. the bounding box around all
    // used kernel elements.
    _calculateRealKernelExt();
  }
 
 
  //-------------------------------------------------------------------------------------------
  // Create kernel coordinate and value tables in separable table format, that means
  // a 2-D kernel where each rows describes the elements of 1-D kernels.
  // At most 6 rows make sense, because the maximum kernel dimension is 6.
  // It is legal to leave out any axis by passing empty vectors to define lower
  // dimensional separable kernels.
  // The separability flag is set to true.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::setSeparableKernel(const std::vector<typename std::vector<KDATATYPE> > &separableRows)
  {
    setSeparable(true);
    _clearTables();
 
    // Get number of rows, permit at most 6 for 6 dimensions.
    const size_t numRows = separableRows.size() > 6 ? 6 : separableRows.size();
 
    // Sum up total number of entries in all rows.
    size_t numEntries = 0;
    for (size_t e=0; e < numRows; ++e){ numEntries += separableRows[e].size(); }
 
    // Create a value table and a coordinate table.
    std::vector<KDATATYPE> valueTab; valueTab.reserve(numEntries);
    std::vector<ImageVector>    coordTab; coordTab.reserve(numEntries);
 
    // Compose coordinate table.
    for (size_t r=0; r < numRows; ++r){
 
      // Get row and its number of elements.
      const std::vector<KDATATYPE> &row = separableRows[r];
      const size_t numRElements = row.size();
 
      // Coordinate in the separable table, only x and y coordinates are set.
      // y defines the dimension, x defines the position in that dimension.
      ImageVector entryCoord(0);
      entryCoord.y = r;
 
      // Count entries with row.
      for (size_t re=0; re < numRElements; ++re){
        // Define x-coordinate of entry in table row.
        entryCoord.x = re;
        coordTab.push_back(entryCoord);
        valueTab.push_back(row[re]);
      }
    }
 
    // Define the kernel table. The extent is 2-D, because it's the
    // separable table format.
    setPartialKernel(coordTab.size(), &(coordTab[0]), &(valueTab[0]));
 
    // Make separable information up to date.
    _validateSeparabilityInfos();
  }
 
 
  //-------------------------------------------------------------------------------------------
  // Resizes the kernel to a new state and tries to maintain the previous elements
  // if possible. Note that new kernel size may differ from desired one if empty
  // areas are at border of the resized kernel so that only a smaller real kernel remains.
  // Regions which are created newly get kernel elements with value \c newVal.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::resizeKernel(const ImageVector &ext, KDATATYPE newVal)
  {
    // Create a buffer kernel.
    TKernel<KDATATYPE> bufferKern = *this;
 
    // Clear all dynamic tables.
    _clearTables();
 
    // Scan new kernel extents.
    for (MLint u=0; u < ext.u; u++){
      for (MLint t=0; t < ext.t; t++){
        for (MLint c=0; c < ext.c; c++){
          for (MLint z=0; z < ext.z; z++){
            for (MLint y=0; y < ext.y; y++){
              for (MLint x=0; x < ext.x; x++){
                // Create vector from new coordinates.
                ImageVector addCoord(x,y,z,c,t,u);
 
                // If new coordinate is within old kernel range the add it only if it exists.
                if ((u < getExtent().u) &&
                    (t < getExtent().t) &&
                    (c < getExtent().c) &&
                    (z < getExtent().z) &&
                    (y < getExtent().y) &&
                    (x < getExtent().x)){
                  // Within old extents. Is it defined?
                  MLint idx = bufferKern.findIndex(addCoord);
                  if (idx >=0){
                    // Old coordinate existed. So add it to new kernel without extent updating.
                    _addCoordinate(addCoord, bufferKern.getValueTab()[idx], false);
                  }
                }
                else{
                  // New coordinate is outside old range. Add newVal's in new area.
                  _addCoordinate(addCoord, newVal, false);
                } // else
 
              } // x
            }
          }
        }
      }
    } // u
 
    // Calculate the real kernel size, i.e. the bounding box around all
    // used kernel elements.
    _calculateRealKernelExt();
  }
 
 
  //-------------------------------------------------------------------------------------------
  // Fill all undefined kernel elements with \c newVal.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::fillUndefinedElements(KDATATYPE newVal)
  {
    // Invalidate current information for separable kernels.
    _invalidateSeparableInfos();
 
    // Scan new kernel extents.
    for (MLint u=0; u < getExtent().u; u++){
      for (MLint t=0; t < getExtent().t; t++){
        for (MLint c=0; c < getExtent().c; c++){
          for (MLint z=0; z < getExtent().z; z++){
            for (MLint y=0; y < getExtent().y; y++){
              for (MLint x=0; x < getExtent().x; x++){
                // Create vector from new coordinates.
                ImageVector addCoord(x,y,z,c,t,u);
 
                // Does the coordinate exist in kernel?
                if (findIndex(addCoord) < 0){
                  // Coordinate does not exits. So add it without extent updating.
                  _addCoordinate(addCoord, newVal, false);
                } // if
 
              } // x
            }
          }
        }
      }
    } // u
 
    // Calculate the real kernel size, i.e. the bounding box around all
    // used kernel elements.
    _calculateRealKernelExt();
  }
 
 
 
  //-------------------------------------------------------------------------------------------
  // Takes the current kernel, computes radii from the extents of the kernel
  // and removes all kernel elements which are outside the ellipsoid defined by
  // these radii.
  // Note that filter properties of kernel changes by this operation.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::makeCircular()
  {
    // Invalidate current information for separable kernels.
    _invalidateSeparableInfos();
 
    // Compute theoretical center of kernel. Subtract 0.5 to shift
    // kernel center from voxel center to voxel origin. Examples in 3D:
    // E.g. : (1,1,1) sized kernels are addressed from 0 to 0 and so they have
    //         their center at (0,0,0)\n
    // E.g. : (2,2,2) sized kernels are addressed from 0 to 1 and so they have
    //         their center at (0.5, 0.5, 0.5).\n
    // E.g. : (3,3,3) sized kernels are addressed from 0 to 2 and so they have
    //         their center at (1.0, 1.0, 1.0).\n
    MLdouble cx = static_cast<MLdouble>(_ext.x) / 2.0 - 0.5;
    MLdouble cy = static_cast<MLdouble>(_ext.y) / 2.0 - 0.5;
    MLdouble cz = static_cast<MLdouble>(_ext.z) / 2.0 - 0.5;
    MLdouble cc = static_cast<MLdouble>(_ext.c) / 2.0 - 0.5;
    MLdouble ct = static_cast<MLdouble>(_ext.t) / 2.0 - 0.5;
    MLdouble cu = static_cast<MLdouble>(_ext.u) / 2.0 - 0.5;
 
    // Compute kernel radii.
    MLdouble rx = static_cast<MLdouble>(_ext.x) / 2.0;
    MLdouble ry = static_cast<MLdouble>(_ext.y) / 2.0;
    MLdouble rz = static_cast<MLdouble>(_ext.z) / 2.0;
    MLdouble rc = static_cast<MLdouble>(_ext.c) / 2.0;
    MLdouble rt = static_cast<MLdouble>(_ext.t) / 2.0;
    MLdouble ru = static_cast<MLdouble>(_ext.u) / 2.0;
 
    // Save local tables and set default table states.
    ImageVector    *coordTab  = _coordTab;
    _coordTab = nullptr;
 
    KDATATYPE *valueTab  = _valueTab;
    _valueTab = nullptr;
 
    size_t tabSize = _tabSize;
    _tabSize = 0;
 
    // Compute distance for all kernel voxels and normalize them.
    // If the coordinate is within ellipsoid distance then
    // add it to the kernel elements, otherwise forget it.
    for (size_t e=0; e < tabSize; e++){
      ImageVector v = coordTab[e];
      MLdouble pu = (static_cast<MLdouble>(v.u) - cu)/ru;
      MLdouble pt = (static_cast<MLdouble>(v.t) - ct)/rt;
      MLdouble pc = (static_cast<MLdouble>(v.c) - cc)/rc;
      MLdouble pz = (static_cast<MLdouble>(v.z) - cz)/rz;
      MLdouble py = (static_cast<MLdouble>(v.y) - cy)/ry;
      MLdouble px = (static_cast<MLdouble>(v.x) - cx)/rx;
 
      // Distance smaller or equal 1? Then we are within kernel radius and
      // we set kernel value 'true'; otherwise we set kernel value 'false'.
      if (sqrt(px*px + py*py + pz*pz + pc*pc + pt*pt + pu*pu) <= 1){
        _addCoordinate(v, valueTab[e], false);
      }
    }
 
    // Remove saved tables.
    delete [] coordTab;
    coordTab = nullptr;
    delete [] valueTab;
    valueTab = nullptr;
 
    // Recalculate the kernelSize.
    _calculateRealKernelExt();
  }
 
 
 
  //-------------------------------------------------------------------------------------------
  // Applies a point symmetric mirroring of all kernel elements.
  // The dimension param specifies a mirroring dimension if in range [0,5],
  // otherwise mirroring is applied in all dimensions. Default is -1.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::mirror(int dimension)
  {
    // Invalidate current information for separable kernels.
    _invalidateSeparableInfos();
 
    // Save local tables and set default table states.
    ImageVector    *coordTabBuf  = _coordTab;
    _coordTab = nullptr;
 
    KDATATYPE *valueTabBuf  = _valueTab;
    _valueTab = nullptr;
 
    size_t tabSizeBuf       = _tabSize;
    _tabSize = 0;
 
    // Calculate maximum kernel positions in all dimensions, that is one
    // less than extent.
    ImageVector maxP(_ext-ImageVector(1));
 
    if ((dimension<0) || (dimension>5)){
      // Subtract each coordinate of kernel elements from the
      // maximum extent position in all 6 dimensions and
      // add it to the kernel elements (which have been made
      // empty above).
      for (size_t e=0; e < tabSizeBuf; e++){
        _addCoordinate(maxP - coordTabBuf[e], valueTabBuf[e], false);
      }
    }
    else{
      ImageVector mirVec;
      const size_t dimensionSize_t = mlrange_cast<size_t>(dimension);
      for (size_t e=0; e < tabSizeBuf; e++){
        // Take normal coordinate.
        mirVec = coordTabBuf[e];
 
        // mirror the selected coordinate component.
        mirVec[dimensionSize_t] = maxP[dimensionSize_t] - coordTabBuf[e][dimensionSize_t];
 
        // Add new coordinate with value to kernel.
        _addCoordinate(mirVec, valueTabBuf[e], false);
      }
    }
 
    // Remove old saved tables.
    delete [] coordTabBuf;
    coordTabBuf = nullptr;
    delete [] valueTabBuf;
    valueTabBuf = nullptr;
 
    // Recalculate the kernelSize.
    _calculateRealKernelExt();
  }
 
  //-------------------------------------------------------------------------------
  // Calculate binomial coefficients for (n)
  //                                     (k)
  //-------------------------------------------------------------------------------
  template <typename KDATATYPE>
    MLldouble TKernel<KDATATYPE>::binomialcoeff(MLint n, MLint k)
  {
    MLint     i=0;
    MLldouble v=1.0;
 
    if (k>n){ v=0; }
    else{
      if((k!=0) && (k!=n)){
        for (i=1; i<(n+1);   i++){ v*=i; }
        for (i=1; i<(k+1);   i++){ v/=i; }
        for (i=1; i<(n-k+1); i++){ v/=i; }
      }
    }
 
    return(v);
  }
 
  //-------------------------------------------------------------------------------------------
  // Returns a vector with \c numSample values binomial coefficients.
  // If \c numSamples is passed as 0 then an empty vector is returned.
  // if \c normalize is passed true (the default) then all values are
  // normalized to the sum of absolute values 1.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    std::vector<KDATATYPE> TKernel<KDATATYPE>::get1DGauss(size_t numSamples, bool normalize)
  {
    // The return value;
    std::vector<KDATATYPE> v;
 
    if (numSamples > 0){
      v.reserve(numSamples);
 
      // Fill the axis with binomial coefficients and sum up all values.
      KDATATYPE sum = 0;
      for (size_t u=0; u<numSamples; ++u){
        // Calculate the entry, write it into the array and sum its absolute value up.
        // As documented in class description we cast in case of integer KDATATYPES.
        KDATATYPE val = static_cast<KDATATYPE>(binomialcoeff(numSamples-1, u));
        v.push_back(val);
        sum += static_cast<KDATATYPE>(mlAbs(val));
      }
 
      // If normalization is desired then divide all elements by the sum.
      if (normalize && (sum > 0)){
        for (size_t u=0; u<numSamples; u++){ v[u] /= sum; }
      }
    }
 
    // Return the vector.
    return v;
  }
 
  //-------------------------------------------------------------------------------------------
  // Replaces the current kernel by a normalized gauss kernel.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::setGauss(const ImageVector &ext)
  {
    // Invalidate current information for separable kernels.
    _invalidateSeparableInfos();
 
    // Build gauss kernel by using separability property of gauss kernel.
    // So calculate a discrete binomial function for each axis.
    // The values of the kernel elements are products of the corresponding
    // values along the six binomial function axis.
    std::vector<KDATATYPE> vx=get1DGauss(ext.x, false);
    std::vector<KDATATYPE> vy=get1DGauss(ext.y, false);
    std::vector<KDATATYPE> vz=get1DGauss(ext.z, false);
    std::vector<KDATATYPE> vc=get1DGauss(ext.c, false);
    std::vector<KDATATYPE> vt=get1DGauss(ext.t, false);
    std::vector<KDATATYPE> vu=get1DGauss(ext.u, false);
 
    // Build kernel by using the setKernel method for separable kernels.
    // Normalize kernel values after building the product to remain in normal grey value areas.
    setKernel(ext, &(vx[0]), &(vy[0]), &(vz[0]), &(vc[0]), &(vt[0]), &(vu[0]), true);
  }
 
 
  //-------------------------------------------------------------------------------------------
  // Initialization. Should be called only by constructors.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::_init()
  {
    _ext.set(0);
    _negExt.set(0);
    _posExt.set(0);
 
    _tabSize  = 0;
 
    _valueTab = nullptr;
    _coordTab = nullptr;
 
    // Separability information.
    _isSeparable                  = false;
    _areSeparabilityInfosUpToDate = false;
    _separableDimEntries.set(0);
    _separableOneDimExt.set(0);
    _separableExt.set(0);
    _separableNegExt.set(0);
    _separablePosExt.set(0);
    _separableCoordTab.clear();
  }
 
 
  //-------------------------------------------------------------------------------------------
  // Clear all internal (dynamic) tables.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::_clearTables()
  {
    // Set table size to 0.
    _tabSize = 0;
 
    // Invalidate separability infos which are calculated on demand.
    _areSeparabilityInfosUpToDate = false;
 
    delete [] _valueTab;
 
    _valueTab = nullptr;
    delete [] _coordTab;
    _coordTab = nullptr;
  }
 
 
 
  //-------------------------------------------------------------------------------------------
  // Add a new coordinate with value to the coordinate and value table.
  // Still not tested! If update is true (default) then the kernel extents
  // are updated.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::_addCoordinate(const ImageVector &pos, KDATATYPE value, bool update)
  {
    // Create new data tables.
    ImageVector    *newCoordTab = nullptr;
    KDATATYPE *newValueTab = nullptr;
 
    // Invalidate separability infos which are calculated on demand.
    _areSeparabilityInfosUpToDate = false;
 
    // Compute new table size.
    size_t newTabSize = _tabSize+1;
 
    newCoordTab = new ImageVector [_tabSize+1];
    newValueTab = new KDATATYPE[_tabSize+1];
 
    // Handle not enough memory.
    if (!newCoordTab || !newValueTab){
      delete [] newCoordTab;
      newCoordTab = nullptr;
      delete [] newValueTab;
      newValueTab = nullptr;
      ML_PRINT_FATAL_ERROR("Kernel<KDATATYPE>::_addCoordinate(...):",
                           ML_NO_MEMORY,
                           "Cannot create kernel tables. Returning with invalid "
                           "reset tables. This will probably cause other errors.");
      return;
    }
 
    // Copy old contents into new tables.
    if (_tabSize > 0){
#if defined(__GNUC__)
      // Suppress warning:
      // "'void* memcpy(void*, const void*, size_t)' writing to an object of type 'ml::ImageVector'
      // {aka 'class ml::TImageVector<long long int>'} with no trivial copy-assignment;
      // use copy-assignment or copy-initialization instead [-Wclass-memaccess]"
      // We are sure this is ok in this case!
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wclass-memaccess"
#endif
      memcpy(newCoordTab, _coordTab, _tabSize * sizeof(ImageVector));
#if defined(__GNUC__)
#pragma GCC diagnostic pop
#endif
      memcpy(newValueTab, _valueTab, _tabSize * sizeof(KDATATYPE));
    }
 
    // Set new entry.
    newCoordTab[_tabSize] = pos;
    newValueTab[_tabSize] = value;
 
    // Remove old tables.
    _clearTables();
 
    // Increase table size.
    _tabSize  = newTabSize;
    _coordTab = newCoordTab;
    _valueTab = newValueTab;
 
    // Update the real kernel size.
    if (update){ _calculateRealKernelExt(); }
  }
 
 
  //-------------------------------------------------------------------------------------------
  // Calculate the correct maximum kernel extent for all used kernel elements.
  // The maximum extent in all dimensions is reduced as much that all kernel coordinates just
  // remain in valid extent. Leading empty entries in the kernel are not removed.
  // Empty kernels are considered as kernels with extent (1,1,1,1,1,1).
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::_calculateRealKernelExt()
  {
    // Search minimum and maximum extent of the kernel.
    // For the case that no coordinate elements are available set
    // values to (0,0,0,0,0,0) so that resulting kernel ext is
    // calculated to (1,1,1,1,1,1) below.
    ImageVector maximum(0,0,0,0,0,0);
    if (getTabSize() > 0){
      for (size_t c=0; c < getTabSize(); c++){
        maximum = ImageVector::compMax(_coordTab[c], maximum);
      }
    }
 
    // Define kernel ext as difference between corners of bounding box.
    _ext = maximum+ImageVector(1,1,1,1,1,1);
 
    // Compute a negative and positive extend of the kernel. So even and odd kernels
    // can be handled without really having a kernel center. Subtract one to avoid a
    // to avoid increment of extents by one.
    _negExt = calculateNegativeExtentFromExtent(_ext);
    _posExt = calculatePositiveExtentFromExtent(_ext);
 
    // The kernel must never have extents smaller than 1 after this call.
    if (!_ext.allBiggerZero()){
      ML_PRINT_FATAL_ERROR("void TKernel<KDATATYPE>::_calculateRealKernelExt()",
                           ML_PROGRAMMING_ERROR,
                           "Kernel has invalid extents now. Continuing anyway.");
    }
  }
 
 
 
  //-------------------------------------------------------------------------------------------
  //
  //        Support for an interpretation as separable kernel.
  //
  //-------------------------------------------------------------------------------------------
 
  //-------------------------------------------------------------------------------------------
  // Set/unset a flag to indicate that the first 6 rows of the first kernel slice
  // is considered as 1-D axes of a separable filter kernel. Note that this flag
  // only denotes how applications shall interpret this kernel; however it does
  // NOT change any behaviour of this class.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::setSeparable(bool isSeparableVal)
  {
    // Invalidate separability information if flag changes.
    if (_isSeparable != isSeparableVal){
      _isSeparable = isSeparableVal;
      _invalidateSeparableInfos();
    }
  }
 
  //-------------------------------------------------------------------------------------------
  // Indicates whether the first 6 rows of the first kernel slice are
  // interpreted as 1-D axes of a separable filter kernel; default is false.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    bool TKernel<KDATATYPE>::isSeparable() const
  {
    return _isSeparable;
  }
 
  //-------------------------------------------------------------------------------------------
  // Returns a ImageVector with the number of entries of the separable kernel for
  // the dimensions 0,...,5. That means the n-th index contains the number
  // of valueTab entries of row n, where n is from [0,...,5].
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    ImageVector TKernel<KDATATYPE>::getSeparableDimEntries() const
  {
    _validateSeparabilityInfos();
    return _separableDimEntries;
  }
 
  //------------------------------------------------------------------------------------------
  //------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    ImageVector TKernel<KDATATYPE>::getSeparableOneDimExtents() const
  {
    _validateSeparabilityInfos();
    return _separableOneDimExt;
  }
 
  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    size_t TKernel<KDATATYPE>::getSeparableDimIndex(size_t dim) const
  {
    _validateSeparabilityInfos();
    if (dim>5){ dim = 5; }
 
    // Sum up elements in lower dimensions to get index to the correct table entry.
    switch (dim){
      case 0: return 0;
 
      case 1: return static_cast<size_t>(_separableDimEntries.x);
 
      case 2: return static_cast<size_t>(_separableDimEntries.x +
                                         _separableDimEntries.y);
 
      case 3: return static_cast<size_t>(_separableDimEntries.x +
                                         _separableDimEntries.y +
                                         _separableDimEntries.z);
 
      case 4: return static_cast<size_t>(_separableDimEntries.x +
                                         _separableDimEntries.y +
                                         _separableDimEntries.z +
                                         _separableDimEntries.c);
 
      case 5: return static_cast<size_t>(_separableDimEntries.x +
                                         _separableDimEntries.y +
                                         _separableDimEntries.z +
                                         _separableDimEntries.c +
                                         _separableDimEntries.t);
      default:{
        ML_PRINT_FATAL_ERROR("TKernel<KDATATYPE>::getSeparableDimIndex",
                             ML_PROGRAMMING_ERROR, "returning x-table entries");
        return 0;
      }
    }
  }
 
  //-------------------------------------------------------------------------------------------
  // Returns the table of all coordinates for filtering with separable kernels.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    const std::vector<ImageVector> &TKernel<KDATATYPE>::getSeparableCoordTab() const
  {
    _validateSeparabilityInfos();
    return _separableCoordTab;
  }
 
  //-------------------------------------------------------------------------------------------
  // Invalidate separability information.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::_invalidateSeparableInfos()
  {
    _areSeparabilityInfosUpToDate = false;
  }
 
  //-------------------------------------------------------------------------------------------
  // This method updates all members which are needed in query methods for
  // separability properties of the kernel.
  // Note:
  // This method changes mutable members since it is required for get-methods
  // on information about separability information.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::_validateSeparabilityInfos() const
  {
    // Speed up performance since it's called first in all separable info methods.
    if (_areSeparabilityInfosUpToDate){ return; }
    if (!_areSeparabilityInfosUpToDate){
 
      // Scan all coordinate components of the kernel until the y-coordinate is >= 6,
      // because all other information after that is not relevant for first 6 separable
      // filter rows. Expand the extent of the separable kernel and determine coordinate
      // of the axis components.
      _separableDimEntries.set(0);
      _separableOneDimExt.set(0);
      _separableExt.set(0);
      _separableCoordTab.clear();
      ImageVector sepCoord(-1);
      size_t idx=0;
      while ((idx<_tabSize) && (_coordTab[idx].y < 6)){
        // Get the current entry.
        const ImageVector &entry = _coordTab[idx];
 
        const size_t entryY = mlrange_cast<size_t>(entry.y);
 
        // Increment the number of entries of a certain dimension of the
        // 6 1-D kernels if the kernel is interpreted as a separable kernel.
        // Note that each of the first six rows contains the separable 1-D kernel
        // for dimensions 1-6.
        _separableDimEntries[entryY]++;
 
        // Increment the extent of the separable kernel dimensions if the
        // coordinate given by entry has a larger x-coordinate than the
        // extent of the corresponding dimension.
        if (entry.x+1 > _separableOneDimExt[entryY]){ _separableOneDimExt[entryY] = entry.x+1; }
 
        // Add the coordinate of the 1-D kernel to the table of coordinates.
        // Note that separable kernels are placed in the "center" of the corresponding
        // "normal" kernel. Hence, all components which are not given by entry.y,
        // need to be set to the position of the kernel "center" which is given by
        // the negative extent of the kernel. However, we do not have this "center"
        // here, because we do not know it before having all extents.
        // So we we mark the missing coordinate components with -1 and replace them
        // in an additional loop afterwards.
        sepCoord.set(-1);
        sepCoord[entryY] = entry.x;
        _separableCoordTab.push_back(sepCoord);
 
        // Count the number of separable entries in the kernel with idx: Go to next entry.
        ++idx;
      }
 
      // Determine second version of separable kernel extent where non filtered
      // dimensions are set to 1.
      _separableExt = _separableOneDimExt;
      _separableExt.fillSmallerComps(1,1);
 
      // Determine the center of the kernel which is considered to be the negative kernel extent.
      const ImageVector negExt =_separableExt / ImageVector(2);
 
      // Replace all -1 component values with the kernel center; thus the 1-D kernel is
      // centered in the region where the corresponding normal kernel would also have its
      // "center".
      for (size_t e=0; e < idx; ++e){
        // Get coordinate.
        ImageVector &coord = _separableCoordTab[e];
        for (size_t d=0; d < 6; ++d){
          // Replace all -1 component values with the kernel center.
          if (coord[d] == -1){ coord[d] = negExt[d]; }
        }
      }
 
      // Calculate negative and positive kernel extent for a separable kernel.
      // If a dimension is empty or zero sized then a zero component is returned.
      _separableNegExt = _separableExt/ImageVector(2);
      _separablePosExt = (_separableExt - (_separableExt/ImageVector(2))) - ImageVector(1);
 
      // Mark infos as valid.
      _areSeparabilityInfosUpToDate = true;
    }
  }
 
  //-------------------------------------------------------------------------------------------
  // End of support for the interpretation of the kernel as separable kernel.
  //-------------------------------------------------------------------------------------------
 
 
 
  typedef TKernel<MLfloat>        FloatKernel;
 
  typedef TKernel<MLdouble>       DoubleKernel;
 
  typedef TKernel<MLldouble>      LDoubleKernel;
 
 
  typedef MLdouble                KernelDataType;
 
  typedef TKernel<KernelDataType> Kernel;
 
 
ML_END_NAMESPACE
 
#endif // __mlKernel_H