MeVisLab/Standard/Sources/ML/MLKernel/mlKernel.h

Go to the documentation of this file.

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

//-------------------------------------------------------------------------
// Prevent multiple including of this file.
#if !defined(__mlKernel_H)
#define __mlKernel_H

// ML-includes
#ifndef __mlInitSystemKernel_H
#include "mlInitSystemKernel.h"
#endif

#include <errno.h>

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=NULL);

    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=NULL, bool * const mask=NULL);

    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)
    {
      ML_TRACE_IN( "TKernel<KDATATYPE::setKernel( )" );
      ML_TRY
      {
        // 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=NULL;

        // If desired then normalize all kernel elements so that their sum becomes 1.
        if (normalize){ manipulateKernelElements(KERN_NORMALIZE, 0); }
      }
      ML_CATCH_RETHROW;
    }

    void setSeparableKernel(const std::vector<ML_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;


#ifdef ML_DEPRECATED



  public:

    inline ML_DEPRECATED KDATATYPE getNegValueSum() const { return getNegativeValueSum(); }

    inline ML_DEPRECATED KDATATYPE getPosValueSum() const { return getPositiveValueSum(); }

    inline ML_DEPRECATED const ImageVector &getExt() const { return getExtent(); }

    inline ML_DEPRECATED const ImageVector &getNegExt() const { return getNegativeExtent(); }

    inline ML_DEPRECATED const ImageVector &getPosExt() const { return getPositiveExtent(); }

    static inline ML_DEPRECATED ImageVector calcNegExtFromExt(const ImageVector &ext)  { return calculateNegativeExtentFromExtent(ext); }

    static inline ML_DEPRECATED ImageVector calcPosExtFromExt(const ImageVector &ext)  { return calculatePositiveExtentFromExtent(ext); }


#endif // ML_DEPRECATED


  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()
  {
    ML_TRACE_IN( "TKernel::TKernel() " );

    _init();
  }

  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    TKernel<KDATATYPE>::TKernel(const TKernel &kern)
  {
    ML_TRACE_IN( "TKernel::TKernel()" );

    _init();
    operator=(kern);
  }

  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    TKernel<KDATATYPE>::~TKernel()
  {
    ML_TRACE_IN( "TKernel::~TKernel( )" );
    ML_TRY
    {
      _clearTables();
    }
    ML_CATCH; // Do not rethrow in destructors.
  }

  //-------------------------------------------------------------------------------------------
  // Assignment operator. Sets current instance to contents of \c kern.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    const TKernel<KDATATYPE>& TKernel<KDATATYPE>::operator=(const TKernel &kern)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::operator=" );
    ML_TRY
    {
      // 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;
    }
    ML_CATCH_RETHROW;
  }

  //-------------------------------------------------------------------------------------------
  // Reset kernel to construction state.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::reset()
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::reset()" );
    ML_TRY
    {
      // Remove dynamic structures.
      _clearTables();

      // Call constructor initialization. We can do that here, because
      // all internal tables have been cleared.
      _init();
    }
    ML_CATCH_RETHROW;
  }

  //-------------------------------------------------------------------------------------------
  // Returns the current number of kernel elements.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    size_t TKernel<KDATATYPE>::getTabSize() const
  {
    ML_TRACE_IN_TIME_CRITICAL("size_t TKernel<KDATATYPE>::getTabSize()");

    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
  {
    ML_TRACE_IN_TIME_CRITICAL("const ImageVector* TKernel<KDATATYPE>::getCoordTab()");

    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
  {
    ML_TRACE_IN_TIME_CRITICAL("const KDATATYPE* TKernel<KDATATYPE>::getValueTab()");

    return _isSeparable ? (_valueTab + getSeparableDimIndex(dim)) : _valueTab;
  }

  //-------------------------------------------------------------------------------------------
  // Return the sum of all kernel element values
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    KDATATYPE TKernel<KDATATYPE>::getValueTabSum() const
  {
    ML_TRACE_IN( "KDATATYPE TKernel<( )::getValueTabSum()" );

    // 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
  {
    ML_TRACE_IN( "KDATATYPE TKernel<KDATATYPE>::getMinValue()" );

    // 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
  {
    ML_TRACE_IN( "KDATATYPE TKernel<KDATATYPE>::getMaxValue() " );

    // 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
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::getNegativeValueSum()" );

    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
  {
    ML_TRACE_IN( "KDATATYPE TKernel<KDATATYPE>::getPositiveValueSum()" );

    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)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::manipulateKernelElements()" );
    ML_TRY
    {
      // 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(_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, 1.0))){
          MLldouble invBase=1.0/log(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(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;

      }
    }
    ML_CATCH_RETHROW;
  }

  //-------------------------------------------------------------------------------------------
  // 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
  {
    ML_TRACE_IN_TIME_CRITICAL("const ImageVector &TKernel<KDATATYPE>::getExtent()");

    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
  {
    ML_TRACE_IN_TIME_CRITICAL("const ImageVector &TKernel<KDATATYPE>::getNegativeExtent()");

    if (_isSeparable){
      _validateSeparabilityInfos();
      return _separableNegExt;
    }
    else{
      return _negExt;
    }
  }

  //-------------------------------------------------------------------------------------------
  // See \c getNegativeExtent().
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    const ImageVector &TKernel<KDATATYPE>::getPositiveExtent() const
  {
    ML_TRACE_IN_TIME_CRITICAL("const ImageVector &TKernel<KDATATYPE>::getPositiveExtent()");

    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)
  {
    ML_TRACE_IN( "MLint TKernel<KDATATYPE>::coordToIndex()" );

    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)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::coordToIndex()" );

    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)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::indexToCoord()" );

    return SubImage::indexToCoord(idx, ext);
  }

  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    ImageVector TKernel<KDATATYPE>::calculateNegativeExtentFromExtent(const ImageVector &ext)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::calculateNegativeExtentFromExtent()" );
    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)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::calculatePositiveExtentFromExtent()" );
    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
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::findIndex()" );

    // 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)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::setPartialKernel()" );
    ML_TRY
    {
      // Clear all dynamic tables.
      _clearTables();

      // Define the new number of kernel voxels.
      _tabSize = numElems;

      // Recreate the coordinate and value tables
      if (_tabSize > 0)      {
        ML_CHECK_NEW(_coordTab, ImageVector[_tabSize]);

        ML_CHECK_NEW(_valueTab, 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();
    }
    ML_CATCH_RETHROW;
  }


  //-------------------------------------------------------------------------------------------
  // 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
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::getKernel()" );
    ML_TRY
    {
      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.
                  MLsnprintf(strLine, 255, "(*,%I64Ld,%I64Ld,%I64Ld,%I64Ld,%I64Ld):", 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 long double because string reading and parsing of carrier
                      // types (which is theoretically possible) 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;
    }
    ML_CATCH_RETHROW;
  }



  //-------------------------------------------------------------------------------------------
  // 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)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::setKerne()" );

    // Initialize an error string to return error line numbers.
    std::string errStr = "";
    char *copy = NULL;

    ML_TRY
    {
      // 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 = MLsscanf(str.c_str(), "(%I64Ld,%I64Ld,%I64Ld,%I64Ld,%I64Ld,%I64Ld):%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 = MLsscanf(str.c_str(), "(*,%I64Ld,%I64Ld,%I64Ld,%I64Ld,%I64Ld):", &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=NULL;

                    // 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]="";
                        MLsnprintf(err, 29, " %I64Ld", lineNumber);
                        errStr += err;
                      }
                    }

                    // Remove copy of line.
                    if (copy){
                      Memory::freeMemory(copy);
                      copy = NULL;
                    }
                  } // 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();
    }
    ML_CATCH_BLOCK(...)
    {
      Memory::freeMemory(copy);
      copy = NULL;
      ML_CHECK(0); // Log the error.
      throw;
    }

    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)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::setKernel()" );
    ML_TRY
    {
      // 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){
        ML_CHECK_NEW(_coordTab, ImageVector[_tabSize]);
        ML_CHECK_NEW(_valueTab, 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();
    }
    ML_CATCH_RETHROW;
  }


  //-------------------------------------------------------------------------------------------
  // 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<ML_TYPENAME std::vector<KDATATYPE> > &separableRows)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE::setSeparableKernel( )" );
    ML_TRY
    {
      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();
    }
    ML_CATCH_RETHROW;
  }


  //-------------------------------------------------------------------------------------------
  // 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)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::resizeKernel()" );
    ML_TRY
    {
      // 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();
    }
    ML_CATCH_RETHROW;
  }


  //-------------------------------------------------------------------------------------------
  // Fill all undefined kernel elements with \c newVal.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::fillUndefinedElements(KDATATYPE newVal)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::fillUndefinedElements()" );
    ML_TRY
    {
      // 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();
    }
    ML_CATCH_RETHROW;
  }



  //-------------------------------------------------------------------------------------------
  // 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()
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::makeCircular()" );
    ML_TRY
    {
      // 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 = NULL;

      KDATATYPE *valueTab  = _valueTab;
      _valueTab = NULL;

      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 (sqrtf(px*px + py*py + pz*pz + pc*pc + pt*pt + pu*pu) <= 1){
          _addCoordinate(v, valueTab[e], false);
        }
      }

      // Remove saved tables.
      ML_DELETE_ARRAY(coordTab);
      ML_DELETE_ARRAY(valueTab);

      // Recalculate the kernelSize.
      _calculateRealKernelExt();
    }
    ML_CATCH_RETHROW;
  }



  //-------------------------------------------------------------------------------------------
  // 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)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::mirror()" );
    ML_TRY
    {
      // Invalidate current information for separable kernels.
      _invalidateSeparableInfos();

      // Save local tables and set default table states.
      ImageVector    *coordTabBuf  = _coordTab;
      _coordTab = NULL;

      KDATATYPE *valueTabBuf  = _valueTab;
      _valueTab = NULL;

      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;
        for (size_t e=0; e < tabSizeBuf; e++){
          // Take normal coordinate.
          mirVec = coordTabBuf[e];

          // mirror the selected coordinate component.
          mirVec[dimension] = maxP[dimension] - coordTabBuf[e][dimension];

          // Add new coordinate with value to kernel.
          _addCoordinate(mirVec, valueTabBuf[e], false);
        }
      }

      // Remove old saved tables.
      ML_DELETE_ARRAY(coordTabBuf);
      ML_DELETE_ARRAY(valueTabBuf);

      // Recalculate the kernelSize.
      _calculateRealKernelExt();
    }
    ML_CATCH_RETHROW;
  }




  // TODO !!! investigate !!!: (Guenter Ott - 08.04.2004) OK:
  //  Maybe this method should be placed at a more central point?
  //  (the template is also not used)
  //  There is easy a danger of overflow.
  //-------------------------------------------------------------------------------
  // Calculate binomial coefficients for (n)
  //                                     (k)
  //-------------------------------------------------------------------------------
  template <typename KDATATYPE>
    MLldouble TKernel<KDATATYPE>::binomialcoeff(MLint n, MLint k)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::binomialcoeff()" );
    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);
  }



  // TODO !!! investigate !!!: (Guenter Ott - 08.04.2004) OK : an incomplete suggestion
  //  to calculate binomialcoeff.
  //  Maybe interesting. If not just delete.
