IntersectionTarget.h

#pragma once

#if !defined IntersectionTarget_h

#define IntersectionTarget_h

#include "DGtal/helpers/StdDefs.h"
#include "DGtal/base/ConstAlias.h"
#include "DGtal/shapes/Mesh.h"
#include "DGtal/shapes/IntersectionTarget.h"
#include "DGtal/kernel/SpaceND.h"
#include "DGtal/kernel/sets/CDigitalSet.h"
#include "DGtal/geometry/tools/determinant/PredicateFromOrientationFunctor2.h"
#include "DGtal/geometry/tools/determinant/InHalfPlaneBySimple3x3Matrix.h"

namespace DGtal
{


  template<typename TSpace, size_t TSeparation, size_t TDimension = 1>
  struct IntersectionTargetTrait
  {
    BOOST_STATIC_ASSERT_MSG( TSpace::dimension == 3, "DigitalSet dimension must be 3");

    BOOST_STATIC_ASSERT_MSG( TSeparation == 6 || TSeparation == 26, "Separation must be 6 or 26");

    BOOST_STATIC_ASSERT_MSG( TDimension == 1, "Currently only 1D intersection target is implemented");

    IntersectionTargetTrait() = delete;

    template<typename Space, size_t Separation, size_t Dimension = 1>
    struct IntersectionTarget;

    /****** Associated types *********************/
    using Space2D  = SpaceND<2, typename TSpace::Integer>;
    using PointR3  = typename TSpace::RealPoint;
    using VectorR3 = typename TSpace::RealPoint;
    using PointR2  = typename Space2D::RealPoint;
    using Type = IntersectionTarget<TSpace, TSeparation, TDimension>;
    /*********************************************/

    struct Edge
    {
      PointR3 myFirst, mySecond;
      Edge() = default;
      Edge(const Edge&) = default;
      Edge(const PointR3& mf, const PointR3& ms) : myFirst(mf), mySecond(ms) {}
    };

    template<typename Space>
    struct IntersectionTarget<Space, 6, 1>
    {
      // ----------------------- Members ------------------------------

      const std::array<Edge, 3> myTarget {{
          { { 0.5, 0.0, 0.0 }, { -0.5,  0.0,  0.0} } , // AB
          { { 0.0, 0.0, 0.5 }, {  0.0,  0.0, -0.5} } , // CD
          { { 0.0, 0.5, 0.0 }, {  0.0, -0.5,  0.0} }   // EF
        }};

      // ----------------------- Standard services --------------------

      const std::array<Edge, 3>& operator()() {
        return myTarget;
      }

      const Edge& operator()(int i) {
        return myTarget[i];
      }

      PointR2 project(int i, const PointR3& p) const {
        ASSERT( 0 <= i && i <= 2 );

        if(myTarget[i].myFirst[0] == 0.5)
          return { p[1], p[2] }; // ignore x
        else if(myTarget[i].myFirst[1] == 0.5)
          return { p[0], p[2] }; // ignore y
        else
          return { p[0], p[1] }; // ignore z
      }
    };

    template<typename Space>
    struct IntersectionTarget<Space, 26, 1>
    {
      // ----------------------- Members ------------------------------

      const std::array<Edge, 4> myTarget {{
          { { -0.5, 0.5,  0.5}, {  0.5, -0.5, -0.5} } , // AG
          { {  0.5, 0.5,  0.5}, { -0.5, -0.5, -0.5} } , // BH
          { {  0.5, 0.5, -0.5}, { -0.5, -0.5,  0.5} } , // CE
          { { -0.5, 0.5, -0.5}, {  0.5, -0.5,  0.5} }   // DF
        }};

      std::array<VectorR3, 4> myE1base;

      std::array<VectorR3, 4> myE2base;

      IntersectionTarget<Space, 26, 1>()
      {
        const double coef = 1. / std::sqrt(3.);

        std::array<VectorR3, 4> myNormal {{
            {  coef, -coef, -coef } ,
            { -coef, -coef, -coef } ,
            { -coef, -coef,  coef } ,
            {  coef, -coef,  coef }
          }};

        // computation of each e1/e2 bases foreach normal
        for(int i = 0; i < 4; i++)
        {
          auto& x = myNormal[i][0];
          auto& y = myNormal[i][1];
          auto& z = myNormal[i][2];

          myE1base[i] = std::move( VectorR3(-y, x, 0).getNormalized() );
          myE2base[i] = std::move( VectorR3(-x*z, -y*z, x*x + y*y).getNormalized() );
        }
      }

      // ----------------------- Standard services --------------------

      const std::array<Edge, 4>& operator()() const {
        return myTarget;
      }

      const Edge& operator()(int i) const {
        return myTarget[i];
      }

      PointR2 project(int i, const PointR3& p) const {
        ASSERT( 0 <= i && i <= 3 );

        return { myE1base[i].dot(p), myE2base[i].dot(p) };
      }
    };
  };
}

#endif // !defined MeshVoxelizer_h