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 <array>
#include <cmath>
 
#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"
 
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, bool keepInside=true):
      myPosDimAdded(dimAdded), mySliceIndex(sliceIndex), myDomain(aDomain3DImg),
      myDimRotated(dimRotated), myRotationAngle(rotationAngle), myDefaultPoint (aDomain3DImg.lowerBound()), myKeepInsideDom(keepInside)
    {
      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, bool keepInside=true):
      myPosDimAdded(dimAdded), mySliceIndex(sliceIndex), myDomain(aDomain3DImg),
      myDimRotated(dimRotated), myRotationAngle(rotationAngle), myDefaultPoint (defaultPoint), myKeepInsideDom(keepInside)
    {
      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,
      bool keepInside=true):
      myPosDimAdded(dimAdded), mySliceIndex(sliceIndex), myDomain(aDomain3DImg),
      myDimRotated(dimRotated), myRotationAngle(rotationAngle), myCenter(ptCenter), myDefaultPoint (defaultPoint), myKeepInsideDom(keepInside)
    {
    };
 
    SliceRotator2D( const Dimension &dimAdded, const TDomain3D &aDomain3DImg, const Integer &sliceIndex,
                    const Dimension &dimRotated,  const Point &ptCenter, double rotationAngle, bool keepInside=true):
      myPosDimAdded(dimAdded), mySliceIndex(sliceIndex), myDomain(aDomain3DImg),
      myDimRotated(dimRotated), myRotationAngle(rotationAngle), myCenter(ptCenter), myDefaultPoint (aDomain3DImg.lowerBound()),
      myKeepInsideDom(keepInside)
    {
    };
 
    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(!myKeepInsideDom || myDomain.isInside(pt))
         return pt;
       else
          return  myDefaultPoint;
    }
  public:
    // 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;
    bool myKeepInsideDom = true;
  };
 
 
  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 & aNormalVector,
                         const typename Point::Component  &aWidth,
                         const Point &aDefautPoint = Point(0,0,0)): myDomain(aDomain3DImg),
                                                                    myDefaultPoint (aDefautPoint)
    {
      double d = -aNormalVector[0]*anOriginPoint[0] - aNormalVector[1]*anOriginPoint[1] - aNormalVector[2]*anOriginPoint[2];
      typename Space::RealPoint pRefOrigin;
      if(aNormalVector[0]!=0){
        pRefOrigin [0]= -d/aNormalVector[0];
        pRefOrigin [1]= 0.0;
        pRefOrigin [2]= 0.0;
        if(pRefOrigin==anOriginPoint){
          pRefOrigin[1]=-1.0;
        }
      }else if (aNormalVector[1]!=0){
        pRefOrigin [0]= 0.0;
        pRefOrigin [1]= -d/aNormalVector[1];
        pRefOrigin [2]= 0.0;
        if(pRefOrigin==anOriginPoint){
          pRefOrigin[0]=-1.0;
        }
      }else if (aNormalVector[2]!=0){
        pRefOrigin [0]= 0.0;
        pRefOrigin [1]= 0.0;
        pRefOrigin [2]= -d/aNormalVector[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]*aNormalVector[2]-uDir1[2]*aNormalVector[1];
      uDir2[1] = uDir1[2]*aNormalVector[0]-uDir1[0]*aNormalVector[2];
      uDir2[2] = uDir1[0]*aNormalVector[1]-uDir1[1]*aNormalVector[0];
 
      uDir2/=uDir2.norm();
      
      myOriginPointEmbeddedIn3D = anOriginPoint + Point(uDir1*aWidth/2) + Point(uDir2*aWidth/2);
      myFirstAxisEmbeddedDirection = -uDir1;
      mySecondAxisEmbeddedDirection = -uDir2;
    }
 
 
    Point2DEmbedderIn3D( const TDomain3D &aDomain3DImg,
                         const Point &anOriginPoint, const typename Space::RealPoint & aNormalVector,
                         const typename Space::RealPoint & orientXaxisVector,
                         const typename Point::Component  &aWidth,
                         const Point &aDefautPoint = Point(0,0,0)): myDomain(aDomain3DImg),
                                                                    myDefaultPoint (aDefautPoint)
    {
     
      typename Space::RealPoint uDir1;
      uDir1 = orientXaxisVector/orientXaxisVector.norm();
      typename Space::RealPoint uDir2;
      uDir2[0] = uDir1[1]*aNormalVector[2]-uDir1[2]*aNormalVector[1];
      uDir2[1] = uDir1[2]*aNormalVector[0]-uDir1[0]*aNormalVector[2];
      uDir2[2] = uDir1[0]*aNormalVector[1]-uDir1[1]*aNormalVector[0];
      uDir2/=uDir2.norm();
      myOriginPointEmbeddedIn3D = anOriginPoint + Point(-uDir1*aWidth/2) + Point(-uDir2*aWidth/2);
      myFirstAxisEmbeddedDirection = uDir1;
      mySecondAxisEmbeddedDirection = uDir2;
    }
 
 
    template <typename TPoint2D>
    inline
    Point  operator()(const TPoint2D& aPoint, bool checkInsideDomain=true) const
    {
      Point pt ;
      for( Dimension i=0; i<pt.size(); i++){
 
        pt[i] = static_cast<Integer>(floor(NumberTraits<Integer>::castToDouble(aPoint[0])
                                                 *myFirstAxisEmbeddedDirection[i]+myOriginPointEmbeddedIn3D[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;
        }
    }
 
 
    inline
    void shiftOriginPoint(const typename Space::RealPoint& shift)
      {
       for( Dimension i=0; i<myOriginPointEmbeddedIn3D.size(); i++){
         myOriginPointEmbeddedIn3D[i] = myOriginPointEmbeddedIn3D[i] + shift[i];
       }
      }
    
    
  private:
    TDomain3D myDomain;
 
    // Origin (or lower point) of the 2D image embedded in the 3D domain
    typename Space::RealPoint  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 = typename TDomain::Size >
 class BasicDomainSubSampler
  {
  public:
    typedef typename TDomain::Space  Space;
    typedef typename TDomain::Size Size;
    typedef typename TDomain::Integer IntergerDom;
    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] = static_cast<IntergerDom>(floor(NumberTraits<IntergerDom>::castToDouble(domainLowerBound[dim]) /
                                                               NumberTraits<TValue>::castToDouble( aGridSize[dim] )));        
        domainUpperBound[dim] = static_cast<IntergerDom>(floor(NumberTraits<IntergerDom>::castToDouble(domainUpperBound[dim]) /
                                                               NumberTraits<TValue>::castToDouble( 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<Dimension> & 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<Dimension> myAxisFlipped;
 };
 
 
 template <typename TRealVector, typename TVector>
 struct VectorRounding
 {
     BOOST_STATIC_ASSERT(( TRealVector::dimension ==  TVector::dimension ));
     inline
     TVector operator () ( const TRealVector & point ) const
     {
         TVector out;
         for ( unsigned int i = 0; i < TVector::dimension; i++ )
             out[i] = std::round ( point[i] );
         return out;
     }
 };
 
  }//namespace functors
} // namespace dgtal
 
 
// Includes inline functions.
#include "DGtal/kernel/BasicPointFunctors.ih"
//                                                                           //
 
#endif // !defined BasicPointFunctors_h
#undef BasicPointFunctors_RECURSES
#endif // else defined(BasicPointFunctors_RECURSES)