/*
  template <typename KDATATYPE>
  MLldouble TKernel<KDATATYPE>::binomialcoeff(MLint n, MLint k) const
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::binomialcoeff()" );
    MLint     i=0;
    MLldouble v=0.0;


    // bei der Bedingung bin ich mir nicht klaren; wuerde rueckfragen, was mit n<0 oder mit k<0 ist.
    // auch das richtige resultat fuer n=0 und k=0 ist mir nicht bekannt (tippe aber auf 1)
    if( k<=n )
    {
      v=1.0;
      MLint kGrenze = n-k > k ? n-k : k;
      for ( MLint i=n; i>kGrenze ;--i)
      {
        v*=i;
        v/=(n+1-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)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::get1DGauss()" );

    // The return value;
    std::vector<KDATATYPE> v;

    ML_TRY
    {
      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>(fabs(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; }
        }
      }
    }
    ML_CATCH_RETHROW;

    // Return the vector.
    return v;
  }

  //-------------------------------------------------------------------------------------------
  // Replaces the current kernel by a normalized gauss kernel.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::setGauss(const ImageVector &ext)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::setGauss()" );

    ML_TRY
    {
      // 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);
    }
    ML_CATCH_BLOCK(...)
    {
      // Value for an 1x1x1x1x1x1 kernel.
      KDATATYPE KerneldataIdKernel[]={ static_cast<KDATATYPE>(1.0f) };

      setKernel(ImageVector(1,1,1,1,1,1), KerneldataIdKernel);
      ML_CHECK(0); // Log the error.
      throw;
    }
  }


  //-------------------------------------------------------------------------------------------
  // Initialization. Should be called only by constructors.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::_init()
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::_init()" );

    _ext.set(0);
    _negExt.set(0);
    _posExt.set(0);

    _tabSize  = 0;

    _valueTab = NULL;
    _coordTab = NULL;

    // 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()
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::_clearTables()" );
    ML_TRY
    {
      // Set table size to 0.
      _tabSize = 0;

      // Invalidate separability infos which are calculated on demand.
      _areSeparabilityInfosUpToDate = false;

      ML_DELETE_ARRAY(_valueTab);
      ML_DELETE_ARRAY(_coordTab);

    }
    ML_CATCH_RETHROW;
  }



  //-------------------------------------------------------------------------------------------
  // 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)
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::_addCoordinate()" );

    // Create new data tables.
    ImageVector    *newCoordTab = NULL;
    KDATATYPE *newValueTab = NULL;

    ML_TRY
    {
      // Invalidate separability infos which are calculated on demand.
      _areSeparabilityInfosUpToDate = false;

      // Compute new table size.
      size_t newTabSize = _tabSize+1;

      ML_CHECK_NEW(newCoordTab, ImageVector [_tabSize+1]);
      ML_CHECK_NEW(newValueTab, KDATATYPE[_tabSize+1]);

      // Handle not enough memory.
      if (!newCoordTab || !newValueTab){
        ML_DELETE_ARRAY(newCoordTab);
        ML_DELETE_ARRAY(newValueTab);
        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){
        memcpy(newCoordTab, _coordTab, _tabSize*sizeof(ImageVector));
        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(); }
    }
    ML_CATCH_BLOCK(...)
    {
      ML_DELETE_ARRAY(newCoordTab);
      ML_DELETE_ARRAY(newValueTab);
      ML_CHECK(0);
      throw;
    }
  }


  //-------------------------------------------------------------------------------------------
  // 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()
  {
    ML_TRACE_IN( "TKernel<KDATATYPE>::_calculateRealKernelExt()" );
    ML_TRY
    {
      // 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.");
      }
    }
    ML_CATCH_RETHROW;
  }



  //-------------------------------------------------------------------------------------------
  //
  //        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)
  {
    ML_TRACE_IN_TIME_CRITICAL("TKernel<KDATATYPE>::setSeparable()");
    ML_TRY
    {
      // Invalidate separability information if flag changes.
      if (_isSeparable != isSeparableVal){
        _isSeparable = isSeparableVal;
        _invalidateSeparableInfos();
      }
    }
    ML_CATCH_RETHROW;
  }

  //-------------------------------------------------------------------------------------------
  // 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
  {
    ML_TRACE_IN_TIME_CRITICAL("TKernel<KDATATYPE>::isSeparable()");

    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
  {
    ML_TRACE_IN("TKernel<KDATATYPE>::getSeparableDimEntries()");

    _validateSeparabilityInfos();
    return _separableDimEntries;
  }

  //------------------------------------------------------------------------------------------
  //------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    ImageVector TKernel<KDATATYPE>::getSeparableOneDimExtents() const
  {
    ML_TRACE_IN("TKernel<KDATATYPE>::getSeparableOneDimExtents()");

    _validateSeparabilityInfos();
    return _separableOneDimExt;
  }

  //-------------------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    size_t TKernel<KDATATYPE>::getSeparableDimIndex(size_t dim) const
  {
    ML_TRACE_IN("TKernel<KDATATYPE>::getSeparableDimIndex()");

    ML_TRY
    {
      _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;
        }
      }
    }
    ML_CATCH_RETHROW;
  }

  //-------------------------------------------------------------------------------------------
  // Returns the table of all coordinates for filtering with separable kernels.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    const std::vector<ImageVector> &TKernel<KDATATYPE>::getSeparableCoordTab() const
  {
    ML_TRACE_IN("TKernel<KDATATYPE>::getSeparableCoordTab()");

    ML_TRY
    {
      _validateSeparabilityInfos();
    }
    ML_CATCH_RETHROW;

    return _separableCoordTab;
  }

  //-------------------------------------------------------------------------------------------
  // Invalidate separability information.
  //-------------------------------------------------------------------------------------------
  template <typename KDATATYPE>
    void TKernel<KDATATYPE>::_invalidateSeparableInfos()
  {
    ML_TRACE_IN_TIME_CRITICAL("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; }

    ML_TRACE_IN("TKernel<KDATATYPE>::_validateSeparabilityInfos()");
    ML_TRY
    {
      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];

          // 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[entry.y]++;

          // 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[entry.y]){ _separableOneDimExt[entry.y] = 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[entry.y] = 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 (int 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;
      }
    }
    ML_CATCH_BLOCK(...)
    {
      // Something went wrong, invalidate separability information.
      _areSeparabilityInfosUpToDate=false;
    }
  }

  //-------------------------------------------------------------------------------------------
  // 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