Rev 468 | Rev 507 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed
#include "triangulate.h"
#include <queue>
#include <vector>
#include <cassert>
#include <CGLA/Vec3f.h>
#include "Manifold.h"
#include "VertexCirculator.h"
#include "FaceCirculator.h"
#include "VertexHandle.h"
#include "HalfEdgeHandle.h"
#include "FaceHandle.h"
namespace HMesh
{
using std::vector;
using std::pair;
using std::cout;
using std::endl;
using CGLA::Vec3f;
using CGLA::sqr;
void get_candidates(const VertexHandle& v, vector<HalfEdgeHandle>& candidates)
{
vector<VertexHandle> vertices;
vector<HalfEdgeHandle> hedges;
VertexCirculator vc(v);
for(;!vc.end();++vc)
{
vertices.push_back(vc.vertex());
hedges.push_back(vc.halfedge());
}
size_t N = vc.no_steps();
vector<VertexHandle> vertices_check(vertices);
assert(N == vertices.size());
for(int i = N - 1; i >= 0; --i)
{
HalfEdgeHandle 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();
// TODO : does not compile in debug mode in visual studio
//assert(h!=v);
}
}
}
float curv(const Vec3f& p, 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 VertexHandle& v, const VertexHandle& n)
{
vector<HalfEdgeHandle> hedges;
VertexCirculator vc(v);
for(; !vc.end(); ++vc)
hedges.push_back(vc.halfedge());
vector<Vec3f> one_ring;
vector<Vec3f> one_ring_n;
int N = vc.no_steps();
for(int i = N - 1; i >= 0; --i)
{
HalfEdgeHandle 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 VertexHandle& v)
{
vector<HalfEdgeHandle> hedges;
VertexCirculator vc(v);
for(; !vc.end(); ++vc)
hedges.push_back(vc.halfedge());
vector<Vec3f> one_ring;
int N = vc.no_steps();
for(int i = N - 1; i >= 0; --i){
HalfEdgeHandle 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 triangulate_by_vertex_face_split(Manifold& m)
{
vector<FaceHandle> fv;
for(FaceHandle f = m.faces_begin(); f != m.faces_end(); ++f){
fv.push_back(f);
}
for(size_t i = 0; i < fv.size(); ++i){
m.split_face_by_vertex(fv[i]);
}
}
void triangulate_by_edge_face_split(Manifold& m)
{
vector<FaceHandle> fv;
for(FaceHandle f = m.faces_begin(); f != m.faces_end(); ++f)
fv.push_back(f);
for(size_t i = 0; i < fv.size(); ++i)
triangulate_face_by_edge_split(m, fv[i].get_idx());
}
struct PotentialEdge
{
int time_tag;
float badness;
FIDX f;
VIDX 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 VertexHandle& v, PotentialEdgeQueue& pot_edges)
{
vector<Vec3f> one_ring;
vector<HalfEdgeHandle> candidates;
get_candidates(v, candidates);
Vec3f p0 = v.get().pos;
Vec3f norm = normal(v);
int n = candidates.size();
for(int i=0;i<n;++i)
{
VertexHandle v_n = candidates[i].vert();
Vec3f edir = normalize(v_n.get().pos - p0);
Vec3f norm_n = 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.get().touched;
pot.v0 = v.get_idx();
pot.v1 = candidates[i].vert().get_idx();
pot.f = candidates[i].face().get_idx();
pot_edges.push(pot);
}
}
void create_candidates(Manifold& m, PotentialEdgeQueue& pot_edges)
{
for(VertexHandle v = m.vertices_begin(); v!= m.vertices_end(); ++v){
v.get().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<VertexHandle> reeval_vec;
FaceHandle f = m.face_handle(pot_edge.f);
for(FaceCirculator fc(f); !fc.end(); ++fc){
reeval_vec.push_back(fc.vertex());
}
VertexHandle v0 = m.vertex_handle(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.get().touched)
{
VertexHandle v1 = m.vertex_handle(pot_edge.v1);
vector<Vec3f> one_ring;
vector<HalfEdgeHandle> candidates;
m.split_face_by_edge(f, v0, v1);
// Recompute priorities.
for(size_t i = 0; i < reeval_vec.size(); ++i)
{
VertexHandle& v = reeval_vec[i];
v.get().touched++;
insert_potential_edges(v, pot_edges);
}
}
pot_edges.pop();
}
}
void shortest_edge_triangulate(Manifold& m)
{
int work;
do{
// For every face.
work = 0;
for(FaceHandle f = m.faces_begin(); f != m.faces_end(); ++f){
// Create a vector of vertices.
vector<VertexHandle> verts;
for(FaceCirculator fc(f); !fc.end(); ++fc)
verts.push_back(fc.vertex());
// If there are just three we are done.
if(verts.size() == 3) continue;
// Find vertex pairs that may be connected.
vector<pair<int,int> > vpairs;
const int N = verts.size();
for(int i = 0; i < N - 2; ++i){
for(int j = i + 2; j < N; ++j){
if(!is_connected(verts[i], verts[j]))
vpairs.push_back(pair<int,int>(i, j));
}
}
if(vpairs.empty()){
cout << "Warning: could not triangulate a face."
<< "Probably a vertex appears more than one time in other vertex's one-ring" << endl;
continue;
}
/* For all vertex pairs, find the edge lengths. Combine the
vertices forming the shortest edge. */
float min_len=FLT_MAX;
int min_k = -1;
for(size_t k = 0; k < vpairs.size(); ++k){
int i = vpairs[k].first;
int j = vpairs[k].second;
float len = sqr_length(verts[i].get().pos - verts[j].get().pos);
if(len<min_len){
min_len = len;
min_k = k;
}
}
assert(min_k != -1);
{
// Split faces along edge whose midpoint is closest to isovalue
int i = vpairs[min_k].first;
int j = vpairs[min_k].second;
m.split_face_by_edge(f, verts[i], verts[j]);
}
++work;
}
}
while(work);
}
void triangulate_face_by_edge_split(Manifold& m, VIDX idx)
{
FaceHandle f = m.face_handle(idx);
assert(f.last().next().next() != f.last());
assert(f.last().next() != f.last());
// as long as f is not a triangle
while(f.last().next().next().next() != f.last()){
assert(no_edges(f) > 3);
// f is split into triangle from first three vertices, and becomes remaining polygon in next iteration
assert(f.last().next().next() != f.last());
VertexHandle v0 = f.last().vert();
VertexHandle v1 = f.last().next().next().vert();
assert(v0 != v1);
f = m.split_face_by_edge(f, v0, v1);
}
}
}