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/RegularTopology.hpp>
 
#include <bundle/Slice.hpp>
 
namespace Fiber
{
 
template <class ArrayType>
struct  FieldExplorer;
 
template <Dims_t 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)
                        if ( dP == 0.0)
                                return 0.5;
 
        ResultType r = (Isolevel - P1) / dP;
                if (r<0.0) return 0.0;
                if (r>1.0) return 1.0;
 
                return r;
        }
 
};
 
template <class ValueType, Dims_t 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_t { 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);
 
        try
        {
                for(MultiIndex<3> CellIndex : NumberOfCells )
                {
                        Polygonize(IsoVertexEdge,
                                   IsoVertexWeight,
                                   IsoSurfaceTriangles,
 
                                   IsosurfaceVerticesAsEdges,
 
                                   FE,
                                   CellIndex,
                                   precision);
                }
        }
        catch(...)
        {
                printf("ComputeIsosurfaceT(): INTERNAL ERROR\n"); fflush(stdout);
                throw;
        }
 
RefPtr<MemBase>
        IsoVertexEdgeM = makeMemArray1D( IsoVertexEdge  ),
        IsoVertexEdgeW = makeMemArray1D( IsoVertexWeight ),
 
        IsoTriangleM   = makeMemArray1D( IsoSurfaceTriangles );
 
 
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 ),
                                  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