SparseGridTopology.h

Go to the documentation of this file.

/******************************************************************************
*       SOFA, Simulation Open-Framework Architecture, version 1.0 beta 3      *
*                (c) 2006-2008 MGH, INRIA, USTL, UJF, CNRS                    *
*                                                                             *
* This library is free software; you can redistribute it and/or modify it     *
* under the terms of the GNU Lesser General Public License as published by    *
* the Free Software Foundation; either version 2.1 of the License, or (at     *
* your option) any later version.                                             *
*                                                                             *
* This library is distributed in the hope that it will be useful, but WITHOUT *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
* for more details.                                                           *
*                                                                             *
* You should have received a copy of the GNU Lesser General Public License    *
* along with this library; if not, write to the Free Software Foundation,     *
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.          *
*******************************************************************************
*                               SOFA :: Modules                               *
*                                                                             *
* Authors: The SOFA Team and external contributors (see Authors.txt)          *
*                                                                             *
* Contact information: contact@sofa-framework.org                             *
******************************************************************************/
#ifndef SOFA_COMPONENT_TOPOLOGY_SPARSEGRIDTOPOLOGY_H
#define SOFA_COMPONENT_TOPOLOGY_SPARSEGRIDTOPOLOGY_H

#include <string>


#include <sofa/component/topology/MeshTopology.h>
#include <sofa/helper/MarchingCubeUtility.h>
#include <sofa/defaulttype/Vec.h>
#include <sofa/component/topology/RegularGridTopology.h>

#include <sofa/helper/io/Mesh.h>
namespace sofa
{

    namespace component
    {

        namespace topology
        {

            using namespace sofa::defaulttype;
            using sofa::helper::MarchingCubeUtility;
            
            class SparseGridTopology : public MeshTopology
            {
                public:
                    
                    typedef fixed_array<Vector3,8> CubeCorners;
                    typedef enum{OUTSIDE,INSIDE,BOUNDARY} Type; 
                    
                    typedef Vec<3, int> Vec3i;
                    
                    SparseGridTopology(bool _isVirtual=false);
                    
//                  virtual void reinit() {updateMesh();};
                    
//                  static const float WEIGHT[8][8];
                    static const float WEIGHT27[8][27];
                    static const int cornerIndicesFromFineToCoarse[8][8];
                    
                    bool load(const char* filename);
                    virtual void init();
                    void buildAsFinest(); 
                    void buildFromFiner(); 
                    void buildVirtualFinerLevels(); 
                    
                    typedef std::map<Vector3,int> MapBetweenCornerPositionAndIndice;
                    
                    typedef std::vector<fixed_array<int,8> > HierarchicalCubeMap; 
                    HierarchicalCubeMap _hierarchicalCubeMap;
                    typedef helper::vector<int> InverseHierarchicalCubeMap; 
                    InverseHierarchicalCubeMap _inverseHierarchicalCubeMap;
                    
                    typedef std::map<int,float> AHierarchicalPointMap;
//                  typedef helper::vector< std::pair<int,float> >  AHierarchicalPointMap;
                    typedef helper::vector< AHierarchicalPointMap > HierarchicalPointMap; 
                    HierarchicalPointMap _hierarchicalPointMap;
                    typedef helper::vector< AHierarchicalPointMap > InverseHierarchicalPointMap; 
                    InverseHierarchicalPointMap _inverseHierarchicalPointMap;
                    typedef helper::vector< int > PointMap;
                    PointMap _pointMap; 
                    PointMap _inversePointMap;  
                    
                    
                    enum{UP,DOWN,RIGHT,LEFT,BEFORE,BEHIND,NUM_CONNECTED_NODES};
                    typedef helper::vector< helper::fixed_array<int,NUM_CONNECTED_NODES> > NodeAdjacency; 
                    NodeAdjacency _nodeAdjacency;
                    typedef helper::vector< helper::vector<int> >NodeCubesAdjacency; 
                    NodeCubesAdjacency _nodeCubesAdjacency;
                    typedef helper::vector< helper::vector<int> >NodeCornersAdjacency; 
                    NodeCornersAdjacency _nodeCornersAdjacency;
                    
                    
                    helper::vector<SparseGridTopology*> _virtualFinerLevels; 
                    
                    
                    Vec<3, int> getN() const { return n.getValue();}
                    int getNx() const { return n.getValue()[0]; }
                    int getNy() const { return n.getValue()[1]; }
                    int getNz() const { return n.getValue()[2]; }
                    
