mlIntegerVector.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_INTEGER_VECTOR_H
#define ML_INTEGER_VECTOR_H

 
// ML-includes
#include "mlLinearAlgebra.h"
#include <mlTypeTraits.h>
#include <mlUtilsSystemC.h>
 
 
ML_START_NAMESPACE
 
  //-------------------------------------------------------------------------
  //-------------------------------------------------------------------------
  typedef int VectorDimIdxType;
 
  //-------------------------------------------------------------------------
  //-------------------------------------------------------------------------
  template <typename CompIntType, VectorDimIdxType NumDim>
  class TVectorNDBase
  {
    public:

    typedef CompIntType ComponentType;

    enum { NumberOfDimensions = NumDim };

    ComponentType array[NumberOfDimensions];
 
    protected:

    TVectorNDBase(){};
  };
 
  //-------------------------------------------------------------------------

  //-------------------------------------------------------------------------
  template<class TVectorBase>
  class TVector : public TVectorBase
  {
  public:

    typedef typename TVectorBase::ComponentType ComponentType;

    typedef TVectorBase ParentClass;

    enum { NumberOfDimensions = ParentClass::NumberOfDimensions };
 
 
    //------------------------------------------------------
    //------------------------------------------------------
    TVector() : ParentClass()
    {
      memset(this->array, 0, static_cast<size_t>(NumberOfDimensions)*sizeof(ComponentType));
    }

    explicit TVector(const ComponentType i) : ParentClass()
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d)
      {
        this->array[d]=i;
      }
    }
 
    // Even though the compiler should generate the copy constructor for us, gcc9+10 have a bug here, so implementation by hand again
    TVector(const TVector<TVectorBase>& v2) : ParentClass()
    {
      *this = v2;
    }
 
    TVector<TVectorBase>& operator=(const TVector<TVectorBase>& v2)
    {
      for (VectorDimIdxType d = 0; d < NumberOfDimensions; ++d)
      {
        this->array[d] = v2.array[d];
      }
      return *this;
    }

    template <typename TVectorBase2>
    explicit TVector(const TVector<TVectorBase2> &v2) : ParentClass()
    {
      VectorDimIdxType numDims = NumberOfDimensions;
      if (v2.NumberOfDimensions < numDims)
      {
        numDims = v2.NumberOfDimensions;
      }
      VectorDimIdxType d=0;
      for (; d < numDims; ++d)
      {
        this->array[d] = static_cast<ComponentType>(v2.array[d]);
      }
      for (; d < NumberOfDimensions; ++d)
      {
        this->array[d] = 0;
      }
    }

    explicit TVector(const ComponentType xp, const TVector &p) : ParentClass()
    {
      this->array[0] = xp;
      for (VectorDimIdxType d=1; d < NumberOfDimensions; ++d){ this->array[d] = p.array[d]; }
    }

    explicit TVector(const MLint num, const ComponentType* const arr, const ComponentType deflt);

    void copy(const MLint num, ComponentType* const arr) const;

    void set(const ComponentType v=0)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d]=v; }
    }

    void set(const TVector &v)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] = v.array[d]; }
    }

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

    TVector<TVectorBase>  operator+ (const TVector &b) const
    {
      return TVector<TVectorBase>(*this) += b;
    }

    TVector<TVectorBase> &operator+= (const TVector<TVectorBase> &b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] += b.array[d]; }
      return *this;
    }

    TVector<TVectorBase>  operator- (const TVector<TVectorBase> &b) const
    {
      return TVector<TVectorBase>(*this) -= b;
    }

    TVector<TVectorBase> &operator-= (const TVector<TVectorBase> &b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] -= b.array[d]; }
      return *this;
    }

    TVector<TVectorBase>  operator* (const TVector<TVectorBase> &b) const
    {
      return TVector<TVectorBase>(*this) *= b;
    }

    TVector<TVectorBase> &operator*= (const TVector<TVectorBase> &b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] *= b.array[d]; }
      return *this;
    }

    TVector<TVectorBase>  operator+ (const ComponentType b) const
    {
      TVector<TVectorBase> retVec(*this);
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ retVec.array[d] += b; }
      return retVec;
    }

    TVector<TVectorBase>  operator- (const ComponentType b) const
    {
      TVector<TVectorBase> retVec(*this);
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ retVec.array[d] -= b; }
      return retVec;
    }

    TVector<TVectorBase>  operator* (const ComponentType b) const
    {
      TVector<TVectorBase> retVec(*this);
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ retVec.array[d] *= b; }
      return retVec;
    }

    TVector<TVectorBase>  operator/ (const TVector<TVectorBase> &b) const
    {
      return TVector<TVectorBase>(*this) /= b;
    }

    TVector<TVectorBase> &operator/= (const TVector<TVectorBase> &b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){
        ML_CHECK_THROW(b.array[d]);
        this->array[d] /= b.array[d];
      }
      return *this;
    }

    bool operator== (const TVector<TVectorBase> &b) const
    {
      bool retFlag = true;
      for (VectorDimIdxType d=0; retFlag && (d < NumberOfDimensions); ++d){ retFlag = retFlag && (this->array[d] == b.array[d]); }
      return retFlag;
    }

    bool operator!= (const TVector<TVectorBase> &b) const
    {
      return !(*this == b);
    }

    bool operator< (const TVector<TVectorBase> &b) const
    {
      bool retFlag = true;
      for (VectorDimIdxType d=0; retFlag && (d < NumberOfDimensions); ++d){ retFlag = retFlag && (this->array[d] < b.array[d]); }
      return retFlag;
    }

    bool operator<= (const TVector<TVectorBase> &b) const
    {
      bool retFlag = true;
      for (VectorDimIdxType d=0; retFlag && (d < NumberOfDimensions); ++d){ retFlag = retFlag && (this->array[d] <= b.array[d]); }
      return retFlag;
    }

    bool operator> (const TVector<TVectorBase> &b) const
    {
      bool retFlag = true;
      for (VectorDimIdxType d=0; retFlag && (d < NumberOfDimensions); ++d){ retFlag = retFlag && (this->array[d] > b.array[d]); }
      return retFlag;
    }

    bool operator>= (const TVector<TVectorBase> &b) const
    {
      bool retFlag = true;
      for (VectorDimIdxType d=0; retFlag && (d < NumberOfDimensions); ++d){ retFlag = retFlag && (this->array[d] >= b.array[d]); }
      return retFlag;
    }

    TVector<TVectorBase>  operator& (const TVector<TVectorBase> &b) const
    {
      return TVector<TVectorBase>(*this) &= b;
    }

    TVector<TVectorBase>  operator| (const TVector<TVectorBase> &b) const
    {
      return TVector<TVectorBase>(*this) |= b;
    }

    TVector<TVectorBase>  operator^  (const TVector<TVectorBase> &b) const
    {
      return TVector<TVectorBase>(*this) ^= b;
    }

    TVector<TVectorBase> &operator &= (const TVector<TVectorBase> &b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] &= b.array[d]; }
      return *this;
    }

    TVector<TVectorBase> &operator |= (const TVector<TVectorBase> &b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] |= b.array[d]; }
      return *this;
    }

    TVector<TVectorBase> &operator ^= (const TVector<TVectorBase> &b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] ^= b.array[d]; }
      return *this;
    }

    TVector<TVectorBase>  operator& (const ComponentType b) const
    {
      return TVector<TVectorBase>(*this) &= b;
    }

    TVector<TVectorBase>  operator| (const ComponentType b) const
    {
      return TVector<TVectorBase>(*this) |= b;
    }

    TVector<TVectorBase>  operator^ (const ComponentType b) const
    {
      return TVector<TVectorBase>(*this) ^= b;
    }

    TVector<TVectorBase> &operator&= (const ComponentType b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] &= b; }
      return *this;
    }

    TVector<TVectorBase> &operator|= (const ComponentType b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] |= b; }
      return *this;
    }

    TVector<TVectorBase> &operator^= (const ComponentType b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] ^= b; }
      return *this;
    }

    TVector<TVectorBase>  operator>> (const int b) const
    {
      return TVector<TVectorBase>(*this) >>= b;
    }
    TVector<TVectorBase>  operator<< (const int b) const
    {
      return TVector<TVectorBase>(*this) <<= b;
    }

    TVector<TVectorBase> &operator>>= (const int b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] >>= b; }
      return *this;
    }
    TVector<TVectorBase> &operator<<= (const int b)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] <<= b; }
      return *this;
    }

    TVector<TVectorBase>  operator- () const
    {
      TVector<TVectorBase> retVec(*this);
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ retVec.array[d] = -retVec.array[d]; }
      return retVec;
    }

    // Note: No tracing! No try, No catch! Do nothing here that could degrade performance!
    ComponentType operator[](const size_t i) const { return this->array[i]; }

    // Note: No tracing! No try, No catch! Do nothing here that could degrade performance!
    ComponentType &operator[](const size_t i)      { return this->array[i]; }

    MLint getMaxIdx() const
    {
      VectorDimIdxType idx = NumberOfDimensions-1;
      ComponentType      val = this->array[idx];
      for (VectorDimIdxType d=NumberOfDimensions-2; d >=0; d--){
        if (this->array[d] >= val){ idx=d; val=this->array[d]; }
      }
      return idx;
    }

    ComponentType getMax() const
    {
      return this->array[getMaxIdx()];
    }

    MLint getMinIdx() const
    {
      VectorDimIdxType idx = NumberOfDimensions-1;
      ComponentType      val = this->array[idx];
      for (VectorDimIdxType d=NumberOfDimensions-2; d >= 0; d--){
        if (this->array[d] <= val){ idx=d; val=this->array[d]; }
      }
      return idx;
    }

    ComponentType getMin() const
    {
      return this->array[getMinIdx()];
    }

    TVector<TVectorBase> getStrides(const ComponentType offset=1) const;

    TVector<TVectorBase> getVectorPosition(ComponentType offsetPos) const;

    MLint getExtDimension() const;

    bool hasNegativeComp() const
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ if (this->array[d]<0){ return true; } }
      return false;
    }

    bool allBiggerZero() const
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ if (this->array[d]<=0){ return false; } }
      return true;
    }

    bool sanityCheckNumberOfVoxels() const
    {
      double product = 1.0;
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d)
      { 
        product *= static_cast<double>(this->array[d]);
      }
      return _ml_numeric_limits<ComponentType>::MaxValue >= product;
    }

 
 
    //------------------------------------------------------
    //------------------------------------------------------
    // Methods:
    TVector<TVectorBase> divCeil(const TVector<TVectorBase> &b) const { return TVector<TVectorBase>::divCeil(*this, b); }
    TVector<TVectorBase> compMin(const TVector<TVectorBase> &b) const { return TVector<TVectorBase>::compMin(*this, b); }
    TVector<TVectorBase> compMax(const TVector<TVectorBase> &b) const { return TVector<TVectorBase>::compMax(*this, b); }
    TVector<TVectorBase> compMod(const TVector<TVectorBase> &b) const { return TVector<TVectorBase>::compMod(*this, b); }
    TVector<TVectorBase> compAbs() const { return TVector<TVectorBase>::compAbs(*this);    }

    ComponentType compMul() const { return TVector<TVectorBase>::compMul(*this); }
    ComponentType compSum() const { return TVector<TVectorBase>::compSum(*this); }

    ComponentType dot(const TVector<TVectorBase> &b) const { return TVector<TVectorBase>::dot(*this, b);     }

    static TVector<TVectorBase> divCeil(const TVector<TVectorBase> &a, const TVector<TVectorBase> &b);
    static TVector<TVectorBase> compMin(const TVector<TVectorBase> &a, const TVector<TVectorBase> &b);
    static TVector<TVectorBase> compMax(const TVector<TVectorBase> &a, const TVector<TVectorBase> &b);
    static TVector<TVectorBase> compMod(const TVector<TVectorBase> &a, const TVector<TVectorBase> &b);
    static TVector<TVectorBase> compAbs(const TVector<TVectorBase> &a);

    static ComponentType compMul(const TVector<TVectorBase> &a)
    {
      ComponentType retVal = a.array[0];
      for (VectorDimIdxType d=1; d < NumberOfDimensions; ++d)
      {
        retVal *= a.array[d];
      }
      return retVal;
    }

    static ComponentType   compSum(const TVector<TVectorBase> &a)
    {
      ComponentType retVal = a.array[0];
      for (VectorDimIdxType d=1; d < NumberOfDimensions; ++d)
      {
        retVal += a.array[d];
      }
      return retVal;
    }
 

    static ComponentType dot(const TVector<TVectorBase> &a, const TVector<TVectorBase> &b)
    {
      ComponentType retVal = a.array[0]*b.array[0];
      for (VectorDimIdxType d=1; d < NumberOfDimensions; ++d)
      {
        retVal += a.array[d]*b.array[d];
      }
      return retVal;
    }

    void fillSmallerComps(const ComponentType threshold=0, const ComponentType fillVal=0)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d)
      {
        if (this->array[d] < threshold)
        {
          this->array[d] = fillVal;
        }
      }
    }

    void fillGreaterComps(const ComponentType threshold=0, const ComponentType fillVal=0)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ if (this->array[d] > threshold){ this->array[d] = fillVal; } }
    }

    void fillEqualComps(const ComponentType threshold=0, const ComponentType fillVal=0)
    {
      for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ if (this->array[d] == threshold){ this->array[d] = fillVal; } }
    }


    std::string print(const std::string &openStr="",
                      const std::string &sepStr=" ",
                      const std::string &termStr="",
                      const MLint16 numDigits=-1) const;
 
  };
 
