LocalFromWorldPoint.hpp

#ifndef __LOCALFROMWORLDPOINT_HPP
#define __LOCALFROMWORLDPOINT_HPP
 
#include <fish/fiber/field/DirectProductArray.hpp>
#include <fish/fiber/vector/Interpolate.hpp>
#include <fish/fiber/vector/LinearIpol.hpp>
#include <fish/fiber/baseop/UniGridMapper.hpp>
#include <fish/fiber/baseop/TreeCellMapper.hpp>
#include <fish/fiber/baseop/ExpandBBox.hpp>
 
#include <fish/fiber/field/Fragment.hpp>
#include <fish/fiber/grid/Chart.hpp>
 
#include <fish/fiber/baseop/BundleBaseSpaceFunctions.hpp>
 
#include <eagle/PhysicalSpace.hpp>
#include <aerie/KDTree.hpp>
#include <eagle/FixedArray.hpp>
 
#include <memcore/RefPtr.hpp>
#include "FragmentSkeleton.hpp"
#include "PointSearch.hpp"
#include <fish/fiber/field/DirectProductArrayTypes.hpp>
 
#include "gridopDllApi.h"
 
#ifdef  USE_OPENMP
#include <omp.h>
#endif
 
namespace Fiber
{
 
using namespace Eagle;
using namespace Eagle::PhysicalSpace;
 
 
struct LocalPoint
{
Eagle::PhysicalSpace::point local_point;
const Eagle::PhysicalSpace::point world_point;
RefPtr<FragmentID>          frag_id;
MultiIndex<1>               cell; // required for hexahedral and pointcloud data, thus 1 dimensional
tvector                     cell_size;
 
        LocalPoint()
        : local_point(0, 0, 0)
        , world_point(0, 0, 0)
        , frag_id(NullPtr())
        , cell( MultiIndex<1>(0) )
        , cell_size(0, 0, 0)
        {}
 
        LocalPoint( const Eagle::PhysicalSpace::point wp )
        : local_point(0, 0, 0)
        , world_point( wp )
        , frag_id(NullPtr())
        , cell( MultiIndex<1>(0) )
        , cell_size(0, 0, 0)
        {}
 
};
 
class   gridop_API LocalFromWorldPoint : public MemCore::ReferenceBase<LocalFromWorldPoint >
{
        // info retrieved from iterator, which is called once
//      bool curvi_coords;
//      bool fragmented_coords;
 
//      bool inited;
 
//      RefPtr<Slice> slice;
        RefPtr<Field> coords;
 
        CartesianFragments      FragmentProperties;
        PointSearch             UniformPointSearch;
 
public:
        RefPtr<Grid>&grid() { return FragmentProperties.TheGrid; }
private:
 
        double TreeQueryFac;
        RefPtr<KDTree<3, int> > MyTree;
 
        static  inline MultiIndex<3> getFragmentOffset(const RefPtr<FragmentID>&f)
        {
         MultiIndex<3> Offset;
                 if ( RefPtr<FragmentLocation<3, FragmentProperty::Offset> >
                      F3 = interface_cast<FragmentLocation<3, FragmentProperty::Offset> > (f) )
                 {
                         Offset = *F3;
                 }
                 else
                 {
                         Offset = 0, 0, 0;
                 }
 
                 return Offset;
        }
 
        typedef MemArray<3, Eagle::tvector3>           VectorArray_t;
        typedef MemArray<3, double>                    ScalarArray_t;
        typedef MemArray<3, Eagle::point3>             CoordsArray_t;
        typedef MemArray<1, Eagle::point3>             CoordsArray1D_t;
        typedef MemArray<1, std::string>               FragmentIDs_t;
        typedef CartesianRectilinearArray3D     RectiLinearArray_t;
 
        typedef FixedArray<point,8>                    Vertexlist_t;
 
        typedef MemArray<1, RefPtr<BoundingBox> >      FragmentBounds_t;
        typedef DirectProductMemArray<Eagle::point3>   ProcArray_t;
 
        // Remember Fragment IDs and Bounding Boxes, and UniGRid object in class after initialization
//      RefPtr<FragmentIDs_t>    MyFragIds;
//      RefPtr<FragmentBounds_t> MyFragBBs;
        RefPtr<Field>&UniMapperField()
        {
                return UniformPointSearch.UniMapperField;
        }
 
        struct  MyIterator : Field::Iterator
        {
                int             c;
 
                RefPtr<Field> FragmentBoxes;
 
                RefPtr<Field>           Coords;
                RefPtr<KDTree<3, int> > Tree;
                double last_radius;
 
                bool build_tree;
 
//              RefPtr<FragmentIDs_t>    MyFragIds;
//              RefPtr<FragmentBounds_t> MyFragBBs;
 
                MyIterator(const RefPtr<Field>&FragmentBoxesPositions,
                           const RefPtr<Field>&CoordsP,
                           const RefPtr<KDTree<3,int> >&TreeP,
                           const bool build_treeP)
//, const RefPtr<FragmentIDs_t>&FragIdsP, const RefPtr<FragmentBounds_t>&FragBBsP )
                :c(0)
                , FragmentBoxes( FragmentBoxesPositions )
                , Coords(CoordsP)
                , Tree(TreeP)
                , last_radius(0)
                , build_tree(build_treeP)
//              , MyFragIds(FragIdsP)
//              , MyFragBBs(FragBBsP)
                {}
 
