BasicPointFunctors.h

#pragma once

#if defined(BasicPointFunctors_RECURSES)
#error Recursive header files inclusion detected in BasicPointFunctors.h
#else // defined(BasicPointFunctors_RECURSES)

#define BasicPointFunctors_RECURSES

#if !defined BasicPointFunctors_h

#define BasicPointFunctors_h

// Inclusions
#include <iostream>
#include <iterator>
#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/kernel/SpaceND.h"
#include "DGtal/kernel/NumberTraits.h"
#include "DGtal/base/BasicBoolFunctors.h"
#include "DGtal/kernel/CPointPredicate.h"
#include "DGtal/base/CQuantity.h"
#include "DGtal/kernel/domains/CDomain.h"
#include "DGtal/base/ConstAlias.h"
#include <array>

namespace DGtal
{
namespace functors
  {
  // template class Projector
  template <typename S = SpaceND< 2, DGtal::Z2i::Integer > >
  struct Projector
  {
    typedef S Space;
    typedef typename Space::Dimension Dimension;
    BOOST_STATIC_CONSTANT( Dimension, dimension = Space::dimension );
    typedef typename Space::Integer Integer;
    typedef typename Space::Point Point;

    Projector(const Integer& aDefaultInteger = NumberTraits<Integer>::zero());


    template<typename TIterator>
    void init ( const TIterator& itb, const TIterator& ite );


    void initRemoveOneDim ( const Dimension  &dimRemoved );

    void initAddOneDim ( const Dimension & newDim );


    template<typename TInputPoint>
    Point operator()( const TInputPoint& aPoint ) const;

  private:
    std::array<Dimension, dimension> myDims;

    Integer myDefaultInteger;

  }; // end of class ConstantPointFunctors





  template <typename TDomain3D, typename TInteger =  DGtal::Z3i::Integer >
  class SliceRotator2D
  {
  public:

    typedef SpaceND< 3, TInteger>  Space;
    typedef typename Space::Dimension Dimension;
    typedef typename Space::Point Point;
    typedef typename Space::Integer Integer;
    SliceRotator2D( const Dimension &dimAdded, const TDomain3D &aDomain3DImg,
                    const Integer &sliceIndex,  const Dimension &dimRotated,
                    double rotationAngle):
      myPosDimAdded(dimAdded), mySliceIndex(sliceIndex), myDomain(aDomain3DImg),
      myDimRotated(dimRotated), myRotationAngle(rotationAngle), myDefaultPoint (aDomain3DImg.lowerBound())
    {
      myCenter[0] = aDomain3DImg.lowerBound()[0]+((aDomain3DImg.upperBound())[0]-(aDomain3DImg.lowerBound())[0])/2.0;
      myCenter[1] = aDomain3DImg.lowerBound()[1]+((aDomain3DImg.upperBound())[1]-(aDomain3DImg.lowerBound())[1])/2.0;
      myCenter[2] = aDomain3DImg.lowerBound()[2]+((aDomain3DImg.upperBound())[2]-(aDomain3DImg.lowerBound())[2])/2.0;
      myCenter[dimAdded]=sliceIndex;
    };

    SliceRotator2D( const Dimension &dimAdded, const TDomain3D &aDomain3DImg,
                    const Integer &sliceIndex,  const Dimension &dimRotated,
                    double rotationAngle,  const Point &defaultPoint):
      myPosDimAdded(dimAdded), mySliceIndex(sliceIndex), myDomain(aDomain3DImg),
      myDimRotated(dimRotated), myRotationAngle(rotationAngle), myDefaultPoint (defaultPoint)
    {
      myCenter[0] = aDomain3DImg.lowerBound()[0]+((aDomain3DImg.upperBound())[0]-(aDomain3DImg.lowerBound())[0])/2.0;
      myCenter[1] = aDomain3DImg.lowerBound()[1]+((aDomain3DImg.upperBound())[1]-(aDomain3DImg.lowerBound())[1])/2.0;
      myCenter[2] = aDomain3DImg.lowerBound()[2]+((aDomain3DImg.upperBound())[2]-(aDomain3DImg.lowerBound())[2])/2.0;
      myCenter[dimAdded]=sliceIndex;
    };

