geodesic_cpp_mod/geodesic_mesh_elements.h

//Copyright (C) 2008 Danil Kirsanov, MIT License
#ifndef GEODESIC_MESH_ELEMENTS_20071231
#define GEODESIC_MESH_ELEMENTS_20071231

// here we define the building elements of the mesh: 
// 3D-points, vertices, edges, faces, and surface points

#include <cstddef>

namespace geodesic{

class Vertex; 
class Edge; 
class Face; 
class Mesh; 
class MeshElementBase;

typedef Vertex* vertex_pointer;
typedef Edge* edge_pointer;
typedef Face* face_pointer;
typedef Mesh* mesh_pointer;
typedef MeshElementBase* base_pointer;

template <class Data>       //simple vector that stores info about mesh references 
class SimpleVector          //for efficiency, it uses an outside memory allocator
{
public:
    SimpleVector():
      m_size(0),
      m_begin(NULL)
    {};

    typedef Data* iterator;

    unsigned size(){return m_size;};
    iterator begin(){return m_begin;};
    iterator end(){return m_begin + m_size;};

    template<class DataPointer>
    void set_allocation(DataPointer begin, unsigned size)
    {
        FEASSERT(begin != NULL || size == 0);
        m_size = size;
        m_begin = (iterator)begin;
    }

    Data& operator[](unsigned i)
    {
        FEASSERT(i < m_size);
        return *(m_begin + i);
    }

    void clear()
    {
        m_size = 0;
        m_begin = NULL;
    }

private:
    unsigned m_size;
    Data* m_begin;
};

enum PointType
{
    VERTEX,
    EDGE,
    FACE,
    UNDEFINED_POINT
};

class MeshElementBase   //prototype of vertices, edges and faces
{
public:
    typedef SimpleVector<vertex_pointer> vertex_pointer_vector; 
    typedef SimpleVector<edge_pointer> edge_pointer_vector; 
    typedef SimpleVector<face_pointer> face_pointer_vector; 

    MeshElementBase():
        m_id(0),
        m_type(UNDEFINED_POINT)
    {};

    vertex_pointer_vector& adjacent_vertices(){return m_adjacent_vertices;};
    edge_pointer_vector& adjacent_edges(){return m_adjacent_edges;};
    face_pointer_vector& adjacent_faces(){return m_adjacent_faces;};

    unsigned& id(){return m_id;}; 
    PointType type(){return m_type;};

protected:
    vertex_pointer_vector m_adjacent_vertices;      //list of the adjacent vertices
    edge_pointer_vector m_adjacent_edges;           //list of the adjacent edges
    face_pointer_vector m_adjacent_faces;           //list of the adjacent faces

    unsigned m_id;                          //unique id
    PointType m_type;                           //vertex, edge or face
};

class Point3D           //point in 3D and corresponding operations
{
public:
    Point3D(){};
    Point3D(Point3D* p)
    {
        x() = p->x();
        y() = p->y();
        z() = p->z();
    };

    double* xyz(){return m_coordinates;};
    double& x(){return *m_coordinates;};
    double& y(){return *(m_coordinates+1);};
    double& z(){return *(m_coordinates+2);};

    void set(double new_x, double new_y, double new_z)
    {
        x() = new_x;
        y() = new_y;
        z() = new_z;
    }

    void set(double* data)
    {
        x() = *data;
        y() = *(data+1);
        z() = *(data+2);
    }

    double distance(double* v)
    {
        double dx = m_coordinates[0] - v[0];
        double dy = m_coordinates[1] - v[1];
        double dz = m_coordinates[2] - v[2];

        return sqrt(dx*dx + dy*dy + dz*dz);
    };

    double distance(Point3D* v)
    {
        return distance(v->xyz());
    };

    void add(Point3D* v)
    {
        x() += v->x();
        y() += v->y();
        z() += v->z();
    };

    void multiply(double v)
    {
        x() *= v;
        y() *= v;
        z() *= v;
    };

private:
    double m_coordinates[3];                    //xyz
};

class Vertex: public MeshElementBase, public Point3D
{
public:
    Vertex()
    {
        m_type = VERTEX;
    };

    ~Vertex(){};

    bool& saddle_or_boundary(){return m_saddle_or_boundary;};
private:                    
                                    //this flag speeds up exact geodesic algorithm
    bool m_saddle_or_boundary;      //it is true if total adjacent angle is larger than 2*PI or this vertex belongs to the mesh boundary
};


class Face: public MeshElementBase
{
public:
    Face()
    {
        m_type = FACE;
    };

    ~Face(){};

    edge_pointer opposite_edge(vertex_pointer v);
    vertex_pointer opposite_vertex(edge_pointer e);
    edge_pointer next_edge(edge_pointer e, vertex_pointer v);

    double vertex_angle(vertex_pointer v)
    {
        for(unsigned i=0; i<3; ++i)
        {
            if(adjacent_vertices()[i]->id() == v->id())
            {
                return m_corner_angles[i];
            }
        }
        FEASSERT(0);
        return 0;
    }

    double* corner_angles(){return m_corner_angles;};

private:
    double m_corner_angles[3];      //triangle angles in radians; angles correspond to vertices in m_adjacent_vertices
};

class Edge: public MeshElementBase
{
public:
    Edge()
    {
        m_type = EDGE;
    };