                bool apply(const RefPtr<FragmentID>&f, const MemCore::StrongPtr<Fiber::CreativeArrayBase>&DC);
 
        } ;
 
 
 
        class gridop_API BBSelectorBase
        {
        protected:
 
 
        public:
                BBSelectorBase() {}
                virtual ~BBSelectorBase(){}
 
                virtual unsigned int getNext() = 0;
                virtual bool hasNext() = 0;
                virtual void reset(const int _last, const point& pos) = 0;
        };
 
 public:
        class gridop_API BBSelectorMemoryTree : public BBSelectorBase
        {
                int last, current, size;
                point reset_pos;
                map<double, int>m;
                map<double, int>::iterator it;
                double TreeQueryFac;
                RefPtr<KDTree<3,int> > Tree;
        public:
                BBSelectorMemoryTree(const RefPtr<KDTree<3, int> >TreeP, const double TreeQueryFacP)
                : reset_pos(point(0,0,0))
                , TreeQueryFac(TreeQueryFacP)
                , Tree(TreeP)
                {
                        m.clear();
                        last=0;
                        current=0;
                        size=0;
                }
 
                unsigned int getNext();
                bool hasNext();
                void reset( const int _last, const point& pos );
        };
 
 
        class FindCell
        {
//      protected:
//              RefPtr<CoordsArray_t>          CoordsArr;
//              RefPtr<VectorArray_t>          VecsArr;
//              MultiIndex<3> iAx, iAy, iAz;
        public:
                FindCell()
                        {
//                              iAx = 1, 0, 0;
//                              iAy = 0, 1, 0;
//                              iAz = 0, 0, 1;
                        }
 
                virtual ~FindCell(){}
 
//              void setCoords(const RefPtr<CoordsArray_t>  CA)  { CoordsArr = CA; }
//              void setVecs  (const RefPtr<VectorArray_t>  VA)  { VecsArr   = VA; }
 
//              inline tvector cross(const tvector&t1, const tvector&t2)
//                      {
//                      tvector res;
//                              res =  t1[1]*t2[2] - t1[2]*t2[1],
//                                     -( t1[0]*t2[2] - t1[2]*t2[0] ),
//                                     t1[0]*t2[1] - t1[1]*t2[0];
//
//                              return res;
//                      }
 
                virtual int    pointInCell( const point&P ) = 0;
                virtual point  localUVW   () = 0;
        };
 
 
        class HexaHedralCellNewton : public FindCell
        {
                RefPtr<UniGridMapper> UniMap;
 
                point uvw;
                int point_type;
                double block_epsilon;
 
        public:
                HexaHedralCellNewton():FindCell()
                {}
 
                virtual ~HexaHedralCellNewton(){}
 
                void  setUniMap  ( RefPtr<UniGridMapper> unimap ) { UniMap = unimap; }
                int   pointInCell( const point&P );
                point localUVW   ();
 
                void setEpsilon(const double e) { block_epsilon = e; }
        };
 
        class ProjectiveCell// : public FindCell
        {
                RefPtr<TreeCellMapper> m_tree_map;
 
                // use LocalPoint here?? Reference??
                point uvw;
                tvector cell_size;
                int point_type;
 
        public:
                ProjectiveCell()//:FindCell()
                {}
 
                virtual ~ProjectiveCell(){}
 
                void  setTreeMap ( RefPtr<TreeCellMapper> tree_map ) { m_tree_map = tree_map; }
                inline point localUVW() { return uvw; }
                inline tvector cellSize() { return cell_size; }
 
                inline int pointInCell( const point&P )
                {
                        point_type = m_tree_map->localCoordinates( uvw, cell_size, P );
                        return point_type;
                }
 
 
        };
 
        // to be removed later
        double block_epsilon;
 
                // parameters for unigrid building
        double  res_scale, prec_scale;
 
 
 
public:
        int     MaxListSize;
private:
        tvector m_tree_query_scale;
public:
        LocalFromWorldPoint( Slice&SourceSliceP, const RefPtr<Grid>&SourceGridP,
                             const string&Gridname,
                             double     res_scale = 1.0, double prec_scale = 0.1,
                             tvector tree_query_scale = tvector(1.,1.,1.) );
 
        ~LocalFromWorldPoint();
 
//      bool get(const point&position, pair<point, RefPtr<FragmentID> >&result );
        bool get(const point&position, LocalPoint&result );
 
        bool getViaTree( LocalPoint&result );
 
        /* Set an epsilon value for widening the block boundary test. This is a hack solution and should be done internally via
         * neighbourhood information avoiding espsilons. Some tests showed that a value of 0.05 seems to do pretty ok.
         */
        void setBlockBoundaryEpsilon(const double block_epsilonP) { block_epsilon = block_epsilonP; }
 
 
        int numberOfBordersOfCell( const pair<point, string> localdata, int&hint1, int&hint2 );
 
 
        bool getBorderIntersectionNormal( const point worldposition,
                                          const pair<point, string> localpoint,
                                          const tvector direction,
                                          point&intersection, bivector&normal );
};
 
} // namespace
 
#endif