Subversion Repositories gelsvn

Rev

Rev 417 | Rev 601 | Go to most recent revision | Show entire file | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed

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