    ~Edge(){};

    double& length(){return m_length;};

    face_pointer opposite_face(face_pointer f)
    {
        if(adjacent_faces().size() == 1)
        {
            FEASSERT(adjacent_faces()[0]->id() == f->id());
            return NULL;
        }

        FEASSERT(adjacent_faces()[0]->id() == f->id() || 
               adjacent_faces()[1]->id() == f->id());

        return adjacent_faces()[0]->id() == f->id() ? 
               adjacent_faces()[1] : adjacent_faces()[0];
    };

    vertex_pointer opposite_vertex(vertex_pointer v)
    {
        FEASSERT(belongs(v));

        return adjacent_vertices()[0]->id() == v->id() ?
               adjacent_vertices()[1] : adjacent_vertices()[0];
    };

    bool belongs(vertex_pointer v)
    {
        return adjacent_vertices()[0]->id() == v->id() || 
               adjacent_vertices()[1]->id() == v->id();
    }

    bool is_boundary(){return adjacent_faces().size() == 1;};

    vertex_pointer v0(){return adjacent_vertices()[0];};
    vertex_pointer v1(){return adjacent_vertices()[1];};

    void local_coordinates(Point3D* point, 
                           double& x, 
                           double& y)
    {
        double d0 = point->distance(v0());
        if(d0 < 1e-50)
        {
            x = 0.0;
            y = 0.0;
            return;
        }

        double d1 = point->distance(v1());
        if(d1 < 1e-50)
        {
            x = m_length;
            y = 0.0;
            return;
        }

        x = m_length/2.0 + (d0*d0 - d1*d1)/(2.0*m_length);
        y = sqrt(std::max(0.0, d0*d0 - x*x));
        return;
    }

private:
    double m_length;                            //length of the edge
};

class SurfacePoint:public Point3D  //point on the surface of the mesh
{
public:
    SurfacePoint():
        m_p(NULL)
    {};

    SurfacePoint(vertex_pointer v):     //set the surface point in the vertex
        SurfacePoint::Point3D(v),
        m_p(v)
    {};

    SurfacePoint(face_pointer f):       //set the surface point in the center of the face
        m_p(f)
    {
        set(0,0,0);
        add(f->adjacent_vertices()[0]);
        add(f->adjacent_vertices()[1]);
        add(f->adjacent_vertices()[2]);
        multiply(1./3.);
    };

    SurfacePoint(edge_pointer e,        //set the surface point in the middle of the edge
                 double a = 0.5):       
        m_p(e)
    {
        double b = 1 - a;

        vertex_pointer v0 = e->adjacent_vertices()[0];
        vertex_pointer v1 = e->adjacent_vertices()[1];

        x() = b*v0->x() + a*v1->x();
        y() = b*v0->y() + a*v1->y();
        z() = b*v0->z() + a*v1->z();
    };

    SurfacePoint(base_pointer g, 
                 double x,
                 double y,
                 double z,
                 PointType t = UNDEFINED_POINT):
        m_p(g)
    {
        set(x,y,z);
    };

    void initialize(SurfacePoint const& p)
    {
        *this = p;
    }

    ~SurfacePoint(){};

    PointType type(){return m_p ? m_p->type() : UNDEFINED_POINT;};
    base_pointer& base_element(){return m_p;};
protected:
    base_pointer m_p;           //could be face, vertex or edge pointer
};

inline edge_pointer Face::opposite_edge(vertex_pointer v)
{
    for(unsigned i=0; i<3; ++i)
    {
        edge_pointer e = adjacent_edges()[i];
        if(!e->belongs(v))
        {
            return e;
        }
    }
    FEASSERT(0);
    return NULL;
}

inline vertex_pointer Face::opposite_vertex(edge_pointer e)
{
    for(unsigned i=0; i<3; ++i)
    {
        vertex_pointer v = adjacent_vertices()[i];
        if(!e->belongs(v))
        {
            return v;
        }
    }
    FEASSERT(0);
    return NULL;
}

inline edge_pointer Face::next_edge(edge_pointer e, vertex_pointer v)
{
    FEASSERT(e->belongs(v));

    for(unsigned i=0; i<3; ++i)
    {
        edge_pointer next = adjacent_edges()[i];
        if(e->id() != next->id() && next->belongs(v))
        {
            return next;
        }
    }
    FEASSERT(0);
    return NULL;
}

struct HalfEdge         //prototype of the edge; used for mesh construction
{
    unsigned face_id;
    unsigned vertex_0;      //adjacent vertices sorted by id value
    unsigned vertex_1;      //they are sorted, vertex_0 < vertex_1
};

inline bool operator < (const HalfEdge &x, const HalfEdge &y)
{
    if(x.vertex_0 == y.vertex_0)
    {
        return x.vertex_1 < y.vertex_1;
    }
    else
    {
        return x.vertex_0 < y.vertex_0;
    }
}

inline bool operator != (const HalfEdge &x, const HalfEdge &y)
{
    return x.vertex_0 != y.vertex_0 || x.vertex_1 != y.vertex_1;
}

inline bool operator == (const HalfEdge &x, const HalfEdge &y)
{
    return x.vertex_0 == y.vertex_0 && x.vertex_1 == y.vertex_1;
}

} //geodesic

#endif