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


#include <cfloat>
#include "CGLA/Vec3d.h"
#include <queue>
#include "mesh_optimization.h"
#include "HMesh/FaceCirculator.h"
#include "HMesh/VertexCirculator.h"
 
using namespace std;
using namespace CGLA;
using namespace HMesh;
 
namespace HMesh
{
	// Small utility functions
	namespace 
	{
		class LineSeg
		{
			const CGLA::Vec3d p0;
			const CGLA::Vec3d p1;
			const float l;
			const CGLA::Vec3d dir;
		public:
 
			LineSeg(const CGLA::Vec3d& _p0, const CGLA::Vec3d& _p1): 
				p0(_p0), p1(_p1), l(length(p1-p0)), dir((p1-p0)/l) {}
			
			bool inseg(const CGLA::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(HalfEdgeIter h)
		{
			HalfEdgeIter ho = h->opp;
			Vec3d v[4];
			v[0]=Vec3d(h->vert->pos);
			v[1]=Vec3d(ho->vert->pos);
			v[2]=Vec3d(h->next->vert->pos);
			v[3]=Vec3d(ho->next->vert->pos);
 
			Vec3d dir(0,0,-1);//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;
		}
	}
 
	// Energy functionals
	namespace
	{
 
 
		class ValencyEnergy: public EnergyFun
		{
		public:
			double delta_energy(HalfEdgeIter he) const
			{
				FaceIter f = he->face;
				
				VertexIter v1 = he->opp->vert;
				VertexIter v2 = he->vert;
				VertexIter vo1 = he->next->vert;
				VertexIter vo2 = he->opp->next->vert;
 
				int val1  = valency(v1);
				int val2  = valency(v2);
				int valo1 = valency(vo1);
				int valo2 = valency(vo2);
				
				// The optimal valency is four for a boundary vertex
				// and six elsewhere.
				int t1 = is_boundary(v1) ? 4 : 6;
				int t2 = is_boundary(v2) ? 4 : 6;
				int to1 = is_boundary(vo1) ? 4 : 6;
				int to2 = is_boundary(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(HalfEdgeIter he) const
			{
				return static_cast<double>(gel_rand()/static_cast<float>(GEL_RAND_MAX));
			}
		};
 
		class MinAngleEnergy: public EnergyFun
		{
			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);
				
				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;
				FaceIter f2 = he->opp->face;
 
				Vec3d v0(he->vert->pos);
				Vec3d v1(he->next->vert->pos);
				Vec3d v2(he->opp->vert->pos);
				Vec3d v3(he->opp->next->vert->pos);
 
				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;
			const bool use_alpha;
 
			double cos_ang(const Vec3d& n1, const Vec3d& n2) const
			{
				return s_max(-1.0, s_min(1.0, dot(n1, n2)));
			}
 
			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); 
			}
 
			mutable double ab_12;
			mutable double ab_a1;
			mutable double ab_b1;
			mutable double ab_2c;
			mutable double ab_2d;
			
			mutable double aa_12;
			mutable double aa_b1;
			mutable double aa_c1;
			mutable double aa_2a;
			mutable double aa_2d;
 
		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));
				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);
				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 energy(HalfEdgeIter he) const
			{
				double a = cos_ang(Vec3d(normal(he->face)),
													 Vec3d(normal(he->opp->face)));
 
				Vec3d va(he->vert->pos);
				Vec3d vb(he->opp->vert->pos);
				if(use_alpha)
					return edge_alpha_energy(va,vb,a);
				
				return edge_c_energy(va,vb,a);
			}
				
 
			double delta_energy(HalfEdgeIter he) const
			{
				compute_angles(he);
 
				Vec3d va(he->vert->pos);
				Vec3d vb(he->opp->vert->pos);
				Vec3d vc(he->next->vert->pos);
				Vec3d vd(he->opp->next->vert->pos);
 
				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)
							+edge_alpha_energy(va,vc,aa_2a)
							+edge_alpha_energy(vd,va,aa_2d);
								
						return after-before;
					}
				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;
			}
 
			double min_angle(HalfEdgeIter he) const
			{
				compute_angles(he);
				
				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);
			// 						if(ln> 1e-20)
			// 							norms[i] = n/ln;
			// 					}
			// 				double H = 0.0f;
			// 				for(int i=0;i<N; ++i)
			// 					{
			// 						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;
				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;
						
						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;
					}
				H /= 4;
 
				return H;
			}
 
		public:
			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;
				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));
					}
				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)
							vd_ring_aft.push_back(Vec3d(vc->pos));
					}
				double before =
					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);
 
				double after =
					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;
				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));
					}
				return 3*abs(2 * M_PI - asum)/area_sum;
			}
 
		public:
			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;
				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));
					}
				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)
							vd_ring_aft.push_back(Vec3d(vc->pos));
					}
				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;
			}
		};
 
		struct PQElem
		{
			double pri;
			HalfEdgeIter he;
			int time;
 
			PQElem() {}
			PQElem(double _pri, HalfEdgeIter _he, int _time): 
				pri(_pri), he(_he), time(_time) {}
		};
		
		bool operator<(const PQElem& e0, const PQElem& e1)
		{
			return e0.pri > e1.pri;
		}
	
 
			void add_to_queue(priority_queue<PQElem>& Q,
												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)
				{
					add_to_queue(Q, vc.get_halfedge(), efun);
				}
		}
	}
 
	void priority_queue_optimization(Manifold& m, const EnergyFun& efun)
	{
		priority_queue<PQElem> Q;
		{
			HalfEdgeIter he = m.halfedges_begin();
			for(;he != m.halfedges_end(); ++he)
				he->touched = 0;
			
			he = m.halfedges_begin();
			for(;he != m.halfedges_end(); ++he)
				if(he->touched == 0)
					add_to_queue(Q, he, efun);
		}
 
		while(!Q.empty())
			{
				PQElem elem = Q.top();
				Q.pop();
				HalfEdgeIter he = elem.he;
				if(he->touched == elem.time)
					{
						if(m.flip(he))
							{
								add_one_ring_to_queue(Q, he->vert, efun);
								add_one_ring_to_queue(Q, he->next->vert, efun);
								add_one_ring_to_queue(Q, he->opp->vert, efun);
								add_one_ring_to_queue(Q, he->opp->next->vert, efun);
							}
					}
			}
	}
 
	void simulated_annealing_optimization(Manifold& m, 
																				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(;he != m.halfedges_end(); ++he)
				{
					if(!is_boundary(he))
						{
							double e = efun.delta_energy(he);
							tmin = s_min(e, tmin);
							tsum += efun.energy(he);
						}
				}
		}
		cout << "Initial energy:" << tsum << endl;
		T = -2*tmin;
			
		if(max_iter>0) do
			{
				vector<HalfEdgeIter>  halfedges;
				HalfEdgeIter he = m.halfedges_begin();
				for(;he != m.halfedges_end(); ++he)
					{
						if(!is_boundary(he))
							halfedges.push_back(he);
					}
				random_shuffle(halfedges.begin(), halfedges.end());
				swaps = 0;
				for(size_t i=0;i<halfedges.size();++i)
					{
						he = halfedges[i];
 
						// 							if(!would_flip(he))
						DihedralEnergy dih_en;
						if(dih_en.min_angle(he) > -0.4)
							{
								double delta = efun.delta_energy(he);
								if(delta < -1e-8)
									{
										if(m.flip(he))
											++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))
													++swaps;
											}
									}
							}
					}
				cout << "Swaps = " << swaps << " T = " << T << endl;
				if(iter % 5 == 0 && iter > 0)
					T *= 0.9;
			}
		while(++iter < max_iter && swaps);
		cout << "Iterations "  << iter << endl; 
	}
	
 
	void minimize_dihedral_angle(Manifold& m, 
															 int iter,
															 bool anneal, 
															 bool alpha,
															 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);
	}
 
 
 
 
}
 

