mirror of https://github.com/CGAL/cgal
This commit is contained in:
Marc Pouget
parent
8d97d550a0
commit
cb1977cedf
|
|
@ -1254,7 +1254,10 @@ Ridges_3/doc_tex/Ridges_3/ellipsoid_ridges.pdf -text
|
|||
Ridges_3/doc_tex/Ridges_3/ellipsoid_ridges_mesh.eps -text
|
||||
Ridges_3/doc_tex/Ridges_3/ellipsoid_ridges_mesh.jpg -text
|
||||
Ridges_3/doc_tex/Ridges_3/ellipsoid_ridges_mesh.pdf -text
|
||||
Ridges_3/include/CGAL/algo.C -text
|
||||
Ridges_3/examples/Ridges_3/PolyhedralSurf.C -text
|
||||
Ridges_3/examples/Ridges_3/README -text
|
||||
Ridges_3/examples/Ridges_3/blind.C -text
|
||||
Ridges_3/examples/Ridges_3/options.C -text
|
||||
Robustness/demo/Robustness/robustness.vcproj -text
|
||||
SearchStructures/doc_tex/SearchStructures/d-range.eps -text
|
||||
SearchStructures/doc_tex/SearchStructures/d-range.gif -text svneol=unset#unset
|
||||
|
|
|
|||
|
|
@ -0,0 +1,109 @@
|
|||
#ifndef _GSL_H_
|
||||
#define _GSL_H_
|
||||
|
||||
#include <gsl/gsl_linalg.h>
|
||||
#include <gsl/gsl_eigen.h>
|
||||
|
||||
////////////////////////class gsl_Vector/////////////////////
|
||||
class gsl_Vector{
|
||||
protected:
|
||||
gsl_vector* m_vector;
|
||||
public:
|
||||
//contructor
|
||||
// initializes all the elements of the vector to zero
|
||||
gsl_Vector(size_t n) { m_vector = gsl_vector_calloc(n); }
|
||||
//access
|
||||
double const operator[](size_t i) const { return gsl_vector_get(m_vector,
|
||||
i); }
|
||||
//use as left value v[i]=10
|
||||
double& operator[](size_t i) { return *gsl_vector_ptr(m_vector, i);}
|
||||
const gsl_vector* vector() const { return m_vector; }
|
||||
gsl_vector* vector() { return m_vector; }
|
||||
};
|
||||
|
||||
////////////////////////class gsl_Matrix/////////////////////
|
||||
class gsl_Matrix{
|
||||
protected:
|
||||
gsl_matrix* m_matrix;
|
||||
public:
|
||||
//contructor
|
||||
// initializes all the elements of the matrix to zero.
|
||||
gsl_Matrix(size_t n1, size_t n2) { m_matrix = gsl_matrix_calloc (n1, n2); }
|
||||
|
||||
//class Row, to define the usual double operator [][] for matrices
|
||||
class Row{
|
||||
gsl_matrix* matrix;
|
||||
size_t row;
|
||||
public:
|
||||
Row(gsl_matrix* _matrix, size_t _row) :
|
||||
matrix(_matrix), row(_row) {}
|
||||
double const operator[](size_t column) const
|
||||
{ return gsl_matrix_get(this->matrix, row, column);}
|
||||
double& operator[](size_t column)
|
||||
{ return *gsl_matrix_ptr(this->matrix, row, column);}
|
||||
};//END class Row
|
||||
|
||||
Row operator[](size_t _row) { return Row(m_matrix, _row);}
|
||||
|
||||
//access
|
||||
const gsl_matrix* matrix() const { return m_matrix; }
|
||||
gsl_matrix* matrix() { return m_matrix; }
|
||||
};
|
||||
|
||||
////////////////////////class GSL/////////////////////
|
||||
class GSL{
|
||||
public:
|
||||
typedef double FT;
|
||||
typedef gsl_Vector Vector;
|
||||
typedef gsl_Matrix Matrix;
|
||||
// solve for eigenvalues and eigenvectors.
|
||||
// eigen values are sorted in descending order,
|
||||
// eigen vectors are sorted in accordance.
|
||||
static
|
||||
void eigen_symm_algo(Matrix& S, Vector& eval, Matrix& evec);
|
||||
//solve MX=B using SVD and give the condition number of M
|
||||
//M replaced by U after gsl_linalg_SV_decomp()
|
||||
//The diagonal and lower triangular part of M are destroyed during the
|
||||
//computation, but the strict upper triangular part is not referenced.
|
||||
static
|
||||
void solve_ls_svd_algo(Matrix& M, Vector& X, const Vector& B,
|
||||
double &cond_nb);
|
||||
};
|
||||
|
||||
void GSL::eigen_symm_algo(Matrix& S, Vector& eval, Matrix& evec)
|
||||
{
|
||||
const size_t common_size = S.matrix()->size1;
|
||||
assert( S.matrix()->size2 == common_size
|
||||
&& eval.vector()->size == common_size
|
||||
&& evec.matrix()->size1 == common_size
|
||||
&& evec.matrix()->size2 == common_size
|
||||
);
|
||||
gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (common_size);
|
||||
gsl_eigen_symmv (S.matrix(), eval.vector(), evec.matrix(), w);
|
||||
gsl_eigen_symmv_free (w);
|
||||
gsl_eigen_symmv_sort (eval.vector(), evec.matrix(),
|
||||
GSL_EIGEN_SORT_VAL_DESC);
|
||||
}
|
||||
|
||||
void GSL::solve_ls_svd_algo(Matrix& M, Vector& X, const Vector& B,
|
||||
double &cond_nb)
|
||||
{
|
||||
const size_t nb_lines = M.matrix()->size1;
|
||||
assert( B.vector()->size == nb_lines );
|
||||
const size_t nb_columns = M.matrix()->size2;
|
||||
assert( X.vector()->size == nb_columns );
|
||||
|
||||
gsl_matrix * V = gsl_matrix_alloc(nb_columns,nb_columns);
|
||||
gsl_vector * S = gsl_vector_alloc(nb_columns);
|
||||
gsl_vector * work = gsl_vector_alloc(nb_columns);
|
||||
|
||||
gsl_linalg_SV_decomp (M.matrix(), V, S, work);
|
||||
gsl_linalg_SV_solve (M.matrix(), V, S, B.vector(), X.vector());
|
||||
|
||||
cond_nb = gsl_vector_get(S,0)/gsl_vector_get(S,nb_columns-1);
|
||||
gsl_matrix_free(V);
|
||||
gsl_vector_free(S);
|
||||
gsl_vector_free(work);
|
||||
}
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,63 @@
|
|||
#include "PolyhedralSurf.h"
|
||||
|
||||
Vector_3 PolyhedralSurf::getHalfedge_vector(Halfedge * h)
|
||||
{
|
||||
Vector_3 v = h->opposite()->vertex()->point() - h->vertex()->point();
|
||||
return v;
|
||||
}
|
||||
|
||||
///???????myterious behavior????????
|
||||
void PolyhedralSurf::compute_facets_normals()
|
||||
{
|
||||
std::for_each(this->facets_begin(), this->facets_end(),
|
||||
Facet_unit_normal());
|
||||
}
|
||||
|
||||
void PolyhedralSurf::compute_edges_length()
|
||||
{
|
||||
std::for_each(this->halfedges_begin(), this->halfedges_end(),
|
||||
Edge_length());
|
||||
}
|
||||
|
||||
double PolyhedralSurf::
|
||||
compute_mean_edges_length_around_vertex(Vertex* v)
|
||||
{
|
||||
Halfedge_around_vertex_circulator
|
||||
hedgeb = v->vertex_begin(), hedgee = hedgeb;
|
||||
int count_he = 0;
|
||||
double sum = 0.;
|
||||
CGAL_For_all(hedgeb, hedgee)
|
||||
{
|
||||
sum += hedgeb->getLength();
|
||||
count_he++;
|
||||
}
|
||||
return sum/count_he;
|
||||
}
|
||||
|
||||
Vector_3 PolyhedralSurf::computeFacetsAverageUnitNormal(Vertex * v)
|
||||
{
|
||||
Halfedge *h;
|
||||
Facet *f;
|
||||
Vector_3 sum(0., 0., 0.), n;
|
||||
|
||||
Halfedge_around_vertex_circulator
|
||||
hedgeb = v->vertex_begin(), hedgee = hedgeb;
|
||||
|
||||
do
|
||||
{
|
||||
h = &(*hedgeb);
|
||||
if (h->is_border_edge())
|
||||
{
|
||||
hedgeb++;
|
||||
continue;
|
||||
}
|
||||
|
||||
f = &(*h->facet());
|
||||
n = f->getUnitNormal();
|
||||
sum = (sum + n);
|
||||
hedgeb++;
|
||||
}
|
||||
while (hedgeb != hedgee);
|
||||
sum = sum / std::sqrt(sum * sum);
|
||||
return sum;
|
||||
}
|
||||
|
|
@ -26,55 +26,49 @@ template < class Refs, class Tag, class Pt, class FGeomTraits >
|
|||
class My_vertex:public CGAL::HalfedgeDS_vertex_base < Refs, Tag, Pt >
|
||||
{
|
||||
protected:
|
||||
typedef typename Refs::Vertex Vertex;
|
||||
typedef typename Refs::Halfedge Halfedge;
|
||||
typedef typename Refs::Face Facet;
|
||||
|
||||
//to be def in Vertex_wrapper too from the Traits
|
||||
typedef typename FGeomTraits::FT FT;
|
||||
typedef typename FGeomTraits::Point_3 Point_3;
|
||||
typedef typename FGeomTraits::Vector_3 Vector_3;
|
||||
|
||||
|
||||
protected:
|
||||
//monge info
|
||||
Vector_3 d1; //max ppal dir
|
||||
Vector_3 d2; //min ppal dir; monge normal is then n_m=d1^d2, should be so that n_m.n_mesh>0
|
||||
FT k1, k2; //max/min ppal curv
|
||||
FT b0, b3; //blue/red extremalities
|
||||
FT P1, P2; //if fourth order quantities
|
||||
//monge data
|
||||
Vector_3 m_d1; //max ppal dir
|
||||
Vector_3 m_d2; //min ppal dir; monge normal is then n_m=d1^d2, should be so that n_m.n_mesh>0
|
||||
FT m_k1, m_k2; //max/min ppal curv
|
||||
FT m_b0, m_b3; //blue/red extremalities
|
||||
FT m_P1, m_P2; //if fourth order quantities
|
||||
|
||||
//this is for collecting i-th ring neighbours
|
||||
char ring_tag;
|
||||
int ring_index;
|
||||
|
||||
public:
|
||||
//monge info
|
||||
const Vector _3 d1() const { return d1; }
|
||||
Vector_3& d1() { return d1; }
|
||||
const Vector _3 d2() const { return d2; }
|
||||
Vector_3& d2() { return d2; }
|
||||
const FT k1() const { return k1; }
|
||||
FT& k1() { return k1; }
|
||||
const FT k2() const { return k2; }
|
||||
FT& k2() { return k2; }
|
||||
const FT b0() const { return b0; }
|
||||
FT& b0() { return b0; }
|
||||
const FT b3() const { return b3; }
|
||||
FT& b3() { return b3; }
|
||||
const FT P1() const { return P1; }
|
||||
FT& P1() { return P1; }
|
||||
const FT P2() const { return P2; }
|
||||
FT& P2() { return P2; }
|
||||
//monge data
|
||||
const Vector_3 d1() const { return m_d1; }
|
||||
Vector_3& d1() { return m_d1; }
|
||||
const Vector_3 d2() const { return d2; }
|
||||
Vector_3& d2() { return m_d2; }
|
||||
const FT k1() const { return m_k1; }
|
||||
FT& k1() { return m_k1; }
|
||||
const FT k2() const { return m_k2; }
|
||||
FT& k2() { return m_k2; }
|
||||
const FT b0() const { return m_b0; }
|
||||
FT& b0() { return m_b0; }
|
||||
const FT b3() const { return m_b3; }
|
||||
FT& b3() { return m_b3; }
|
||||
const FT P1() const { return m_P1; }
|
||||
FT& P1() { return m_P1; }
|
||||
const FT P2() const { return m_P2; }
|
||||
FT& P2() { return m_P2; }
|
||||
|
||||
//this is for collecting i-th ring neighbours
|
||||
void setRingIndex(int i) { ring_index = i; }
|
||||
int getRingIndex() { return ring_index; }
|
||||
void resetRingIndex() { ring_index = -1; }
|
||||
void setRingTag() { ring_tag = 1; }
|
||||
char getRingTag() { return ring_tag; }
|
||||
|
||||
|
||||
My_vertex(const Point_3 & pt):
|
||||
CGAL::HalfedgeDS_vertex_base < Refs, Tag, Point_3 > (pt),
|
||||
ring_tag(0), ring_index(-1) {}
|
||||
ring_index(-1) {}
|
||||
My_vertex() {}
|
||||
};
|
||||
|
||||
|
|
@ -85,52 +79,42 @@ class My_vertex:public CGAL::HalfedgeDS_vertex_base < Refs, Tag, Pt >
|
|||
template < class Refs, class Tag, class FGeomTraits >
|
||||
class My_facet:public CGAL::HalfedgeDS_face_base < Refs, Tag >
|
||||
{
|
||||
//protected:
|
||||
public:
|
||||
typedef typename Refs::Vertex Vertex;
|
||||
typedef typename Refs::Vertex_handle Vertex_handle;
|
||||
typedef typename Refs::Halfedge Halfedge;
|
||||
typedef typename Refs::Halfedge_handle Halfedge_handle;
|
||||
typedef typename FGeomTraits::Vector_3 Vector_3;
|
||||
|
||||
typedef typename FGeomTraits::Vector_3 Vector_3;
|
||||
typedef typename FGeomTraits::Point_3 Point_3;
|
||||
|
||||
public:
|
||||
protected:
|
||||
Vector_3 normal;
|
||||
|
||||
int ring_index;
|
||||
bool m_is_visited;
|
||||
public:
|
||||
My_facet(): ring_index(-1) {}
|
||||
Vector_3 & getUnitNormal() { return normal; }
|
||||
void setNormal(Vector_3 & n) { normal = n; }
|
||||
void setNormal(Vector_3 n) { normal = n; }
|
||||
|
||||
//this is for collecting i-th ring neighbours
|
||||
protected:
|
||||
int ring_index;
|
||||
public:
|
||||
void setRingIndex(int i) { ring_index = i; }
|
||||
int getRingIndex() { return ring_index; }
|
||||
void resetRingIndex() { ring_index = -1; }
|
||||
//this is for following ridge lines
|
||||
void set_visited(bool b) { m_is_visited = b; }
|
||||
const bool is_visited() const { return m_is_visited;}
|
||||
void reset_is_visited() { m_is_visited = false; }
|
||||
};
|
||||
|
||||
//----------------------------------------------------------------
|
||||
// Halfedge
|
||||
//----------------------------------------------------------------
|
||||
template < class Refs, class Tprev, class Tvertex, class Tface, class FGeomTraits >
|
||||
template < class Refs, class Tprev, class Tvertex, class Tface >
|
||||
class My_halfedge:public CGAL::HalfedgeDS_halfedge_base < Refs, Tprev, Tvertex, Tface >
|
||||
{
|
||||
protected:
|
||||
typedef typename FGeomTraits::Point_3 Point_3;
|
||||
typedef typename FGeomTraits::Vector_3 Vector_3;
|
||||
|
||||
protected:
|
||||
int ring_index;
|
||||
double len;
|
||||
public:
|
||||
My_halfedge(): ring_index(-1) {}
|
||||
void setRingIndex(int i) { ring_index = i; }
|
||||
int getRingIndex() {return ring_index; }
|
||||
void resetRingIndex() {ring_index = -1; }
|
||||
public:
|
||||
My_halfedge(): ring_index(-1) {}
|
||||
void setLength(double l) { len = l; }
|
||||
double getLength() { return len; }
|
||||
};
|
||||
|
|
@ -143,9 +127,11 @@ struct Wrappers_VFH:public CGAL::Polyhedron_items_3 {
|
|||
template < class Refs, class Traits > struct Vertex_wrapper {
|
||||
typedef struct {
|
||||
public:
|
||||
//typedef typename Traits::Vector_3 Vector_3;
|
||||
//all types needed by the vertex...
|
||||
typedef typename Traits::FT FT;
|
||||
typedef typename Traits::Point_3 Point_3;
|
||||
typedef typename Traits::Vector_3 Vector_3;
|
||||
typedef typename Traits::Plane_3 Plane_3;
|
||||
} FGeomTraits;
|
||||
typedef typename Traits::Point_3 Point_3;
|
||||
typedef My_vertex < Refs, CGAL::Tag_true, Point_3, FGeomTraits > Vertex;
|
||||
|
|
@ -158,65 +144,34 @@ struct Wrappers_VFH:public CGAL::Polyhedron_items_3 {
|
|||
//all types needed by the facet...
|
||||
typedef struct {
|
||||
public:
|
||||
typedef typename Traits::Point_3 Point_3;
|
||||
// typedef typename Traits::Point_3 Point_3;
|
||||
typedef typename Traits::Vector_3 Vector_3;
|
||||
typedef typename Traits::Plane_3 Plane_3;
|
||||
} FGeomTraits;
|
||||
} FGeomTraits;
|
||||
//custom type instantiated...
|
||||
typedef My_facet < Refs, CGAL::Tag_true, FGeomTraits > Face;
|
||||
};
|
||||
|
||||
// wrap halfedge
|
||||
template < class Refs, class Traits > struct Halfedge_wrapper {
|
||||
typedef struct {
|
||||
public:
|
||||
typedef typename Traits::Point_3 Point_3;
|
||||
typedef typename Traits::Vector_3 Vector_3;
|
||||
typedef typename Traits::Plane_3 Plane_3;
|
||||
} FGeomTraits;
|
||||
typedef My_halfedge < Refs,
|
||||
typedef My_halfedge < Refs,
|
||||
CGAL::Tag_true,
|
||||
CGAL::Tag_true, CGAL::Tag_true, FGeomTraits > Halfedge;
|
||||
CGAL::Tag_true, CGAL::Tag_true > Halfedge;
|
||||
};
|
||||
};
|
||||
|
||||
//------------------------------------------------
|
||||
// Polyhedron
|
||||
//PolyhedralSurf
|
||||
//------------------------------------------------
|
||||
//using standard Cartesian kernel
|
||||
typedef double DFT;
|
||||
typedef CGAL::Cartesian<DFT> Data_Kernel;
|
||||
|
||||
typedef CGAL::Polyhedron_3 < Data_Kernel, Wrappers_VFH > Polyhedron;
|
||||
typedef Polyhedron::Vertex Vertex;
|
||||
typedef Polyhedron::Vertex_handle Vertex_handle;
|
||||
typedef Polyhedron::Halfedge Halfedge;
|
||||
typedef Polyhedron::Halfedge_handle Halfedge_handle;
|
||||
|
||||
typedef Polyhedron::Vertex_iterator Vertex_iterator;
|
||||
typedef Polyhedron::Halfedge_iterator Halfedge_iterator;
|
||||
typedef Polyhedron::
|
||||
Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator;
|
||||
typedef Polyhedron::
|
||||
Halfedge_around_facet_circulator Halfedge_around_facet_circulator;
|
||||
typedef Polyhedron::Facet_iterator Facet_iterator;
|
||||
|
||||
typedef Data_Kernel::Point_3 DPoint;
|
||||
typedef Data_Kernel::Vector_3 Vector_3;
|
||||
typedef Data_Kernel::Plane_3 Plane_3;
|
||||
typedef Data_Kernel::Sphere_3 Sphere_3;
|
||||
|
||||
/////////////////////class PolyhedralSurf///////////////////////////
|
||||
class PolyhedralSurf:public Polyhedron {
|
||||
public:
|
||||
typedef Data_Kernel::Point_3 DPoint;
|
||||
typedef Data_Kernel::Vector_3 Vector_3;
|
||||
typedef Data_Kernel::Plane_3 Plane_3;
|
||||
|
||||
public:
|
||||
static Vector_3 getHalfedge_vector(Halfedge * h);
|
||||
|
||||
PolyhedralSurf() {}
|
||||
|
||||
static Vector_3 getHalfedge_vector(Halfedge * h);
|
||||
void compute_edges_length();
|
||||
double compute_mean_edges_length_around_vertex(Vertex * v);
|
||||
void compute_facets_normals();
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
#ifndef _POLY_OP_H_
|
||||
#define _POLY_OP_H_
|
||||
|
||||
//?????????mysterious , used in Poly.C with mysterious
|
||||
//std::for_all??????? why arent they member fct of polysurf?
|
||||
|
||||
struct Edge_length {
|
||||
template < class HalfEdge >
|
||||
void operator() (HalfEdge & h)
|
||||
{
|
||||
double d =
|
||||
CGAL::squared_distance(h.prev()->vertex()->point(),
|
||||
h.vertex()->point());
|
||||
h.setLength(CGAL::sqrt(d));
|
||||
}
|
||||
};
|
||||
|
||||
//the facet stores the normal
|
||||
struct Facet_unit_normal {
|
||||
template < class Facet >
|
||||
void operator() (Facet & f)
|
||||
{
|
||||
typename Facet::Halfedge_handle h = f.halfedge();
|
||||
typename Facet::Vector_3 normal =
|
||||
CGAL::cross_product(h->vertex()->point() -
|
||||
h->opposite()->vertex()->point(),
|
||||
h->next()->vertex()->point() -
|
||||
h->opposite()->vertex()->point());
|
||||
f.setNormal( normal / CGAL::sqrt(normal * normal));
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,640 @@
|
|||
#ifndef _SURPOLYRINGS_H_
|
||||
#define _SURPOLYRINGS_H_
|
||||
|
||||
using namespace std;
|
||||
|
||||
//--------------------------------------------------
|
||||
//we store hedges. exple: 1ring neighbours of
|
||||
//v=hedges anchored at v
|
||||
//--------------------------------------------------
|
||||
// when storing hedges to a ith ring: do we store hedges pointing to
|
||||
// the ring or opposite hedges?
|
||||
enum HTO_HOPPOSITE_ONERING { HTO_ONERING, HOPPOSITE_ONERING };
|
||||
|
||||
template < class TPoly > class T_PolyhedralSurf_rings
|
||||
{
|
||||
public:
|
||||
// typedef typename TPoly::Point_3 Point_3;
|
||||
typedef typename TPoly::Vertex Vertex;
|
||||
// typedef typename TPoly::Vertex_handle Vertex_handle;
|
||||
typedef typename TPoly::Halfedge Halfedge;
|
||||
typedef typename TPoly::Facet Facet;
|
||||
typedef typename TPoly::Halfedge_around_vertex_circulator
|
||||
Halfedge_around_vertex_circulator;
|
||||
typedef typename TPoly::Vertex_iterator Vertex_iterator;
|
||||
//debug
|
||||
static void check_vertex_indices(Vertex* v);
|
||||
static void check_ring_indices(TPoly& poly);
|
||||
|
||||
// static void from_into(ptrSet & a, ptrSet & b, char erase_a = 0);
|
||||
|
||||
|
||||
//---------------VERTEX---------------------------------
|
||||
//neighbours using a tag on each vertex ---default is -1, vertex is 0,
|
||||
// then 1,2,.. for ith rings
|
||||
//---------------------------------------------------------
|
||||
static int
|
||||
push_neighbours_of(Vertex * start, int ith,
|
||||
std::vector < Vertex * >&nextRing,
|
||||
std::vector < Vertex * >&all);
|
||||
|
||||
static int
|
||||
collect_ith_ring_neighbours(int ith,
|
||||
std::vector < Vertex * >¤tRing,
|
||||
std::vector < Vertex * >&nextRing,
|
||||
std::vector < Vertex * >&all);
|
||||
|
||||
|
||||
static void reset_ring_indices(std::vector < Vertex * >&vces);
|
||||
// static void reset_ring_indices(std::vector< Vertex_handle >&vces);
|
||||
//Neighbours collected in ccw order on each ring
|
||||
static int collect_ith_ring_ccw(int ith,
|
||||
vector<int>& sizes,
|
||||
vector<Vertex*>& vces);
|
||||
static int i_rings_ccw(int i,
|
||||
Vertex* start,
|
||||
vector<int>& sizes,
|
||||
vector<Vertex*>& vces);
|
||||
|
||||
//----------------------------------------------------------------
|
||||
// Facet
|
||||
//----------------------------------------------------------------
|
||||
|
||||
|
||||
|
||||
//------------------HEDGES--------------------------------
|
||||
//hdges incident to endpoint of h
|
||||
//we store opposite hedges, ie those
|
||||
//to the next ith ring neighbours
|
||||
//---------------------------------------------------------
|
||||
static int
|
||||
push_hedges_of(Vertex * start, int ith,
|
||||
HTO_HOPPOSITE_ONERING to_opp,
|
||||
std::vector < Halfedge * >&nextRing,
|
||||
std::vector < Halfedge * >&recordedOnce);
|
||||
//, char& multiple);
|
||||
|
||||
//grabbing the next ring from currentRing
|
||||
static int
|
||||
next_ring_hedges(int ith, HTO_HOPPOSITE_ONERING to_opp,
|
||||
std::vector < Halfedge * >¤tRing,
|
||||
std::vector < Halfedge * >&nextRing,
|
||||
std::vector < Halfedge * >&recordedOnce);
|
||||
//, char& multiple);
|
||||
|
||||
//ith ring from scratch ie the center vertex
|
||||
static int
|
||||
collect_ith_ring_hedges(Vertex * start, int ith,
|
||||
HTO_HOPPOSITE_ONERING to_opp,
|
||||
std::vector < Halfedge * >&recordedOnce,
|
||||
char cleanup);
|
||||
|
||||
//Collect contour edges of the ith ring, these edges oriented CCW
|
||||
// i.e. their incident facet is inside.
|
||||
static void
|
||||
contour_ith_ring_hedges(Vertex * start, int ith,
|
||||
std::vector < Halfedge * > &contour);
|
||||
|
||||
static int
|
||||
contour_ccw_ith_ring_hedges(Vertex * start, int ith,
|
||||
std::vector < Halfedge * > &contourCCW);
|
||||
|
||||
static void reset_ring_indices(std::vector < Halfedge * >&hedges);
|
||||
|
||||
static void reset_ring_indices(std::vector < Facet * >&facets);
|
||||
|
||||
// to report 1st crossing edge with a sphere
|
||||
// static Point_3 orig_point;
|
||||
// static void set_orig_point(Point_3 & p);
|
||||
// static void
|
||||
// getFarthestNeighbourOfV_ofOrigPoint(Vertex * v,
|
||||
// Halfedge * &h_farthest);
|
||||
|
||||
};
|
||||
|
||||
|
||||
//////////////////////////////////////////////
|
||||
//////////IMLPEMENTATION///////////////////////
|
||||
//////////////////////////////////////////////
|
||||
//--------------------------------------------------
|
||||
//-- collecting ith ring neighbours
|
||||
//--------------------------------------------------
|
||||
|
||||
template < class TPoly >
|
||||
void T_PolyhedralSurf_rings <TPoly >::
|
||||
check_vertex_indices(Vertex* v)
|
||||
{
|
||||
Halfedge_around_vertex_circulator itb = v->vertex_begin(),
|
||||
ite = v->vertex_begin();
|
||||
assert( ((&(*itb))->opposite()->vertex()->getRingIndex()==-1) );
|
||||
itb++;
|
||||
for( ; itb != ite ; itb++)
|
||||
assert(((&(*itb))->opposite()->vertex()->getRingIndex()==-1) );
|
||||
}
|
||||
|
||||
|
||||
template < class TPoly >
|
||||
void T_PolyhedralSurf_rings < TPoly >::
|
||||
check_ring_indices(TPoly& poly)
|
||||
{
|
||||
Vertex_iterator itb = poly.vertices_begin (),
|
||||
ite = poly.vertices_end();
|
||||
for( ; itb != ite ; itb++)
|
||||
T_PolyhedralSurf_rings < TPoly >::check_vertex_indices( &(*itb));
|
||||
}
|
||||
|
||||
// //neighbours using a map to make sure we do not collect twice
|
||||
// //---------------------------------------------------------
|
||||
// template < class TPoly >
|
||||
// void T_PolyhedralSurf_rings < TPoly >::
|
||||
// from_into(ptrSet & a, ptrSet & b, char erase_a)
|
||||
// {
|
||||
// ptrSet_iterator itb = a.begin(), ite = a.end();
|
||||
// for (; itb != ite; itb++)
|
||||
// b.insert(*itb);
|
||||
|
||||
// if (erase_a)
|
||||
// a.erase(a.begin(), ite);
|
||||
// }
|
||||
|
||||
//neighbours using a tag on each vertex ---default is -1, vertex is 0,
|
||||
// then 1,2,.. for ith rings
|
||||
//-------------------------------------------------------------------
|
||||
|
||||
template < class TPoly >
|
||||
int T_PolyhedralSurf_rings < TPoly >::
|
||||
push_neighbours_of(Vertex * start, int ith,
|
||||
std::vector < Vertex * >&nextRing,
|
||||
std::vector < Vertex * >&all)
|
||||
{
|
||||
int nb = 0;
|
||||
Vertex *v;
|
||||
Halfedge *h;
|
||||
|
||||
Halfedge_around_vertex_circulator
|
||||
hedgeb = start->vertex_begin(), hedgee = hedgeb;
|
||||
|
||||
do
|
||||
{
|
||||
h = &(*hedgeb);
|
||||
hedgeb++;
|
||||
v = &(*h->opposite()->vertex());
|
||||
|
||||
//make sure vertex not visited yet
|
||||
if (v->getRingIndex() != -1)
|
||||
continue;
|
||||
//tag the vertex as belonging to the its ring
|
||||
v->setRingIndex(ith);
|
||||
nb++;
|
||||
nextRing.push_back(v);
|
||||
all.push_back(v);
|
||||
}
|
||||
while (hedgeb != hedgee);
|
||||
return nb;
|
||||
}
|
||||
|
||||
template < class TPoly >
|
||||
int T_PolyhedralSurf_rings < TPoly >::
|
||||
collect_ith_ring_neighbours(int ith, std::vector < Vertex * >¤tRing,
|
||||
std::vector < Vertex * >&nextRing,
|
||||
std::vector < Vertex * >&all)
|
||||
{
|
||||
Vertex *v;
|
||||
typename std::vector < Vertex * >::iterator itb =
|
||||
currentRing.begin(), ite = currentRing.end();
|
||||
|
||||
for (; itb != ite; itb++)
|
||||
{
|
||||
v = *itb;
|
||||
T_PolyhedralSurf_rings < TPoly >::
|
||||
push_neighbours_of(v, ith, nextRing, all);
|
||||
}
|
||||
return nextRing.size();
|
||||
}
|
||||
|
||||
template < class TPoly >
|
||||
void T_PolyhedralSurf_rings < TPoly >::
|
||||
reset_ring_indices(std::vector < Vertex * >&vces)
|
||||
{
|
||||
typename std::vector < Vertex * >::iterator itb = vces.begin(), ite =
|
||||
vces.end();
|
||||
for (; itb != ite; itb++)
|
||||
(*itb)->resetRingIndex();
|
||||
}
|
||||
|
||||
// //surcharge pour les handles
|
||||
// template < class TPoly >
|
||||
// void T_PolyhedralSurf_rings < TPoly >::
|
||||
// reset_ring_indices(std::vector < Vertex_handle >&vces)
|
||||
// {
|
||||
// typename std::vector < Vertex_handle >::iterator itb = vces.begin(), ite =
|
||||
// vces.end();
|
||||
// for (; itb != ite; itb++)
|
||||
// (&(*(*itb)))->resetRingIndex();
|
||||
// }
|
||||
|
||||
//see also the other method below--------------
|
||||
//Neighbours collected in ccw order on each ring
|
||||
//if a border edge is encountered, return 0
|
||||
//sizes[j] = nb of vertives of the jth ring
|
||||
// template <class TPoly>
|
||||
// int T_PolyhedralSurf_rings<TPoly>::
|
||||
// collect_ith_ring_ccw(int ith,
|
||||
// vector<int>& sizes,
|
||||
// vector<Vertex*>& vces)
|
||||
// {
|
||||
// //find a vertex of the ith ring from a vertex of the (i-1)th ring
|
||||
// int size=0;
|
||||
// Vertex* v=vces[vces.size()-1];
|
||||
// Halfedge* h;
|
||||
// Halfedge_around_vertex_circulator
|
||||
// hedgeb = v->vertex_begin(), hedgee = hedgeb;
|
||||
// h = &(*hedgeb);
|
||||
// if (h->is_border_edge()) return 0;
|
||||
// while (h->opposite()->vertex()->getRingIndex() != -1)
|
||||
// {
|
||||
// hedgeb++;
|
||||
// h = &(*hedgeb);
|
||||
// if (h->is_border_edge()) return 0;
|
||||
// }
|
||||
// v = &(*(h->opposite()->vertex()));
|
||||
// vces.push_back(v);
|
||||
// size++;
|
||||
// v->setRingIndex(ith);
|
||||
// h=&(*h->opposite()->next());
|
||||
// if (h->is_border_edge()) return 0;
|
||||
// v=&(*h->vertex());
|
||||
// //follow ccw the ith ring
|
||||
// while (v->getRingIndex() != ith)
|
||||
// {
|
||||
// while ( (v->getRingIndex() != -1) && (v->getRingIndex() != ith) )
|
||||
// {
|
||||
// h=&(*h->opposite()->next());
|
||||
// if (h->is_border_edge()) return 0;
|
||||
// v=&(*h->vertex());
|
||||
// }
|
||||
// if (v->getRingIndex() == -1)
|
||||
// {
|
||||
// vces.push_back(v);
|
||||
// size++;
|
||||
// v->setRingIndex(ith);
|
||||
// h=&(*h->next());
|
||||
// if (h->is_border_edge()) return 0;
|
||||
// v=&(*h->vertex());
|
||||
// }
|
||||
// }
|
||||
// sizes.push_back(size);
|
||||
// return 1;
|
||||
// }
|
||||
//Other method using contour_ccw_ith_ring_hedges--------------------------
|
||||
template <class TPoly>
|
||||
int T_PolyhedralSurf_rings<TPoly>::
|
||||
collect_ith_ring_ccw(int ith,
|
||||
vector<int>& sizes,
|
||||
vector<Vertex*>& vces)
|
||||
{
|
||||
Vertex* v=vces[0];
|
||||
std::vector < Halfedge * > contourCCW;
|
||||
//check if no border edge is encountered and set contourCCW of halfedges
|
||||
if (T_PolyhedralSurf_rings<TPoly>::
|
||||
contour_ccw_ith_ring_hedges(v, ith, contourCCW) == 0) return 0;
|
||||
|
||||
typename std::vector < Halfedge * >::iterator itb =
|
||||
contourCCW.begin(), ite = contourCCW.end();
|
||||
|
||||
//retrieve vertices from hedge contour
|
||||
for (; itb != ite; itb++) vces.push_back(&(*((*itb)->vertex())));
|
||||
sizes.push_back(contourCCW.size());
|
||||
return 1;
|
||||
}
|
||||
|
||||
|
||||
template <class TPoly>
|
||||
int T_PolyhedralSurf_rings<TPoly>::
|
||||
i_rings_ccw(int i,
|
||||
Vertex* start,
|
||||
vector<int>& sizes,
|
||||
vector<Vertex*>& vces)
|
||||
{
|
||||
//set the 0-ring
|
||||
start->setRingIndex(0);
|
||||
sizes.push_back(1);
|
||||
vces.push_back(start);
|
||||
int ith;
|
||||
//set ith ring, ith from 1 to i
|
||||
for(ith=1;ith<=i;ith++)
|
||||
if (T_PolyhedralSurf_rings<TPoly>::
|
||||
collect_ith_ring_ccw(ith, sizes, vces) == 0)
|
||||
{
|
||||
T_PolyhedralSurf_rings<TPoly>::reset_ring_indices(vces);
|
||||
return 0;
|
||||
}
|
||||
T_PolyhedralSurf_rings<TPoly>::reset_ring_indices(vces);
|
||||
return 1;
|
||||
}
|
||||
|
||||
//-----------------------------------------------------
|
||||
//---------- hedges.
|
||||
//-----------------------------------------------------
|
||||
template < class TPoly >
|
||||
int T_PolyhedralSurf_rings < TPoly >::
|
||||
push_hedges_of(Vertex * start, int ith, HTO_HOPPOSITE_ONERING to_opp,
|
||||
std::vector < Halfedge * >&nextRing,
|
||||
std::vector < Halfedge * >&recordedOnce)
|
||||
{
|
||||
// char ok;
|
||||
int n = 0;
|
||||
Halfedge *h, *hopp, *hpushed;
|
||||
Halfedge_around_vertex_circulator
|
||||
hedgeb = start->vertex_begin(), hedgee = hedgeb;
|
||||
do
|
||||
{
|
||||
h = &(*hedgeb);
|
||||
hopp = &(*h->opposite());
|
||||
hedgeb++;
|
||||
//store either h or hopp
|
||||
hpushed = (to_opp == HTO_ONERING) ? hopp : h;
|
||||
|
||||
//in any case
|
||||
nextRing.push_back(hpushed);
|
||||
|
||||
//record or not?
|
||||
if (hpushed->getRingIndex() == -1)
|
||||
{
|
||||
//tag the vertex as belonging to the its ring
|
||||
//marc: isn't it: tag the edge and its opposite?
|
||||
h->setRingIndex(ith);
|
||||
hopp->setRingIndex(ith);
|
||||
n++; //xfc
|
||||
recordedOnce.push_back(hpushed);
|
||||
}
|
||||
//hedge already recorded
|
||||
//else multiple=1;
|
||||
}
|
||||
while (hedgeb != hedgee);
|
||||
return n;
|
||||
}
|
||||
|
||||
//notice that for each helfedge, we need to recover the endpoint, ie
|
||||
//the vertex on the current ring being explored
|
||||
template < class TPoly >
|
||||
int T_PolyhedralSurf_rings < TPoly >::
|
||||
next_ring_hedges(int ith, HTO_HOPPOSITE_ONERING to_opp,
|
||||
std::vector < Halfedge * >¤tRing,
|
||||
std::vector < Halfedge * >&nextRing,
|
||||
std::vector < Halfedge * >&recordedOnce)
|
||||
{
|
||||
int s = 0;
|
||||
Vertex *v;
|
||||
Halfedge *h;
|
||||
typename std::vector < Halfedge * >::iterator itb =
|
||||
currentRing.begin(), ite = currentRing.end();
|
||||
for (; itb != ite; itb++)
|
||||
{
|
||||
h = (Halfedge *) (*itb);
|
||||
|
||||
//grab vertex
|
||||
if (to_opp == HTO_ONERING)
|
||||
v = &*(h->vertex());
|
||||
else
|
||||
v = &*(h->opposite()->vertex());
|
||||
|
||||
s += T_PolyhedralSurf_rings < TPoly >::
|
||||
push_hedges_of(v, ith, to_opp, nextRing, recordedOnce);
|
||||
}
|
||||
return s;
|
||||
}
|
||||
|
||||
//Warning collect hedges from rings 1 to ith, better to call it
|
||||
//collect_i_rings_hedges ???????
|
||||
template < class TPoly >
|
||||
int T_PolyhedralSurf_rings < TPoly >::
|
||||
collect_ith_ring_hedges(Vertex * start, int ith,
|
||||
HTO_HOPPOSITE_ONERING to_opp,
|
||||
std::vector < Halfedge * >&recordedOnce,
|
||||
char cleanup)
|
||||
{
|
||||
// char multiple=0;
|
||||
int i = 1, np;
|
||||
std::vector < Halfedge * >baseRing, nextRing;
|
||||
// typename std::vector<Halfedge*>::iterator itb, ite;
|
||||
std::vector < Halfedge * >*p_baseRing = &baseRing, *p_nextRing =
|
||||
&nextRing;
|
||||
|
||||
//collect the one ring
|
||||
np = T_PolyhedralSurf_rings < TPoly >::
|
||||
push_hedges_of(start, int (1), to_opp, nextRing, recordedOnce);
|
||||
|
||||
while (1)
|
||||
{
|
||||
//exit?
|
||||
i++;
|
||||
if (i > ith)
|
||||
break;
|
||||
|
||||
//find next ring
|
||||
p_baseRing->erase(p_baseRing->begin(), p_baseRing->end());
|
||||
std::swap(p_baseRing, p_nextRing);
|
||||
|
||||
np = T_PolyhedralSurf_rings < TPoly >::
|
||||
next_ring_hedges(i, HTO_ONERING, *p_baseRing, *p_nextRing,
|
||||
recordedOnce);
|
||||
}
|
||||
//clean up hedges tags
|
||||
np = recordedOnce.size(); //#triangles
|
||||
if (cleanup)
|
||||
T_PolyhedralSurf_rings < TPoly >::reset_ring_indices(recordedOnce);
|
||||
return np;
|
||||
}
|
||||
|
||||
//Warning contour edges are not given in a cyclic order
|
||||
template < class TPoly >
|
||||
void T_PolyhedralSurf_rings < TPoly >::
|
||||
contour_ith_ring_hedges(Vertex * start, int ith,
|
||||
std::vector < Halfedge * > &contour)
|
||||
{
|
||||
std::vector < Halfedge * > ringHedges;
|
||||
typename std::vector < Halfedge * >::iterator itb =
|
||||
ringHedges.begin(), ite = ringHedges.end();
|
||||
T_PolyhedralSurf_rings < TPoly >::
|
||||
collect_ith_ring_hedges(start, ith, HTO_ONERING, ringHedges, 0);
|
||||
//debug
|
||||
// fprintf(stderr, "\n");
|
||||
// int iii;
|
||||
// for( iii=0;iii<ringHedges.size();iii++)
|
||||
// fprintf(stderr, "[ ringedge %d , edgeindex %d, vertex %d ,opposite vertex %d] \n",
|
||||
// ringHedges[iii],
|
||||
// ringHedges[iii]->getRingIndex(),
|
||||
// &(*( ringHedges[iii]->vertex() )),
|
||||
// &(*(ringHedges[iii]->opposite()->vertex() )) );
|
||||
//debug
|
||||
|
||||
|
||||
//set ringHedges index to i
|
||||
// for (; itb != ite; itb++) (*itb)->setRingIndex(ith);non! garder
|
||||
// les index!!
|
||||
|
||||
//for a hedge h in ringHedges and h index = ith, h->next is on the
|
||||
// contour if its opposite is not visited
|
||||
//Note that contour edges are then oriented CCW
|
||||
itb = ringHedges.begin(), ite = ringHedges.end();
|
||||
for (; itb != ite; itb++)
|
||||
{
|
||||
if ( ( (*itb)->getRingIndex()==ith )&&
|
||||
((*itb)->next()->opposite()->getRingIndex() == -1) )
|
||||
contour.push_back( &(*((*itb)->next())) );
|
||||
}
|
||||
//debug
|
||||
// fprintf(stderr, "\n\n");
|
||||
// for( iii=0;iii<contour.size();iii++)
|
||||
// fprintf(stderr, "[ contour edge %d , edgeindex %d, vertex %d ,opposite vertex %d] \n",
|
||||
// contour[iii], contour[iii]->getRingIndex(),
|
||||
// &(*( contour[iii]->vertex() )),
|
||||
// &(*(contour[iii]->opposite()->vertex() )) );
|
||||
// //debug
|
||||
|
||||
//reset indices
|
||||
T_PolyhedralSurf_rings < TPoly >::reset_ring_indices(ringHedges);
|
||||
}
|
||||
|
||||
//if the contour is a circle, then put edges head to tail, ccw
|
||||
template < class TPoly >
|
||||
int T_PolyhedralSurf_rings < TPoly >::
|
||||
contour_ccw_ith_ring_hedges(Vertex * start, int ith,
|
||||
std::vector < Halfedge * > &contourCCW)
|
||||
{
|
||||
std::vector < Halfedge * > contour;
|
||||
T_PolyhedralSurf_rings < TPoly >::
|
||||
contour_ith_ring_hedges(start, ith, contour);
|
||||
|
||||
//check if no border edge is encountered
|
||||
typename std::vector < Halfedge * >::iterator itb =
|
||||
contour.begin(), ite = contour.end();
|
||||
for (; itb != ite; itb++) if ( (*itb)->is_border_edge()) return 0;
|
||||
|
||||
int size = contour.size();
|
||||
Halfedge* h_cur, *h;
|
||||
h_cur = contour[0];
|
||||
// fprintf(stderr, "\n\n");
|
||||
// int ii;
|
||||
// for( ii=0, itb=contour.begin();ii<contour.size();ii++, itb++){
|
||||
// fprintf(stderr, "[%d ", contour[ii]->vertex());
|
||||
// fprintf(stderr, " %d] ", (*itb)->vertex());}
|
||||
|
||||
// fprintf(stderr, "\n\n");
|
||||
contourCCW.push_back(h_cur);
|
||||
do
|
||||
{
|
||||
for (itb=contour.begin(); itb != ite; itb++)
|
||||
{//debug
|
||||
h = &(*( ((*itb)->opposite()) ));
|
||||
// cerr << "h_cur->vertex() " << h_cur->vertex() << endl;
|
||||
// << "h->vertex() " << h->vertex() << endl;
|
||||
if ( (h_cur->vertex()) == (h->vertex()) )
|
||||
{//the edge following the current edge is found
|
||||
contourCCW.push_back(*itb);
|
||||
h_cur = *itb;
|
||||
break;
|
||||
}
|
||||
}//debug
|
||||
}
|
||||
while (contourCCW.size() != size);
|
||||
|
||||
//debug
|
||||
// int iii;
|
||||
|
||||
// fprintf(stderr, "\n");
|
||||
|
||||
// for( iii=0;iii<contour.size();iii++)
|
||||
// fprintf(stderr, "[ contour edge %d , edgeindex %d, vertex %d ,opp vertex %d] \n",
|
||||
// contour[iii],
|
||||
// contour[iii]->getRingIndex(),
|
||||
// &(*( contour[iii]->vertex() )),
|
||||
// &(*(contour[iii]->opposite()->vertex() )) );
|
||||
// //debug
|
||||
// //check the contour debug! //debug
|
||||
// fprintf(stderr, "\n");
|
||||
// for( iii=0;iii<contourCCW.size();iii++)
|
||||
// fprintf(stderr, "[ contourCCW edge %d , edgeindex %d, vertex %d ,opp vertex %d] \n",
|
||||
// contourCCW[iii],
|
||||
// contourCCW[iii]->getRingIndex(),
|
||||
// &(*( contourCCW[iii]->vertex() )),
|
||||
// &(*(contourCCW[iii]->opposite()->vertex() )) );
|
||||
// //debug
|
||||
|
||||
for (int i=0; i<contourCCW.size(); i++)
|
||||
{
|
||||
assert ( &(*(contourCCW[i]->vertex())) ==
|
||||
&(*(contourCCW[(i+1) % size]->opposite()->vertex())) );
|
||||
}
|
||||
return 1;
|
||||
}
|
||||
|
||||
|
||||
template < class TPoly >
|
||||
void T_PolyhedralSurf_rings < TPoly >::
|
||||
reset_ring_indices(std::vector < Halfedge * >&hedges)
|
||||
{
|
||||
typename std::vector < Halfedge * >::iterator itb =
|
||||
hedges.begin(), ite = hedges.end();
|
||||
for (; itb != ite; itb++)
|
||||
{
|
||||
(*itb)->resetRingIndex();
|
||||
(*itb)->opposite()->resetRingIndex();
|
||||
}
|
||||
}
|
||||
|
||||
template < class TPoly >
|
||||
void T_PolyhedralSurf_rings < TPoly >::
|
||||
reset_ring_indices(std::vector < Facet * >&facets)
|
||||
{
|
||||
typename std::vector < Facet * >::iterator itb = facets.begin(), ite =
|
||||
facets.end();
|
||||
for (; itb != ite; itb++)
|
||||
(*itb)->resetRingIndex();
|
||||
}
|
||||
|
||||
// // to report 1st crossing edge with a sphere
|
||||
// template < class TPoly >
|
||||
// typename TPoly::Point_3 T_PolyhedralSurf_rings < TPoly >::orig_point;
|
||||
|
||||
// template < class TPoly >
|
||||
// void T_PolyhedralSurf_rings < TPoly >::set_orig_point(Point_3 & p)
|
||||
// {
|
||||
// T_PolyhedralSurf_rings < TPoly >::orig_point = p;
|
||||
// }
|
||||
|
||||
// template < class TPoly >
|
||||
// void T_PolyhedralSurf_rings < TPoly >::
|
||||
// getFarthestNeighbourOfV_ofOrigPoint(Vertex * v, Halfedge * &h_farthest)
|
||||
// {
|
||||
// double d, dref;
|
||||
// Halfedge *h;
|
||||
// Halfedge_around_vertex_circulator
|
||||
// hedgeb = v->vertex_begin(), hedgee = hedgeb;
|
||||
|
||||
// //init
|
||||
// h = h_farthest = &(*hedgeb);
|
||||
// hedgeb++;
|
||||
// dref = squared_distance(h->opposite()->vertex()->point(),
|
||||
// T_PolyhedralSurf_rings < TPoly >::orig_point);
|
||||
// //Visit neighbours
|
||||
// do
|
||||
// {
|
||||
// h = &(*hedgeb);
|
||||
// hedgeb++;
|
||||
|
||||
// assert(v == (&(*h->vertex())));
|
||||
|
||||
// d = squared_distance(h->opposite()->vertex()->point(),
|
||||
// T_PolyhedralSurf_rings < TPoly >::orig_point);
|
||||
// if (d > dref)
|
||||
// {
|
||||
// dref = d;
|
||||
// h_farthest = h;
|
||||
// }
|
||||
// }
|
||||
// while (hedgeb != hedgee);
|
||||
// }
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,44 @@
|
|||
Program blind_1pt.exe
|
||||
---------------------
|
||||
takes as input a <inputPoints.txt> which is a xyz text file
|
||||
it compute the fitting for this single entry and
|
||||
it outputs results in the <output.txt> file and on the standard std::cout
|
||||
|
||||
Usage is : main <inputPoints.txt> <output.txt> <d_fitting>, <d_monge>"
|
||||
|
||||
./blind_1pt.exe data/in_points_file.txt output.txt 2 2
|
||||
|
||||
|
||||
|
||||
Program blind.exe
|
||||
-----------------
|
||||
takes an filename.off file as input,
|
||||
it computes a fitting for each vertex
|
||||
|
||||
it outputs the results in :
|
||||
1. filename.off.4ogl.txt which can be visualized with
|
||||
the demo program visu.exe
|
||||
2. if option -v, filename.off.verb.txt contains human readable results
|
||||
|
||||
Usage is : blind with options
|
||||
"f:fName <string>", // filename.off
|
||||
"d:deg <int>", // degree of the jet
|
||||
"m:mdegree <int>", // degree of the Monge rep
|
||||
"a:nrings <int>", // # rings
|
||||
0 means collect enough rings to make appro possible
|
||||
k>=1 fixes the nb of rings to be collected
|
||||
"p:npoints <int>", //# points
|
||||
0 means this option is not considered, this is the default
|
||||
n>=1 fixes the nb of points to be used
|
||||
"v|", // verbose?
|
||||
|
||||
Note : if the nb of collected points is less than the required min number of
|
||||
points to make the approxiamtion possible (which is constrained by the deg)
|
||||
then the vertex is skipped.
|
||||
|
||||
./blind.exe -f data/ellipe0.003.off -d2 -m2 -a2
|
||||
./blind.exe -f data/poly2x\^2+y\^2-0.062500-off -d2 -m2 -a2
|
||||
|
||||
visu with:
|
||||
|
||||
./visu.exe ../../examples/Jet_fitting_3/data/poly2x\^2+y\^2-0.062500-off ../../examples/Jet_fitting_3/data_poly2x\^2+y\^2-0.062500-off.4ogl.txt
|
||||
|
|
@ -0,0 +1,313 @@
|
|||
#include <CGAL/basic.h>
|
||||
#include <CGAL/Cartesian.h>
|
||||
|
||||
#include <stdio.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <stdlib.h>
|
||||
#include <vector>
|
||||
#include <list>
|
||||
|
||||
#include "../../include/CGAL/Ridges.h"
|
||||
#include "../../../Jet_fitting_3/include/CGAL/Monge_via_jet_fitting.h"
|
||||
#include "GSL.h"
|
||||
|
||||
#include "PolyhedralSurf.h"
|
||||
#include "PolyhedralSurf_rings.h"
|
||||
#include "options.h"//parsing command line
|
||||
|
||||
//Kernel of the PolyhedralSurf
|
||||
typedef double DFT;
|
||||
typedef CGAL::Cartesian<DFT> Data_Kernel;
|
||||
typedef Data_Kernel::Point_3 DPoint;
|
||||
typedef Data_Kernel::Vector_3 DVector;
|
||||
typedef PolyhedralSurf::Vertex Vertex;
|
||||
typedef PolyhedralSurf::Vertex_iterator Vertex_iterator;
|
||||
typedef T_PolyhedralSurf_rings<PolyhedralSurf> Poly_rings;
|
||||
//Kernel for local computations
|
||||
typedef double LFT;
|
||||
typedef CGAL::Cartesian<LFT> Local_Kernel;
|
||||
typedef CGAL::Monge_via_jet_fitting<Data_Kernel, Local_Kernel, GSL> My_Monge_via_jet_fitting;
|
||||
typedef CGAL::Monge_rep<Data_Kernel> My_Monge_rep;
|
||||
typedef CGAL::Monge_info<Local_Kernel> My_Monge_info;
|
||||
|
||||
//RIDGES
|
||||
typedef CGAL::Ridge_line<PolyhedralSurf> Ridge_line;
|
||||
typedef CGAL::Ridge_approximation<PolyhedralSurf,
|
||||
std::list<Ridge_line>::iterator > Ridge_approximation;
|
||||
|
||||
//Syntax requirred by Options
|
||||
static const char *const optv[] = {
|
||||
"?|?",
|
||||
"f:fName <string>", //name of the input off file
|
||||
"d:deg <int>", //degree of the jet
|
||||
"m:mdegree <int>", //degree of the Monge rep
|
||||
"a:nrings <int>", //# rings
|
||||
"p:npoints <int>", //# points
|
||||
"v|",//verbose?
|
||||
NULL
|
||||
};
|
||||
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// fct parameters
|
||||
unsigned int d_fitting = 2;
|
||||
unsigned int d_monge = 2;
|
||||
unsigned int nb_rings = 0;//seek min # of rings to get the required #pts
|
||||
unsigned int nb_points_to_use = 0;//
|
||||
bool verbose = false;
|
||||
|
||||
char *if_name = NULL, //input file name
|
||||
*w_if_name = NULL; //as above, but / replaced by _
|
||||
char* res4openGL_fname;
|
||||
char* verbose_fname;
|
||||
std::ofstream *out_4ogl = NULL, *out_verbose = NULL;
|
||||
|
||||
int optchar;
|
||||
char *optarg;
|
||||
Options opts(*argv, optv);
|
||||
OptArgvIter iter(--argc, ++argv);
|
||||
|
||||
while ((optchar = opts(iter, (const char *&) optarg))){
|
||||
switch (optchar){
|
||||
case 'f': if_name = optarg; break;
|
||||
case 'd': d_fitting = atoi(optarg); break;
|
||||
case 'm': d_monge = atoi(optarg); break;
|
||||
case 'a': nb_rings = atoi(optarg); break;
|
||||
case 'p': nb_points_to_use = atoi(optarg); break;
|
||||
case 'v': verbose=true; break;
|
||||
default:
|
||||
cerr << "Unknown command line option " << optarg;
|
||||
exit(0);
|
||||
}
|
||||
}
|
||||
|
||||
//prepare output file names
|
||||
assert(if_name != NULL);
|
||||
w_if_name = new char[strlen(if_name)];
|
||||
strcpy(w_if_name, if_name);
|
||||
for(unsigned int i=0; i<strlen(w_if_name); i++)
|
||||
if (w_if_name[i] == '/') w_if_name[i]='_';
|
||||
cerr << if_name << '\n';
|
||||
cerr << w_if_name << '\n';
|
||||
|
||||
res4openGL_fname = new char[strlen(w_if_name) + 9];// append .4ogl.txt
|
||||
sprintf(res4openGL_fname, "%s.4ogl.txt", w_if_name);
|
||||
out_4ogl = new std::ofstream(res4openGL_fname, std::ios::out);
|
||||
assert(out_4ogl!=NULL);
|
||||
//if verbose only...
|
||||
if(verbose){
|
||||
verbose_fname = new char[strlen(w_if_name) + 9];// append .verb.txt
|
||||
sprintf(verbose_fname, "%s.verb.txt", w_if_name);
|
||||
out_verbose = new std::ofstream( verbose_fname, std::ios::out);
|
||||
assert(out_verbose != NULL);
|
||||
CGAL::set_pretty_mode(*out_verbose);
|
||||
}
|
||||
unsigned int nb_vertices_considered = 0;//count vertices for verbose
|
||||
// output
|
||||
|
||||
//load the model from <mesh.off>
|
||||
PolyhedralSurf P;
|
||||
std::ifstream stream(if_name);
|
||||
stream >> P;
|
||||
fprintf(stderr, "loadMesh %d Ves %d Facets\n",
|
||||
P.size_of_vertices(), P.size_of_facets());
|
||||
if(verbose)
|
||||
(*out_verbose) << "Polysurf with " << P.size_of_vertices()
|
||||
<< " vertices and " << P.size_of_facets()
|
||||
<< " facets. " << std::endl;
|
||||
|
||||
//exit if not enough points in the model
|
||||
unsigned int min_nb_points = (d_fitting + 1) * (d_fitting + 2) / 2;
|
||||
if (min_nb_points > P.size_of_vertices()) exit(0);
|
||||
//initialize Polyhedral data : length of edges, normal of facets
|
||||
P.compute_edges_length();
|
||||
P.compute_facets_normals();
|
||||
|
||||
//container for approximation points
|
||||
std::vector<DPoint> in_points;
|
||||
//container to collect vertices of the PolyhedralSurf
|
||||
std::vector<Vertex*> current_ring, next_ring, gathered;
|
||||
std::vector<Vertex*> *p_current_ring, *p_next_ring;
|
||||
|
||||
//MAIN LOOP
|
||||
Vertex_iterator vitb = P.vertices_begin(), vite = P.vertices_end();
|
||||
for (; vitb != vite; vitb++) {
|
||||
//initialize
|
||||
Vertex* v = &(*vitb);//handle to *
|
||||
unsigned int nbp = 0, //current nb of collected points
|
||||
ith = 0; //i-th ring index
|
||||
current_ring.clear();
|
||||
next_ring.clear();
|
||||
p_current_ring = ¤t_ring;
|
||||
p_next_ring = &next_ring;
|
||||
gathered.clear();
|
||||
in_points.clear();
|
||||
My_Monge_rep monge_rep;
|
||||
My_Monge_info monge_info;
|
||||
|
||||
//DO NOT FORGET TO UNTAG AT THE END!
|
||||
v->setRingIndex(ith);
|
||||
//collect 0th ring : the vertex v!
|
||||
gathered.push_back(v);
|
||||
nbp = 1;
|
||||
//collect 1-ring
|
||||
ith = 1;
|
||||
nbp = Poly_rings::push_neighbours_of(v, ith, current_ring, gathered);
|
||||
//collect more neighbors depending on options...
|
||||
|
||||
//OPTION -p nb, with nb != 0
|
||||
//for approximation/interpolation with a fixed nb of points, collect
|
||||
// enough rings and discard some points of the last collected ring
|
||||
// to get the exact "nb_points_to_use"
|
||||
if ( nb_points_to_use != 0 ) {
|
||||
while( gathered.size() < nb_points_to_use ) {
|
||||
ith++;
|
||||
//using tags
|
||||
nbp += Poly_rings::
|
||||
collect_ith_ring_neighbours(ith, *p_current_ring,
|
||||
*p_next_ring, gathered);
|
||||
//next round must be launched from p_nextRing...
|
||||
p_current_ring->clear();
|
||||
std::swap(p_current_ring, p_next_ring);
|
||||
}
|
||||
//clean up
|
||||
Poly_rings::reset_ring_indices(gathered);
|
||||
//discard non-required collected points of the last ring
|
||||
gathered.resize(nb_points_to_use, NULL);
|
||||
assert(gathered.size() == nb_points_to_use );
|
||||
}
|
||||
else{
|
||||
//OPTION -a nb, with nb = 0
|
||||
// select the mini nb of rings needed to make approx possible
|
||||
if (nb_rings == 0) {
|
||||
while ( gathered.size() < min_nb_points ) {
|
||||
ith++;
|
||||
nbp += Poly_rings::
|
||||
collect_ith_ring_neighbours(ith, *p_current_ring,
|
||||
*p_next_ring, gathered);
|
||||
//next round must be launched from p_nextRing...
|
||||
p_current_ring->clear();
|
||||
std::swap(p_current_ring, p_next_ring);
|
||||
}
|
||||
}
|
||||
//OPTION -a nb, with nb = 1, nothing to do! we have already
|
||||
// collected the 1 ring
|
||||
//OPTION -a nb, with nb > 1
|
||||
//for approximation with a fixed nb of rings, collect
|
||||
// neighbors up to the "nb_rings"th ring
|
||||
if (nb_rings > 1)
|
||||
while (ith < nb_rings) {
|
||||
ith++;
|
||||
nbp += Poly_rings::
|
||||
collect_ith_ring_neighbours(ith, *p_current_ring,
|
||||
*p_next_ring, gathered);
|
||||
//next round must be launched from p_nextRing...
|
||||
p_current_ring->clear();
|
||||
std::swap(p_current_ring, p_next_ring);
|
||||
}
|
||||
|
||||
//clean up
|
||||
Poly_rings::reset_ring_indices(gathered);
|
||||
} //END ELSE
|
||||
|
||||
//skip if the nb of gathered points is to small
|
||||
if ( gathered.size() < min_nb_points ) continue;
|
||||
|
||||
//store the gathered points
|
||||
std::vector<Vertex*>::iterator itb = gathered.begin(),
|
||||
ite = gathered.end();
|
||||
CGAL_For_all(itb,ite) in_points.push_back((*itb)->point());
|
||||
|
||||
// run the main fct : perform the fitting
|
||||
My_Monge_via_jet_fitting do_it(in_points.begin(), in_points.end(),
|
||||
d_fitting, d_monge,
|
||||
monge_rep, monge_info);
|
||||
|
||||
//switch min-max ppal curv/dir wrt the mesh orientation
|
||||
DVector normal_mesh = P.computeFacetsAverageUnitNormal(v);
|
||||
DVector normal_monge = monge_rep.n();
|
||||
if ( normal_mesh*normal_monge < 0 )
|
||||
{
|
||||
monge_rep.n() = -monge_rep.n();
|
||||
std::swap(monge_rep.d1(), monge_rep.d2());
|
||||
if ( d_monge >= 2)
|
||||
std::swap(monge_rep.coefficients()[0],monge_rep.coefficients()[1]);
|
||||
if ( d_monge >= 3) {
|
||||
std::swap(monge_rep.coefficients()[2],monge_rep.coefficients()[5]);
|
||||
std::swap(monge_rep.coefficients()[3],monge_rep.coefficients()[4]);}
|
||||
if ( d_monge >= 4) {
|
||||
std::swap(monge_rep.coefficients()[6],monge_rep.coefficients()[10]);
|
||||
std::swap(monge_rep.coefficients()[7],monge_rep.coefficients()[9]);}
|
||||
std::vector<DFT>::iterator itb = monge_rep.coefficients().begin(),
|
||||
ite = monge_rep.coefficients().end();
|
||||
for (;itb!=ite;itb++) { *itb = -(*itb); }
|
||||
}
|
||||
|
||||
//OpenGL output
|
||||
//scaling for ppal dir,
|
||||
// may be optimized with a global mean edges length computed only once
|
||||
// on all edges of P
|
||||
DFT scale_ppal_dir = P.compute_mean_edges_length_around_vertex(v)/2;
|
||||
(*out_4ogl) << v->point() << " "
|
||||
<< monge_rep.origin_pt() << " "
|
||||
<< monge_rep.d1() * scale_ppal_dir << " "
|
||||
<< monge_rep.d2() * scale_ppal_dir << " "
|
||||
<< monge_rep.coefficients()[0] << " "
|
||||
<< monge_rep.coefficients()[1] << " "
|
||||
<< std::endl;
|
||||
|
||||
//verbose txt output
|
||||
if (verbose){
|
||||
(*out_verbose) << "--- vertex " << ++nb_vertices_considered
|
||||
<< " : " << v->point() << std::endl
|
||||
<< "number of points used : " << in_points.size() << std::endl
|
||||
<< "number of rings used : " << ith << std::endl
|
||||
<< "origin : " << monge_rep.origin_pt() << std::endl
|
||||
<< "d1 : " << monge_rep.d1() << std::endl
|
||||
<< "d2 : " << monge_rep.d2() << std::endl
|
||||
<< "n : " << monge_rep.n() << std::endl
|
||||
<< "cond_nb : " << monge_info.cond_nb() << std::endl
|
||||
<< "pca_eigen_vals " << monge_info.pca_eigen_vals()[0]
|
||||
<< " " << monge_info.pca_eigen_vals()[2]
|
||||
<< " " << monge_info.pca_eigen_vals()[3] << std::endl
|
||||
<< "pca_eigen_vecs : " << std::endl
|
||||
<< monge_info.pca_eigen_vecs()[0] << std::endl
|
||||
<< monge_info.pca_eigen_vecs()[1] << std::endl
|
||||
<< monge_info.pca_eigen_vecs()[2] << std::endl;
|
||||
if ( d_monge >= 2)
|
||||
(*out_verbose) << "k1 : " << monge_rep.coefficients()[0] << std::endl
|
||||
<< "k2 : " << monge_rep.coefficients()[1] << std::endl
|
||||
<< std::endl ;
|
||||
} //END IF
|
||||
} //END FOR LOOP
|
||||
|
||||
///////////////////////////////////////////
|
||||
//RIDGES!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
///////////////////////////////////////////
|
||||
Ridge_approximation ridge_approximation;
|
||||
std::list<Ridge_line> ridge_lines;
|
||||
// std::back_insert_iterator<std::list<Ridge_line> > iter_lines(ridge_lines), iter_end;
|
||||
std::list<Ridge_line>::iterator iter_lines = ridge_lines.begin(), iter_end;
|
||||
iter_end = ridge_approximation.do_it(P, iter_lines);
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// std::cout << "ok 1 " << std::endl;
|
||||
//cleanup filenames
|
||||
delete res4openGL_fname;
|
||||
out_4ogl->close();
|
||||
delete out_4ogl;
|
||||
if(verbose) {
|
||||
delete verbose_fname;
|
||||
out_verbose->close();
|
||||
delete out_verbose;
|
||||
}
|
||||
// std::cout << "ok 2 " << std::endl ;
|
||||
return 1;
|
||||
}
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,492 @@
|
|||
// ****************************************************************************
|
||||
// ^FILE: options.h - option parsing classes
|
||||
//
|
||||
// ^DESCRIPTION:
|
||||
// This file defines classes used to parse command-line options.
|
||||
// Options may be parsed from an array of strings, or from any structure
|
||||
// for which a corresponding option-iterator exists.
|
||||
//
|
||||
// ^HISTORY:
|
||||
// 03/06/92 Brad Appleton <bradapp@enteract.com> Created
|
||||
//
|
||||
// 03/23/93 Brad Appleton <bradapp@enteract.com>
|
||||
// - Added OptIstreamIter class
|
||||
//
|
||||
// 03/08/94 Brad Appleton <bradapp@enteract.com>
|
||||
// - Added Options::reset() member function
|
||||
//
|
||||
// 07/31/97 Brad Appleton <bradapp@enteract.com>
|
||||
// - Added PARSE_POS control flag and POSITIONAL return value
|
||||
// ^^**************************************************************************
|
||||
|
||||
#ifndef _options_h
|
||||
#define _options_h
|
||||
|
||||
using namespace std;
|
||||
#include <iostream>
|
||||
|
||||
|
||||
// Abstract class to iterate through options/arguments
|
||||
//
|
||||
class OptIter {
|
||||
public:
|
||||
OptIter(void) {}
|
||||
virtual ~OptIter(void);
|
||||
|
||||
// curr() returns the current item in the iterator without
|
||||
// advancing on to the next item. If we are at the end of items
|
||||
// then NULL is returned.
|
||||
virtual const char *
|
||||
curr(void) = 0;
|
||||
|
||||
// next() advances to the next item.
|
||||
virtual void
|
||||
next(void) = 0;
|
||||
|
||||
// operator() returns the current item in the iterator and then
|
||||
// advances on to the next item. If we are at the end of items
|
||||
// then NULL is returned.
|
||||
virtual const char *
|
||||
operator()(void);
|
||||
} ;
|
||||
|
||||
// Abstract class for a rewindable OptIter
|
||||
//
|
||||
class OptIterRwd : public OptIter {
|
||||
public:
|
||||
OptIterRwd(void);
|
||||
|
||||
virtual ~OptIterRwd(void);
|
||||
|
||||
virtual const char *
|
||||
curr(void) = 0;
|
||||
|
||||
virtual void
|
||||
next(void) = 0;
|
||||
|
||||
virtual const char *
|
||||
operator()(void) = 0;
|
||||
|
||||
// rewind() resets the "current-element" to the first one in the "list"
|
||||
virtual void
|
||||
rewind(void) = 0;
|
||||
} ;
|
||||
|
||||
// Class to iterate through an array of tokens. The array may be terminated
|
||||
// by NULL or a count containing the number of tokens may be given.
|
||||
//
|
||||
class OptArgvIter : public OptIterRwd {
|
||||
private:
|
||||
const char * const * av; // arg vector
|
||||
int ac; // arg count
|
||||
int ndx; // index of current arg
|
||||
|
||||
public:
|
||||
OptArgvIter(const char * const argv[])
|
||||
: av(argv), ac(-1), ndx(0) {}
|
||||
|
||||
OptArgvIter(int argc, const char * const argv[])
|
||||
: av(argv), ac(argc), ndx(0) {init(argc, argv); }
|
||||
|
||||
void init(int argc, const char * const argv[]);// {v=argv; ac=argc; ndx=0;}
|
||||
|
||||
virtual
|
||||
~OptArgvIter(void);
|
||||
|
||||
virtual const char *
|
||||
curr(void);
|
||||
|
||||
virtual void
|
||||
next(void);
|
||||
|
||||
virtual const char *
|
||||
operator()(void);
|
||||
|
||||
virtual void
|
||||
rewind(void);
|
||||
|
||||
// index returns the current index to use for argv[]
|
||||
int index(void) { return ndx; }
|
||||
} ;
|
||||
|
||||
|
||||
// Class to iterate through a string containing delimiter-separated tokens
|
||||
//
|
||||
class OptStrTokIter : public OptIterRwd {
|
||||
private:
|
||||
unsigned len; // length of token-string
|
||||
const char * str; // the token-string
|
||||
const char * seps; // delimiter-set (separator-characters)
|
||||
const char * cur; // current token
|
||||
char * tokstr; // our copy of the token-string
|
||||
|
||||
static const char * default_delims; // default delimiters = whitespace
|
||||
|
||||
public:
|
||||
OptStrTokIter(const char * tokens, const char * delimiters =0);
|
||||
|
||||
virtual
|
||||
~OptStrTokIter(void);
|
||||
|
||||
virtual const char *
|
||||
curr(void);
|
||||
|
||||
virtual void
|
||||
next(void);
|
||||
|
||||
virtual const char *
|
||||
operator()(void);
|
||||
|
||||
virtual void
|
||||
rewind(void);
|
||||
|
||||
// delimiters() with NO arguments returns the current set of delimiters,
|
||||
// If an argument is given then it is used as the new set of delimiters.
|
||||
const char *
|
||||
delimiters(void) { return seps; }
|
||||
|
||||
void
|
||||
delimiters(const char * delims) {
|
||||
seps = (delims) ? delims : default_delims ;
|
||||
}
|
||||
} ;
|
||||
|
||||
|
||||
// OptIstreamIter is a class for iterating over arguments that come
|
||||
// from an input stream. Each line of the input stream is considered
|
||||
// to be a set of white-space separated tokens. If the the first
|
||||
// non-white character on a line is '#' ('!' for VMS systems) then
|
||||
// the line is considered a comment and is ignored.
|
||||
//
|
||||
// *Note:: If a line is more than 1022 characters in length then we
|
||||
// treat it as if it were several lines of length 1022 or less.
|
||||
//
|
||||
// *Note:: The string tokens returned by this iterator are pointers
|
||||
// to temporary buffers which may not necessarily stick around
|
||||
// for too long after the call to curr() or operator(), hence
|
||||
// if you need the string value to persist - you will need to
|
||||
// make a copy.
|
||||
//
|
||||
class OptIstreamIter : public OptIter {
|
||||
private:
|
||||
istream & is ;
|
||||
OptStrTokIter * tok_iter ;
|
||||
|
||||
void
|
||||
fill(void);
|
||||
|
||||
public:
|
||||
static const unsigned MAX_LINE_LEN ;
|
||||
|
||||
OptIstreamIter(istream & input);
|
||||
|
||||
virtual
|
||||
~OptIstreamIter(void);
|
||||
|
||||
virtual const char *
|
||||
curr(void);
|
||||
|
||||
virtual void
|
||||
next(void);
|
||||
|
||||
virtual const char *
|
||||
operator()(void);
|
||||
} ;
|
||||
|
||||
|
||||
// Now we are ready to define a class to declare and parse command-options
|
||||
//
|
||||
//
|
||||
// CLASS
|
||||
// =====
|
||||
// Options -- parse command-line options
|
||||
//
|
||||
// SYNOPSIS
|
||||
// ========
|
||||
// #include <options.h>
|
||||
//
|
||||
// Options opts(cmdname, optv);
|
||||
// char cmdname[], *optv[];
|
||||
//
|
||||
// DESCRIPTION
|
||||
// ===========
|
||||
// The Options constructor expects a command-name (usually argv[0]) and
|
||||
// a pointer to an array of strings. The last element in this array MUST
|
||||
// be NULL. Each non-NULL string in the array must have the following format:
|
||||
//
|
||||
// The 1st character must be the option-name ('c' for a -c option).
|
||||
//
|
||||
// The 2nd character must be one of '|', '?', ':', '*', or '+'.
|
||||
// '|' -- indicates that the option takes NO argument;
|
||||
// '?' -- indicates that the option takes an OPTIONAL argument;
|
||||
// ':' -- indicates that the option takes a REQUIRED argument;
|
||||
// '*' -- indicates that the option takes 0 or more arguments;
|
||||
// '+' -- indicates that the option takes 1 or more arguments;
|
||||
//
|
||||
// The remainder of the string must be the long-option name.
|
||||
//
|
||||
// If desired, the long-option name may be followed by one or more
|
||||
// spaces and then by the name of the option value. This name will
|
||||
// be used when printing usage messages. If the option-value-name
|
||||
// is not given then the string "<value>" will be used in usage
|
||||
// messages.
|
||||
//
|
||||
// One may use a space to indicate that a particular option does not
|
||||
// have a corresponding long-option. For example, "c: " (or "c:")
|
||||
// means the -c option takes a value & has NO corresponding long-option.
|
||||
//
|
||||
// To specify a long-option that has no corresponding single-character
|
||||
// option is a bit trickier: Options::operator() still needs an "option-
|
||||
// character" to return when that option is matched. One may use a whitespace
|
||||
// character or a non-printable character as the single-character option
|
||||
// in such a case. (hence " |hello" would only match "--hello").
|
||||
//
|
||||
// EXCEPTIONS TO THE ABOVE:
|
||||
// ------------------------
|
||||
// If the 1st character of the string is '-', then the rest of the
|
||||
// string must correspond to the above format, and the option is
|
||||
// considered to be a hidden-option. This means it will be parsed
|
||||
// when actually matching options from the command-line, but will
|
||||
// NOT show-up if a usage message is printed using the usage() member
|
||||
// function. Such an example might be "-h|hidden". If you want to
|
||||
// use any "dummy" options (options that are not parsed, but that
|
||||
// to show up in the usage message), you can specify them along with
|
||||
// any positional parameters to the usage() member function.
|
||||
//
|
||||
// If the 2nd character of the string is '\0' then it is assumed
|
||||
// that there is no corresponding long-option and that the option
|
||||
// takes no argument (hence "f", and "f| " are equivalent).
|
||||
//
|
||||
// Examples:
|
||||
// const char * optv[] = {
|
||||
// "c:count <number>",
|
||||
// "s?str <string>",
|
||||
// "x",
|
||||
// " |hello",
|
||||
// "g+groups <newsgroup>",
|
||||
// NULL
|
||||
// } ;
|
||||
//
|
||||
// optv[] now corresponds to the following:
|
||||
//
|
||||
// usage: cmdname [-c|--count <number>] [-s|--str [<string>]]
|
||||
// [-x] [--hello] [-g|--groups <newsgroup> ...]
|
||||
//
|
||||
// Long-option names are matched case-insensitive and only a unique prefix
|
||||
// of the name needs to be specified.
|
||||
//
|
||||
// Option-name characters are case-sensitive!
|
||||
//
|
||||
// CAVEAT
|
||||
// ======
|
||||
// Because of the way in which multi-valued options and options with optional
|
||||
// values are handled, it is NOT possible to supply a value to an option in
|
||||
// a separate argument (different argv[] element) if the value is OPTIONAL
|
||||
// and begins with a '-'. What this means is that if an option "-s" takes an
|
||||
// optional value value and you wish to supply a value of "-foo" then you must
|
||||
// specify this on the command-line as "-s-foo" instead of "-s -foo" because
|
||||
// "-s -foo" will be considered to be two separate sets of options.
|
||||
//
|
||||
// A multi-valued option is terminated by another option or by the end-of
|
||||
// options. The following are all equivalent (if "-l" is a multi-valued
|
||||
// option and "-x" is an option that takes no value):
|
||||
//
|
||||
// cmdname -x -l item1 item2 item3 -- arg1 arg2 arg3
|
||||
// cmdname -x -litem1 -litem2 -litem3 -- arg1 arg2 arg3
|
||||
// cmdname -l item1 item2 item3 -x arg1 arg2 arg3
|
||||
//
|
||||
//
|
||||
// EXAMPLE
|
||||
// =======
|
||||
// #include <options.h>
|
||||
//
|
||||
// static const char * optv[] = {
|
||||
// "H|help",
|
||||
// "c:count <number>",
|
||||
// "s?str <string>",
|
||||
// "x",
|
||||
// " |hello",
|
||||
// "g+groups <newsgroup>",
|
||||
// NULL
|
||||
// } ;
|
||||
//
|
||||
// main(int argc, char * argv[]) {
|
||||
// int optchar;
|
||||
// const char * optarg;
|
||||
// const char * str = "default_string";
|
||||
// int count = 0, xflag = 0, hello = 0;
|
||||
// int errors = 0, ngroups = 0;
|
||||
//
|
||||
// Options opts(*argv, optv);
|
||||
// OptArgvIter iter(--argc, ++argv);
|
||||
//
|
||||
// while( optchar = opts(iter, optarg) ) {
|
||||
// switch (optchar) {
|
||||
// case 'H' :
|
||||
// opts.usage(cout, "files ...");
|
||||
// exit(0);
|
||||
// break;
|
||||
// case 'g' :
|
||||
// ++ngroups; break; // the groupname is in "optarg"
|
||||
// case 's' :
|
||||
// str = optarg; break;
|
||||
// case 'x' :
|
||||
// ++xflag; break;
|
||||
// case ' ' :
|
||||
// ++hello; break;
|
||||
// case 'c' :
|
||||
// if (optarg == NULL) ++errors;
|
||||
// else count = (int) atol(optarg);
|
||||
// break;
|
||||
// default : ++errors; break;
|
||||
// } //switch
|
||||
// }
|
||||
//
|
||||
// if (errors || (iter.index() == argc)) {
|
||||
// if (! errors) {
|
||||
// cerr << opts.name() << ": no filenames given." << endl ;
|
||||
// }
|
||||
// opts.usage(cerr, "files ...");
|
||||
// exit(1);
|
||||
// }
|
||||
//
|
||||
// cout << "xflag=" << ((xflag) ? "ON" : "OFF") << endl
|
||||
// << "hello=" << ((hello) ? "YES" : "NO") << endl
|
||||
// << "count=" << count << endl
|
||||
// << "str=\"" << ((str) ? str : "No value given!") << "\"" << endl
|
||||
// << "ngroups=" << ngroups << endl ;
|
||||
//
|
||||
// if (iter.index() < argc) {
|
||||
// cout << "files=" ;
|
||||
// for (int i = iter.index() ; i < argc ; i++) {
|
||||
// cout << "\"" << argv[i] << "\" " ;
|
||||
// }
|
||||
// cout << endl ;
|
||||
// }
|
||||
// }
|
||||
//
|
||||
class Options {
|
||||
private:
|
||||
const char * cmdname; // name of the command
|
||||
const char * const * optvec; // vector of option-specifications (last=NULL)
|
||||
unsigned explicit_end : 1; // were we terminated because of "--"?
|
||||
unsigned optctrls : 7; // control settings (a set of OptCtrl masks)
|
||||
const char * nextchar; // next option-character to process
|
||||
const char * listopt; // last list-option we matched
|
||||
|
||||
void
|
||||
check_syntax(void) const;
|
||||
|
||||
const char *
|
||||
match_opt(char opt, int ignore_case =0) const;
|
||||
|
||||
const char *
|
||||
match_longopt(const char * opt, int len, int & ambiguous) const;
|
||||
|
||||
int
|
||||
parse_opt(OptIter & iter, const char * & optarg);
|
||||
|
||||
int
|
||||
parse_longopt(OptIter & iter, const char * & optarg);
|
||||
|
||||
public:
|
||||
enum OptCtrl {
|
||||
DEFAULT = 0x00, // Default setting
|
||||
ANYCASE = 0x01, // Ignore case when matching short-options
|
||||
QUIET = 0x02, // Dont print error messages
|
||||
PLUS = 0x04, // Allow "+" as a long-option prefix
|
||||
SHORT_ONLY = 0x08, // Dont accept long-options
|
||||
LONG_ONLY = 0x10, // Dont accept short-options
|
||||
// (also allows "-" as a long-option prefix).
|
||||
NOGUESSING = 0x20, // Normally, when we see a short (long) option
|
||||
// on the command line that doesnt match any
|
||||
// known short (long) options, then we try to
|
||||
// "guess" by seeing if it will match any known
|
||||
// long (short) option. Setting this mask prevents
|
||||
// this "guessing" from occurring.
|
||||
PARSE_POS = 0x40 // By default, Options will not present positional
|
||||
// command-line arguments to the user and will
|
||||
// instead stop parsing when the first positonal
|
||||
// argument has been encountered. If this flag
|
||||
// is given, Options will present positional
|
||||
// arguments to the user with a return code of
|
||||
// POSITIONAL; ENDOPTS will be returned only
|
||||
// when the end of the argument list is reached.
|
||||
} ;
|
||||
|
||||
// Error return values for operator()
|
||||
//
|
||||
enum OptRC {
|
||||
ENDOPTS = 0,
|
||||
BADCHAR = -1,
|
||||
BADKWD = -2,
|
||||
AMBIGUOUS = -3,
|
||||
POSITIONAL = -4
|
||||
} ;
|
||||
|
||||
Options(const char * name, const char * const optv[]);
|
||||
|
||||
virtual
|
||||
~Options(void);
|
||||
|
||||
// name() returns the command name
|
||||
const char *
|
||||
name(void) const { return cmdname; }
|
||||
|
||||
// ctrls() (with no arguments) returns the existing control settings
|
||||
unsigned
|
||||
ctrls(void) const { return optctrls; }
|
||||
|
||||
// ctrls() (with 1 argument) sets new control settings
|
||||
void
|
||||
ctrls(unsigned newctrls) { optctrls = newctrls; }
|
||||
|
||||
// reset for another pass to parse for options
|
||||
void
|
||||
reset(void) { nextchar = listopt = NULL; }
|
||||
|
||||
// usage() prints options usage (followed by any positional arguments
|
||||
// listed in the parameter "positionals") on the given outstream
|
||||
void
|
||||
usage(ostream & os, const char * positionals) const ;
|
||||
|
||||
// operator() iterates through the arguments as necessary (using the
|
||||
// given iterator) and returns the character value of the option
|
||||
// (or long-option) that it matched. If the option has a value
|
||||
// then the value given may be found in optarg (otherwise optarg
|
||||
// will be NULL).
|
||||
//
|
||||
// 0 is returned upon end-of-options. At this point, "iter" may
|
||||
// be used to process any remaining positional parameters. If the
|
||||
// PARSE_POS control-flag is set then 0 is returned only when all
|
||||
// arguments in "iter" have been exhausted.
|
||||
//
|
||||
// If an invalid option is found then BADCHAR is returned and *optarg
|
||||
// is the unrecognized option character.
|
||||
//
|
||||
// If an invalid long-option is found then BADKWD is returned and optarg
|
||||
// points to the bad long-option.
|
||||
//
|
||||
// If an ambiguous long-option is found then AMBIGUOUS is returned and
|
||||
// optarg points to the ambiguous long-option.
|
||||
//
|
||||
// If the PARSE_POS control-flag is set then POSITIONAL is returned
|
||||
// when a positional argument is encountered and optarg points to
|
||||
// the positonal argument (and "iter" is advanced to the next argument
|
||||
// in the iterator).
|
||||
//
|
||||
// Unless Options::QUIET is used, missing option-arguments and
|
||||
// invalid options (and the like) will automatically cause error
|
||||
// messages to be issued to cerr.
|
||||
int
|
||||
operator()(OptIter & iter, const char * & optarg) ;
|
||||
|
||||
// Call this member function after operator() has returned 0
|
||||
// if you want to know whether or not options were explicitly
|
||||
// terminated because "--" appeared on the command-line.
|
||||
//
|
||||
int
|
||||
explicit_endopts() const { return explicit_end; }
|
||||
} ;
|
||||
|
||||
#endif /* _options_h */
|
||||
|
|
@ -1,379 +0,0 @@
|
|||
#ifndef _RIDGE_3_H_
|
||||
#define _RIDGE_3_H_
|
||||
|
||||
#include <pair>
|
||||
|
||||
|
||||
//Orient monge normal according to mesh normals.
|
||||
|
||||
CGAL_BEGIN_NAMESPACE
|
||||
|
||||
|
||||
enum Ridge_type {NONE=0, BLUE, RED, CREST, BE, BH, BC, RE, RH, RC};
|
||||
|
||||
|
||||
//---------------------------------------
|
||||
//Ridge_line : a connected sequence of edges crossed by a ridge, with type and weigths
|
||||
//---------------------------------------
|
||||
template < class Poly > class Ridge_line
|
||||
{
|
||||
public:
|
||||
typedef typename Poly::Traits::FT FT;
|
||||
typedef typename Poly::Vertex_handle Vertex_handle;
|
||||
typedef typename Poly::Halfedge_handle Halfedge_handle;
|
||||
typedef typename Poly::Facet_handle Facet_handle;
|
||||
typedef std::pair<Halfedge_handle, FT> ridge_he;
|
||||
|
||||
protected:
|
||||
Ridge_type line_type
|
||||
std::list<ridge_he> line;
|
||||
FT strength;
|
||||
FT sharpness;
|
||||
|
||||
public:
|
||||
const FT weight() const {return weight;}
|
||||
const FT sharpness() const {return sharpness;}
|
||||
std::list<ridge_he>* line() { return &line;}
|
||||
|
||||
//constructor
|
||||
Ridge_line( Halfedge_handle h1, Halfedge_handle h2, Ridge_type r_type) :
|
||||
line_type(r_type), strength(0.)
|
||||
{};
|
||||
|
||||
//compute the barycentric coordinate of the xing point (blue or red)
|
||||
//for he: p->q coord is st xing_point = coord*p + (1-coord)*q
|
||||
FT bary_coord(Halfedge_handle he);
|
||||
void compute_weight(char color);
|
||||
void compute_sharpness(char color);
|
||||
//When the line is extended with a he, the bary coord of the
|
||||
//crossing point is computed, the pair (he,coord) is added and the
|
||||
//weigths are updated
|
||||
void addback( Halfedge_handle he);
|
||||
void addfront( Halfedge_handle he);
|
||||
|
||||
};
|
||||
|
||||
// IMPLEMENTATION OF Ridge_line members
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
//constructor
|
||||
template < class Poly >
|
||||
Ridge_line( Halfedge_handle h1, Halfedge_handle h2, Ridge_type r_type) :
|
||||
line_type(r_type), strength(0.)
|
||||
{
|
||||
line.push_back(h1);
|
||||
addback(h2);
|
||||
}
|
||||
|
||||
template < class Poly >
|
||||
void Ridge_line<Poly>::
|
||||
addback( Halfedge_handle he)
|
||||
{
|
||||
Halfedge_handle he_cur = ( --(line.end()) )->first;
|
||||
FT coord = bary_coord(he);
|
||||
FT coord_cur = bary_coord(he_cur);
|
||||
Vertex_handle v_p = he->opposite()->vertex(), v_q = he->vertex(),
|
||||
v_p_cur = he_cur->opposite()->vertex(), v_q_cur = he->vertex(); // he: p->q
|
||||
FT k;
|
||||
if ( (line_type == BE) || (line_type == BH) || (line_type == BC) ) {
|
||||
k =( std::fabs(v_p->k1) * coord + std::fabs(v_q->k1) * (1-coord) )/2;
|
||||
}
|
||||
if ( (line_type == RE) || (line_type == RH) || (line_type == RC) ) {
|
||||
k =( std::fabs(v_p->k2) * coord + std::fabs(v_q->k2) * (1-coord) )/2;
|
||||
}
|
||||
Vector_3 segment = (v_p->point()-ORIGIN)*coord + (v_q->point()-ORIGIN)*(1-coord) -
|
||||
((v_p_cur->point()-ORIGIN)*coord_cur + (v_q_cur->point()-ORIGIN)*(1-coord_cur));
|
||||
strength += k * CGAL::sqrt(segment * segment);
|
||||
|
||||
line.push_back( pair(he, coord));
|
||||
}
|
||||
|
||||
template < class Poly >
|
||||
void Ridge_line<Poly>::
|
||||
addfront( Halfedge_handle he)
|
||||
{
|
||||
Halfedge_handle he_cur = ( line.begin() )->first;
|
||||
FT coord = bary_coord(he);
|
||||
FT coord_cur = bary_coord(he_cur);
|
||||
Vertex_handle v_p = he->opposite()->vertex(), v_q = he->vertex(),
|
||||
v_p_cur = he_cur->opposite()->vertex(), v_q_cur = he->vertex(); // he: p->q
|
||||
FT k;
|
||||
if ( (line_type == BE) || (line_type == BH) || (line_type == BC) ) {
|
||||
k =( std::fabs(v_p->k1) * coord + std::fabs(v_q->k1) * (1-coord) )/2;
|
||||
}
|
||||
if ( (line_type == RE) || (line_type == RH) || (line_type == RC) ) {
|
||||
k =( std::fabs(v_p->k2) * coord + std::fabs(v_q->k2) * (1-coord) )/2;
|
||||
}
|
||||
Vector_3 segment = (v_p->point()-ORIGIN)*coord + (v_q->point()-ORIGIN)*(1-coord) -
|
||||
((v_p_cur->point()-ORIGIN)*coord_cur + (v_q_cur->point()-ORIGIN)*(1-coord_cur));
|
||||
strength += k * CGAL::sqrt(segment * segment);
|
||||
|
||||
line.push_front( pair(he, coord));
|
||||
}
|
||||
|
||||
|
||||
template < class Poly >
|
||||
FT Ridge_line<Poly>::
|
||||
bary_coord(Halfedge_handle he)
|
||||
{
|
||||
FT b_p, b_q; // extremalities at p and q for he: p->q
|
||||
if ( (line_type == BE) || (line_type == BH) || (line_type == BC) ) {
|
||||
b_p = he->opposite()->vertex()->b0();
|
||||
b_q = he->vertex()->b0();
|
||||
}
|
||||
if ( (line_type == RE) || (line_type == RH) || (line_type == RC) ) {
|
||||
b_p = he->opposite()->vertex()->b3();
|
||||
b_q = he->vertex()->b3();
|
||||
}
|
||||
return std::fabs(b_q) / ( std::fabs(b_q) + std::fabs(b_p) );
|
||||
}
|
||||
|
||||
|
||||
/////////Ridge_approximation//////////////////////////////////////////
|
||||
|
||||
//Find BLUE ridges (Elliptic or Hyperbolic)
|
||||
//do the same for RED and CREST
|
||||
//iterate on P facets, find a non-visited, regular, 2BXing triangle,
|
||||
//follow non-visited, regular, 2BXing triangles in both sens to create
|
||||
//a Ridge line.
|
||||
//Each time a edge is added the strength of the current line is updated
|
||||
// + length(edge)*|k|
|
||||
void
|
||||
compute_ridges(Ridge_type r_type)
|
||||
{
|
||||
//set all facets non visited
|
||||
|
||||
|
||||
Facet_iterator itb = P.facets_begin(), ite = P.facets_end();
|
||||
for(;itb!=ite;itb++)
|
||||
{
|
||||
Facet_handle f= &(*itb);
|
||||
if (f->is_visited()) continue;
|
||||
f->set_visited(true);
|
||||
Halfedge_handle h1, h2, curhe1, curhe2, curhe;
|
||||
|
||||
//h1 h2 are the hedges crossed if any, r_type should be BLUE,
|
||||
//RED or CREST ; cur_ridge_type should be BE, BH, BC, RE, RH, RC or NONE
|
||||
Ridge_type cur_ridge_type = facet_ridge_type(f,h1,h2,r_type)
|
||||
if ( cur_ridge_type == NONE ) continue;
|
||||
|
||||
//When a he is inserted in a Ridge_line, a pair <he,bary_coord>
|
||||
//is created and the weight(s) are updated
|
||||
Ridge_line *cur_ridge_line = new Ridge_line(h1,h2,cur_ridge_type);
|
||||
*ridge_lines_it++ = cur_ridge_line;
|
||||
|
||||
//next triangle adjacent to h1 (push_front)
|
||||
if ( !(h1->is_border_edge()) )
|
||||
{
|
||||
f = h1->opposite()->facet();
|
||||
curhe = h1;
|
||||
while (cur_ridge_type == facet_ridge_type(f,curhe1,curhe2,r_type))
|
||||
{
|
||||
//follow the ridge from curhe
|
||||
if (f->is_visited()) break;
|
||||
f->set_visited(true);
|
||||
if (curhe->opposite() == curhe1) curhe = curhe2;
|
||||
else curhe = curhe1;
|
||||
cur_ridge_line->addfront(curhe);
|
||||
if ( !(curhe->is_border_edge()) ) f =
|
||||
curhe->opposite()->facet();
|
||||
else break;
|
||||
}
|
||||
//exit from the while if
|
||||
//1. border
|
||||
//2. not same type, then do not set visisted cause a BE
|
||||
// follows a BH
|
||||
}
|
||||
|
||||
//next triangle adjacent to h2 (push_back)
|
||||
if ( !(h2->is_border_edge()) )
|
||||
{
|
||||
f = h2->opposite()->facet();
|
||||
curhe = h2;
|
||||
while (cur_ridge_type == facet_ridge_type(f,curhe1,curhe2,r_type))
|
||||
{
|
||||
//follow the ridge from curhe
|
||||
if (f->is_visitedB()) break;
|
||||
f->set_visitedB(true);
|
||||
if (curhe->opposite() == curhe1) curhe = curhe2;
|
||||
else curhe = curhe1;
|
||||
cur_ridge_line->addback(curhe);
|
||||
if ( !(curhe->is_border_edge()) ) f =
|
||||
curhe->opposite()->facet();
|
||||
else break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Ridge_type
|
||||
facet_ridge_type(Facet_handle f, Halfedge_handle& he1, Halfedge_handle&
|
||||
he2, Ridge_type r_type)
|
||||
{
|
||||
//polyhedral data
|
||||
//we have v1--h1-->v2--h2-->v3--h3-->v1
|
||||
Halfedge_handle h1 = f->halfedge(), h2, h3;
|
||||
Vertex_handle v1, v2, v3;
|
||||
v2 = h1->vertex();
|
||||
h2 = h1->next();
|
||||
v3 = h2->vertex();
|
||||
h3 = h2->next();
|
||||
v1 = h3->vertex();
|
||||
|
||||
//check for regular facet
|
||||
if ( v1->d1()*v2->d1() * v1->d1()*v3->d1() * v2->d1()*v3->d1() < 0 ) return NONE;
|
||||
|
||||
//determine potential crest color
|
||||
Ridge_type crest_color = NONE;
|
||||
if (r_type == CREST)
|
||||
{
|
||||
if ( std::fabs(v1->k1()+v2->k1()+v3->k1()) > std::fabs(v1->k2()+v2->k2()+v3->k2()) )
|
||||
crest_color = BC;
|
||||
if ( std::fabs(v1->k1()+v2->k1()+v3->k1()) < std::fabs(v1->k2()+v2->k2()+v3->k2()) )
|
||||
crest_color = RC;
|
||||
if ( std::fabs(v1->k1()+v2->k1()+v3->k1()) = std::fabs(v1->k2()+v2->k2()+v3->k2()) )
|
||||
return NONE;
|
||||
}
|
||||
|
||||
//compute Xing on the 3 edges
|
||||
bool h1_is_crossed, h2_is_crossed, h3_is_crossed;
|
||||
if ( r_type == BLUE || crest_color == BC )
|
||||
{
|
||||
xing_on_edge(h1, h1_is_crossed, BLUE);
|
||||
xing_on_edge(h2, h2_is_crossed, BLUE);
|
||||
xing_on_edge(h3, h3_is_crossed, BLUE);
|
||||
}
|
||||
if ( r_type == RED || crest_color == RC )
|
||||
{
|
||||
xing_on_edge(h1, h1_is_crossed, RED);
|
||||
xing_on_edge(h2, h2_is_crossed, RED);
|
||||
xing_on_edge(h3, h3_is_crossed, RED);
|
||||
}
|
||||
|
||||
//test of 2Xing
|
||||
if ( !(h1_is_crossed && h2_is_crossed && !h3_is_crossed) &&
|
||||
!(h1_is_crossed && !h2_is_crossed && h3_is_crossed) &&
|
||||
(!h1_is_crossed && h2_is_crossed && h3_is_crossed) )
|
||||
return NONE;
|
||||
if (h1_is_crossed && h2_is_crossed && !h3_is_crossed)
|
||||
{
|
||||
he1 = h1;
|
||||
he2 = h2;
|
||||
}
|
||||
if (h1_is_crossed && !h2_is_crossed && h3_is_crossed)
|
||||
{
|
||||
he1 = h1;
|
||||
he2 = h3;
|
||||
}
|
||||
if (!h1_is_crossed && h2_is_crossed && h3_is_crossed)
|
||||
{
|
||||
he1 = h2;
|
||||
he2 = h3;
|
||||
}
|
||||
|
||||
Vertex_handle v_p1 = he1->opposite()->vertex(), v_q1 = he1->vertex(),
|
||||
v_p2 = he2->opposite()->vertex(), v_q2 = he2->vertex(); // he1: p1->q1
|
||||
|
||||
if ( r_type == BLUE || crest_color == BC )
|
||||
{
|
||||
FT coord1 = std::fabs(v_q1->b0()) / ( std::fabs(v_p1->b0()) + std::fabs(v_q1->b0()) ),
|
||||
coord2 = std::fabs(v_q2->b0()) / ( std::fabs(v_p2->b0()) + std::fabs(v_q2->b0()) );
|
||||
Vector_3 r1 = (v_p1->point()-ORIGIN)*coord1 + (v_q1->point()-ORIGIN)*(1-coord1),
|
||||
r2 = (v_p2->point()-ORIGIN)*coord2 + (v_q2->point()-ORIGIN)*(1-coord2);
|
||||
int b_sign = b_sign_pointing_to_ridge(v1, v2, v3, r1, r2, BLUE);
|
||||
if (b_sign == 1) { if (r_type == BLUE) return BE; else return BC;}
|
||||
if (b_sign == -1) return BH;
|
||||
}
|
||||
if ( r_type == RED || crest_color == RC )
|
||||
{
|
||||
FT coord1 = std::fabs(v_q1->b3()) / ( std::fabs(v_p1->b3()) + std::fabs(v_q1->b3()) ),
|
||||
coord2 = std::fabs(v_q2->b3()) / ( std::fabs(v_p2->b3()) + std::fabs(v_q2->b3()) );
|
||||
Vector_3 r1 = (v_p1->point()-ORIGIN)*coord1 + (v_q1->point()-ORIGIN)*(1-coord1),
|
||||
r2 = (v_p2->point()-ORIGIN)*coord2 + (v_q2->point()-ORIGIN)*(1-coord2);
|
||||
int b_sign = b_sign_pointing_to_ridge(v1, v2, v3, r1, r2, RED);
|
||||
if (b_sign == -1) { if (r_type == RED) return RE; else return RC;}
|
||||
if (b_sign == 1) return BH;
|
||||
}
|
||||
assert(0);
|
||||
}
|
||||
|
||||
void
|
||||
xing_on_edge(Halfedge_handle he, bool& is_crossed, Ridge_type color)
|
||||
{
|
||||
is_crossed = false;
|
||||
FT sign;
|
||||
FT b_p, b_q; // extremalities at p and q for he: p->q
|
||||
Vector_3 d_p = he->opposite()->vertex()->d1(),
|
||||
d_q = he->vertex()->d1(); //ppal dir
|
||||
if ( color == BLUE ) {
|
||||
b_p = he->opposite()->vertex()->b0();
|
||||
b_q = he->vertex()->b0();
|
||||
}
|
||||
else {
|
||||
b_p = he->opposite()->vertex()->b3();
|
||||
b_q = he->vertex()->b3();
|
||||
}
|
||||
if ( b_p == 0 && b_q == 0 ) return;
|
||||
if ( b_p == 0 && b_q !=0 ) sign = d_p*d_q * b_q;
|
||||
if ( b_p != 0 && b_q ==0 ) sign = d_p*d_q * b_p;
|
||||
if ( b_p != 0 && b_q !=0 ) sign = d_p*d_q * b_p * b_q;
|
||||
if ( sign < 0 ) is_crossed = true;
|
||||
}
|
||||
|
||||
|
||||
|
||||
//for a ridge segment in a triangle [r1,r2], let r = r2 - r1 and normalize,
|
||||
//the projection of a point p on the line (r1,r2) is pp=r1+tr, with t=(p-r1)*r
|
||||
//then the vector v starting at p is pointing to the ridge line (r1,r2) if
|
||||
//(pp-p)*v >0
|
||||
int
|
||||
b_sign_pointing_to_ridge(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3,
|
||||
Vector_3 r1, Vector_3 r2, Ridge_type color)
|
||||
{
|
||||
Vector_3 r = r2 - r1, dv1, dv2, dv3;
|
||||
FT bv1, bv2, bv3;
|
||||
if ( color == BLUE ) {
|
||||
bv1 = v1->b0();
|
||||
bv2 = v2->b0();
|
||||
bv3 = v3->b0();
|
||||
dv1 = v1->d1();
|
||||
dv2 = v2->d1();
|
||||
dv3 = v3->d1();
|
||||
}
|
||||
else {
|
||||
bv1 = v1->b3();
|
||||
bv2 = v2->b3();
|
||||
bv3 = v3->b3();
|
||||
dv1 = v1->d2();
|
||||
dv2 = v2->d2();
|
||||
dv3 = v3->d2();
|
||||
}
|
||||
if ( r !=0 ) r = r/CGAL::sqrt(r*r);
|
||||
FT sign1, sign2, sign3;
|
||||
sign1 = (r1 - (v1->point()-ORIGIN) + (((v1->point()-ORIGIN)-r1)*r)*r )*dv1;
|
||||
sign2 = (r1 - (v2->point()-ORIGIN) + (((v2->point()-ORIGIN)-r1)*r)*r )*dv2;
|
||||
sign3 = (r1 - (v3->point()-ORIGIN) + (((v3->point()-ORIGIN)-r1)*r)*r )*dv3;
|
||||
|
||||
int compt = 0;
|
||||
if ( sign1 > 0 ) compt++; else if (sign1 < 0) compt--;
|
||||
if ( sign2 > 0 ) compt++; else if (sign2 < 0) compt--;
|
||||
if ( sign3 > 0 ) compt++; else if (sign3 < 0) compt--;
|
||||
|
||||
if (compt > 0) return 1; else return -1;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
CGAL_END_NAMESPACE
|
||||
|
||||
#endif
|
||||
Loading…
Reference in New Issue