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

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

Subversion Repositories gelsvn

(root)/trunk/src/HMesh/mesh_optimization.cpp – Rev 599