Subversion Repositories gelsvn

#include <queue>

#include <vector>

#include "triangulate.h"

#include "HMesh/VertexCirculator.h"

#include "HMesh/FaceCirculator.h"

using namespace std;

using namespace CGLA;

using namespace HMesh;

namespace HMeshUtil

{

  void get_candidates(const VertexIter& v,

											vector<HalfEdgeIter>& candidates)

  {

    vector<VertexIter> vertices;

    vector<HalfEdgeIter> hedges;

    VertexCirculator vc(v);

    for(;!vc.end();++vc)

      {

				vertices.push_back(vc.get_vertex());

				hedges.push_back(vc.get_halfedge());

      }

    int N = vc.no_steps();

    vector<VertexIter> vertices_check(vertices);

    assert(N==vertices.size());

    for(int i=N-1;i>=0;--i)

      {

				HalfEdgeIter h = hedges[i]->next;

				while(h->vert != vertices[(i+N-1)%N])

					{

						if(find(vertices_check.begin(),vertices_check.end(), 

										h->vert) == vertices_check.end())

							candidates.push_back(h);

						h = h->next;

						assert(h!=v);

					}

      }

  }

  float curv(const Vec3f& p, std::vector<Vec3f>& vec)

  {

    int N = vec.size();

    std::vector<Vec3f> normals;

    for(int i=0;i<N;++i)

      {

				Vec3f norm = normalize(cross(vec[i]-p, vec[(i+1)%N]-p));

				normals.push_back(norm);

      }

    float alpha = 0;

    for(int i=0;i<N;++i)

      {

				alpha += (vec[i]-p).length()*acos(dot(normals[i],normals[(i+1)%N]));

      }

    return alpha;

  }

  float get_badness(const VertexIter& v, const VertexIter& n)

  {

    vector<HalfEdgeIter> hedges;

    VertexCirculator vc(v);

    for(;!vc.end();++vc)

      hedges.push_back(vc.get_halfedge());

    vector<Vec3f> one_ring;

    vector<Vec3f> one_ring_n;

    int N = vc.no_steps();

    for(int i=N-1;i>=0;--i)

      {

				HalfEdgeIter h = hedges[i];

				h = h->next;

				while(h->vert != v)

					{

						one_ring.push_back(h->vert->get_pos());

						if(h->vert != n)

							one_ring_n.push_back(h->vert->get_pos());

						h = h->next;

					}

      }

    return curv(v->get_pos(),one_ring)-curv(v->get_pos(),one_ring_n);

  }

  const CGLA::Vec3f get_normal(const VertexIter& v)

  {

    vector<HalfEdgeIter> hedges;

    VertexCirculator vc(v);

    for(;!vc.end();++vc)

      hedges.push_back(vc.get_halfedge());

    vector<Vec3f> one_ring;

    int N = vc.no_steps();

    for(int i=N-1;i>=0;--i)

      {

				HalfEdgeIter h = hedges[i];

				h = h->next;

				while(h->vert != v)

					{

						one_ring.push_back(h->vert->get_pos());

						h = h->next;

					}

      }

    Vec3f norm(0);

    N= one_ring.size();

    Vec3f p0 = v->get_pos();

    for(int i=0;i<N;++i)

      {

				Vec3f p1 = one_ring[i];

				Vec3f p2 = one_ring[(i+1)%N];

				Vec3f e0 = normalize(p1-p0);

				Vec3f e1 = normalize(p2-p0);

				norm += cross(e0,e1)*acos(dot(e0,e1));

      }

    return normalize(norm);

  }

  void safe_triangulate(Manifold& m)

  {

    vector<FaceIter> fv;

    for(FaceIter  fi=m.faces_begin(); fi != m.faces_end(); ++fi)

      fv.push_back(fi);

    for(int i=0;i<fv.size();++i)

      m.safe_triangulate(fv[i]);

  }

  void triangulate(Manifold& m)

  {

    vector<FaceIter> fv;

    for(FaceIter  fi=m.faces_begin(); fi != m.faces_end(); ++fi)

      fv.push_back(fi);

    for(int i=0;i<fv.size();++i)

      m.triangulate(fv[i]);

  }

  struct PotentialEdge

  {

    int time_tag;

    float badness;

    FaceIter f;

    VertexIter v0, v1;

  };

  bool operator>(const PotentialEdge& e0, const PotentialEdge& e1)

  {

    return e0.badness>e1.badness;

  }

  typedef std::priority_queue<PotentialEdge,

															std::vector<PotentialEdge>,

															std::greater<PotentialEdge> > 

  PotentialEdgeQueue;

  void insert_potential_edges(const VertexIter& v,

															PotentialEdgeQueue& pot_edges)

  {

    vector<Vec3f> one_ring;

    vector<HalfEdgeIter> candidates;

    get_candidates(v, candidates);

    Vec3f p0 = v->get_pos();

    Vec3f norm = get_normal(v);

    int n = candidates.size();

    for(int i=0;i<n;++i)

      {

				VertexIter v_n = candidates[i]->vert;

				Vec3f edir = normalize(v_n->get_pos()-p0);

				Vec3f norm_n = get_normal(v_n);

				float bad = sqr(dot(edir, norm));

				float bad_n = sqr(dot(edir, norm_n));

				PotentialEdge pot;

				/* So if the edge between two vertices is not orthogonal to 

					 their normals, the badness is increased. Badness is also

					 increased if the normals are very different. */

				pot.badness = bad+bad_n - dot(norm_n,norm);

				pot.time_tag = v->touched;

				pot.v0 = v;

				pot.v1 = candidates[i]->vert;

				pot.f = candidates[i]->face;

				pot_edges.push(pot);

      }

  }

  void create_candidates(Manifold& m, PotentialEdgeQueue& pot_edges)

  {

    for(VertexIter v= m.vertices_begin(); v!= m.vertices_end();++v)

      {

				v->touched = 0;

				insert_potential_edges(v, pot_edges);

      }

  }

  void curvature_triangulate(Manifold& m)

  {

    PotentialEdgeQueue pot_edges;

    // Create candidates

    create_candidates(m,pot_edges);

    while(!pot_edges.empty())

      {

				const PotentialEdge& pot_edge = pot_edges.top();

				// Record all the vertices of the face. We need to 

				// recompute the candidates.

				std::vector<VertexIter> reeval_vec;

				FaceCirculator fc(pot_edge.f);

				while(!fc.end())

					{

						reeval_vec.push_back(fc.get_vertex());

						++fc;

					}

				VertexIter v0 = pot_edge.v0;

				// Make sure that the vertex has not been touched since 

				// we created the potential edge. If the vertex has been

				// touched the candidate edge has either (a) been processed,

				// (b) received a lower priority or (c) become invalid.

				if(pot_edge.time_tag == v0->touched)

					{

						VertexIter v1 = pot_edge.v1;

						vector<Vec3f> one_ring;

						vector<HalfEdgeIter> candidates;

						m.split_face(pot_edge.f, v0, v1);

						// Recompute priorities.

						for(int i=0;i<reeval_vec.size();++i)

							{

								VertexIter& v = reeval_vec[i];

								v->touched++;

								insert_potential_edges(v, pot_edges);

							}

					}

				pot_edges.pop();

      }

  }

}