marchingCube.h

00001 /* This file is part of the Palabos library.
00002  *
00003  * Copyright (C) 2011 FlowKit Sarl
00004  * Avenue de Chailly 23
00005  * 1012 Lausanne, Switzerland
00006  * E-mail contact: contact@flowkit.com
00007  *
00008  * The most recent release of Palabos can be downloaded at 
00009  * <http://www.palabos.org/>
00010  *
00011  * The library Palabos is free software: you can redistribute it and/or
00012  * modify it under the terms of the GNU Affero General Public License as
00013  * published by the Free Software Foundation, either version 3 of the
00014  * License, or (at your option) any later version.
00015  *
00016  * The library is distributed in the hope that it will be useful,
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  * GNU Affero General Public License for more details.
00020  *
00021  * You should have received a copy of the GNU Affero General Public License
00022  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00023 */
00024 
00025 #ifndef MARCHING_CUBE_H
00026 #define MARCHING_CUBE_H
00027 
00028 #include "core/globalDefs.h"
00029 #include "offLattice/boundaryShapes3D.h"
00030 #include "offLattice/triangleSet.h"
00031 #include "offLattice/triangularSurfaceMesh.h"
00032 #include "offLattice/offLatticeBoundaryProfiles3D.h"
00033 #include "offLattice/triangleBoundary3D.h"
00034 #include <vector>
00035 
00036 namespace plb {
00037 
00038 template<typename T> class IsoSurfaceDefinition3D;
00039 
00040 
00042 
00049 template<typename T>
00050 void isoSurfaceMarchingCube (
00051         std::vector<typename TriangleSet<T>::Triangle>& triangles,
00052         std::vector<MultiBlock3D*> surfDefinitionArgs,
00053         IsoSurfaceDefinition3D<T>* isoSurfaceDefinition,
00054         Box3D const& domain, std::vector<plint> surfaceIds = std::vector<plint>() );
00055 
00057 template<typename T>
00058 void isoSurfaceMarchingCube (
00059         std::vector<typename TriangleSet<T>::Triangle>& triangles,
00060         VoxelizedDomain3D<T>& voxelizedDomain, Box3D const& domain );
00061 
00063 template<typename T>
00064 void isoSurfaceMarchingCube (
00065         std::vector<typename TriangleSet<T>::Triangle>& triangles,
00066         MultiScalarField3D<T>& scalarField, std::vector<T> const& isoLevels, Box3D const& domain );
00067 
00068 
00069 template<typename T>
00070 class IsoSurfaceDefinition3D {
00071 public:
00072     virtual ~IsoSurfaceDefinition3D() { }
00073     virtual bool isInside (
00074             plint surfaceId, Array<plint,3> const& position ) const =0;
00075     virtual bool isValid(Array<plint,3> const& position) const { return true; }
00076     virtual Array<T,3> getSurfacePosition (
00077              plint surfaceId, Array<plint,3> const& p1, Array<plint,3> const& p2 ) const =0;
00078     virtual void setArguments(std::vector<AtomicBlock3D*> const& arguments) =0;
00079     virtual IsoSurfaceDefinition3D<T>* clone() const =0;
00080     virtual plint getNumArgs() const =0;
00081     virtual std::vector<plint> getSurfaceIds() const=0;
00082 public:
00083     bool edgeIsValid(plint iX, plint iY, plint iZ, int edge) const;
00084 };
00085 
00086 template<typename T>
00087 class ScalarFieldIsoSurface3D : public IsoSurfaceDefinition3D<T> {
00088 public:
00089     ScalarFieldIsoSurface3D(std::vector<T> const& isoValues_);
00090     virtual bool isInside (
00091             plint surfaceId, Array<plint,3> const& position ) const;
00092     virtual Array<T,3> getSurfacePosition (
00093             plint surfaceId, Array<plint,3> const& p1, Array<plint,3> const& p2 ) const;
00094     virtual void setArguments(std::vector<AtomicBlock3D*> const& arguments);
00095     virtual ScalarFieldIsoSurface3D<T>* clone() const;
00096     virtual plint getNumArgs() const { return 1; }
00097     virtual std::vector<plint> getSurfaceIds() const;
00098 private:
00099     std::vector<T> isoValues;
00100     ScalarField3D<T>* scalar;
00101     Dot3D location;
00102 };
00103 
00104 
00105 
00106 template<typename T, class SurfaceData>
00107 class BoundaryShapeIsoSurface3D : public IsoSurfaceDefinition3D<T> {
00108 public:
00109     BoundaryShapeIsoSurface3D(BoundaryShape3D<T,SurfaceData>* shape_);
00110     ~BoundaryShapeIsoSurface3D();
00111     BoundaryShapeIsoSurface3D(BoundaryShapeIsoSurface3D<T,SurfaceData> const& rhs);
00112     BoundaryShapeIsoSurface3D<T,SurfaceData>& operator=(BoundaryShapeIsoSurface3D<T,SurfaceData> const& rhs);
00113     void swap(BoundaryShapeIsoSurface3D<T,SurfaceData>& rhs);
00114     virtual bool isInside (
00115             plint surfaceId, Array<plint,3> const& position ) const;
00116     virtual Array<T,3> getSurfacePosition (
00117             plint surfaceId, Array<plint,3> const& p1, Array<plint,3> const& p2 ) const;
00122     virtual void setArguments(std::vector<AtomicBlock3D*> const& arguments);
00123     virtual BoundaryShapeIsoSurface3D<T,SurfaceData>* clone() const;
00124     virtual plint getNumArgs() const { return 3; }
00125     virtual std::vector<plint> getSurfaceIds() const;
00126 private:
00127     BoundaryShape3D<T,SurfaceData>* shape;
00128 };
00129 
00130 
00131 template<typename T>
00132 class MarchingCubeSurfaces3D : public BoxProcessingFunctional3D {
00133 public:
00134     typedef typename TriangleSet<T>::Triangle Triangle;
00135 public:
00136     MarchingCubeSurfaces3D( std::vector<plint> surfaceIds_,
00137                             IsoSurfaceDefinition3D<T>* isoSurface_,
00138                             bool edgeOrientedData_=false );
00139     ~MarchingCubeSurfaces3D();
00140     MarchingCubeSurfaces3D(MarchingCubeSurfaces3D<T> const& rhs);
00141     MarchingCubeSurfaces3D<T>& operator=(MarchingCubeSurfaces3D<T> const& rhs);
00142     void swap(MarchingCubeSurfaces3D<T>& rhs);
00143     virtual void processGenericBlocks (
00144                 Box3D domain, std::vector<AtomicBlock3D*> fields );
00145     virtual void defaultImplementation (
00146                 Box3D domain, AtomicContainerBlock3D* triangleContainer );
00147     virtual void edgeOriented (
00148                 Box3D domain, AtomicContainerBlock3D* triangleContainer );
00149     virtual MarchingCubeSurfaces3D<T>* clone() const;
00150     virtual void getTypeOfModification(std::vector<modif::ModifT>& modified) const;
00151     void setEdgeOrientedEnvelope(plint edgeOrientedEnvelope_) {
00152         edgeOrientedEnvelope = edgeOrientedEnvelope_;
00153     }
00154 public:
00155     class TriangleSetData : public ContainerBlockData {
00156     public:
00157         std::vector<Triangle> triangles;
00158         virtual TriangleSetData* clone() const {
00159             return new TriangleSetData(*this);
00160         }
00161     };
00162     class EdgeOrientedTriangleSetData : public ContainerBlockData {
00163     public:
00166         struct OnEdgeTriangle {
00169             Array<plint,4> vertex1, vertex2;
00170         };
00171     public:
00172         EdgeOrientedTriangleSetData(plint nx, plint ny, plint nz)
00173             : data(nx,ny,nz)
00174         { }
00175         virtual EdgeOrientedTriangleSetData* clone() const {
00176             return new EdgeOrientedTriangleSetData(*this);
00177         }
00178         void addTriangle( plint iX, plint iY, plint iZ, int iEdge,
00179                           OnEdgeTriangle const& triangle )
00180         {
00181             switch(iEdge) {
00182                 case 0: data.get(iX,iY,iZ).edge1triangles.push_back(triangle); break;
00183                 case 1: data.get(iX,iY,iZ).edge2triangles.push_back(triangle); break;
00184                 case 2: data.get(iX,iY,iZ).edge3triangles.push_back(triangle); break;
00185                 default: PLB_ASSERT( false );
00186             }
00187         }
00188         std::vector<OnEdgeTriangle> const& getTriangles(plint iX, plint iY, plint iZ, int iEdge) const {
00189             switch(iEdge) {
00190                 case 0: return data.get(iX,iY,iZ).edge1triangles;
00191                 case 1: return data.get(iX,iY,iZ).edge2triangles;
00192                 case 2: return data.get(iX,iY,iZ).edge3triangles;
00193                 default: PLB_ASSERT( false );
00194             }
00195         }
00196         void setVertex(plint iX, plint iY, plint iZ, int iEdge, Array<T,3> const& vertex) {
00197             switch(iEdge) {
00198                 case 0:
00199                     data.get(iX,iY,iZ).edge1Vertex = vertex;
00200                     data.get(iX,iY,iZ).edge1VertexDefined = true;
00201                     break;
00202                 case 1:
00203                     data.get(iX,iY,iZ).edge2Vertex = vertex;
00204                     data.get(iX,iY,iZ).edge2VertexDefined = true;
00205                     break;
00206                 case 2:
00207                     data.get(iX,iY,iZ).edge3Vertex = vertex;
00208                     data.get(iX,iY,iZ).edge3VertexDefined = true;
00209                     break;
00210                 default: PLB_ASSERT( false );
00211             }
00212         }
00213         bool getVertex(plint iX, plint iY, plint iZ, int iEdge, Array<T,3>& vertex) const {
00214             switch(iEdge) {
00215                 case 0:
00216                     vertex = data.get(iX,iY,iZ).edge1Vertex;
00217                     return data.get(iX,iY,iZ).edge1VertexDefined;
00218                 case 1:
00219                     vertex = data.get(iX,iY,iZ).edge2Vertex;
00220                     return data.get(iX,iY,iZ).edge2VertexDefined;
00221                 case 2:
00222                     vertex = data.get(iX,iY,iZ).edge3Vertex;
00223                     return data.get(iX,iY,iZ).edge3VertexDefined;
00224                 default: PLB_ASSERT( false );
00225             }
00226         }
00227         std::vector<T> const& getScalars(plint iX, plint iY, plint iZ, int iEdge) const {
00228             switch(iEdge) {
00229                 case 0: return data.get(iX,iY,iZ).scalars1;
00230                 case 1: return data.get(iX,iY,iZ).scalars2;
00231                 case 2: return data.get(iX,iY,iZ).scalars3;
00232                 default: PLB_ASSERT( false );
00233             }
00234         }
00235         std::vector<T>& getScalars(plint iX, plint iY, plint iZ, int iEdge) {
00236             switch(iEdge) {
00237                 case 0: return data.get(iX,iY,iZ).scalars1;
00238                 case 1: return data.get(iX,iY,iZ).scalars2;
00239                 case 2: return data.get(iX,iY,iZ).scalars3;
00240                 default: PLB_ASSERT( false );
00241             }
00242         }
00243         bool isEdgeVertexDefined(plint iX, plint iY, plint iZ, int iEdge) {
00244             switch(iEdge) {
00245                 case 0: return data.get(iX,iY,iZ).edge1VertexDefined;
00246                 case 1: return data.get(iX,iY,iZ).edge2VertexDefined;
00247                 case 2: return data.get(iX,iY,iZ).edge3VertexDefined;
00248                 default: PLB_ASSERT( false );
00249             }
00250         }
00251         plint getNx() const { return data.getNx(); }
00252         plint getNy() const { return data.getNy(); }
00253         plint getNz() const { return data.getNz(); }
00254         Box3D getBoundingBox() const { return data.getBoundingBox(); }
00255     private:
00256         struct EdgeData {
00257             EdgeData()
00258                 : edge1VertexDefined(false),
00259                   edge2VertexDefined(false),
00260                   edge3VertexDefined(false)
00261             { }
00262             // Each element in the three following vectors defines the topology of
00263             // a triangle.
00264             std::vector<OnEdgeTriangle> edge1triangles;
00265             std::vector<OnEdgeTriangle> edge2triangles;
00266             std::vector<OnEdgeTriangle> edge3triangles;
00267             // Every edge has at most one vertex which is shared by all
00268             // triangles that have a vertex on this edge.
00269             Array<T,3> edge1Vertex, edge2Vertex, edge3Vertex;
00270             std::vector<T> scalars1, scalars2, scalars3;
00271             bool edge1VertexDefined, edge2VertexDefined, edge3VertexDefined;
00272         };
00273     private:
00274         ScalarField3D<EdgeData> data;
00275     };
00276 private:
00277     void marchingCubeImpl (
00278              plint iX, plint iY, plint iZ, plint surfaceId,
00279              std::vector<Triangle>& triangles,
00280              int& cubeIndex, std::vector<Array<T,3> >& vertlist );
00281     void polygonize (
00282              plint iX, plint iY, plint iZ, plint surfaceId,
00283              std::vector<Triangle>& triangles );
00286     void polygonize (
00287              plint iX, plint iY, plint iZ, plint surfaceId,
00288              std::vector<Triangle>& triangles,
00289              std::vector<Array<plint,4> >& edgeAttributions );
00290 private:
00291     std::vector<plint> surfaceIds;
00292     IsoSurfaceDefinition3D<T>* isoSurface;
00293     bool edgeOrientedData;
00294     plint edgeOrientedEnvelope;
00295 };
00296 
00297 struct MarchingCubeConstants {
00298     static const int edgeTable[256];
00299     static const int triTable[256][16];
00300     static const int edgeNeighb[12][3];
00301     static const int edgeOrient[12];
00302 };
00303 
00304 }  // namespace plb
00305 
00306 #endif  // MARCHING_CUBE_H