    SliceRotator2D( const Dimension &dimAdded, const TDomain3D &aDomain3DImg, const Integer &sliceIndex,
                    const Dimension &dimRotated,  const Point &ptCenter, double rotationAngle, const Point &defaultPoint):
      myPosDimAdded(dimAdded), mySliceIndex(sliceIndex), myDomain(aDomain3DImg),
      myDimRotated(dimRotated), myRotationAngle(rotationAngle), myCenter(ptCenter), myDefaultPoint (defaultPoint)
    {
    };

    SliceRotator2D( const Dimension &dimAdded, const TDomain3D &aDomain3DImg, const Integer &sliceIndex,
                    const Dimension &dimRotated,  const Point &ptCenter, double rotationAngle):
      myPosDimAdded(dimAdded), mySliceIndex(sliceIndex), myDomain(aDomain3DImg),
      myDimRotated(dimRotated), myRotationAngle(rotationAngle), myCenter(ptCenter), myDefaultPoint (aDomain3DImg.lowerBound())
    {
    };

    template <typename TPointDimMinus>
    inline
    Point  operator()(const TPointDimMinus& aPoint) const
    {
      Point pt;
      Dimension pos=0;
      std::vector<Dimension> indexesRotate;
      for( Dimension i=0; i<pt.size(); i++)
        {
          if(i!=myPosDimAdded)
            {
              pt[i]= aPoint[pos];
              pos++;
            }else
            {
              pt[i]=mySliceIndex;
            }
        }
      for( Dimension i=0; i<pt.size(); i++)
        {
          if(i!=myDimRotated)
            indexesRotate.push_back(i);
        }
      double d1 = pt[indexesRotate[0]] - myCenter[indexesRotate[0]];
      double d2 = pt[indexesRotate[1]] - myCenter[indexesRotate[1]];

      pt[indexesRotate[0]] = myCenter[indexesRotate[0]] + static_cast<Integer>(round(d1*cos(myRotationAngle)-d2*sin(myRotationAngle) ));
      pt[indexesRotate[1]] = myCenter[indexesRotate[1]] + static_cast<Integer>(round(d1*sin(myRotationAngle)+d2*cos(myRotationAngle) ));

      if(myDomain.isInside(pt))
        return pt;
      else
        return  myDefaultPoint;
    }
  private:
    // position of insertion of the new dimension
    Dimension myPosDimAdded;
    // the index of the slice associated to the new domain.
    Integer mySliceIndex;
    TDomain3D myDomain;
    Dimension myDimRotated;
    double myRotationAngle;
    PointVector<3, double> myCenter;
    Point myDefaultPoint;
  };


  template <typename TDomain3D, typename TInteger =  DGtal::Z3i::Integer >
  class Point2DEmbedderIn3D
  {
  public:

    typedef SpaceND< 3, TInteger>  Space;
    typedef typename Space::Point Point;
    typedef typename Space::Integer Integer;

    Point2DEmbedderIn3D( const TDomain3D &aDomain3DImg,
                         const Point &anOriginPoint, const Point &anUpperPointOnAxis1,
                         const Point &anUpperPointOnAxis2,
                         const Point &aDefautPoint = Point(0,0,0)): myDomain(aDomain3DImg),
                                                                    myOriginPointEmbeddedIn3D(anOriginPoint),
                                                                    myDefaultPoint (aDefautPoint),
                                                                    myFirstAxisEmbeddedDirection(Point(anUpperPointOnAxis1[0]-anOriginPoint[0],
                                                                                                       anUpperPointOnAxis1[1]-anOriginPoint[1],
                                                                                                       anUpperPointOnAxis1[2]-anOriginPoint[2])),
      mySecondAxisEmbeddedDirection(Point(anUpperPointOnAxis2[0]-anOriginPoint[0],
                                          anUpperPointOnAxis2[1]-anOriginPoint[1],
                                          anUpperPointOnAxis2[2]-anOriginPoint[2]))


    {
      myFirstAxisEmbeddedDirection /= myFirstAxisEmbeddedDirection.norm();
      mySecondAxisEmbeddedDirection /= mySecondAxisEmbeddedDirection.norm();
    }


