WebSVN – gelsvn – Diff – /trunk/src/HMesh/mesh_optimization.cpp

+/* ----------------------------------------------------------------------- *
+ * This file is part of GEL, http://www.imm.dtu.dk/GEL
+ * Copyright (C) the authors and DTU Informatics
+ * For license and list of authors, see ../../doc/intro.pdf
+ * ----------------------------------------------------------------------- */
+#include "mesh_optimization.h"
 #include <cfloat>
-#include "CGLA/Vec3d.h"
 #include <queue>
-#include "mesh_optimization.h"
-#include "HMesh/FaceCirculator.h"
+#include <vector>
-#include "HMesh/VertexCirculator.h"
-using namespace std;
+#include <CGLA/Vec3d.h>
-using namespace CGLA;
+//#include "Manifold.h"
-using namespace HMesh;
+#include "AttributeVector.h"
 namespace HMesh
+{
+    using namespace std;
+    using namespace CGLA;
-	// Small utility functions
+    // Small utility functions
-	namespace
+    namespace
-	{
+    {
-		class LineSeg
+        class LineSeg
-		{
+        {
-			const CGLA::Vec3d p0;
+            const Vec3d p0;
-			const CGLA::Vec3d p1;
+            const Vec3d p1;
-			const float l;
+            const float l;
-			const CGLA::Vec3d dir;
+            const CGLA::Vec3d dir;
-		public:
+        public:
-			LineSeg(const CGLA::Vec3d& _p0, const CGLA::Vec3d& _p1):
+            LineSeg(const Vec3d& _p0, const Vec3d& _p1):
-				p0(_p0), p1(_p1), l(length(p1-p0)), dir((p1-p0)/l) {}
+			p0(_p0), p1(_p1), l(length(p1-p0)), dir((p1-p0)/l) {}
-			bool inseg(const CGLA::Vec3d& p) const
+			bool inseg(const Vec3d& p) const
+			{
 				double t = dot(dir, p-p0);
 				if(t<0)
 					return false;
 				if(t>l)
 					return false;
 				return true;
+			}
+        };
+        Vec3d compute_normal(Vec3d* v)
+        {
+            Vec3d norm;
+            for(int i = 0; i < 4; ++i)
+            {
+                norm[0] += (v[i][1]-v[(i+1)%4][1])*(v[i][2]+v[(i+1)%4][2]);
+                norm[1] += (v[i][2]-v[(i+1)%4][2])*(v[i][0]+v[(i+1)%4][0]);
+                norm[2] += (v[i][0]-v[(i+1)%4][0])*(v[i][1]+v[(i+1)%4][1]);
+            }
+            float l = norm.length();
+            if(l>0.0f)
+                norm /= l;
+            return norm;
+        }
+        bool would_flip(const Manifold& m, HalfEdgeID h)
+        {
+            Walker w = m.walker(h);
+            VertexID hv = w.vertex();
+            VertexID hov = w.opp().vertex();
+            VertexID hnv = w.next().vertex();
+            VertexID honv = w.opp().next().vertex();
+            Vec3d v[4];
+            v[0] = Vec3d(m.pos(hv));
+            v[1] = Vec3d(m.pos(hov));
+            v[2] = Vec3d(m.pos(hnv));
+            v[3] = Vec3d(m.pos(honv));
+            Vec3d dir = compute_normal(v);
+            Vec3d n1a = cross(v[3]-v[0], v[2]-v[0]);
+            Vec3d n2a = cross(v[2]-v[1], v[3]-v[1]);
-		};
+            if(dot(normalize(n1a), dir) < 0)
+                return true;
+            if(dot(normalize(n2a), dir) < 0)
+                return true;
+            return false;
+        }
+    }
+    double ValencyEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const
+    {
+        Walker w = m.walker(h);
+        VertexID v1 = w.opp().vertex();
+        VertexID v2 = w.vertex();
+        VertexID vo1 = w.next().vertex();
+        VertexID vo2 = w.opp().next().vertex();
+        int val1  = valency(m, v1);
-		Vec3d compute_normal(Vec3d* v)
+        int val2  = valency(m, v2);
+        int valo1 = valency(m, vo1);
+        int valo2 = valency(m, vo2);
+        // The optimal valency is four for a boundary vertex
+        // and six elsewhere.
+        int t1 = boundary(m, v1) ? 4 : 6;
+        int t2 = boundary(m, v2) ? 4 : 6;
+        int to1 = boundary(m, vo1) ? 4 : 6;
+        int to2 = boundary(m, vo2) ? 4 : 6;
+        int before =
+        sqr(val1-t1)+sqr(val2-t2)+
+        sqr(valo1-to1)+sqr(valo2-to2);
+        int after =
+        sqr(valo1+1-to1)+sqr(val1-1-t1)+
+        sqr(val2-1-t2)+sqr(valo2+1-to2);
+        return static_cast<double>(after-before);
+    }
+	class RandomEnergy: public EnergyFun
+	{
+	public:
+		double delta_energy(const Manifold& m, HalfEdgeID he) const
+		{
-			Vec3d norm;
-			for(int i=0;i<4;++i)
-				{
-					norm[0] += (v[i][1]-v[(i+1)%4][1])*(v[i][2]+v[(i+1)%4][2]);
-					norm[1] += (v[i][2]-v[(i+1)%4][2])*(v[i][0]+v[(i+1)%4][0]);
+			return static_cast<double>(gel_rand()/static_cast<float>(GEL_RAND_MAX));
-					norm[2] += (v[i][0]-v[(i+1)%4][0])*(v[i][1]+v[(i+1)%4][1]);
-				}
-			float l = norm.length();
-			if(l>0.0f)
-				norm /= l;
-			return norm;
+		}
+	};
+	double MinAngleEnergy::min_angle(const Vec3d& v0, const Vec3d& v1, const Vec3d& v2) const
+	{
+		Vec3d a = normalize(v1-v0);
+		Vec3d b = normalize(v2-v1);
+		Vec3d c = normalize(v0-v2);
+		return min(dot(a,-c), min(dot(b,-a), dot(c,-b)));
+	}
+	double MinAngleEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const
+	{
-		bool would_flip(HalfEdgeIter h)
+		Walker w = m.walker(h);
-		{
-			HalfEdgeIter ho = h->opp;
+		VertexID hv = w.vertex();
+		VertexID hnv = w.next().vertex();
+		VertexID hov= w.opp().vertex();
+		VertexID honv = w.opp().next().vertex();
-			Vec3d v[4];
+		Vec3d v0(m.pos(hv));
+		Vec3d v1(m.pos(hnv));
-			v[0]=Vec3d(h->vert->pos);
+		Vec3d v2(m.pos(hov));
+		Vec3d v3(m.pos(honv));
+		Vec3d n1a = normalize(cross(v1-v0,v3-v0));
+		Vec3d n2a = normalize(cross(v3-v2,v1-v2));
+		if(dot(n1a, n2a) > thresh){
+			double before = min(min_angle(v0,v1,v2), min_angle(v0,v2,v3));
+			double after = min(min_angle(v0,v1,v3), min_angle(v1,v2,v3));
-			v[1]=Vec3d(ho->vert->pos);
+			return -(after-before);
+		}
+		return DBL_MAX;
+	}
+	void DihedralEnergy::compute_angles(const Manifold & m, HalfEdgeID h) const
+	{
+		Walker w = m.walker(h);
-			v[2]=Vec3d(h->next->vert->pos);
+		VertexID hv = w.vertex();
+		VertexID hov = w.opp().vertex();
-			v[3]=Vec3d(ho->next->vert->pos);
+		VertexID hnv = w.next().vertex();
+		VertexID honv = w.opp().next().vertex();
+		Vec3d va(m.pos(hv));
+		Vec3d vb(m.pos(hov));
+		Vec3d vc(m.pos(hnv));
+		Vec3d vd(m.pos(honv));
+		FaceID fa = w.next().opp().face();
+		FaceID fb = w.next().next().opp().face();
+		FaceID fc = w.opp().next().opp().face();
+		FaceID fd = w.opp().next().next().opp().face();
+		Vec3d n1 = normalize(cross(vc-va, vb-va));
+		Vec3d n2 = normalize(cross(vb-va, vd-va));
+		Vec3d na = fa == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fa));
+		Vec3d nb = fb == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fb));
+		Vec3d nc = fc == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fc));
+		Vec3d nd = fd == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fd));
+		Vec3d fn1 = normalize(cross(vb-vc, vd-vc));
+		Vec3d fn2 = normalize(cross(vd-vc, va-vc));
+		ab_12 = cos_ang(n1,n2);
+		ab_a1 = cos_ang(na,n1);
+		ab_b1 = cos_ang(nb,n1);
+		ab_2c = cos_ang(n2,nc);
+		ab_2d = cos_ang(n2,nd);
+		aa_12 = cos_ang(fn1,fn2);
+		aa_b1 = cos_ang(nb,fn1);
+		aa_c1 = cos_ang(nc, fn1);
+		aa_2a = cos_ang(fn2, na);
+		aa_2d = cos_ang(fn2,nd);
+	}
+	double DihedralEnergy::energy(const Manifold& m, HalfEdgeID h) const
+	{
+		Walker w = m.walker(h);
+		FaceID hf = w.face();
+		FaceID hof = w.opp().face();
+		double a = cos_ang(Vec3d(normal(m, hf)), Vec3d(normal(m, hof)));
+		VertexID hv = w.vertex();
-			Vec3d dir(0,0,-1);//compute_normal(v);
+		VertexID hov = w.opp().vertex();
+		Vec3d va(m.pos(hv));
+		Vec3d vb(m.pos(hov));
+		if(use_alpha)
+			return edge_alpha_energy(va,vb,a);
+		return edge_c_energy(va,vb,a);
+	}
+	double DihedralEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const
+	{
+		compute_angles(m, h);
+		Walker w = m.walker(h);
+		VertexID hv = w.vertex();
+		VertexID hov = w.opp().vertex();
+		VertexID hnv = w.next().vertex();
+		VertexID honv = w.opp().next().vertex();
+		Vec3d va(m.pos(hv));
+		Vec3d vb(m.pos(hov));
+		Vec3d vc(m.pos(hnv));
+		Vec3d vd(m.pos(honv));
+		if(use_alpha){
+			double before =
+			edge_alpha_energy(va,vb,ab_12)
+			+edge_alpha_energy(va,vc,ab_a1)
+			+edge_alpha_energy(vc,vb,ab_b1)
+			+edge_alpha_energy(vd,vb,ab_2c)
+			+edge_alpha_energy(vd,va,ab_2d);
+			double after =
+			edge_alpha_energy(vd,vc,aa_12)
+			+edge_alpha_energy(vb,vc,aa_b1)
+			+edge_alpha_energy(vd,vb,aa_c1)
-			Vec3d n1a = cross(v[3]-v[0], v[2]-v[0]);
+			+edge_alpha_energy(va,vc,aa_2a)
-			Vec3d n2a = cross(v[2]-v[1], v[3]-v[1]);
+			+edge_alpha_energy(vd,va,aa_2d);
-			if(dot(normalize(n1a),dir) < 0)
-				return true;
-			if(dot(normalize(n2a),dir) < 0)
-				return true;
+			return (after-before);
-			return false;
+		}
+		double before =
+		edge_c_energy(va,vb,ab_12)
+		+edge_c_energy(va,vc,ab_a1)
+		+edge_c_energy(vc,vb,ab_b1)
+		+edge_c_energy(vd,vb,ab_2c)
+		+edge_c_energy(vd,va,ab_2d);
+		double after =
+		edge_c_energy(vd,vc,aa_12)
+		+edge_c_energy(vb,vc,aa_b1)
+		+edge_c_energy(vd,vb,aa_c1)
+		+edge_c_energy(va,vc,aa_2a)
+		+edge_c_energy(vd,va,aa_2d);
+		return after-before;
+	}
-	// Energy functionals
-	namespace
+	double CurvatureEnergy::abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const
+	{
+		const size_t N = ring.size();
+		double H = 0;
+		for(int i = 0; i < N; ++i){
+			Vec3d vnim1 = ring[(i+N-1)%N] - v;
+			Vec3d vni   = ring[i] - v;
+			Vec3d vnip1 = ring[(i+1)%N] - v;
-		class ValencyEnergy: public EnergyFun
+			Vec3d Nm = normalize(cross(vni, vnim1));
+			Vec3d Np = normalize(cross(vnip1, vni));
+			double beta = acos(max(-1.0, min(1.0, dot(Nm, Np))));
+			H += vni.length() * beta;
-		{
+		}
-		public:
+		H /= 4;
+		return H;
+	}
-			double delta_energy(HalfEdgeIter he) const
+	double CurvatureEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const
-			{
+	{
-				FaceIter f = he->face;
+		Walker w = m.walker(h);
-				VertexIter v1 = he->opp->vert;
+		VertexID va = w.vertex();
-				VertexIter v2 = he->vert;
+		VertexID vb = w.opp().vertex();
-				VertexIter vo1 = he->next->vert;
+		VertexID vc = w.next().vertex();
-				VertexIter vo2 = he->opp->next->vert;
+		VertexID vd = w.opp().next().vertex();
-				int val1  = valency(v1);
+		Vec3d va_pos(m.pos(va));
+		Vec3d vb_pos(m.pos(vb));
+		Vec3d vc_pos(m.pos(vc));
+		Vec3d vd_pos(m.pos(vd));
+		vector<Vec3d> va_ring_bef;
+		vector<Vec3d> va_ring_aft;
+		vector<Vec3d> vb_ring_bef;
+		vector<Vec3d> vb_ring_aft;
+		vector<Vec3d> vc_ring_bef;
+		vector<Vec3d> vc_ring_aft;
+		vector<Vec3d> vd_ring_bef;
+		vector<Vec3d> vd_ring_aft;
+		for(Walker wv = m.walker(va); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
-				int val2  = valency(v2);
+			VertexID v = wv.vertex();
+			Vec3d pos(m.pos(v));
+			va_ring_bef.push_back(pos);
+			if(v != vb)
-				int valo1 = valency(vo1);
+				va_ring_aft.push_back(pos);
+		}
+		for(Walker wv = m.walker(vb); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
-				int valo2 = valency(vo2);
+			VertexID v = wv.vertex();
+			Vec3d pos(m.pos(v));
+			vb_ring_bef.push_back(pos);
+			if(v != va)
+				vb_ring_aft.push_back(pos);
+		}
-				// The optimal valency is four for a boundary vertex
+		for(Walker wv = m.walker(vc); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
-				// and six elsewhere.
+			VertexID v = wv.vertex();
+			Vec3d pos(m.pos(v));
-				int t1 = is_boundary(v1) ? 4 : 6;
+			vc_ring_bef.push_back(pos);
-				int t2 = is_boundary(v2) ? 4 : 6;
+			vc_ring_aft.push_back(pos);
+			if(v == va)
-				int to1 = is_boundary(vo1) ? 4 : 6;
+				vc_ring_aft.push_back(vd_pos);
+		}
+		for(Walker wv = m.walker(vd); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
-				int to2 = is_boundary(vo2) ? 4 : 6;
+			VertexID v = wv.vertex();
+			Vec3d pos(m.pos(v));
+			vd_ring_bef.push_back(pos);
+			vd_ring_aft.push_back(pos);
+			if(v == vb)
+				vd_ring_aft.push_back(vc_pos);
+		}
-				int before =
+		double before =
-					sqr(val1-t1)+sqr(val2-t2)+
+		abs_mean_curv(va_pos, va_ring_bef) +
+		abs_mean_curv(vb_pos, vb_ring_bef) +
+		abs_mean_curv(vc_pos, vc_ring_bef) +
-					sqr(valo1-to1)+sqr(valo2-to2);
+		abs_mean_curv(vd_pos, vd_ring_bef);
-				int after =
+		double after =
-					sqr(valo1+1-to1)+sqr(val1-1-t1)+
+		abs_mean_curv(va_pos, va_ring_aft) +
+		abs_mean_curv(vb_pos, vb_ring_aft) +
+		abs_mean_curv(vc_pos, vc_ring_aft) +
-					sqr(val2-1-t2)+sqr(valo2+1-to2);
+		abs_mean_curv(vd_pos, vd_ring_aft);
-				return static_cast<double>(after-before);
+		return after-before;
-			}
+	}
-		};
-		class RandomEnergy: public EnergyFun
+	class GaussCurvatureEnergy: public EnergyFun
+	{
+		double gauss_curv(const Vec3d& v, const vector<Vec3d>& ring) const
+		{
+			const size_t N = ring.size();
-		public:
+			double asum=0.0f;
+			double area_sum=0;
+			for(int i = 0; i < N; ++i){
+				const Vec3d& v1 = ring[i];
-			double delta_energy(HalfEdgeIter he) const
+				const Vec3d& v2 = ring[(i+1)%N];
+				Vec3d a = v1-v;
-			{
+				Vec3d b = v2-v;
-				return static_cast<double>(gel_rand()/static_cast<float>(GEL_RAND_MAX));
+				asum += acos(max(-1.0, min(1.0, dot(a,b)/(length(a)*length(b)))));
+				area_sum += 0.5 * length(cross(a,b));
+			}
+			return 3*abs(2 * M_PI - asum)/area_sum;
-		};
+		}
+	public:
-		class MinAngleEnergy: public EnergyFun
+		double delta_energy(const Manifold& m, HalfEdgeID h) const
+		{
-			double min_angle(const Vec3d& v0,const Vec3d& v1,const Vec3d& v2) const
-			{
-				Vec3d a = normalize(v1-v0);
-				Vec3d b = normalize(v2-v1);
-				Vec3d c = normalize(v0-v2);
+			Walker w = m.walker(h);
-				return min(dot(a,-c), min(dot(b,-a), dot(c,-b)));
-			}
-			double thresh;
-		public:
-			MinAngleEnergy(double _thresh): thresh(_thresh) {}
-			double delta_energy(HalfEdgeIter he) const
-			{
-				FaceIter f = he->face;
+			VertexID va = w.vertex();
-				FaceIter f2 = he->opp->face;
+			VertexID vb = w.opp().vertex();
-				Vec3d v0(he->vert->pos);
-				Vec3d v1(he->next->vert->pos);
+			VertexID vc = w.next().vertex();
-				Vec3d v2(he->opp->vert->pos);
-				Vec3d v3(he->opp->next->vert->pos);
+			VertexID vd = w.opp().next().vertex();
-				Vec3d n1a = normalize(cross(v1-v0,v3-v0));
-				Vec3d n2a = normalize(cross(v3-v2,v1-v2));
-				if(dot(n1a,n2a)> thresh)
-					{
-						double before = min(min_angle(v0,v1,v2), min_angle(v0,v2,v3));
-						double after = min(min_angle(v0,v1,v3), min_angle(v1,v2,v3));
-						return -(after-before);
-					}
-				return DBL_MAX;
-			}
-		};
-		class DihedralEnergy: public EnergyFun
-		{
-			const double gamma;
+			Vec3d va_pos(m.pos(va));
-			const bool use_alpha;
+			Vec3d vb_pos(m.pos(vb));
-			double cos_ang(const Vec3d& n1, const Vec3d& n2) const
-			{
-				return s_max(-1.0, s_min(1.0, dot(n1, n2)));
+			Vec3d vc_pos(m.pos(vc));
-			}
-			double edge_alpha_energy(Vec3d v1, Vec3d v2, double ca) const
-			{
-				return pow(length(v1-v2)*(acos(ca)), 1.0f/gamma);
-			}
-			double edge_c_energy(Vec3d v1, Vec3d v2, double ca) const
-			{
-				return pow(length(v1-v2)*(1-ca), 1.0f/gamma);
+			Vec3d vd_pos(m.pos(vd));
-			}
-			mutable double ab_12;
+			vector<Vec3d> va_ring_bef;
-			mutable double ab_a1;
+			vector<Vec3d> va_ring_aft;
-			mutable double ab_b1;
+			vector<Vec3d> vb_ring_bef;
-			mutable double ab_2c;
-			mutable double ab_2d;
-			mutable double aa_12;
+			vector<Vec3d> vb_ring_aft;
-			mutable double aa_b1;
+			vector<Vec3d> vc_ring_bef;
-			mutable double aa_c1;
+			vector<Vec3d> vc_ring_aft;
-			mutable double aa_2a;
+			vector<Vec3d> vd_ring_bef;
-			mutable double aa_2d;
+			vector<Vec3d> vd_ring_aft;
-		public:
-			DihedralEnergy(double _gamma = 4.0, bool _use_alpha=false):
-				gamma(_gamma), use_alpha(_use_alpha) {}
-			void compute_angles(HalfEdgeIter he) const
-			{
-				Vec3d va(he->vert->pos);
-				Vec3d vb(he->opp->vert->pos);
-				Vec3d vc(he->next->vert->pos);
-				Vec3d vd(he->opp->next->vert->pos);
-				FaceIter fa = he->next->opp->face;
-				FaceIter fb = he->next->next->opp->face;
-				FaceIter fc = he->opp->next->opp->face;
-				FaceIter fd = he->opp->next->next->opp->face;
-				Vec3d n1 = normalize(cross(vc-va, vb-va));
-				Vec3d n2 = normalize(cross(vb-va, vd-va));
-				Vec3d na = fa==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fa));
+			for(Walker wv = m.walker(va); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
-				Vec3d nb = fb==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fb));
-				Vec3d nc = fc==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fc));
-				Vec3d nd = fd==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fd));
-				Vec3d fn1 = normalize(cross(vb-vc, vd-vc));
-				Vec3d fn2 = normalize(cross(vd-vc, va-vc));
-				ab_12 =	cos_ang(n1,n2);
-				ab_a1 =	cos_ang(na,n1);
-				ab_b1 =	cos_ang(nb,n1);
-				ab_2c =	cos_ang(n2,nc);
+				VertexID v = wv.vertex();
-				ab_2d =	cos_ang(n2,nd);
+				Vec3d pos(m.pos(v));
-				aa_12 = cos_ang(fn1,fn2);
+				va_ring_bef.push_back(pos);
-				aa_b1 = cos_ang(nb,fn1);
+				if(v != vb)
-				aa_c1 = cos_ang(nc, fn1);
-				aa_2a = cos_ang(fn2, na);
+					va_ring_aft.push_back(pos);
-				aa_2d = cos_ang(fn2,nd);
 			}
