Isosurface.hpp

#ifndef __BASEOP_ISOSURFACE_HPP
#define __BASEOP_ISOSURFACE_HPP

#include "gridopDllApi.h"

#include <eagle/PhysicalSpace.hpp>

#include <vector/MultiIndex.hpp>

#include <field/Cell.hpp>

#include <grid/Grid.hpp>
#include <grid/CartesianChart.hpp>
#include <grid/Topology.hpp>

#include <bundle/Slice.hpp>

namespace Fiber
{

template <class ArrayType>
struct  FieldExplorer;

template <int Dims, class ValueType>
struct  FieldExplorer<MultiArray<Dims, ValueType> >
{
        const MultiArray<Dims, ValueType>&Field;
        const ValueType                   Isolevel;

        FieldExplorer(const MultiArray<Dims, ValueType>&FieldData, const ValueType&Threshold)
        : Field(FieldData)
        , Isolevel(Threshold)
        {}

        bool    isBelowThreshold(const MultiIndex<Dims>&I) const
        {
                return Field[ I ] < Isolevel;
        }

        template <class ResultType>
        ResultType InterpolationWeight(const MultiIndex<Dims>&FirstVertex, const MultiIndex<Dims>&SecondVertex,
                                       const ResultType&precision) const
        {
        const ValueType &P1 = Field[ FirstVertex  ],
                        &P2 = Field[ SecondVertex ]; 

                if ( fabs(Isolevel - P1) < precision)
                        return 0.0; 

                if ( fabs(Isolevel - P2) < precision)
                        return 1.0; 

                ValueType dP = P2-P1; 
                        if ( fabs(dP) < precision)
                                return 0.0; 

                return (Isolevel - P1) / dP;
        }

};

template <class ValueType, int Dims, class ResultType>
void    ComputeInterpolationEdgeAndWeight(index_t&EdgeID,
                                          int    &EdgeOrientation,
                                          ResultType&Weight,
                                          const FieldExplorer<MultiArray<Dims, ValueType> >&FE,
                                          const MultiIndex<Dims>&FirstVertex, const MultiIndex<Dims>&SecondVertex,
                                          const ResultType&precision)
{
MultiIndex<Dims> EdgeDir = SecondVertex - FirstVertex;
MultiIndex<Dims> IndexDifference = distance(SecondVertex, FirstVertex); 

        EdgeOrientation =  FirstVertex.getSignedOrientation(SecondVertex);
        assert( EdgeOrientation != 0 ); 

        if (EdgeOrientation<0)
        {
                EdgeOrientation = -EdgeOrientation;
                EdgeID = ComputeEdgeIDfromVertexAndOrientation(SecondVertex, EdgeOrientation-1, FE.Field.Size() ); 
                assert( EdgeID < NumberOfEdges( FE.Field.Size() ));
                Weight = 1.0 - FE.InterpolationWeight(FirstVertex, SecondVertex, precision);
        } 
        else
        {
                EdgeID = ComputeEdgeIDfromVertexAndOrientation(FirstVertex, EdgeOrientation -1, FE.Field.Size() ); 
                assert( EdgeID < NumberOfEdges( FE.Field.Size() ));
                Weight = FE.InterpolationWeight(FirstVertex, SecondVertex, precision);
        }
}

extern gridop_API const int edgeTable[256];
extern gridop_API const int triTable[256][16];

template <class ArrayType>
inline int      Polygonize(// data given per iso vertex
                           std::vector<index_t>       &IsoVertexEdge,
                           std::vector<double>        &IsoVertexWeight,

                           // data given per iso triangle
                           std::vector<TriangleCell>    &IsoSurfaceTriangles,

                           // internal collector
                           map<index_t, index_t>        &IsosurfaceVerticesAsEdges,

                           const FieldExplorer<ArrayType>&FE,

                           const MultiIndex<3>&CellIndex,
                           double  precision,
                           index_t EdgeOffset = 0)
{
typedef map<index_t, index_t>   IsosurfaceVerticesAsEdges_t;

        /*
          Determine the index into the edge table which
          tells us which vertices are inside of the surface
        */
int     cubeindex = 0;

MultiIndex<3> VertexIndices[1<<3];

static  const int VertexConvention[] =
{
0, // 0,0,0
1, // 1,0,0
3, // 0,1,0
2, // 1,1,0
4, // 0,0,1
5, // 1,0,1
7, // 1,1,0
6, // 1,1,1
}; 

//
// Build up mapping from vertex number (counted) to vertex index (computed) 
//
        for(int i=0; i< (1<<3); i++)
        {
        MultiIndex<3> I = CellIndex + MultiIndex<3>::BitIndex(i); 
                VertexIndices[ VertexConvention[i] ] = I;
        } 

// 
//  Compute edge bitmap 
//
        for(int i=0; i< (1<<3); i++)
        {
        MultiIndex<3> I = VertexIndices[i]; 

//              if ( ScalarField[ I ] < Isolevel )
                if ( FE.isBelowThreshold( I ) )
                {
                        cubeindex |= 1<<i;
                }
        } 
        /* Cube is entirely in/out of the surface */
        if (edgeTable[cubeindex] == 0)
                return 0; 

index_t Edges[12] = { -1u, -1u, -1u, -1u,
                      -1u, -1u, -1u, -1u,
                      -1u, -1u, -1u, -1u };

int     EdgeOrientations[12];
double  Weights[12];

static  int EdgeVertices[12][2] = {
{0,1},
{1,2},
{2,3},
{3,0},

{4,5},
{5,6},
{6,7},
{7,4},

{0,4},
{1,5},
{2,6},
{3,7}
};

int     Nints = 0; 

        // Compute Edges array 

        // 
        // Find the edges where the surface intersects the cube 
        //
        for(int E=0; E<12; E++)
        {
                if (edgeTable[cubeindex] & (1<<E) )
                {
                int edge = E; 

                MultiIndex<3> V0 = VertexIndices[ EdgeVertices[edge][0] ],
                              V1 = VertexIndices[ EdgeVertices[edge][1] ];

                        Nints++;

                        ComputeInterpolationEdgeAndWeight(Edges[edge], EdgeOrientations[edge], Weights[edge],
                                                          FE,
                                                          V0, V1,
                                                          precision);
                }
        } 

        // 
        // Create and the triangles
        //
int     ntriang = 0;
        for (int i=0;triTable[cubeindex][i]!=-1; i+=3) 
        {
                // Indices in locale Edge/Weights table
        int     I[3] = { triTable[cubeindex][i  ],
                         triTable[cubeindex][i+1],
                         triTable[cubeindex][i+2] };

        TriangleCell Triangle;
                for(int k=0; k<3; k++)
                {
                int     LocalEdge = I[k];
                        assert( LocalEdge <12 );
                index_t edge = Edges[ LocalEdge ] + EdgeOffset;

                index_t SurfaceVertexIndex;
        
                IsosurfaceVerticesAsEdges_t::const_iterator it =
                        IsosurfaceVerticesAsEdges.find( edge ); 
                        if (it == IsosurfaceVerticesAsEdges.end() ) 
                        {
                                SurfaceVertexIndex = IsoVertexEdge.size();
                                IsosurfaceVerticesAsEdges[ edge ] = SurfaceVertexIndex;

                                IsoVertexEdge           .push_back( edge );

                                assert(EdgeOrientations[ LocalEdge ] > 0);
                                IsoVertexWeight         .push_back( Weights[ LocalEdge ] ); 
                        }
                        else
                        {
                                SurfaceVertexIndex = it->second;
                        }
                        Triangle[ k ] = SurfaceVertexIndex;
                }
                IsoSurfaceTriangles.push_back( Triangle );

                ntriang++;
        }
        return ntriang;
}

template <class VertexFieldType>
inline void     EvalVertexFieldOnEdges(std::vector<typename VertexFieldType::value_type>&EvaluatedField,
                                       // data given per vertex
                                       const std::vector<index_t>    &Edges,
                                       const std::vector<double>     &Weights,

                                       // original vertex field
                                       const VertexFieldType&RegularField,
                                       const MultiIndex<3>&FragmentSize,
                                       const MultiIndex<3>&FragmentOffset,

                                       index_t  StartEvaluationAtVertex = 0,
                                       index_t  EdgeOffset = 0)
{
        enum { Dims = VertexFieldType::Dims };
        typedef typename VertexFieldType::value_type value_type;

        EvaluatedField.resize( Edges.size() );

        for(index_t i = StartEvaluationAtVertex; i<Edges.size(); i++)
        {
        index_t EdgeIndex = Edges[i] - EdgeOffset;

        std::pair<MultiIndex<3>, MultiIndex<3> > Edge =
                ComputeEdgeVertices(EdgeIndex, FragmentSize );

        const   value_type&P0 = RegularField[ Edge.first  + FragmentOffset ], 
                          &P1 = RegularField[ Edge.second + FragmentOffset ];

        double  w = Weights[ i ];

                EvaluatedField[ i ] = P0 + w*(P1-P0);
        }
}



typedef std::pair<RefPtr<MemArray<1, ::Eagle::PhysicalSpace::point> >, RefPtr<MemArray<1, TriangleCell> > > TriangleSurface_t;

template <class VertexField_t, class ArrayType>
inline gridop_API TriangleSurface_t
        ComputeIsosurfaceT(const FieldExplorer<ArrayType>&FE,
                           const VertexField_t&VertexField,
                           double precision)
{
        // data given per iso vertex
MemCore::MemVector<index_t>             IsoVertexEdge(0);
MemCore::MemVector<double>              IsoVertexWeight(0);

        // data given per iso triangle 
MemCore::MemVector<TriangleCell>        IsoSurfaceTriangles(0);

        // internal collector
map<index_t, index_t>                   IsosurfaceVerticesAsEdges;

MultiIndex<3> NumberOfCells = FE.Field.Size() - MIndex(1,1,1);
MultiIndex<3> CellIndex;

        try
        {
                do
                {
                        Polygonize(IsoVertexEdge,
                                   IsoVertexWeight,
                                   IsoSurfaceTriangles,
        
                                   IsosurfaceVerticesAsEdges,
        
                                   FE,
                                   CellIndex,
                                   precision);
                }
                while(  CellIndex.inc( NumberOfCells ) );
        }
        catch(...)
        {
                printf("ComputeIsosurfaceT(): INTERNAL ERROR\n"); fflush(stdout);
                throw;
        } 

MemBase::Creator_t Crec;

RefPtr<MemBase>
        IsoVertexEdgeM = makeMemArray1D( IsoVertexEdge  , Crec ),
        IsoVertexEdgeW = makeMemArray1D( IsoVertexWeight, Crec ),

        IsoTriangleM   = makeMemArray1D( IsoSurfaceTriangles, Crec ); 


MemCore::MemVector< ::Eagle::PhysicalSpace::point>
        IsoVertices( IsoVertexEdge.size() ); 

//      ComputeVertices 

        EvalVertexFieldOnEdges(IsoVertices.std_vector(),
                        // data given per iso vertex
                               IsoVertexEdge.std_vector(),
                               IsoVertexWeight.std_vector(),
                               
                        // original vertex field
                               VertexField,
                               VertexField.Size(),
                               MIndex(0,0,0)      );

        return TriangleSurface_t( makeMemArray1D( IsoVertices, Crec ),
                                  IsoTriangleM);
}



extern gridop_API
TriangleSurface_t ComputeIsosurface(double Isolevel,
                                    const MultiArray<3,double>&ScalarField,
                                    const MultiArray<3,::Eagle::PhysicalSpace::point>&VertexField,
                                    double precision);


} // namespace Fiber

#endif