    Point2DEmbedderIn3D( const TDomain3D &aDomain3DImg,
                         const Point &anOriginPoint, const typename Space::RealPoint & anNormalVector,
                         const typename Point::Component  &anWidth,
                         const Point &aDefautPoint = Point(0,0,0)): myDomain(aDomain3DImg),
                                                                    myDefaultPoint (aDefautPoint)
    {
      double d = -anNormalVector[0]*anOriginPoint[0] - anNormalVector[1]*anOriginPoint[1] - anNormalVector[2]*anOriginPoint[2];
      typename Space::RealPoint pRefOrigin;
      if(anNormalVector[0]!=0){
        pRefOrigin [0]= -d/anNormalVector[0];
        pRefOrigin [1]= 0.0;
        pRefOrigin [2]= 0.0;
        if(pRefOrigin==anOriginPoint){
          pRefOrigin[1]=-1.0;
        }
      }else if (anNormalVector[1]!=0){
        pRefOrigin [0]= 0.0;
        pRefOrigin [1]= -d/anNormalVector[1];
        pRefOrigin [2]= 0.0;
        if(pRefOrigin==anOriginPoint){
          pRefOrigin[0]=-1.0;
        }
      }else if (anNormalVector[2]!=0){
        pRefOrigin [0]= 0.0;
        pRefOrigin [1]= 0.0;
        pRefOrigin [2]= -d/anNormalVector[2];
        if(pRefOrigin==anOriginPoint){
          pRefOrigin[0]=-1.0;
        }
      }
      typename Space::RealPoint uDir1;
      uDir1=(pRefOrigin-anOriginPoint)/((pRefOrigin-anOriginPoint).norm());
      typename Space::RealPoint uDir2;
      uDir2[0] = uDir1[1]*anNormalVector[2]-uDir1[2]*anNormalVector[1];
      uDir2[1] = uDir1[2]*anNormalVector[0]-uDir1[0]*anNormalVector[2];
      uDir2[2] = uDir1[0]*anNormalVector[1]-uDir1[1]*anNormalVector[0];

      uDir2/=uDir2.norm();

      myOriginPointEmbeddedIn3D = anOriginPoint + uDir1*anWidth/2 + uDir2*anWidth/2;
      myFirstAxisEmbeddedDirection = -uDir1;
      mySecondAxisEmbeddedDirection = -uDir2;

    }




    template <typename TPoint2D>
    inline
    Point  operator()(const TPoint2D& aPoint, bool checkInsideDomain=true) const
    {
      Point pt = myOriginPointEmbeddedIn3D;
      for( Dimension i=0; i<pt.size(); i++){

        pt[i] = pt[i]+static_cast<Integer>(floor(NumberTraits<Integer>::castToDouble(aPoint[0])
                                                 *myFirstAxisEmbeddedDirection[i]));
        pt[i] = pt[i]+static_cast<Integer>(floor(NumberTraits<Integer>::castToDouble(aPoint[1])
                                                 *mySecondAxisEmbeddedDirection[i]));
      }

      if(myDomain.isInside(pt)|| !checkInsideDomain)
        {
          return pt;
        }
      else
        {
#ifdef DEBUG_VERBOSE
          trace.warning() << "Warning pt outside the 3D domain " << pt << std::endl;
#endif
          return  myDefaultPoint;
        }
    }

  private:
    TDomain3D myDomain;

    // Origin (or lower point) of the 2D image embedded in the 3D domain
    Point  myOriginPointEmbeddedIn3D;

    Point myDefaultPoint;

    // Point giving the direction of the embedded first axis of the 2D image.
    typename Space::RealPoint myFirstAxisEmbeddedDirection;

    // Point giving the direction of the embedded second axis of the 2D image.
    typename Space::RealPoint mySecondAxisEmbeddedDirection;

  };


  template< typename TPointPredicate, typename TDomain, typename TValue=typename TDomain::Integer >
  struct PointFunctorFromPointPredicateAndDomain
  {
      typedef TPointPredicate PointPredicate;
      typedef TDomain Domain;
      typedef TValue Value;
      typedef typename Domain::Point Point;

      BOOST_CONCEPT_ASSERT(( concepts::CPointPredicate< PointPredicate > ));
      BOOST_CONCEPT_ASSERT(( concepts::CDomain< Domain > ));
      BOOST_CONCEPT_ASSERT(( concepts::CQuantity< Value > ));

      PointFunctorFromPointPredicateAndDomain( ConstAlias< PointPredicate > aPtrPredicate, ConstAlias< Domain > aDomain,
                                               const Value aTrueValue, const Value aFalseValue );