-			double energy(HalfEdgeIter he) const
-			{
-				double a = cos_ang(Vec3d(normal(he->face)),
+			for(Walker wv = m.walker(vb); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
-													 Vec3d(normal(he->opp->face)));
-				Vec3d va(he->vert->pos);
+				VertexID v = wv.vertex();
-				Vec3d vb(he->opp->vert->pos);
+				Vec3d pos(m.pos(v));
-				if(use_alpha)
-					return edge_alpha_energy(va,vb,a);
+				vb_ring_bef.push_back(pos);
+				if(v != va)
-				return edge_c_energy(va,vb,a);
+					vb_ring_aft.push_back(pos);
 			}
+			for(Walker wv = m.walker(vc); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
+				VertexID v = wv.vertex();
+				Vec3d pos(m.pos(v));
-			double delta_energy(HalfEdgeIter he) const
-			{
-				compute_angles(he);
-				Vec3d va(he->vert->pos);
-				Vec3d vb(he->opp->vert->pos);
+				vc_ring_bef.push_back(pos);
-				Vec3d vc(he->next->vert->pos);
+				vc_ring_aft.push_back(pos);
-				Vec3d vd(he->opp->next->vert->pos);
-				if(use_alpha)
+				if(v == va)
-					{
-						double before =
-							edge_alpha_energy(va,vb,ab_12)
-							+edge_alpha_energy(va,vc,ab_a1)
-							+edge_alpha_energy(vc,vb,ab_b1)
-							+edge_alpha_energy(vd,vb,ab_2c)
-							+edge_alpha_energy(vd,va,ab_2d);
-						double after =
-							edge_alpha_energy(vd,vc,aa_12)
-							+edge_alpha_energy(vb,vc,aa_b1)
-							+edge_alpha_energy(vd,vb,aa_c1)
-							+edge_alpha_energy(va,vc,aa_2a)
-							+edge_alpha_energy(vd,va,aa_2d);
-						return after-before;
+					vc_ring_aft.push_back(pos);
-					}
+			}
-				double before =
-					edge_c_energy(va,vb,ab_12)
-					+edge_c_energy(va,vc,ab_a1)
-					+edge_c_energy(vc,vb,ab_b1)
-					+edge_c_energy(vd,vb,ab_2c)
-					+edge_c_energy(vd,va,ab_2d);
-				double after =
-					edge_c_energy(vd,vc,aa_12)
-					+edge_c_energy(vb,vc,aa_b1)
-					+edge_c_energy(vd,vb,aa_c1)
-					+edge_c_energy(va,vc,aa_2a)
+			for(Walker wv = m.walker(vd); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
-					+edge_c_energy(vd,va,aa_2d);
+				VertexID v = wv.vertex();
-				return after-before;
-			}
-			double min_angle(HalfEdgeIter he) const
-			{
-				compute_angles(he);
+				Vec3d pos(m.pos(v));
-				return s_min(s_min(s_min(s_min(aa_12, aa_b1), aa_c1), aa_2a), aa_2d);
-			}
-		};
-		class CurvatureEnergy: public EnergyFun
-		{
-			// 			double abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const
-			// 			{
-			// 				const int N = ring.size();
-			// 				vector<Vec3d> norms(N,Vec3d(0));
-			// 				for(int i=0;i<N; ++i)
-			// 					{
-			// 						const Vec3d& vn = ring[i];
-			// 						const Vec3d& vnn = ring[(i+1)%N];
-			// 						Vec3d n = cross(vnn-v,vn-v);
-			// 						double ln = length(vn);
+				vd_ring_bef.push_back(pos);
-			// 						if(ln> 1e-20)
-			// 							norms[i] = n/ln;
+				vd_ring_aft.push_back(pos);
-			// 					}
+				if(v == vb)
-			// 				double H = 0.0f;
-			// 				for(int i=0;i<N; ++i)
+					vd_ring_aft.push_back(pos);
-			// 					{
-			// 						const Vec3d& vn = ring[i];
-			// 						double e_len = length(v-vn);
-			// 						double beta =
-			// 							acos(max(-1.0, min(1.0, dot(norms[i], norms[(i+N-1)%N]))));
-			// 						H += e_len * beta;
-			// 					}
-			// 				return H/4.0f;
-			// 			}
- 			double abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const
-			{
+			}
-				const int N = ring.size();
-				double H = 0;
+			double before =
-				for(int i=0;i<N;++i)
-					{
-						Vec3d vnim1 = ring[(i+N-1)%N] - v;
+			gauss_curv(va_pos,va_ring_bef) +
-						Vec3d vni   = ring[i] - v;
+			gauss_curv(vb_pos,vb_ring_bef) +
-						Vec3d vnip1 = ring[(i+1)%N] - v;
+			gauss_curv(vc_pos,vc_ring_bef) +
-						Vec3d Nm = normalize(cross(vni, vnim1));
-						Vec3d Np = normalize(cross(vnip1, vni));
-						double beta = acos(max(-1.0, min(1.0, dot(Nm, Np))));
-						H += vni.length() * beta;
+			gauss_curv(vd_pos,vd_ring_bef);
-					}
-				H /= 4;
-				return H;
-			}
-		public:
+			double after =
-			double delta_energy(HalfEdgeIter he) const
-			{
-				VertexIter va = he->vert;
-				VertexIter vb = he->opp->vert;
-				VertexIter vc = he->next->vert;
-				VertexIter vd = he->opp->next->vert;
-				Vec3d va_pos(va->pos);
-				Vec3d vb_pos(vb->pos);
-				Vec3d vc_pos(vc->pos);
-				Vec3d vd_pos(vd->pos);
-				vector<Vec3d> va_ring_bef;
-				vector<Vec3d> va_ring_aft;
+			gauss_curv(va_pos,va_ring_aft) +
-				vector<Vec3d> vb_ring_bef;
-				vector<Vec3d> vb_ring_aft;
+			gauss_curv(vb_pos,vb_ring_aft) +
-				vector<Vec3d> vc_ring_bef;
-				vector<Vec3d> vc_ring_aft;
+			gauss_curv(vc_pos,vc_ring_aft) +
-				vector<Vec3d> vd_ring_bef;
-				vector<Vec3d> vd_ring_aft;
+			gauss_curv(vd_pos,vd_ring_aft);
-				for(VertexCirculator c(va); !c.end(); ++c)
-					{
-						va_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
-						if(c.get_vertex() != vb)
+			return after-before;
-							va_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
-					}
+		}
-				for(VertexCirculator c(vb); !c.end(); ++c)
-					{
-						vb_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
-						if(c.get_vertex() != va)
-							vb_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
-					}
+	};
-				for(VertexCirculator c(vc); !c.end(); ++c)
-					{
-						vc_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
-						vc_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
-						if(c.get_vertex() == va)
-							vc_ring_aft.push_back(Vec3d(vd->pos));
-					}
+	struct PQElem
-				for(VertexCirculator c(vd); !c.end(); ++c)
-					{
+	{
-						vd_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
-						vd_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
-						if(c.get_vertex() == vb)
-							vd_ring_aft.push_back(Vec3d(vc->pos));
-					}
-				double before =
+		double pri;
-					abs_mean_curv(va_pos,va_ring_bef) +
-					abs_mean_curv(vb_pos,vb_ring_bef) +
-					abs_mean_curv(vc_pos,vc_ring_bef) +
-					abs_mean_curv(vd_pos,vd_ring_bef);
+		HalfEdgeID h;
-				double after =
+		int time;
-					abs_mean_curv(va_pos,va_ring_aft) +
-					abs_mean_curv(vb_pos,vb_ring_aft) +
-					abs_mean_curv(vc_pos,vc_ring_aft) +
-					abs_mean_curv(vd_pos,vd_ring_aft);
-				return after-before;
-			}
-		};
-		class GaussCurvatureEnergy: public EnergyFun
-		{
-			double gauss_curv(const Vec3d& v, const vector<Vec3d>& ring) const
-			{
-				const int N = ring.size();
-				double asum=0.0f;
+		//PQElem() {}
-				double area_sum=0;
-				for(int i=0;i<N; ++i)
-					{
-						const Vec3d& v1 = ring[i];
-						const Vec3d& v2 = ring[(i+1)%N];
-						Vec3d a = v1-v;
-						Vec3d b = v2-v;
-						asum += acos(max(-1.0, min(1.0, dot(a,b)/(length(a)*length(b)))));
-						area_sum += 0.5 * length(cross(a,b));
+		PQElem(double _pri, HalfEdgeID _h, int _time):
-					}
-				return 3*abs(2 * M_PI - asum)/area_sum;
+		pri(_pri), h(_h), time(_time) {}
-			}
+	};
-		public:
-			double delta_energy(HalfEdgeIter he) const
+	bool operator<(const PQElem& e0, const PQElem& e1)
-			{
+	{
-				VertexIter va = he->vert;
+		return e0.pri > e1.pri;
-				VertexIter vb = he->opp->vert;
-				VertexIter vc = he->next->vert;
-				VertexIter vd = he->opp->next->vert;
+	}
-				Vec3d va_pos(va->pos);
-				Vec3d vb_pos(vb->pos);
-				Vec3d vc_pos(vc->pos);
-				Vec3d vd_pos(vd->pos);
-				vector<Vec3d> va_ring_bef;
-				vector<Vec3d> va_ring_aft;
-				vector<Vec3d> vb_ring_bef;
-				vector<Vec3d> vb_ring_aft;
-				vector<Vec3d> vc_ring_bef;
-				vector<Vec3d> vc_ring_aft;
-				vector<Vec3d> vd_ring_bef;
-				vector<Vec3d> vd_ring_aft;
-				for(VertexCirculator c(va); !c.end(); ++c)
-					{
-						va_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
-						if(c.get_vertex() != vb)
-							va_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
-					}
-				for(VertexCirculator c(vb); !c.end(); ++c)
-					{
-						vb_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
-						if(c.get_vertex() != va)
-							vb_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
+	void add_to_queue(const Manifold& m, HalfEdgeAttributeVector<int>& touched, priority_queue<PQElem>& Q, HalfEdgeID h, const EnergyFun& efun)
-					}
-				for(VertexCirculator c(vc); !c.end(); ++c)
-					{
-						vc_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
-						vc_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
-						if(c.get_vertex() == va)
-							vc_ring_aft.push_back(Vec3d(vd->pos));
-					}
-				for(VertexCirculator c(vd); !c.end(); ++c)
-					{
+	{
-						vd_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
-						vd_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
-						if(c.get_vertex() == vb)
+		if(boundary(m, h))
-							vd_ring_aft.push_back(Vec3d(vc->pos));
-					}
+			return;
-				double before =
-					gauss_curv(va_pos,va_ring_bef) +
-					gauss_curv(vb_pos,vb_ring_bef) +
-					gauss_curv(vc_pos,vc_ring_bef) +
-					gauss_curv(vd_pos,vd_ring_bef);
-				double after =
-					gauss_curv(va_pos,va_ring_aft) +
-					gauss_curv(vb_pos,vb_ring_aft) +
-					gauss_curv(vc_pos,vc_ring_aft) +
-					gauss_curv(vd_pos,vd_ring_aft);
-				return after-before;
+		Walker w = m.walker(h);
-			}
-		};
-		struct PQElem
-		{
-			double pri;
-			HalfEdgeIter he;
+		HalfEdgeID ho = w.opp().halfedge();
-			int time;
-			PQElem() {}
-			PQElem(double _pri, HalfEdgeIter _he, int _time):
-				pri(_pri), he(_he), time(_time) {}
-		};
-		bool operator<(const PQElem& e0, const PQElem& e1)
+		double energy = efun.delta_energy(m, h);
+		int t = touched[h] + 1;
-		{
+		touched[h] = t;
+		touched[ho] = t;
+		if((energy<0) && (t < 10000)){
-			return e0.pri > e1.pri;
+			Q.push(PQElem(energy, h, t));
+		}
+	}
-			void add_to_queue(priority_queue<PQElem>& Q,
+	void add_one_ring_to_queue(const Manifold& m, HalfEdgeAttributeVector<int>& touched, priority_queue<PQElem>& Q, VertexID v, const EnergyFun& efun)
-												HalfEdgeIter he,
-												const EnergyFun& efun)
-			{
+	{
-				if(!is_boundary(he))
-					{
-						double energy = efun.delta_energy(he);
-						int t = he->touched + 1;
-						he->touched = t;
-						he->opp->touched = t;
-						if((energy<0) && (t < 10000))
-							Q.push(PQElem(energy, he, t));
-					}
-			}
-		void add_one_ring_to_queue(priority_queue<PQElem>& Q,
-															 VertexIter vi,
-															 const EnergyFun& efun)
-		{
-			for(VertexCirculator vc(vi); !vc.end(); ++vc)
+		for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
-				{
-					add_to_queue(Q, vc.get_halfedge(), efun);
+			add_to_queue(m, touched, Q, w.halfedge(), efun);
-				}
+		}
+	}
 	void priority_queue_optimization(Manifold& m, const EnergyFun& efun)