                    void setN(Vec3i _n) {n.setValue(_n);}
                    void setNx(int _n) { n.setValue(Vec3i(_n             ,n.getValue()[1],n.getValue()[2])); }
                    void setNy(int _n) { n.setValue(Vec3i(n.getValue()[0],_n             ,n.getValue()[2])); }
                    void setNz(int _n) { n.setValue(Vec3i(n.getValue()[0],n.getValue()[1],_n)             ); }
                    
                    int getNbVirtualFinerLevels() const { return _nbVirtualFinerLevels.getValue();}
                    void setNbVirtualFinerLevels(int n) {_nbVirtualFinerLevels.setValue(n);}
                                        
                    void setMin(Vector3 val) {_min.setValue(val);}                  
                    void setXmin(SReal val) { _min.setValue(Vector3(val             ,_min.getValue()[1],_min.getValue()[2])); }
                    void setYmin(SReal val) { _min.setValue(Vector3(_min.getValue()[0],val             ,_min.getValue()[2])); }
                    void setZmin(SReal val) { _min.setValue(Vector3(_min.getValue()[0],_min.getValue()[1],val)             ); }

                    
                    void setMax(Vector3 val) {_min.setValue(val);}
                    
                    void setXmax(SReal val) { _max.setValue(Vector3(val             ,_max.getValue()[1],_max.getValue()[2])); }
                    void setYmax(SReal val) { _max.setValue(Vector3(_max.getValue()[0],val             ,_max.getValue()[2])); }
                    void setZmax(SReal val) { _max.setValue(Vector3(_max.getValue()[0],_max.getValue()[1],val)             ); }
                    
                    Vector3 getMin() {return _min.getValue();}
                    SReal getXmin() { return _min.getValue()[0]; }
                    SReal getYmin() { return _min.getValue()[1]; }
                    SReal getZmin() { return _min.getValue()[2]; }
                    
                    Vector3 getMax() {return _max.getValue();}
                    SReal getXmax() { return _max.getValue()[0]; }
                    SReal getYmax() { return _max.getValue()[1]; }
                    SReal getZmax() { return _max.getValue()[2]; }
    
                    bool hasPos()  const{ return true; }
                    
                    virtual int findCube(const Vector3& pos, SReal& fx, SReal &fy, SReal &fz);
    
                    virtual int findNearestCube(const Vector3& pos, SReal& fx, SReal &fy, SReal &fz);
                    
                    virtual helper::fixed_array<int,6> findneighboorCubes( int indice );
                    
                    
                    virtual Type getType( int i );
                    
                    SparseGridTopology* getFinerSparseGrid() const {return _finerSparseGrid;}
                    void setFinerSparseGrid( SparseGridTopology* fsp ){_finerSparseGrid=fsp;}
                    SparseGridTopology* getCoarserSparseGrid() const {return _coarserSparseGrid;}
                    void setCoarserSparseGrid( SparseGridTopology* csp ){_coarserSparseGrid=csp;}
                    
                    void updateMesh();
                        
                    RegularGridTopology _regularGrid; 
                    vector< int > _indicesOfRegularCubeInSparseGrid; 
                    vector< int > _indicesOfCubeinRegularGrid; 
                    
                    Vector3 getPointPos( int i ){ return Vector3( seqPoints[i][0],seqPoints[i][1],seqPoints[i][2] ); }
                    
                    void getMesh( sofa::helper::io::Mesh &m);
                    
                    
                    void setDimVoxels( int a, int b, int c){ dataResolution.setValue(Vec3i(a,b,c));}
                    void setSizeVoxel( float a, float b, float c){ voxelSize.setValue(Vector3(a,b,c));}
                    
                    bool getVoxel(unsigned int x, unsigned int y, unsigned int z)
                    {                     
                      return getVoxel(dataResolution.getValue()[0]*dataResolution.getValue()[1]*z + dataResolution.getValue()[0]*y + x);
                    }
                       
                    bool getVoxel(unsigned int index) const
                    {
                      const int i = index%8;
                      unsigned char c = dataVoxels.getValue()[index>>3];
                      return ((c&((int)(pow(2.0f, i)))) >> i) == 1;
                    }; 
                    
                    
                protected:
                    bool isVirtual;
                    Data< Vec<3, int>    > n;
                    Data< Vector3 > _min;
                    Data< Vector3 > _max;
                    Data< int > _nbVirtualFinerLevels; 
                    
                    Data< Vec3i >           dataResolution;
                    Data< Vector3 >         voxelSize;
                    Data< unsigned int >    marchingCubeStep;
                    Data< unsigned int >    convolutionSize;
                    
