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

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

Subversion Repositories gelsvn

(root)/trunk/src/HMesh/mesh_optimization.cpp @ 346 – Rev 182