+	{
+		HalfEdgeAttributeVector<int> touched(m.allocated_halfedges(), 0);
 		priority_queue<PQElem> Q;
-		{
-			HalfEdgeIter he = m.halfedges_begin();
+		cout << "Building priority queue"<< endl;
-			for(;he != m.halfedges_end(); ++he)
-				he->touched = 0;
-			he = m.halfedges_begin();
-			for(;he != m.halfedges_end(); ++he)
+		for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){
-				if(he->touched == 0)
+			if(!touched[*h])
-					add_to_queue(Q, he, efun);
+				add_to_queue(m, touched, Q, *h, efun);
+		}
+		cout << "Emptying priority queue of size: " << Q.size() << " ";
 		while(!Q.empty())
-			{
+		{
+			if(Q.size() % 1000 == 0)
+				cout << ".";
+			if(Q.size() % 10000 == 0)
+				cout << Q.size();
-				PQElem elem = Q.top();
+			PQElem elem = Q.top();
-				Q.pop();
+			Q.pop();
-				HalfEdgeIter he = elem.he;
+			if(touched[elem.h] != elem.time)
+				continue;
-				if(he->touched == elem.time)
+			if(!precond_flip_edge(m, elem.h))
+				continue;
-					{
-						if(m.flip(he))
+			m.flip_edge(elem.h);
-							{
+			Walker w = m.walker(elem.h);
-								add_one_ring_to_queue(Q, he->vert, efun);
+			add_one_ring_to_queue(m, touched, Q, w.vertex(), efun);
-								add_one_ring_to_queue(Q, he->next->vert, efun);
+			add_one_ring_to_queue(m, touched, Q, w.next().vertex(), efun);
-								add_one_ring_to_queue(Q, he->opp->vert, efun);
+			add_one_ring_to_queue(m, touched, Q, w.opp().vertex(), efun);
-								add_one_ring_to_queue(Q, he->opp->next->vert, efun);
+			add_one_ring_to_queue(m, touched, Q, w.opp().next().vertex(), efun);
-							}
-					}
-			}
+		}
+		cout << endl;
+	}
-	void simulated_annealing_optimization(Manifold& m,
+	void simulated_annealing_optimization(Manifold& m, const EnergyFun& efun, int max_iter)
-																				const EnergyFun& efun,
-																				int max_iter)
+	{
 		gel_srand(0);
 		int swaps;
 		int iter = 0;
 		double T = 1;
 		double tmin=0;
-		double tsum=0;
+		{
-			HalfEdgeIter he = m.halfedges_begin();
+			for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){
-			for(;he != m.halfedges_end(); ++he)
-				{
-					if(!is_boundary(he))
+				if(boundary(m, *h))
-						{
+					continue;
-							double e = efun.delta_energy(he);
+				double e = efun.delta_energy(m, *h);
-							tmin = s_min(e, tmin);
+				tmin = s_min(e, tmin);
-							tsum += efun.energy(he);
-						}
-				}
+			}
+		}
-		cout << "Initial energy:" << tsum << endl;
+		if (tmin < 0.0)
-		T = -2*tmin;
+			T = -2*tmin;
-		if(max_iter>0) do
+		if(max_iter>0){
-			{
+			do{
-				vector<HalfEdgeIter>  halfedges;
+				cout << "Temperature : " << T << endl;
-				HalfEdgeIter he = m.halfedges_begin();
+				vector<HalfEdgeID>  halfedges;
-				for(;he != m.halfedges_end(); ++he)
+				for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){
-					{
+					if(boundary(m, *h))
-						if(!is_boundary(he))
+						continue;
-							halfedges.push_back(he);
+					halfedges.push_back(*h);
-					}
+				}
 				random_shuffle(halfedges.begin(), halfedges.end());
 				swaps = 0;
-				for(size_t i=0;i<halfedges.size();++i)
+				for(size_t i = 0; i < halfedges.size(); ++i){
-					{
-						he = halfedges[i];
+					HalfEdgeID h = halfedges[i];
-						// 							if(!would_flip(he))
-						DihedralEnergy dih_en;
+					DihedralEnergy dih_en;
-						if(dih_en.min_angle(he) > -0.4)
+					double dma = dih_en.min_angle(m, h);
-							{
+					if(dma > -0.4){
-								double delta = efun.delta_energy(he);
+						double delta = efun.delta_energy(m, h);
-								if(delta < -1e-8)
+						if(delta < -1e-8){
-									{
-										if(m.flip(he))
+							if(precond_flip_edge(m, h)){
-											++swaps;
-									}
-								else
-									{
-										delta = max(1e-8, delta);
-										double prob = min(0.9999, exp(-delta/T));
-										if(gel_rand()/double(GEL_RAND_MAX) < prob)
-											{
-												if(m.flip(he))
+								m.flip_edge(h);
-													++swaps;
+								++swaps;
-											}
-									}
+							}
+						}
+						else{
+							delta = max(1e-8, delta);
+							double prob = min(0.9999, exp(-delta/T));
+							if(gel_rand()/double(GEL_RAND_MAX) < prob){
+								if(precond_flip_edge(m, h)){
+									m.flip_edge(h);
+									++swaps;
+								}
+							}
+						}
+					}
+				}
 				cout << "Swaps = " << swaps << " T = " << T << endl;
 				if(iter % 5 == 0 && iter > 0)
 					T *= 0.9;
+			}
-		while(++iter < max_iter && swaps);
+			while(++iter < max_iter && swaps);
+		}
 		cout << "Iterations "  << iter << endl;
+	}
 	void minimize_dihedral_angle(Manifold& m,
-															 int iter,
+								 int iter,
-															 bool anneal,
+								 bool anneal,
-															 bool alpha,
+								 bool alpha,
-															 double gamma)
+								 double gamma)
+	{
 		DihedralEnergy energy_fun(gamma, alpha);
 		if(anneal)
 			simulated_annealing_optimization(m, energy_fun, iter);
 		else
 			priority_queue_optimization(m, energy_fun);
+	}
 	void randomize_mesh(Manifold& m, int max_iter)