                    virtual void updateEdges();
                    virtual void updateQuads();
                    virtual void updateHexas();
                    
                    MarchingCubeUtility                 marchingCubes;
                    Data< vector< unsigned char > >     dataVoxels; 
                    bool                                _usingMC;
                
                    sofa::helper::vector<Type> _types; 
                    
                    void propagateFrom( const int i, const int j, const int k,  
                                        RegularGridTopology& regularGrid, 
                                        vector<Type>& regularGridTypes, 
                                        vector<bool>& alreadyTested  ) const;
                    
                    void computeBoundingBox(const helper::vector<Vector3>& vertices,
                                                    SReal& xmin, SReal& xmax,
                                                    SReal& ymin, SReal& ymax,
                                                    SReal& zmin, SReal& zmax) const;

                    void voxelizeTriangleMesh(helper::io::Mesh* mesh,
                                              RegularGridTopology& regularGrid, 
                                              vector<Type>& regularGridTypes) const;
                    
                    void buildFromTriangleMesh(const std::string& filename);

                    void buildFromRegularGridTypes(RegularGridTopology& regularGrid, const vector<Type>& regularGridTypes);

                    void buildFromVoxelFile(const std::string& filename);
                    void buildFromRawVoxelFile(const std::string& filename);
                    
                    template< class T>
                    void constructCollisionModels(const sofa::helper::vector< sofa::core::componentmodel::topology::BaseMeshTopology * > &list_mesh,
                                      const sofa::helper::vector< sofa::helper::vector< Vec<3,T> >* >            &list_X) const;
                    
                    SparseGridTopology* _finerSparseGrid; 
                    SparseGridTopology* _coarserSparseGrid; 
                    
                    void setVoxel(int index, unsigned char value){
                      if (value)
                      {
                        (*dataVoxels.beginEdit())[index>>3] |= (int) pow(2.0f, (int)(index%8));
                      }
                      else
                      {
                        const int i = index%8;
                        const int mask = (int) pow(2.0f, i);
                        if (((*dataVoxels.beginEdit())[index>>3]&mask)>>i) (*dataVoxels.beginEdit())[index>>3] -= mask;
                      }
                       };
                       
                    
                    /*  /// to compute valid cubes (intersection between mesh segments and cubes)
                    typedef struct segmentForIntersection{
                        Vector3 center;
                        Vector3 dir;
                        SReal norm;
                        segmentForIntersection(const Vector3& s0, const Vector3& s1)
                        {
                            center = (s0+s1)*.5;
                            dir = center-s0;
                            norm = dir.norm();
                            dir /= norm;
                        };
                    } SegmentForIntersection;
                    struct ltSegmentForIntersection // for set of SegmentForIntersection
                    {
                        bool operator()(const SegmentForIntersection& s0, const SegmentForIntersection& s1) const
                        {
                            return s0.center < s1.center || s0.norm < s1.norm;
                        }
                    };
                    typedef struct cubeForIntersection{
                        Vector3 center;
                        fixed_array<Vector3,3> dir;
                        Vector3 norm;
                        cubeForIntersection( const CubeCorners&  corners )
                        {
                            center = (corners[7] + corners[0]) * .5;
                            
                            norm[0] = (center[0] - corners[0][0]);
                            dir[0] = Vector3(1,0,0);
                            
                            norm[1] = (center[1] - corners[0][1]);
                            dir[1] = Vector3(0,1,0);
                            
                            norm[2] = (center[2] - corners[0][2]);
                            dir[2] = Vector3(0,0,1);
                        }
                    } CubeForIntersection;
                    bool intersectionSegmentBox( const SegmentForIntersection& seg, const CubeForIntersection& cube  ); */
                    
                    bool _alreadyInit;
                    
                public :
                  
#ifdef SOFA_NEW_HEXA          
                    virtual const SeqHexas& getHexas()
                    {
                      if( !_alreadyInit ) init();
                      return sofa::component::topology::MeshTopology::getHexas();
                    }       
#else
                    virtual const SeqCubes& getHexas()
                    {
                      if( !_alreadyInit ) init();
                      return sofa::component::topology::MeshTopology::getHexas();
                    }       
#endif        
                    virtual int getNbPoints() const
                    {
                      if( !_alreadyInit ) const_cast<SparseGridTopology*>(this)->init();
                      return sofa::component::topology::MeshTopology::getNbPoints();
                    }
                    
                    virtual int getNbHexas() { return this->getHexas().size();}
                    
            };

        } // namespace topology

    } // namespace component

} // namespace sofa

#endif