      PointFunctorFromPointPredicateAndDomain( const PointFunctorFromPointPredicateAndDomain & other  );

      Value operator()( const Point& aPoint ) const;

      PointFunctorFromPointPredicateAndDomain & operator= ( const PointFunctorFromPointPredicateAndDomain & other );

  private:
      const PointPredicate * myPtrPredicate;
      const Domain * myDomain;
      Value myTrueValue;
      Value myFalseValue;

  };


 template <typename TDomain, typename TInteger =  DGtal::int32_t, typename TValue = DGtal::uint32_t >
 class BasicDomainSubSampler
  {
  public:
    typedef typename TDomain::Space  Space;
    typedef typename TDomain::Size Size;
    typedef typename Space::Dimension Dimension;
    typedef typename Space::Point Point;

    BasicDomainSubSampler(const TDomain &aSourceDomain, const std::vector<TValue> &aGridSize,
                          const Point  &aGridShift): mySourceDomain(aSourceDomain),
                                                     myGridShift(aGridShift),
                                                     myGridSize(aGridSize)
    {
      Point domainUpperBound=aSourceDomain.upperBound();
      Point domainLowerBound=aSourceDomain.lowerBound();

      for (Dimension dim=0; dim< Space::dimension; dim++){
        domainLowerBound[dim] /= aGridSize[dim];
        domainUpperBound[dim] /= aGridSize[dim];
      }
      myNewDomain = TDomain(domainLowerBound,
                            domainUpperBound);
      Point upperGrid;
      for (Dimension dim=0; dim < Space::dimension; dim++)
        upperGrid[dim] = myGridSize[dim];
      myGridSampleDomain = TDomain(Point::diagonal(0), upperGrid);
    };


    inline
    Point  operator()(const Point& aPoint) const
    {
      Point ptRes = Point::diagonal(0);
      if(!myNewDomain.isInside(aPoint)){
        trace.error() << " The point is not in the source domain: "<<  aPoint << std::endl;
        return ptRes;
      }

      for (Dimension dim=0; dim< Space::dimension; dim++){
        ptRes[dim] = static_cast<TInteger>(floor(NumberTraits<TInteger>::castToDouble(aPoint[dim])*
                                                 NumberTraits<TValue>::castToDouble(myGridSize[dim])));
      }
      ptRes +=myGridShift;

      if(!mySourceDomain.isInside(ptRes)){
        // we are looking for a point inside the domain
        for(typename TDomain::ConstIterator it = myGridSampleDomain.begin();
            it!= myGridSampleDomain.end(); it++){
          if (mySourceDomain.isInside(ptRes+(*it)))
            return ptRes+(*it);
        }
      }
      return ptRes;
    }

    inline
    const TDomain & getSubSampledDomain(){
      return myNewDomain;
    }

  private:
    TDomain mySourceDomain;
    TDomain myNewDomain;
    // used to search a point when the resulting point is outside the source domain.
    TDomain myGridSampleDomain;
    Point myGridShift;
    std::vector<TValue> myGridSize;
 };


 template <typename TDomain>
 class FlipDomainAxis
  {
  public:
    typedef typename TDomain::Space  Space;
    typedef typename TDomain::Size Size;
    typedef typename Space::Dimension Dimension;
    typedef typename Space::Point Point;

    FlipDomainAxis(const TDomain &aSourceDomain, const std::vector<Size> & axisFlipped): mySourceDomain(aSourceDomain),
                                                                                         myAxisFlipped(axisFlipped){
    };


    inline
    Point  operator()(const Point& aPoint) const
    {
      Point ptRes;
      for (Dimension dim=0; dim< Space::dimension; dim++){
        ptRes[dim] = aPoint[dim];
      }
      for(Dimension i = 0; i< myAxisFlipped.size(); i++){
        ptRes[myAxisFlipped[i]] = mySourceDomain.upperBound()[myAxisFlipped[i]]-aPoint[myAxisFlipped[i]];
      }
      return ptRes;
    }

  private:
    TDomain mySourceDomain;
    std::vector<Size> myAxisFlipped;
 };









  }//namespace functors
} // namespace dgtal


// Includes inline functions.
#include "DGtal/kernel/BasicPointFunctors.ih"
//                                                                           //

#endif // !defined BasicPointFunctors_h
#undef BasicPointFunctors_RECURSES
#endif // else defined(BasicPointFunctors_RECURSES)