Subversion Repositories gelsvn

Rev

Rev 59 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

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