WebSVN – gelsvn – Blame – /trunk/src/HMeshUtil/mesh_optimization.cpp

Rev	Author	Line No.	Line
39	bj	1	`#include <cfloat>`
		2	`#include "CGLA/Vec3d.h"`
		3	`#include <queue>`
		4	`#include "mesh_optimization.h"`
		5	`#include "HMesh/FaceCirculator.h"`
		6	`#include "HMesh/VertexCirculator.h"`
		7
		8	`using namespace std;`
		9	`using namespace CGLA;`
		10	`using namespace HMesh;`
		11
		12	`namespace HMeshUtil`
		13	`{`
		14	`// Small utility functions`
		15	`namespace`
		16	`{`
		17	`class LineSeg`
		18	`{`
		19	`const CGLA::Vec3d p0;`
		20	`const CGLA::Vec3d p1;`
		21	`const float l;`
		22	`const CGLA::Vec3d dir;`
		23	`public:`
		24
		25	`LineSeg(const CGLA::Vec3d& _p0, const CGLA::Vec3d& _p1):`
		26	`p0(_p0), p1(_p1), l(length(p1-p0)), dir((p1-p0)/l) {}`
		27
		28	`bool inseg(const CGLA::Vec3d& p) const`
		29	`{`
		30	`double t = dot(dir, p-p0);`
		31	`if(t<0)`
		32	`return false;`
		33	`if(t>l)`
		34	`return false;`
		35	`return true;`
		36	`}`
		37	`};`
		38
		39	`Vec3d compute_normal(Vec3d* v)`
		40	`{`
		41	`Vec3d norm;`
		42	`for(int i=0;i<4;++i)`
		43	`{`
		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]);`
		46	`norm[2] += (v[i][0]-v[(i+1)%4][0])*(v[i][1]+v[(i+1)%4][1]);`
		47	`}`
		48	`float l = norm.length();`
		49	`if(l>0.0f)`
		50	`norm /= l;`
		51	`return norm;`
		52	`}`
		53
		54	`bool would_flip(HalfEdgeIter h)`
		55	`{`
		56	`HalfEdgeIter ho = h->opp;`
		57	`Vec3d v[4];`
		58	`v[0]=Vec3d(h->vert->get_pos());`
		59	`v[1]=Vec3d(ho->vert->get_pos());`
		60	`v[2]=Vec3d(h->next->vert->get_pos());`
		61	`v[3]=Vec3d(ho->next->vert->get_pos());`
		62
		63	`Vec3d dir(0,0,-1);//compute_normal(v);`
		64
		65	`Vec3d n1a = cross(v[3]-v[0], v[2]-v[0]);`
		66	`Vec3d n2a = cross(v[2]-v[1], v[3]-v[1]);`
		67
		68	`if(dot(normalize(n1a),dir) < 0)`
		69	`return true;`
		70	`if(dot(normalize(n2a),dir) < 0)`
		71	`return true;`
		72	`return false;`
		73	`}`
		74	`}`
		75
		76	`// Energy functionals`
		77	`namespace`
		78	`{`
		79	`class EnergyFun`
		80	`{`
		81	`public:`
		82	`virtual double delta_energy(HalfEdgeIter) const = 0;`
		83	`virtual double energy(HalfEdgeIter) const {return 0;}`
		84	`};`
		85
		86
		87	`class ValencyEnergy: public EnergyFun`
		88	`{`
		89	`public:`
		90	`double delta_energy(HalfEdgeIter he) const`
		91	`{`
		92	`FaceIter f = he->face;`
		93
		94	`VertexIter v1 = he->opp->vert;`
		95	`VertexIter v2 = he->vert;`
		96	`VertexIter vo1 = he->next->vert;`
		97	`VertexIter vo2 = he->opp->next->vert;`
		98
		99	`int val1 = valency(v1);`
		100	`int val2 = valency(v2);`
		101	`int valo1 = valency(vo1);`
		102	`int valo2 = valency(vo2);`
		103
		104	`// The optimal valency is four for a boundary vertex`
		105	`// and six elsewhere.`
		106	`int t1 = v1->is_boundary() ? 4 : 6;`
		107	`int t2 = v2->is_boundary() ? 4 : 6;`
		108	`int to1 = vo1->is_boundary() ? 4 : 6;`
		109	`int to2 = vo2->is_boundary() ? 4 : 6;`
		110
		111	`int before =`
		112	`sqr(val1-t1)+sqr(val2-t2)+`
		113	`sqr(valo1-to1)+sqr(valo2-to2);`
		114	`int after =`
		115	`sqr(valo1+1-to1)+sqr(val1-1-t1)+`
		116	`sqr(val2-1-t2)+sqr(valo2+1-to2);`
		117
		118	`return static_cast<double>(after-before);`
		119	`}`
		120	`};`
		121
		122	`class RandomEnergy: public EnergyFun`
		123	`{`
		124	`public:`
		125	`double delta_energy(HalfEdgeIter he) const`
		126	`{`
		127	`return static_cast<double>(rand()/static_cast<float>(RAND_MAX));`
		128	`}`
		129	`};`
		130
		131	`class MinAngleEnergy: public EnergyFun`
		132	`{`
		133	`double min_angle(const Vec3d& v0,const Vec3d& v1,const Vec3d& v2) const`
		134	`{`
		135	`Vec3d a = normalize(v1-v0);`
		136	`Vec3d b = normalize(v2-v1);`
		137	`Vec3d c = normalize(v0-v2);`
		138
		139	`return min(dot(a,-c), min(dot(b,-a), dot(c,-b)));`
		140	`}`
		141
		142	`double thresh;`
		143
		144	`public:`
		145
		146	`MinAngleEnergy(double _thresh): thresh(_thresh) {}`
		147
		148	`double delta_energy(HalfEdgeIter he) const`
		149	`{`
		150	`FaceIter f = he->face;`
		151	`FaceIter f2 = he->opp->face;`
		152	`double t = dot(get_normal(f), get_normal(f2));`
		153	`if(t> thresh \|\| t < -thresh)`
		154	`{`
		155	`Vec3d v0(he->vert->get_pos());`
		156	`Vec3d v1(he->next->vert->get_pos());`
		157	`Vec3d v2(he->opp->vert->get_pos());`
		158	`Vec3d v3(he->opp->next->vert->get_pos());`
		159
		160	`double before = min(min_angle(v0,v1,v2), min_angle(v0,v2,v3));`
		161	`double after = min(min_angle(v0,v1,v3), min_angle(v1,v2,v3));`
		162	`return -(after-before);`
		163	`}`
		164	`return DBL_MAX;`
		165	`}`
		166	`};`
		167
		168	`class DihedralEnergy: public EnergyFun`
		169	`{`
		170	`const double gamma;`
		171	`const bool use_alpha;`
		172
		173	`double cos_ang(const Vec3d& n1, const Vec3d& n2) const`
		174	`{`
		175	`return s_max(-1.0, s_min(1.0, dot(n1, n2)));`
		176	`}`
		177
		178	`double edge_alpha_energy(Vec3d v1, Vec3d v2, double ca) const`
		179	`{`
		180	`return pow(length(v1-v2)*(acos(ca)), 1.0f/gamma);`
		181	`}`
		182
		183	`double edge_c_energy(Vec3d v1, Vec3d v2, double ca) const`
		184	`{`
		185	`return pow(length(v1-v2)*(1-ca), 1.0f/gamma);`
		186	`}`
		187
		188	`mutable double ab_12;`
		189	`mutable double ab_a1;`
		190	`mutable double ab_b1;`
		191	`mutable double ab_2c;`
		192	`mutable double ab_2d;`
		193
		194	`mutable double aa_12;`
		195	`mutable double aa_b1;`
		196	`mutable double aa_c1;`
		197	`mutable double aa_2a;`
		198	`mutable double aa_2d;`
		199
		200	`public:`
		201
		202	`DihedralEnergy(double _gamma = 4.0, bool _use_alpha=false):`
		203	`gamma(_gamma), use_alpha(_use_alpha) {}`
		204
		205	`void compute_angles(HalfEdgeIter he) const`
		206	`{`
		207	`Vec3d va(he->vert->get_pos());`
		208	`Vec3d vb(he->opp->vert->get_pos());`
		209	`Vec3d vc(he->next->vert->get_pos());`
		210	`Vec3d vd(he->opp->next->vert->get_pos());`
		211
		212	`FaceIter fa = he->next->opp->face;`
		213	`FaceIter fb = he->next->next->opp->face;`
		214	`FaceIter fc = he->opp->next->opp->face;`
		215	`FaceIter fd = he->opp->next->next->opp->face;`
		216
		217	`Vec3d n1 = normalize(cross(vc-va, vb-va));`
		218	`Vec3d n2 = normalize(cross(vb-va, vd-va));`
		219
		220	`Vec3d na = fa==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fa));`
		221	`Vec3d nb = fb==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fb));`
		222	`Vec3d nc = fc==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fc));`
		223	`Vec3d nd = fd==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fd));`
		224
		225
		226	`Vec3d fn1 = normalize(cross(vb-vc, vd-vc));`
		227	`Vec3d fn2 = normalize(cross(vd-vc, va-vc));`
		228
		229	`ab_12 = cos_ang(n1,n2);`
		230	`ab_a1 = cos_ang(na,n1);`
		231	`ab_b1 = cos_ang(nb,n1);`
		232	`ab_2c = cos_ang(n2,nc);`
		233	`ab_2d = cos_ang(n2,nd);`
		234
		235	`aa_12 = cos_ang(fn1,fn2);`
		236	`aa_b1 = cos_ang(nb,fn1);`
		237	`aa_c1 = cos_ang(nc, fn1);`
		238	`aa_2a = cos_ang(fn2, na);`
		239	`aa_2d = cos_ang(fn2,nd);`
		240	`}`
		241
		242	`double energy(HalfEdgeIter he) const`
		243	`{`
		244	`double a = cos_ang(Vec3d(get_normal(he->face)),`
		245	`Vec3d(get_normal(he->opp->face)));`
		246
		247	`Vec3d va(he->vert->get_pos());`
		248	`Vec3d vb(he->opp->vert->get_pos());`
		249	`if(use_alpha)`
		250	`return edge_alpha_energy(va,vb,a);`
		251
		252	`return edge_c_energy(va,vb,a);`
		253	`}`
		254
		255
		256	`double delta_energy(HalfEdgeIter he) const`
		257	`{`
		258	`compute_angles(he);`
		259
		260	`Vec3d va(he->vert->get_pos());`
		261	`Vec3d vb(he->opp->vert->get_pos());`
		262	`Vec3d vc(he->next->vert->get_pos());`
		263	`Vec3d vd(he->opp->next->vert->get_pos());`
		264
		265	`if(use_alpha)`
		266	`{`
		267	`double before =`
		268	`edge_alpha_energy(va,vb,ab_12)`
		269	`+edge_alpha_energy(va,vc,ab_a1)`
		270	`+edge_alpha_energy(vc,vb,ab_b1)`
		271	`+edge_alpha_energy(vd,vb,ab_2c)`
		272	`+edge_alpha_energy(vd,va,ab_2d);`
		273
		274	`double after =`
		275	`edge_alpha_energy(vd,vc,aa_12)`
		276	`+edge_alpha_energy(vb,vc,aa_b1)`
		277	`+edge_alpha_energy(vd,vb,aa_c1)`
		278	`+edge_alpha_energy(va,vc,aa_2a)`
		279	`+edge_alpha_energy(vd,va,aa_2d);`
		280
		281	`return after-before;`
		282	`}`
		283	`double before =`
		284	`edge_c_energy(va,vb,ab_12)`
		285	`+edge_c_energy(va,vc,ab_a1)`
		286	`+edge_c_energy(vc,vb,ab_b1)`
		287	`+edge_c_energy(vd,vb,ab_2c)`
		288	`+edge_c_energy(vd,va,ab_2d);`
		289
		290	`double after =`
		291	`edge_c_energy(vd,vc,aa_12)`
		292	`+edge_c_energy(vb,vc,aa_b1)`
		293	`+edge_c_energy(vd,vb,aa_c1)`
		294	`+edge_c_energy(va,vc,aa_2a)`
		295	`+edge_c_energy(vd,va,aa_2d);`
		296
		297	`return after-before;`
		298	`}`
		299
		300	`double min_angle(HalfEdgeIter he) const`
		301	`{`
		302	`compute_angles(he);`
		303
		304	`return s_min(s_min(s_min(s_min(aa_12, aa_b1), aa_c1), aa_2a), aa_2d);`
		305	`}`
		306	`};`
		307
		308	`class CurvatureEnergy: public EnergyFun`
		309	`{`
		310
		311	`// double abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const`
		312	`// {`
		313	`// const int N = ring.size();`
		314	`// vector<Vec3d> norms(N,Vec3d(0));`
		315	`// for(int i=0;i<N; ++i)`
		316	`// {`
		317	`// const Vec3d& vn = ring[i];`
		318	`// const Vec3d& vnn = ring[(i+1)%N];`
		319	`// Vec3d n = cross(vnn-v,vn-v);`
		320	`// double ln = length(vn);`
		321	`// if(ln> 1e-20)`
		322	`// norms[i] = n/ln;`
		323	`// }`
		324	`// double H = 0.0f;`
		325	`// for(int i=0;i<N; ++i)`
		326	`// {`
		327	`// const Vec3d& vn = ring[i];`
		328	`// double e_len = length(v-vn);`
		329	`// double beta =`
		330	`// acos(max(-1.0, min(1.0, dot(norms[i], norms[(i+N-1)%N]))));`
		331	`// H += e_len * beta;`
		332	`// }`
		333	`// return H/4.0f;`
		334	`// }`
		335
		336	`double abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const`
		337	`{`
		338	`const int N = ring.size();`
		339
		340	`double H = 0;`
		341	`for(int i=0;i<N;++i)`
		342	`{`
		343	`Vec3d vnim1 = ring[(i+N-1)%N] - v;`
		344	`Vec3d vni = ring[i] - v;`
		345	`Vec3d vnip1 = ring[(i+1)%N] - v;`
		346
		347	`Vec3d Nm = normalize(cross(vni, vnim1));`
		348	`Vec3d Np = normalize(cross(vnip1, vni));`
		349
		350	`double beta = acos(max(-1.0, min(1.0, dot(Nm, Np))));`
		351	`H += vni.length() * beta;`
		352	`}`
		353	`H /= 4;`
		354
		355	`return H;`
		356	`}`
		357
		358	`public:`
		359	`double delta_energy(HalfEdgeIter he) const`
		360	`{`
		361	`VertexIter va = he->vert;`
		362	`VertexIter vb = he->opp->vert;`
		363	`VertexIter vc = he->next->vert;`
		364	`VertexIter vd = he->opp->next->vert;`
		365
		366	`Vec3d va_pos(va->get_pos());`
		367	`Vec3d vb_pos(vb->get_pos());`
		368	`Vec3d vc_pos(vc->get_pos());`
		369	`Vec3d vd_pos(vd->get_pos());`
		370
		371	`vector<Vec3d> va_ring_bef;`
		372	`vector<Vec3d> va_ring_aft;`
		373	`vector<Vec3d> vb_ring_bef;`
		374	`vector<Vec3d> vb_ring_aft;`
		375	`vector<Vec3d> vc_ring_bef;`
		376	`vector<Vec3d> vc_ring_aft;`
		377	`vector<Vec3d> vd_ring_bef;`
		378	`vector<Vec3d> vd_ring_aft;`
		379
		380	`for(VertexCirculator c(va); !c.end(); ++c)`
		381	`{`
		382	`va_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));`
		383	`if(c.get_vertex() != vb)`
		384	`va_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));`
		385	`}`
		386	`for(VertexCirculator c(vb); !c.end(); ++c)`
		387	`{`
		388	`vb_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));`
		389	`if(c.get_vertex() != va)`
		390	`vb_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));`
		391	`}`
		392	`for(VertexCirculator c(vc); !c.end(); ++c)`
		393	`{`
		394	`vc_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));`
		395	`vc_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));`
		396	`if(c.get_vertex() == va)`
		397	`vc_ring_aft.push_back(Vec3d(vd->get_pos()));`
		398	`}`
		399	`for(VertexCirculator c(vd); !c.end(); ++c)`
		400	`{`
		401	`vd_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));`
		402	`vd_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));`
		403	`if(c.get_vertex() == vb)`
		404	`vd_ring_aft.push_back(Vec3d(vc->get_pos()));`
		405	`}`
		406	`double before =`
		407	`abs_mean_curv(va_pos,va_ring_bef) +`
		408	`abs_mean_curv(vb_pos,vb_ring_bef) +`
		409	`abs_mean_curv(vc_pos,vc_ring_bef) +`
		410	`abs_mean_curv(vd_pos,vd_ring_bef);`
		411
		412	`double after =`
		413	`abs_mean_curv(va_pos,va_ring_aft) +`
		414	`abs_mean_curv(vb_pos,vb_ring_aft) +`
		415	`abs_mean_curv(vc_pos,vc_ring_aft) +`
		416	`abs_mean_curv(vd_pos,vd_ring_aft);`
		417
		418	`return after-before;`
		419	`}`
		420	`};`
		421
		422	`class GaussCurvatureEnergy: public EnergyFun`
		423	`{`
		424
		425	`double gauss_curv(const Vec3d& v, const vector<Vec3d>& ring) const`
		426	`{`
		427	`const int N = ring.size();`
		428	`double asum=0.0f;`
		429	`double area_sum=0;`
		430	`for(int i=0;i<N; ++i)`
		431	`{`
		432	`const Vec3d& v1 = ring[i];`
		433	`const Vec3d& v2 = ring[(i+1)%N];`
		434	`Vec3d a = v1-v;`
		435	`Vec3d b = v2-v;`
		436	`asum += acos(max(-1.0, min(1.0, dot(a,b)/(length(a)*length(b)))));`
		437	`area_sum += 0.5 * length(cross(a,b));`
		438	`}`
		439	`return 3abs(2 M_PI - asum)/area_sum;`
		440	`}`
		441
		442	`public:`
		443	`double delta_energy(HalfEdgeIter he) const`
		444	`{`
		445	`VertexIter va = he->vert;`
		446	`VertexIter vb = he->opp->vert;`
		447	`VertexIter vc = he->next->vert;`
		448	`VertexIter vd = he->opp->next->vert;`
		449
		450	`Vec3d va_pos(va->get_pos());`
		451	`Vec3d vb_pos(vb->get_pos());`
		452	`Vec3d vc_pos(vc->get_pos());`
		453	`Vec3d vd_pos(vd->get_pos());`
		454
		455	`vector<Vec3d> va_ring_bef;`
		456	`vector<Vec3d> va_ring_aft;`
		457	`vector<Vec3d> vb_ring_bef;`
		458	`vector<Vec3d> vb_ring_aft;`
		459	`vector<Vec3d> vc_ring_bef;`
		460	`vector<Vec3d> vc_ring_aft;`
		461	`vector<Vec3d> vd_ring_bef;`
		462	`vector<Vec3d> vd_ring_aft;`
		463
		464	`for(VertexCirculator c(va); !c.end(); ++c)`
		465	`{`
		466	`va_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));`
		467	`if(c.get_vertex() != vb)`
		468	`va_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));`
		469	`}`
		470	`for(VertexCirculator c(vb); !c.end(); ++c)`
		471	`{`
		472	`vb_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));`
		473	`if(c.get_vertex() != va)`
		474	`vb_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));`
		475	`}`
		476	`for(VertexCirculator c(vc); !c.end(); ++c)`
		477	`{`
		478	`vc_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));`
		479	`vc_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));`
		480	`if(c.get_vertex() == va)`
		481	`vc_ring_aft.push_back(Vec3d(vd->get_pos()));`
		482	`}`
		483	`for(VertexCirculator c(vd); !c.end(); ++c)`
		484	`{`
		485	`vd_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));`
		486	`vd_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));`
		487	`if(c.get_vertex() == vb)`
		488	`vd_ring_aft.push_back(Vec3d(vc->get_pos()));`
		489	`}`
		490	`double before =`
		491	`gauss_curv(va_pos,va_ring_bef) +`
		492	`gauss_curv(vb_pos,vb_ring_bef) +`
		493	`gauss_curv(vc_pos,vc_ring_bef) +`
		494	`gauss_curv(vd_pos,vd_ring_bef);`
		495
		496	`double after =`
		497	`gauss_curv(va_pos,va_ring_aft) +`
		498	`gauss_curv(vb_pos,vb_ring_aft) +`
		499	`gauss_curv(vc_pos,vc_ring_aft) +`
		500	`gauss_curv(vd_pos,vd_ring_aft);`
		501
		502	`return after-before;`
		503	`}`
		504	`};`
		505
		506	`struct PQElem`
		507	`{`
		508	`double pri;`
		509	`HalfEdgeIter he;`
		510	`int time;`
		511
		512	`PQElem() {}`
		513	`PQElem(double _pri, HalfEdgeIter _he, int _time):`
		514	`pri(_pri), he(_he), time(_time) {}`
		515	`};`
		516
		517	`bool operator<(const PQElem& e0, const PQElem& e1)`
		518	`{`
		519	`return e0.pri > e1.pri;`
		520	`}`
		521
		522
		523	`void add_to_queue(priority_queue<PQElem>& Q,`
		524	`HalfEdgeIter he,`
		525	`const EnergyFun& efun)`
		526	`{`
		527	`if(!he->is_boundary())`
		528	`{`
		529	`double energy = efun.delta_energy(he);`
		530	`int t = he->touched + 1;`
		531	`he->touched = t;`
		532	`he->opp->touched = t;`
		533	`if((energy<0) && (t < 10000))`
		534	`Q.push(PQElem(energy, he, t));`
		535	`}`
		536	`}`
		537
		538	`void add_one_ring_to_queue(priority_queue<PQElem>& Q,`
		539	`VertexIter vi,`
		540	`const EnergyFun& efun)`
		541	`{`
		542	`for(VertexCirculator vc(vi); !vc.end(); ++vc)`
		543	`{`
		544	`add_to_queue(Q, vc.get_halfedge(), efun);`
		545	`}`
		546	`}`
		547
		548	`void priority_queue_optimization(Manifold& m, const EnergyFun& efun)`
		549	`{`
		550	`priority_queue<PQElem> Q;`
		551	`{`
		552	`HalfEdgeIter he = m.halfedges_begin();`
		553	`for(;he != m.halfedges_end(); ++he)`
		554	`he->touched = 0;`
		555
		556	`he = m.halfedges_begin();`
		557	`for(;he != m.halfedges_end(); ++he)`
		558	`if(he->touched == 0)`
		559	`add_to_queue(Q, he, efun);`
		560	`}`
		561	`double pri_old = -DBL_MAX;`
		562	`while(!Q.empty())`
		563	`{`
		564	`PQElem elem = Q.top();`
		565	`Q.pop();`
		566	`HalfEdgeIter he = elem.he;`
		567	`if(he->touched == elem.time)`
		568	`{`
		569	`if(m.flip(he))`
		570	`{`
		571	`add_one_ring_to_queue(Q, he->vert, efun);`
		572	`add_one_ring_to_queue(Q, he->next->vert, efun);`
		573	`add_one_ring_to_queue(Q, he->opp->vert, efun);`
		574	`add_one_ring_to_queue(Q, he->opp->next->vert, efun);`
		575	`}`
		576	`}`
		577	`}`
		578	`}`
		579
		580	`void simulated_annealing_optimization(Manifold& m,`
		581	`const EnergyFun& efun,`
		582	`int max_iter=10000)`
		583	`{`
		584	`srand(0);`
		585	`int swaps;`
		586	`int iter = 0;`
		587	`double T = 1;`
		588
		589	`double tmin=0;`
		590	`double tsum=0;`
		591	`{`
		592	`HalfEdgeIter he = m.halfedges_begin();`
		593	`for(;he != m.halfedges_end(); ++he)`
		594	`{`
		595	`if(!he->is_boundary())`
		596	`{`
		597	`double e = efun.delta_energy(he);`
		598	`tmin = s_min(e, tmin);`
		599	`tsum += efun.energy(he);`
		600	`}`
		601	`}`
		602	`}`
		603	`cout << "Initial energy:" << tsum << endl;`
		604	`T = -2*tmin;`
		605
		606	`if(max_iter>0) do`
		607	`{`
		608	`vector<HalfEdgeIter> halfedges;`
		609	`HalfEdgeIter he = m.halfedges_begin();`
		610	`for(;he != m.halfedges_end(); ++he)`
		611	`{`
		612	`if(!he->is_boundary())`
		613	`halfedges.push_back(he);`
		614	`}`
		615	`random_shuffle(halfedges.begin(), halfedges.end());`
		616	`swaps = 0;`
		617	`for(int i=0;i<halfedges.size();++i)`
		618	`{`
		619	`he = halfedges[i];`
		620
		621	`// if(!would_flip(he))`
		622	`DihedralEnergy dih_en;`
		623	`if(dih_en.min_angle(he) > -0.4)`
		624	`{`
		625	`double delta = efun.delta_energy(he);`
		626	`if(delta < -1e-8)`
		627	`{`
		628	`if(m.flip(he))`
		629	`++swaps;`
		630	`}`
		631	`else`
		632	`{`
		633	`delta = max(1e-8, delta);`
		634	`double prob = min(0.9999, exp(-delta/T));`
		635	`if(rand()/double(RAND_MAX) < prob)`
		636	`{`
		637	`if(m.flip(he))`
		638	`++swaps;`
		639	`}`
		640	`}`
		641	`}`
		642	`}`
		643	`cout << "Swaps = " << swaps << " T = " << T << endl;`
		644	`if(iter % 5 == 0 && iter > 0)`
		645	`T *= 0.9;`
		646	`}`
		647	`while(++iter < max_iter && swaps);`
		648	`cout << "Iterations " << iter << endl;`
		649	`}`
		650	`}`
		651
		652	`void minimize_dihedral_angle(Manifold& m,`
		653	`int iter,`
		654	`bool anneal,`
		655	`bool alpha,`
		656	`double gamma)`
		657	`{`
		658	`DihedralEnergy energy_fun(gamma, alpha);`
		659	`if(anneal)`
		660	`simulated_annealing_optimization(m, energy_fun, iter);`
		661	`else`
		662	`priority_queue_optimization(m, energy_fun);`
		663	`}`
		664
		665	`void randomize_mesh(Manifold& m, int max_iter)`
		666	`{`
		667	`RandomEnergy energy_fun;`
		668	`simulated_annealing_optimization(m, energy_fun, max_iter);`
		669	`}`
		670
		671	`void minimize_curvature(Manifold& m, bool anneal)`
		672	`{`
		673	`CurvatureEnergy energy_fun;`
		674	`if(anneal)`
		675	`simulated_annealing_optimization(m, energy_fun);`
		676	`else`
		677	`priority_queue_optimization(m, energy_fun);`
		678	`}`
		679
		680	`void minimize_gauss_curvature(Manifold& m, bool anneal)`
		681	`{`
		682	`GaussCurvatureEnergy energy_fun;`
		683	`if(anneal)`
		684	`simulated_annealing_optimization(m, energy_fun);`
		685	`else`
		686	`priority_queue_optimization(m, energy_fun);`
		687	`}`
		688
		689	`void maximize_min_angle(Manifold& m, float thresh, bool anneal)`
		690	`{`
		691	`MinAngleEnergy energy_fun(thresh);`
		692	`if(anneal)`
		693	`simulated_annealing_optimization(m, energy_fun);`
		694	`else`
		695	`priority_queue_optimization(m, energy_fun);`
		696	`}`
		697
		698	`void optimize_valency(Manifold& m, bool anneal)`
		699	`{`
		700	`ValencyEnergy energy_fun;`
		701	`if(anneal)`
		702	`simulated_annealing_optimization(m, energy_fun);`
		703	`else`
		704	`priority_queue_optimization(m, energy_fun);`
		705	`}`
		706
		707
		708
		709
		710	`}`

Subversion Repositories gelsvn

(root)/trunk/src/HMeshUtil/mesh_optimization.cpp @ 136 – Rev 39