ML_END_NAMESPACE
 
 
//-----------------------------------------------------------------------------------
//   Stream output for std::ostream
//-----------------------------------------------------------------------------------
namespace std {

  template <class TVectorBase>
  ostream& operator<<(ostream& s, const ML_NAMESPACE::TVector<TVectorBase> &v){
    s << "(";
    for (ML_NAMESPACE::VectorDimIdxType d=0; d < ML_NAMESPACE::TVector<TVectorBase>::NumberOfDimensions-1; ++d){ s << v.array[d] << ","; }
    s << v.array[ML_NAMESPACE::TVector<TVectorBase>::NumberOfDimensions-1];
    s << ")";
    return s;
  }
 
}
 
 
 
ML_START_NAMESPACE
 
 
  //------------------------------------------------------
  //            TVector<TVectorBase>
  //------------------------------------------------------
  // Constructor that initializes the vector components by copying them from
  // arr[0] ,..., arr[num-1]. All other elements are initialized to
  // the value \p deflt.
  template <class TVectorBase>
  TVector<TVectorBase>::TVector(const MLint num, const ComponentType* const arr, const ComponentType deflt) : ParentClass()
  {
    ML_CHECK_THROW(arr);
    for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ this->array[d] = (num>d ? arr[d] : deflt); }
  }
 
  // Copies the first \p num elements from internal array into \p arr,
  // starting with array[0] to arr[0] to array[num-1] to arr[num-1], respectively.
  // All other components are left unchanged.
  // \p arr must have at least \p num entries and at most NumberOfDimensions elements are copied.
  template <class TVectorBase>
  void TVector<TVectorBase>::copy(const MLint num, ComponentType* const arr) const
  {
    ML_CHECK_THROW(arr);
    for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){ if (num>d){ arr[d] = this->array[d]; } }
  }
 
  // Interprets vector as an image extent and returns a stride vector.
  // The smallest stride (the x stride) is set to offset, i.e., an
  // image voxel is assumed to have size = offset.
  // This way, stride.x = offset, stride.y = offset*x, stride.z = offset*x*y, ...,
  // stride.u=offset*x*y*z*c*t.
  template <class TVectorBase>
  TVector<TVectorBase> TVector<TVectorBase>::getStrides(const ComponentType offset) const
  {
    TVector<TVectorBase> stride;
    stride.array[0] = offset;
    for (VectorDimIdxType d=1; d < NumberOfDimensions; ++d){ stride.array[d] = stride.array[d-1] * this->array[d-1]; }
    return stride;
  }
 
  // Interprets the vector as a stride vector of an image and returns the offset
  // position \p offsetPos from the image origin to the coordinate specified by
  // the vector. For offsetPos, an image with extents given by *this and the return
  // value rv then holds:
  // \code offsetPos = rv.x*x + rv.y*y + ... + rv.u*u \endcode )
  template <class TVectorBase>
  TVector<TVectorBase> TVector<TVectorBase>::getVectorPosition(ComponentType offsetPos) const
  {
    TVector<TVectorBase> vectorPos;
    for (VectorDimIdxType d=NumberOfDimensions-1; d>=1; d--){
      ML_CHECK_THROW(this->array[d]);
      vectorPos.array[d] = offsetPos / this->array[d];
      offsetPos = offsetPos % this->array[d];
    }
    ML_CHECK_THROW(this->array[0]);
    vectorPos.array[0] = offsetPos / this->array[0];
 
    return vectorPos;
  }
 
  // Returns the highest vector component that is not 1. Useful to get the actual dimension of an image.
  template <class TVectorBase>
  MLint TVector<TVectorBase>::getExtDimension() const
  {
    VectorDimIdxType dimen = 1;
    for (VectorDimIdxType d=NumberOfDimensions-1; d >=1; d--){
      if (1 != this->array[d]){
        dimen = d+1;
        break;
      }
    }
    return dimen;
  }
 
  // Returns the component-wise division of the vector this by vector \p b with rounding up.
  // Division by zero caused by zero components of \p b must be avoided by caller.
  template <class TVectorBase>
  TVector<TVectorBase> TVector<TVectorBase>::divCeil(const TVector<TVectorBase> &a, const TVector<TVectorBase> &b)
  {
    TVector<TVectorBase> c;
    for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){
      ML_CHECK_THROW(b.array[d]);
      c.array[d] = static_cast<ComponentType>(ceil(static_cast<MLdouble>(a.array[d])/b.array[d]));
    }
    return c;
  }
 
 
  // Returns the component-wise minimum of TVector<TVectorBase>s this and \p b.
  template <class TVectorBase>
  TVector<TVectorBase> TVector<TVectorBase>::compMin(const TVector<TVectorBase> &a, const TVector<TVectorBase> &b)
  {
    TVector<TVectorBase> c(b);
    for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){
      if (a.array[d] < c.array[d]){ c.array[d] = a.array[d]; }
    }
    return c;
  }
 
  // Returns the component-wise maximum of vectors this and \p b.
  template <class TVectorBase>
  TVector<TVectorBase> TVector<TVectorBase>::compMax(const TVector<TVectorBase> &a, const TVector<TVectorBase> &b)
  {
    TVector<TVectorBase> c(b);
    for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){
      if (a.array[d] > c.array[d]){ c.array[d] = a.array[d]; }
    }
    return c;
  }
 
  // Returns the rest dividing vector this a by \p b component-wise.
  // Division by zero caused by zero components of \p b must be avoided by caller and will
  // be handled as fatal error exception otherwise.
  template <class TVectorBase>
  TVector<TVectorBase> TVector<TVectorBase>::compMod(const TVector<TVectorBase> &a, const TVector<TVectorBase> &b)
  {
    TVector<TVectorBase> c;
    for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){
      ML_CHECK_THROW(b.array[d]);
      c.array[d] = a.array[d] % b.array[d];
    }
    return c;
  }
 
  // Returns a vector with all components of a replaced by their absolute values.
  template <class TVectorBase>
  TVector<TVectorBase> TVector<TVectorBase>::compAbs(const TVector<TVectorBase> &a)
  {
    TVector<TVectorBase> c(a);
    for (VectorDimIdxType d=0; d < NumberOfDimensions; ++d){
      if (c.array[d]<0){ c.array[d] = -c.array[d]; }
    }
    return c;
  }
 
  // Returns the string openStr + v.x + sepStr + ... + v.t + sepStr + v.u + termStr.
  // \p numDigits is the number of digits to be used as print parameter.
  // If numDigits is <=0, then no space specification is used.
  template <class TVectorBase>
  std::string TVector<TVectorBase>::print(const std::string &openStr,
                                                 const std::string &sepStr,
                                                 const std::string &termStr,
                                                 const MLint16      numDigits) const
  {
    std::string result = openStr;
    // 512 chars are sufficient, because MLint64 or even in future MLint
    // as 128 bit version does not use more than 2^128 digits, which are
    // about 10^43.
    char buf [512]="";
    char frmt[100]="";
 
    // Create a format string with or without space specifications
    // dependent on numDigits.
    if (numDigits>=0){
      snprintf(frmt, 99, "%%%dlld", numDigits);
    }
    else{
      snprintf(frmt, 99, "%%lld");
    }
 
    for (VectorDimIdxType d=0; d < NumberOfDimensions-1; ++d){
      snprintf(buf, 511, frmt, this->array[d]);
      result += buf + sepStr;
    }
    snprintf(buf, 511, frmt, this->array[NumberOfDimensions-1]);
    result += buf + termStr;
    return result;
  }
 
ML_END_NAMESPACE
 
#endif //of __mlIntegerVector_H