Subversion Repositories gelsvn

Rev

Rev 39 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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