+	{
 		RandomEnergy energy_fun;
 		simulated_annealing_optimization(m, energy_fun, max_iter);
+	}
 	void minimize_curvature(Manifold& m, bool anneal)
+	{
 		CurvatureEnergy energy_fun;
 		if(anneal)
 			simulated_annealing_optimization(m, energy_fun);
 		else
 			priority_queue_optimization(m, energy_fun);
+	}
 	void minimize_gauss_curvature(Manifold& m, bool anneal)
+	{
 		GaussCurvatureEnergy energy_fun;
 		if(anneal)
 			simulated_annealing_optimization(m, energy_fun);
 		else
 			priority_queue_optimization(m, energy_fun);
+	}
 	void maximize_min_angle(Manifold& m, float thresh, bool anneal)
+	{
 		MinAngleEnergy energy_fun(thresh);
 		if(anneal)
 			simulated_annealing_optimization(m, energy_fun);
 		else
 			priority_queue_optimization(m, energy_fun);
+	}
 	void optimize_valency(Manifold& m, bool anneal)
+	{
 		ValencyEnergy energy_fun;
 		if(anneal)
 			simulated_annealing_optimization(m, energy_fun);
 		else
 			priority_queue_optimization(m, energy_fun);
+	}
+}

Subversion Repositories gelsvn

(root)/trunk/src/HMesh/mesh_optimization.cpp @ 601 – Rev 417 → 595