Rev 182	Rev 417
1	`#include <cfloat>`	1	`#include <cfloat>`
2	`#include "CGLA/Vec3d.h"`	2	`#include "CGLA/Vec3d.h"`
3	`#include <queue>`	3	`#include <queue>`
4	`#include "mesh_optimization.h"`	4	`#include "mesh_optimization.h"`
5	`#include "HMesh/FaceCirculator.h"`	5	`#include "HMesh/FaceCirculator.h"`
6	`#include "HMesh/VertexCirculator.h"`	6	`#include "HMesh/VertexCirculator.h"`
7		7
8	`using namespace std;`	8	`using namespace std;`
9	`using namespace CGLA;`	9	`using namespace CGLA;`
10	`using namespace HMesh;`	10	`using namespace HMesh;`
11		11
12	`namespace HMesh`	12	`namespace HMesh`
13	`{`	13	`{`
14	`// Small utility functions`	14	`// Small utility functions`
15	`namespace`	15	`namespace`
16	`{`	16	`{`
17	`class LineSeg`	17	`class LineSeg`
18	`{`	18	`{`
19	`const CGLA::Vec3d p0;`	19	`const CGLA::Vec3d p0;`
20	`const CGLA::Vec3d p1;`	20	`const CGLA::Vec3d p1;`
21	`const float l;`	21	`const float l;`
22	`const CGLA::Vec3d dir;`	22	`const CGLA::Vec3d dir;`
23	`public:`	23	`public:`
24		24
25	`LineSeg(const CGLA::Vec3d& _p0, const CGLA::Vec3d& _p1):`	25	`LineSeg(const CGLA::Vec3d& _p0, const CGLA::Vec3d& _p1):`
26	`p0(_p0), p1(_p1), l(length(p1-p0)), dir((p1-p0)/l) {}`	26	`p0(_p0), p1(_p1), l(length(p1-p0)), dir((p1-p0)/l) {}`
27		27
28	`bool inseg(const CGLA::Vec3d& p) const`	28	`bool inseg(const CGLA::Vec3d& p) const`
29	`{`	29	`{`
30	`double t = dot(dir, p-p0);`	30	`double t = dot(dir, p-p0);`
31	`if(t<0)`	31	`if(t<0)`
32	`return false;`	32	`return false;`
33	`if(t>l)`	33	`if(t>l)`
34	`return false;`	34	`return false;`
35	`return true;`	35	`return true;`
36	`}`	36	`}`
37	`};`	37	`};`
38		38
39	`Vec3d compute_normal(Vec3d* v)`	39	`Vec3d compute_normal(Vec3d* v)`
40	`{`	40	`{`
41	`Vec3d norm;`	41	`Vec3d norm;`
42	`for(int i=0;i<4;++i)`	42	`for(int i=0;i<4;++i)`
43	`{`	43	`{`
44	`norm[0] += (v[i][1]-v[(i+1)%4][1])*(v[i][2]+v[(i+1)%4][2]);`	44	`norm[0] += (v[i][1]-v[(i+1)%4][1])*(v[i][2]+v[(i+1)%4][2]);`
45	`norm[1] += (v[i][2]-v[(i+1)%4][2])*(v[i][0]+v[(i+1)%4][0]);`	45	`norm[1] += (v[i][2]-v[(i+1)%4][2])*(v[i][0]+v[(i+1)%4][0]);`
46	`norm[2] += (v[i][0]-v[(i+1)%4][0])*(v[i][1]+v[(i+1)%4][1]);`	46	`norm[2] += (v[i][0]-v[(i+1)%4][0])*(v[i][1]+v[(i+1)%4][1]);`
47	`}`	47	`}`
48	`float l = norm.length();`	48	`float l = norm.length();`
49	`if(l>0.0f)`	49	`if(l>0.0f)`
50	`norm /= l;`	50	`norm /= l;`
51	`return norm;`	51	`return norm;`
52	`}`	52	`}`
53		53
54	`bool would_flip(HalfEdgeIter h)`	54	`bool would_flip(HalfEdgeIter h)`
55	`{`	55	`{`
56	`HalfEdgeIter ho = h->opp;`	56	`HalfEdgeIter ho = h->opp;`
57	`Vec3d v[4];`	57	`Vec3d v[4];`
58	`v[0]=Vec3d(h->vert->pos);`	58	`v[0]=Vec3d(h->vert->pos);`
59	`v[1]=Vec3d(ho->vert->pos);`	59	`v[1]=Vec3d(ho->vert->pos);`
60	`v[2]=Vec3d(h->next->vert->pos);`	60	`v[2]=Vec3d(h->next->vert->pos);`
61	`v[3]=Vec3d(ho->next->vert->pos);`	61	`v[3]=Vec3d(ho->next->vert->pos);`
62		62
63	`Vec3d dir(0,0,-1);//compute_normal(v);`	63	`Vec3d dir(0,0,-1);//compute_normal(v);`
64		64
65	`Vec3d n1a = cross(v[3]-v[0], v[2]-v[0]);`	65	`Vec3d n1a = cross(v[3]-v[0], v[2]-v[0]);`
66	`Vec3d n2a = cross(v[2]-v[1], v[3]-v[1]);`	66	`Vec3d n2a = cross(v[2]-v[1], v[3]-v[1]);`
67		67
68	`if(dot(normalize(n1a),dir) < 0)`	68	`if(dot(normalize(n1a),dir) < 0)`
69	`return true;`	69	`return true;`
70	`if(dot(normalize(n2a),dir) < 0)`	70	`if(dot(normalize(n2a),dir) < 0)`
71	`return true;`	71	`return true;`
72	`return false;`	72	`return false;`
73	`}`	73	`}`
74	`}`	74	`}`
75		75
76	`// Energy functionals`	76	`// Energy functionals`
77	`namespace`	77	`namespace`
78	`{`	78	`{`
79		79
80		80
81	`class ValencyEnergy: public EnergyFun`	81	`class ValencyEnergy: public EnergyFun`
82	`{`	82	`{`
83	`public:`	83	`public:`
84	`double delta_energy(HalfEdgeIter he) const`	84	`double delta_energy(HalfEdgeIter he) const`
85	`{`	85	`{`
86	`FaceIter f = he->face;`	86	`FaceIter f = he->face;`
87		87
88	`VertexIter v1 = he->opp->vert;`	88	`VertexIter v1 = he->opp->vert;`
89	`VertexIter v2 = he->vert;`	89	`VertexIter v2 = he->vert;`
90	`VertexIter vo1 = he->next->vert;`	90	`VertexIter vo1 = he->next->vert;`
91	`VertexIter vo2 = he->opp->next->vert;`	91	`VertexIter vo2 = he->opp->next->vert;`
92		92
93	`int val1 = valency(v1);`	93	`int val1 = valency(v1);`
94	`int val2 = valency(v2);`	94	`int val2 = valency(v2);`
95	`int valo1 = valency(vo1);`	95	`int valo1 = valency(vo1);`
96	`int valo2 = valency(vo2);`	96	`int valo2 = valency(vo2);`
97		97
98	`// The optimal valency is four for a boundary vertex`	98	`// The optimal valency is four for a boundary vertex`
99	`// and six elsewhere.`	99	`// and six elsewhere.`
100	`int t1 = is_boundary(v1) ? 4 : 6;`	100	`int t1 = is_boundary(v1) ? 4 : 6;`
101	`int t2 = is_boundary(v2) ? 4 : 6;`	101	`int t2 = is_boundary(v2) ? 4 : 6;`
102	`int to1 = is_boundary(vo1) ? 4 : 6;`	102	`int to1 = is_boundary(vo1) ? 4 : 6;`
103	`int to2 = is_boundary(vo2) ? 4 : 6;`	103	`int to2 = is_boundary(vo2) ? 4 : 6;`
104		104
105	`int before =`	105	`int before =`
106	`sqr(val1-t1)+sqr(val2-t2)+`	106	`sqr(val1-t1)+sqr(val2-t2)+`
107	`sqr(valo1-to1)+sqr(valo2-to2);`	107	`sqr(valo1-to1)+sqr(valo2-to2);`
108	`int after =`	108	`int after =`
109	`sqr(valo1+1-to1)+sqr(val1-1-t1)+`	109	`sqr(valo1+1-to1)+sqr(val1-1-t1)+`
110	`sqr(val2-1-t2)+sqr(valo2+1-to2);`	110	`sqr(val2-1-t2)+sqr(valo2+1-to2);`
111		111
112	`return static_cast<double>(after-before);`	112	`return static_cast<double>(after-before);`
113	`}`	113	`}`
114	`};`	114	`};`
115		115
116	`class RandomEnergy: public EnergyFun`	116	`class RandomEnergy: public EnergyFun`
117	`{`	117	`{`
118	`public:`	118	`public:`
119	`double delta_energy(HalfEdgeIter he) const`	119	`double delta_energy(HalfEdgeIter he) const`
120	`{`	120	`{`
121	`return static_cast<double>(rand()/static_cast<float>(RAND_MAX));`	121	`return static_cast<double>(gel_rand()/static_cast<float>(GEL_RAND_MAX));`
122	`}`	122	`}`
123	`};`	123	`};`
124		124
125	`class MinAngleEnergy: public EnergyFun`	125	`class MinAngleEnergy: public EnergyFun`
126	`{`	126	`{`
127	`double min_angle(const Vec3d& v0,const Vec3d& v1,const Vec3d& v2) const`	127	`double min_angle(const Vec3d& v0,const Vec3d& v1,const Vec3d& v2) const`
128	`{`	128	`{`
129	`Vec3d a = normalize(v1-v0);`	129	`Vec3d a = normalize(v1-v0);`
130	`Vec3d b = normalize(v2-v1);`	130	`Vec3d b = normalize(v2-v1);`
131	`Vec3d c = normalize(v0-v2);`	131	`Vec3d c = normalize(v0-v2);`
132		132
133	`return min(dot(a,-c), min(dot(b,-a), dot(c,-b)));`	133	`return min(dot(a,-c), min(dot(b,-a), dot(c,-b)));`
134	`}`	134	`}`
135		135
136	`double thresh;`	136	`double thresh;`
137		137
138	`public:`	138	`public:`
139		139
140	`MinAngleEnergy(double _thresh): thresh(_thresh) {}`	140	`MinAngleEnergy(double _thresh): thresh(_thresh) {}`
141		141
142	`double delta_energy(HalfEdgeIter he) const`	142	`double delta_energy(HalfEdgeIter he) const`
143	`{`	143	`{`
144	`FaceIter f = he->face;`	144	`FaceIter f = he->face;`
145	`FaceIter f2 = he->opp->face;`	145	`FaceIter f2 = he->opp->face;`
146		146
147	`Vec3d v0(he->vert->pos);`	147	`Vec3d v0(he->vert->pos);`
148	`Vec3d v1(he->next->vert->pos);`	148	`Vec3d v1(he->next->vert->pos);`
149	`Vec3d v2(he->opp->vert->pos);`	149	`Vec3d v2(he->opp->vert->pos);`
150	`Vec3d v3(he->opp->next->vert->pos);`	150	`Vec3d v3(he->opp->next->vert->pos);`
151		151
152	`Vec3d n1a = normalize(cross(v1-v0,v3-v0));`	152	`Vec3d n1a = normalize(cross(v1-v0,v3-v0));`
153	`Vec3d n2a = normalize(cross(v3-v2,v1-v2));`	153	`Vec3d n2a = normalize(cross(v3-v2,v1-v2));`
154	`if(dot(n1a,n2a)> thresh)`	154	`if(dot(n1a,n2a)> thresh)`
155	`{`	155	`{`
156		156
157	`double before = min(min_angle(v0,v1,v2), min_angle(v0,v2,v3));`	157	`double before = min(min_angle(v0,v1,v2), min_angle(v0,v2,v3));`
158	`double after = min(min_angle(v0,v1,v3), min_angle(v1,v2,v3));`	158	`double after = min(min_angle(v0,v1,v3), min_angle(v1,v2,v3));`
159	`return -(after-before);`	159	`return -(after-before);`
160	`}`	160	`}`
161	`return DBL_MAX;`	161	`return DBL_MAX;`
162	`}`	162	`}`
163	`};`	163	`};`
164		164
165	`class DihedralEnergy: public EnergyFun`	165	`class DihedralEnergy: public EnergyFun`
166	`{`	166	`{`
167	`const double gamma;`	167	`const double gamma;`
168	`const bool use_alpha;`	168	`const bool use_alpha;`
169		169
170	`double cos_ang(const Vec3d& n1, const Vec3d& n2) const`	170	`double cos_ang(const Vec3d& n1, const Vec3d& n2) const`
171	`{`	171	`{`
172	`return s_max(-1.0, s_min(1.0, dot(n1, n2)));`	172	`return s_max(-1.0, s_min(1.0, dot(n1, n2)));`
173	`}`	173	`}`
174		174
175	`double edge_alpha_energy(Vec3d v1, Vec3d v2, double ca) const`	175	`double edge_alpha_energy(Vec3d v1, Vec3d v2, double ca) const`
176	`{`	176	`{`
177	`return pow(length(v1-v2)*(acos(ca)), 1.0f/gamma);`	177	`return pow(length(v1-v2)*(acos(ca)), 1.0f/gamma);`
178	`}`	178	`}`
179		179
180	`double edge_c_energy(Vec3d v1, Vec3d v2, double ca) const`	180	`double edge_c_energy(Vec3d v1, Vec3d v2, double ca) const`
181	`{`	181	`{`
182	`return pow(length(v1-v2)*(1-ca), 1.0f/gamma);`	182	`return pow(length(v1-v2)*(1-ca), 1.0f/gamma);`
183	`}`	183	`}`
184		184
185	`mutable double ab_12;`	185	`mutable double ab_12;`
186	`mutable double ab_a1;`	186	`mutable double ab_a1;`
187	`mutable double ab_b1;`	187	`mutable double ab_b1;`
188	`mutable double ab_2c;`	188	`mutable double ab_2c;`
189	`mutable double ab_2d;`	189	`mutable double ab_2d;`
190		190
191	`mutable double aa_12;`	191	`mutable double aa_12;`
192	`mutable double aa_b1;`	192	`mutable double aa_b1;`
193	`mutable double aa_c1;`	193	`mutable double aa_c1;`
194	`mutable double aa_2a;`	194	`mutable double aa_2a;`
195	`mutable double aa_2d;`	195	`mutable double aa_2d;`
196		196
197	`public:`	197	`public:`
198		198
199	`DihedralEnergy(double _gamma = 4.0, bool _use_alpha=false):`	199	`DihedralEnergy(double _gamma = 4.0, bool _use_alpha=false):`
200	`gamma(_gamma), use_alpha(_use_alpha) {}`	200	`gamma(_gamma), use_alpha(_use_alpha) {}`
201		201
202	`void compute_angles(HalfEdgeIter he) const`	202	`void compute_angles(HalfEdgeIter he) const`
203	`{`	203	`{`
204	`Vec3d va(he->vert->pos);`	204	`Vec3d va(he->vert->pos);`
205	`Vec3d vb(he->opp->vert->pos);`	205	`Vec3d vb(he->opp->vert->pos);`
206	`Vec3d vc(he->next->vert->pos);`	206	`Vec3d vc(he->next->vert->pos);`
207	`Vec3d vd(he->opp->next->vert->pos);`	207	`Vec3d vd(he->opp->next->vert->pos);`
208		208
209	`FaceIter fa = he->next->opp->face;`	209	`FaceIter fa = he->next->opp->face;`
210	`FaceIter fb = he->next->next->opp->face;`	210	`FaceIter fb = he->next->next->opp->face;`
211	`FaceIter fc = he->opp->next->opp->face;`	211	`FaceIter fc = he->opp->next->opp->face;`
212	`FaceIter fd = he->opp->next->next->opp->face;`	212	`FaceIter fd = he->opp->next->next->opp->face;`
213		213
214	`Vec3d n1 = normalize(cross(vc-va, vb-va));`	214	`Vec3d n1 = normalize(cross(vc-va, vb-va));`
215	`Vec3d n2 = normalize(cross(vb-va, vd-va));`	215	`Vec3d n2 = normalize(cross(vb-va, vd-va));`
216		216
217	`Vec3d na = fa==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fa));`	217	`Vec3d na = fa==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fa));`
218	`Vec3d nb = fb==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fb));`	218	`Vec3d nb = fb==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fb));`
219	`Vec3d nc = fc==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fc));`	219	`Vec3d nc = fc==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fc));`
220	`Vec3d nd = fd==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fd));`	220	`Vec3d nd = fd==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fd));`
221		221
222		222
223	`Vec3d fn1 = normalize(cross(vb-vc, vd-vc));`	223	`Vec3d fn1 = normalize(cross(vb-vc, vd-vc));`
224	`Vec3d fn2 = normalize(cross(vd-vc, va-vc));`	224	`Vec3d fn2 = normalize(cross(vd-vc, va-vc));`
225		225
226	`ab_12 = cos_ang(n1,n2);`	226	`ab_12 = cos_ang(n1,n2);`
227	`ab_a1 = cos_ang(na,n1);`	227	`ab_a1 = cos_ang(na,n1);`
228	`ab_b1 = cos_ang(nb,n1);`	228	`ab_b1 = cos_ang(nb,n1);`
229	`ab_2c = cos_ang(n2,nc);`	229	`ab_2c = cos_ang(n2,nc);`
230	`ab_2d = cos_ang(n2,nd);`	230	`ab_2d = cos_ang(n2,nd);`
231		231
232	`aa_12 = cos_ang(fn1,fn2);`	232	`aa_12 = cos_ang(fn1,fn2);`
233	`aa_b1 = cos_ang(nb,fn1);`	233	`aa_b1 = cos_ang(nb,fn1);`
234	`aa_c1 = cos_ang(nc, fn1);`	234	`aa_c1 = cos_ang(nc, fn1);`
235	`aa_2a = cos_ang(fn2, na);`	235	`aa_2a = cos_ang(fn2, na);`
236	`aa_2d = cos_ang(fn2,nd);`	236	`aa_2d = cos_ang(fn2,nd);`
237	`}`	237	`}`
238		238
239	`double energy(HalfEdgeIter he) const`	239	`double energy(HalfEdgeIter he) const`
240	`{`	240	`{`
241	`double a = cos_ang(Vec3d(normal(he->face)),`	241	`double a = cos_ang(Vec3d(normal(he->face)),`
242	`Vec3d(normal(he->opp->face)));`	242	`Vec3d(normal(he->opp->face)));`
243		243
244	`Vec3d va(he->vert->pos);`	244	`Vec3d va(he->vert->pos);`
245	`Vec3d vb(he->opp->vert->pos);`	245	`Vec3d vb(he->opp->vert->pos);`
246	`if(use_alpha)`	246	`if(use_alpha)`
247	`return edge_alpha_energy(va,vb,a);`	247	`return edge_alpha_energy(va,vb,a);`
248		248
249	`return edge_c_energy(va,vb,a);`	249	`return edge_c_energy(va,vb,a);`
250	`}`	250	`}`
251		251
252		252
253	`double delta_energy(HalfEdgeIter he) const`	253	`double delta_energy(HalfEdgeIter he) const`
254	`{`	254	`{`
255	`compute_angles(he);`	255	`compute_angles(he);`
256		256
257	`Vec3d va(he->vert->pos);`	257	`Vec3d va(he->vert->pos);`
258	`Vec3d vb(he->opp->vert->pos);`	258	`Vec3d vb(he->opp->vert->pos);`
259	`Vec3d vc(he->next->vert->pos);`	259	`Vec3d vc(he->next->vert->pos);`
260	`Vec3d vd(he->opp->next->vert->pos);`	260	`Vec3d vd(he->opp->next->vert->pos);`
261		261
262	`if(use_alpha)`	262	`if(use_alpha)`
263	`{`	263	`{`
264	`double before =`	264	`double before =`
265	`edge_alpha_energy(va,vb,ab_12)`	265	`edge_alpha_energy(va,vb,ab_12)`
266	`+edge_alpha_energy(va,vc,ab_a1)`	266	`+edge_alpha_energy(va,vc,ab_a1)`
267	`+edge_alpha_energy(vc,vb,ab_b1)`	267	`+edge_alpha_energy(vc,vb,ab_b1)`
268	`+edge_alpha_energy(vd,vb,ab_2c)`	268	`+edge_alpha_energy(vd,vb,ab_2c)`
269	`+edge_alpha_energy(vd,va,ab_2d);`	269	`+edge_alpha_energy(vd,va,ab_2d);`
270		270
271	`double after =`	271	`double after =`
272	`edge_alpha_energy(vd,vc,aa_12)`	272	`edge_alpha_energy(vd,vc,aa_12)`
273	`+edge_alpha_energy(vb,vc,aa_b1)`	273	`+edge_alpha_energy(vb,vc,aa_b1)`
274	`+edge_alpha_energy(vd,vb,aa_c1)`	274	`+edge_alpha_energy(vd,vb,aa_c1)`
275	`+edge_alpha_energy(va,vc,aa_2a)`	275	`+edge_alpha_energy(va,vc,aa_2a)`
276	`+edge_alpha_energy(vd,va,aa_2d);`	276	`+edge_alpha_energy(vd,va,aa_2d);`
277		277
278	`return after-before;`	278	`return after-before;`
279	`}`	279	`}`
280	`double before =`	280	`double before =`
281	`edge_c_energy(va,vb,ab_12)`	281	`edge_c_energy(va,vb,ab_12)`
282	`+edge_c_energy(va,vc,ab_a1)`	282	`+edge_c_energy(va,vc,ab_a1)`
283	`+edge_c_energy(vc,vb,ab_b1)`	283	`+edge_c_energy(vc,vb,ab_b1)`
284	`+edge_c_energy(vd,vb,ab_2c)`	284	`+edge_c_energy(vd,vb,ab_2c)`
285	`+edge_c_energy(vd,va,ab_2d);`	285	`+edge_c_energy(vd,va,ab_2d);`
286		286
287	`double after =`	287	`double after =`
288	`edge_c_energy(vd,vc,aa_12)`	288	`edge_c_energy(vd,vc,aa_12)`
289	`+edge_c_energy(vb,vc,aa_b1)`	289	`+edge_c_energy(vb,vc,aa_b1)`
290	`+edge_c_energy(vd,vb,aa_c1)`	290	`+edge_c_energy(vd,vb,aa_c1)`
291	`+edge_c_energy(va,vc,aa_2a)`	291	`+edge_c_energy(va,vc,aa_2a)`
292	`+edge_c_energy(vd,va,aa_2d);`	292	`+edge_c_energy(vd,va,aa_2d);`
293		293
294	`return after-before;`	294	`return after-before;`
295	`}`	295	`}`
296		296
297	`double min_angle(HalfEdgeIter he) const`	297	`double min_angle(HalfEdgeIter he) const`
298	`{`	298	`{`
299	`compute_angles(he);`	299	`compute_angles(he);`
300		300
301	`return s_min(s_min(s_min(s_min(aa_12, aa_b1), aa_c1), aa_2a), aa_2d);`	301	`return s_min(s_min(s_min(s_min(aa_12, aa_b1), aa_c1), aa_2a), aa_2d);`
302	`}`	302	`}`
303	`};`	303	`};`
304		304
305	`class CurvatureEnergy: public EnergyFun`	305	`class CurvatureEnergy: public EnergyFun`
306	`{`	306	`{`
307		307
308	`// double abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const`	308	`// double abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const`
309	`// {`	309	`// {`
310	`// const int N = ring.size();`	310	`// const int N = ring.size();`
311	`// vector<Vec3d> norms(N,Vec3d(0));`	311	`// vector<Vec3d> norms(N,Vec3d(0));`
312	`// for(int i=0;i<N; ++i)`	312	`// for(int i=0;i<N; ++i)`
313	`// {`	313	`// {`
314	`// const Vec3d& vn = ring[i];`	314	`// const Vec3d& vn = ring[i];`
315	`// const Vec3d& vnn = ring[(i+1)%N];`	315	`// const Vec3d& vnn = ring[(i+1)%N];`
316	`// Vec3d n = cross(vnn-v,vn-v);`	316	`// Vec3d n = cross(vnn-v,vn-v);`
317	`// double ln = length(vn);`	317	`// double ln = length(vn);`
318	`// if(ln> 1e-20)`	318	`// if(ln> 1e-20)`
319	`// norms[i] = n/ln;`	319	`// norms[i] = n/ln;`
320	`// }`	320	`// }`
321	`// double H = 0.0f;`	321	`// double H = 0.0f;`
322	`// for(int i=0;i<N; ++i)`	322	`// for(int i=0;i<N; ++i)`
323	`// {`	323	`// {`
324	`// const Vec3d& vn = ring[i];`	324	`// const Vec3d& vn = ring[i];`
325	`// double e_len = length(v-vn);`	325	`// double e_len = length(v-vn);`
326	`// double beta =`	326	`// double beta =`
327	`// acos(max(-1.0, min(1.0, dot(norms[i], norms[(i+N-1)%N]))));`	327	`// acos(max(-1.0, min(1.0, dot(norms[i], norms[(i+N-1)%N]))));`
328	`// H += e_len * beta;`	328	`// H += e_len * beta;`
329	`// }`	329	`// }`
330	`// return H/4.0f;`	330	`// return H/4.0f;`
331	`// }`	331	`// }`
332		332
333	`double abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const`	333	`double abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const`
334	`{`	334	`{`
335	`const int N = ring.size();`	335	`const int N = ring.size();`
336		336
337	`double H = 0;`	337	`double H = 0;`
338	`for(int i=0;i<N;++i)`	338	`for(int i=0;i<N;++i)`
339	`{`	339	`{`
340	`Vec3d vnim1 = ring[(i+N-1)%N] - v;`	340	`Vec3d vnim1 = ring[(i+N-1)%N] - v;`
341	`Vec3d vni = ring[i] - v;`	341	`Vec3d vni = ring[i] - v;`
342	`Vec3d vnip1 = ring[(i+1)%N] - v;`	342	`Vec3d vnip1 = ring[(i+1)%N] - v;`
343		343
344	`Vec3d Nm = normalize(cross(vni, vnim1));`	344	`Vec3d Nm = normalize(cross(vni, vnim1));`
345	`Vec3d Np = normalize(cross(vnip1, vni));`	345	`Vec3d Np = normalize(cross(vnip1, vni));`
346		346
347	`double beta = acos(max(-1.0, min(1.0, dot(Nm, Np))));`	347	`double beta = acos(max(-1.0, min(1.0, dot(Nm, Np))));`
348	`H += vni.length() * beta;`	348	`H += vni.length() * beta;`
349	`}`	349	`}`
350	`H /= 4;`	350	`H /= 4;`
351		351
352	`return H;`	352	`return H;`
353	`}`	353	`}`
354		354
355	`public:`	355	`public:`
356	`double delta_energy(HalfEdgeIter he) const`	356	`double delta_energy(HalfEdgeIter he) const`
357	`{`	357	`{`
358	`VertexIter va = he->vert;`	358	`VertexIter va = he->vert;`
359	`VertexIter vb = he->opp->vert;`	359	`VertexIter vb = he->opp->vert;`
360	`VertexIter vc = he->next->vert;`	360	`VertexIter vc = he->next->vert;`
361	`VertexIter vd = he->opp->next->vert;`	361	`VertexIter vd = he->opp->next->vert;`
362		362
363	`Vec3d va_pos(va->pos);`	363	`Vec3d va_pos(va->pos);`
364	`Vec3d vb_pos(vb->pos);`	364	`Vec3d vb_pos(vb->pos);`
365	`Vec3d vc_pos(vc->pos);`	365	`Vec3d vc_pos(vc->pos);`
366	`Vec3d vd_pos(vd->pos);`	366	`Vec3d vd_pos(vd->pos);`
367		367
368	`vector<Vec3d> va_ring_bef;`	368	`vector<Vec3d> va_ring_bef;`
369	`vector<Vec3d> va_ring_aft;`	369	`vector<Vec3d> va_ring_aft;`
370	`vector<Vec3d> vb_ring_bef;`	370	`vector<Vec3d> vb_ring_bef;`
371	`vector<Vec3d> vb_ring_aft;`	371	`vector<Vec3d> vb_ring_aft;`
372	`vector<Vec3d> vc_ring_bef;`	372	`vector<Vec3d> vc_ring_bef;`
373	`vector<Vec3d> vc_ring_aft;`	373	`vector<Vec3d> vc_ring_aft;`
374	`vector<Vec3d> vd_ring_bef;`	374	`vector<Vec3d> vd_ring_bef;`
375	`vector<Vec3d> vd_ring_aft;`	375	`vector<Vec3d> vd_ring_aft;`
376		376
377	`for(VertexCirculator c(va); !c.end(); ++c)`	377	`for(VertexCirculator c(va); !c.end(); ++c)`
378	`{`	378	`{`
379	`va_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`	379	`va_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`
380	`if(c.get_vertex() != vb)`	380	`if(c.get_vertex() != vb)`
381	`va_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`	381	`va_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`
382	`}`	382	`}`
383	`for(VertexCirculator c(vb); !c.end(); ++c)`	383	`for(VertexCirculator c(vb); !c.end(); ++c)`
384	`{`	384	`{`
385	`vb_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`	385	`vb_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`
386	`if(c.get_vertex() != va)`	386	`if(c.get_vertex() != va)`
387	`vb_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`	387	`vb_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`
388	`}`	388	`}`
389	`for(VertexCirculator c(vc); !c.end(); ++c)`	389	`for(VertexCirculator c(vc); !c.end(); ++c)`
390	`{`	390	`{`
391	`vc_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`	391	`vc_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`
392	`vc_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`	392	`vc_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`
393	`if(c.get_vertex() == va)`	393	`if(c.get_vertex() == va)`
394	`vc_ring_aft.push_back(Vec3d(vd->pos));`	394	`vc_ring_aft.push_back(Vec3d(vd->pos));`
395	`}`	395	`}`
396	`for(VertexCirculator c(vd); !c.end(); ++c)`	396	`for(VertexCirculator c(vd); !c.end(); ++c)`
397	`{`	397	`{`
398	`vd_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`	398	`vd_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`
399	`vd_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`	399	`vd_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`
400	`if(c.get_vertex() == vb)`	400	`if(c.get_vertex() == vb)`
401	`vd_ring_aft.push_back(Vec3d(vc->pos));`	401	`vd_ring_aft.push_back(Vec3d(vc->pos));`
402	`}`	402	`}`
403	`double before =`	403	`double before =`
404	`abs_mean_curv(va_pos,va_ring_bef) +`	404	`abs_mean_curv(va_pos,va_ring_bef) +`
405	`abs_mean_curv(vb_pos,vb_ring_bef) +`	405	`abs_mean_curv(vb_pos,vb_ring_bef) +`
406	`abs_mean_curv(vc_pos,vc_ring_bef) +`	406	`abs_mean_curv(vc_pos,vc_ring_bef) +`
407	`abs_mean_curv(vd_pos,vd_ring_bef);`	407	`abs_mean_curv(vd_pos,vd_ring_bef);`
408		408
409	`double after =`	409	`double after =`
410	`abs_mean_curv(va_pos,va_ring_aft) +`	410	`abs_mean_curv(va_pos,va_ring_aft) +`
411	`abs_mean_curv(vb_pos,vb_ring_aft) +`	411	`abs_mean_curv(vb_pos,vb_ring_aft) +`
412	`abs_mean_curv(vc_pos,vc_ring_aft) +`	412	`abs_mean_curv(vc_pos,vc_ring_aft) +`
413	`abs_mean_curv(vd_pos,vd_ring_aft);`	413	`abs_mean_curv(vd_pos,vd_ring_aft);`
414		414
415	`return after-before;`	415	`return after-before;`
416	`}`	416	`}`
417	`};`	417	`};`
418		418
419	`class GaussCurvatureEnergy: public EnergyFun`	419	`class GaussCurvatureEnergy: public EnergyFun`
420	`{`	420	`{`
421		421
422	`double gauss_curv(const Vec3d& v, const vector<Vec3d>& ring) const`	422	`double gauss_curv(const Vec3d& v, const vector<Vec3d>& ring) const`
423	`{`	423	`{`
424	`const int N = ring.size();`	424	`const int N = ring.size();`
425	`double asum=0.0f;`	425	`double asum=0.0f;`
426	`double area_sum=0;`	426	`double area_sum=0;`
427	`for(int i=0;i<N; ++i)`	427	`for(int i=0;i<N; ++i)`
428	`{`	428	`{`
429	`const Vec3d& v1 = ring[i];`	429	`const Vec3d& v1 = ring[i];`
430	`const Vec3d& v2 = ring[(i+1)%N];`	430	`const Vec3d& v2 = ring[(i+1)%N];`
431	`Vec3d a = v1-v;`	431	`Vec3d a = v1-v;`
432	`Vec3d b = v2-v;`	432	`Vec3d b = v2-v;`
433	`asum += acos(max(-1.0, min(1.0, dot(a,b)/(length(a)*length(b)))));`	433	`asum += acos(max(-1.0, min(1.0, dot(a,b)/(length(a)*length(b)))));`
434	`area_sum += 0.5 * length(cross(a,b));`	434	`area_sum += 0.5 * length(cross(a,b));`
435	`}`	435	`}`
436	`return 3abs(2 M_PI - asum)/area_sum;`	436	`return 3abs(2 M_PI - asum)/area_sum;`
437	`}`	437	`}`
438		438
439	`public:`	439	`public:`
440	`double delta_energy(HalfEdgeIter he) const`	440	`double delta_energy(HalfEdgeIter he) const`
441	`{`	441	`{`
442	`VertexIter va = he->vert;`	442	`VertexIter va = he->vert;`
443	`VertexIter vb = he->opp->vert;`	443	`VertexIter vb = he->opp->vert;`
444	`VertexIter vc = he->next->vert;`	444	`VertexIter vc = he->next->vert;`
445	`VertexIter vd = he->opp->next->vert;`	445	`VertexIter vd = he->opp->next->vert;`
446		446
447	`Vec3d va_pos(va->pos);`	447	`Vec3d va_pos(va->pos);`
448	`Vec3d vb_pos(vb->pos);`	448	`Vec3d vb_pos(vb->pos);`
449	`Vec3d vc_pos(vc->pos);`	449	`Vec3d vc_pos(vc->pos);`
450	`Vec3d vd_pos(vd->pos);`	450	`Vec3d vd_pos(vd->pos);`
451		451
452	`vector<Vec3d> va_ring_bef;`	452	`vector<Vec3d> va_ring_bef;`
453	`vector<Vec3d> va_ring_aft;`	453	`vector<Vec3d> va_ring_aft;`
454	`vector<Vec3d> vb_ring_bef;`	454	`vector<Vec3d> vb_ring_bef;`
455	`vector<Vec3d> vb_ring_aft;`	455	`vector<Vec3d> vb_ring_aft;`
456	`vector<Vec3d> vc_ring_bef;`	456	`vector<Vec3d> vc_ring_bef;`
457	`vector<Vec3d> vc_ring_aft;`	457	`vector<Vec3d> vc_ring_aft;`
458	`vector<Vec3d> vd_ring_bef;`	458	`vector<Vec3d> vd_ring_bef;`
459	`vector<Vec3d> vd_ring_aft;`	459	`vector<Vec3d> vd_ring_aft;`
460		460
461	`for(VertexCirculator c(va); !c.end(); ++c)`	461	`for(VertexCirculator c(va); !c.end(); ++c)`
462	`{`	462	`{`
463	`va_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`	463	`va_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`
464	`if(c.get_vertex() != vb)`	464	`if(c.get_vertex() != vb)`
465	`va_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`	465	`va_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`
466	`}`	466	`}`
467	`for(VertexCirculator c(vb); !c.end(); ++c)`	467	`for(VertexCirculator c(vb); !c.end(); ++c)`
468	`{`	468	`{`
469	`vb_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`	469	`vb_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`
470	`if(c.get_vertex() != va)`	470	`if(c.get_vertex() != va)`
471	`vb_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`	471	`vb_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`
472	`}`	472	`}`
473	`for(VertexCirculator c(vc); !c.end(); ++c)`	473	`for(VertexCirculator c(vc); !c.end(); ++c)`
474	`{`	474	`{`
475	`vc_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`	475	`vc_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`
476	`vc_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`	476	`vc_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`
477	`if(c.get_vertex() == va)`	477	`if(c.get_vertex() == va)`
478	`vc_ring_aft.push_back(Vec3d(vd->pos));`	478	`vc_ring_aft.push_back(Vec3d(vd->pos));`
479	`}`	479	`}`
480	`for(VertexCirculator c(vd); !c.end(); ++c)`	480	`for(VertexCirculator c(vd); !c.end(); ++c)`
481	`{`	481	`{`
482	`vd_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`	482	`vd_ring_bef.push_back(Vec3d(c.get_vertex()->pos));`
483	`vd_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`	483	`vd_ring_aft.push_back(Vec3d(c.get_vertex()->pos));`
484	`if(c.get_vertex() == vb)`	484	`if(c.get_vertex() == vb)`
485	`vd_ring_aft.push_back(Vec3d(vc->pos));`	485	`vd_ring_aft.push_back(Vec3d(vc->pos));`
486	`}`	486	`}`
487	`double before =`	487	`double before =`
488	`gauss_curv(va_pos,va_ring_bef) +`	488	`gauss_curv(va_pos,va_ring_bef) +`
489	`gauss_curv(vb_pos,vb_ring_bef) +`	489	`gauss_curv(vb_pos,vb_ring_bef) +`
490	`gauss_curv(vc_pos,vc_ring_bef) +`	490	`gauss_curv(vc_pos,vc_ring_bef) +`
491	`gauss_curv(vd_pos,vd_ring_bef);`	491	`gauss_curv(vd_pos,vd_ring_bef);`
492		492
493	`double after =`	493	`double after =`
494	`gauss_curv(va_pos,va_ring_aft) +`	494	`gauss_curv(va_pos,va_ring_aft) +`
495	`gauss_curv(vb_pos,vb_ring_aft) +`	495	`gauss_curv(vb_pos,vb_ring_aft) +`
496	`gauss_curv(vc_pos,vc_ring_aft) +`	496	`gauss_curv(vc_pos,vc_ring_aft) +`
497	`gauss_curv(vd_pos,vd_ring_aft);`	497	`gauss_curv(vd_pos,vd_ring_aft);`
498		498
499	`return after-before;`	499	`return after-before;`
500	`}`	500	`}`

Subversion Repositories gelsvn

(root)/trunk/src/HMesh/mesh_optimization.cpp – Rev 182 → 417