Subversion Repositories gelsvn

Rev

Rev 59 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 59 Rev 62
Line 53... Line 53...
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->pos);
59 
			v[1]=Vec3d(ho->vert->get_pos());
 59 
			v[1]=Vec3d(ho->vert->pos);
60 
			v[2]=Vec3d(h->next->vert->get_pos());
 60 
			v[2]=Vec3d(h->next->vert->pos);
61 
			v[3]=Vec3d(ho->next->vert->get_pos());
 61 
			v[3]=Vec3d(ho->next->vert->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]);
Line 95... Line 95...
95 
				int valo1 = valency(vo1);
 95 
				int valo1 = valency(vo1);
96 
				int valo2 = valency(vo2);
 96 
				int valo2 = valency(vo2);
97
97
98 
				// The optimal valency is four for a boundary vertex
 98 
				// The optimal valency is four for a boundary vertex
99 
				// and six elsewhere.
 99 
				// and six elsewhere.
100 
				int t1 = v1->is_boundary() ? 4 : 6;
 100 
				int t1 = is_boundary(v1) ? 4 : 6;
101 
				int t2 = v2->is_boundary() ? 4 : 6;
 101 
				int t2 = is_boundary(v2) ? 4 : 6;
102 
				int to1 = vo1->is_boundary() ? 4 : 6;
 102 
				int to1 = is_boundary(vo1) ? 4 : 6;
103 
				int to2 = vo2->is_boundary() ? 4 : 6;
 103 
				int to2 = is_boundary(vo2) ? 4 : 6;
104
104
105 
				int before =
105 
				int before =
106 
					sqr(val1-t1)+sqr(val2-t2)+
 106 
					sqr(val1-t1)+sqr(val2-t2)+
107 
					sqr(valo1-to1)+sqr(valo2-to2);
 107 
					sqr(valo1-to1)+sqr(valo2-to2);
108 
				int after =
108 
				int after =
Line 141... Line 141...
141 
 
 141 
 
142 
			double delta_energy(HalfEdgeIter he) const
 142 
			double delta_energy(HalfEdgeIter he) const
143 
			{
 143 
			{
144 
				FaceIter f = he->face;
 144 
				FaceIter f = he->face;
145 
				FaceIter f2 = he->opp->face;
 145 
				FaceIter f2 = he->opp->face;
146 
				double t = dot(get_normal(f), get_normal(f2));
 146 
				double t = dot(normal(f), normal(f2));
147 
				if(t> thresh || t < -thresh)
 147 
				if(t> thresh || t < -thresh)
148 
					{
 148 
					{
149 
						Vec3d v0(he->vert->get_pos());
 149 
						Vec3d v0(he->vert->pos);
150 
						Vec3d v1(he->next->vert->get_pos());
 150 
						Vec3d v1(he->next->vert->pos);
151 
						Vec3d v2(he->opp->vert->get_pos());
 151 
						Vec3d v2(he->opp->vert->pos);
152 
						Vec3d v3(he->opp->next->vert->get_pos());
 152 
						Vec3d v3(he->opp->next->vert->pos);
153
153
154 
						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));
155 
						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));
156 
						return -(after-before);
 156 
						return -(after-before);
157 
					}
 157 
					}
Line 196... Line 196...
196 
			DihedralEnergy(double _gamma = 4.0, bool _use_alpha=false):
196 
			DihedralEnergy(double _gamma = 4.0, bool _use_alpha=false):
197 
				gamma(_gamma), use_alpha(_use_alpha) {}
 197 
				gamma(_gamma), use_alpha(_use_alpha) {}
198 
 
 198 
 
199 
			void compute_angles(HalfEdgeIter he) const
 199 
			void compute_angles(HalfEdgeIter he) const
200 
			{
 200 
			{
201 
				Vec3d va(he->vert->get_pos());
 201 
				Vec3d va(he->vert->pos);
202 
				Vec3d vb(he->opp->vert->get_pos());
 202 
				Vec3d vb(he->opp->vert->pos);
203 
				Vec3d vc(he->next->vert->get_pos());
 203 
				Vec3d vc(he->next->vert->pos);
204 
				Vec3d vd(he->opp->next->vert->get_pos());
 204 
				Vec3d vd(he->opp->next->vert->pos);
205 
 
 205 
 
206 
				FaceIter fa = he->next->opp->face;
 206 
				FaceIter fa = he->next->opp->face;
207 
				FaceIter fb = he->next->next->opp->face;
 207 
				FaceIter fb = he->next->next->opp->face;
208 
				FaceIter fc = he->opp->next->opp->face;
 208 
				FaceIter fc = he->opp->next->opp->face;
209 
				FaceIter fd = he->opp->next->next->opp->face;
 209 
				FaceIter fd = he->opp->next->next->opp->face;
210 
 
 210 
 
211 
				Vec3d n1 = normalize(cross(vc-va, vb-va));
 211 
				Vec3d n1 = normalize(cross(vc-va, vb-va));
212 
				Vec3d n2 = normalize(cross(vb-va, vd-va));
 212 
				Vec3d n2 = normalize(cross(vb-va, vd-va));
213 
 
 213 
 
214 
				Vec3d na = fa==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fa));
 214 
				Vec3d na = fa==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fa));
215 
				Vec3d nb = fb==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fb));
 215 
				Vec3d nb = fb==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fb));
216 
				Vec3d nc = fc==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fc));
 216 
				Vec3d nc = fc==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fc));
217 
				Vec3d nd = fd==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fd));
 217 
				Vec3d nd = fd==NULL_FACE_ITER? Vec3d(0) : Vec3d(normal(fd));
218 
 
 218 
 
219 
 
 219 
 
220 
				Vec3d fn1 = normalize(cross(vb-vc, vd-vc));
 220 
				Vec3d fn1 = normalize(cross(vb-vc, vd-vc));
221 
				Vec3d fn2 = normalize(cross(vd-vc, va-vc));
 221 
				Vec3d fn2 = normalize(cross(vd-vc, va-vc));
222 
 
 222 
 
Line 233... Line 233...
233 
				aa_2d = cos_ang(fn2,nd);
 233 
				aa_2d = cos_ang(fn2,nd);
234 
			}
 234 
			}
235 
 
 235 
 
236 
			double energy(HalfEdgeIter he) const
 236 
			double energy(HalfEdgeIter he) const
237 
			{
 237 
			{
238 
				double a = cos_ang(Vec3d(get_normal(he->face)),
 238 
				double a = cos_ang(Vec3d(normal(he->face)),
239 
													 Vec3d(get_normal(he->opp->face)));
 239 
													 Vec3d(normal(he->opp->face)));
240 
 
 240 
 
241 
				Vec3d va(he->vert->get_pos());
 241 
				Vec3d va(he->vert->pos);
242 
				Vec3d vb(he->opp->vert->get_pos());
 242 
				Vec3d vb(he->opp->vert->pos);
243 
				if(use_alpha)
 243 
				if(use_alpha)
244 
					return edge_alpha_energy(va,vb,a);
 244 
					return edge_alpha_energy(va,vb,a);
245
245
246 
				return edge_c_energy(va,vb,a);
 246 
				return edge_c_energy(va,vb,a);
247 
			}
 247 
			}
Line 249... Line 249...
249 
 
 249 
 
250 
			double delta_energy(HalfEdgeIter he) const
 250 
			double delta_energy(HalfEdgeIter he) const
251 
			{
 251 
			{
252 
				compute_angles(he);
 252 
				compute_angles(he);
253 
 
 253 
 
254 
				Vec3d va(he->vert->get_pos());
 254 
				Vec3d va(he->vert->pos);
255 
				Vec3d vb(he->opp->vert->get_pos());
 255 
				Vec3d vb(he->opp->vert->pos);
256 
				Vec3d vc(he->next->vert->get_pos());
 256 
				Vec3d vc(he->next->vert->pos);
257 
				Vec3d vd(he->opp->next->vert->get_pos());
 257 
				Vec3d vd(he->opp->next->vert->pos);
258 
 
 258 
 
259 
				if(use_alpha)
 259 
				if(use_alpha)
260 
					{
 260 
					{
261 
						double before =
261 
						double before =
262 
							edge_alpha_energy(va,vb,ab_12)
 262 
							edge_alpha_energy(va,vb,ab_12)
Line 355... Line 355...
355 
				VertexIter va = he->vert;
 355 
				VertexIter va = he->vert;
356 
				VertexIter vb = he->opp->vert;
 356 
				VertexIter vb = he->opp->vert;
357 
				VertexIter vc = he->next->vert;
 357 
				VertexIter vc = he->next->vert;
358 
				VertexIter vd = he->opp->next->vert;
 358 
				VertexIter vd = he->opp->next->vert;
359 
 
 359 
 
360 
				Vec3d va_pos(va->get_pos());
 360 
				Vec3d va_pos(va->pos);
361 
				Vec3d vb_pos(vb->get_pos());
 361 
				Vec3d vb_pos(vb->pos);
362 
				Vec3d vc_pos(vc->get_pos());
 362 
				Vec3d vc_pos(vc->pos);
363 
				Vec3d vd_pos(vd->get_pos());
 363 
				Vec3d vd_pos(vd->pos);
364 
 
 364 
 
365 
				vector<Vec3d> va_ring_bef;
 365 
				vector<Vec3d> va_ring_bef;
366 
				vector<Vec3d> va_ring_aft;
 366 
				vector<Vec3d> va_ring_aft;
367 
				vector<Vec3d> vb_ring_bef;
 367 
				vector<Vec3d> vb_ring_bef;
368 
				vector<Vec3d> vb_ring_aft;
 368 
				vector<Vec3d> vb_ring_aft;
Line 371... Line 371...
371 
				vector<Vec3d> vd_ring_bef;
 371 
				vector<Vec3d> vd_ring_bef;
372 
				vector<Vec3d> vd_ring_aft;
 372 
				vector<Vec3d> vd_ring_aft;
373 
 
 373 
 
374 
				for(VertexCirculator c(va); !c.end(); ++c)
 374 
				for(VertexCirculator c(va); !c.end(); ++c)
375 
					{
 375 
					{
376 
						va_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
 376 
						va_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
377 
						if(c.get_vertex() != vb)
 377 
						if(c.get_vertex() != vb)
378 
							va_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
 378 
							va_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
379 
					}
 379 
					}
380 
				for(VertexCirculator c(vb); !c.end(); ++c)
 380 
				for(VertexCirculator c(vb); !c.end(); ++c)
381 
					{
 381 
					{
382 
						vb_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
 382 
						vb_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
383 
						if(c.get_vertex() != va)
 383 
						if(c.get_vertex() != va)
384 
							vb_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
 384 
							vb_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
385 
					}
 385 
					}
386 
				for(VertexCirculator c(vc); !c.end(); ++c)
 386 
				for(VertexCirculator c(vc); !c.end(); ++c)
387 
					{
 387 
					{
388 
						vc_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
 388 
						vc_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
389 
						vc_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
 389 
						vc_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
390 
						if(c.get_vertex() == va)
 390 
						if(c.get_vertex() == va)
391 
							vc_ring_aft.push_back(Vec3d(vd->get_pos()));
 391 
							vc_ring_aft.push_back(Vec3d(vd->pos));
392 
					}
 392 
					}
393 
				for(VertexCirculator c(vd); !c.end(); ++c)
 393 
				for(VertexCirculator c(vd); !c.end(); ++c)
394 
					{
 394 
					{
395 
						vd_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
 395 
						vd_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
396 
						vd_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
 396 
						vd_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
397 
						if(c.get_vertex() == vb)
 397 
						if(c.get_vertex() == vb)
398 
							vd_ring_aft.push_back(Vec3d(vc->get_pos()));
 398 
							vd_ring_aft.push_back(Vec3d(vc->pos));
399 
					}
 399 
					}
400 
				double before =
 400 
				double before =
401 
					abs_mean_curv(va_pos,va_ring_bef) +
 401 
					abs_mean_curv(va_pos,va_ring_bef) +
402 
					abs_mean_curv(vb_pos,vb_ring_bef) +
 402 
					abs_mean_curv(vb_pos,vb_ring_bef) +
403 
					abs_mean_curv(vc_pos,vc_ring_bef) +
 403 
					abs_mean_curv(vc_pos,vc_ring_bef) +
Line 439... Line 439...
439 
				VertexIter va = he->vert;
 439 
				VertexIter va = he->vert;
440 
				VertexIter vb = he->opp->vert;
 440 
				VertexIter vb = he->opp->vert;
441 
				VertexIter vc = he->next->vert;
 441 
				VertexIter vc = he->next->vert;
442 
				VertexIter vd = he->opp->next->vert;
 442 
				VertexIter vd = he->opp->next->vert;
443 
 
 443 
 
444 
				Vec3d va_pos(va->get_pos());
 444 
				Vec3d va_pos(va->pos);
445 
				Vec3d vb_pos(vb->get_pos());
 445 
				Vec3d vb_pos(vb->pos);
446 
				Vec3d vc_pos(vc->get_pos());
 446 
				Vec3d vc_pos(vc->pos);
447 
				Vec3d vd_pos(vd->get_pos());
 447 
				Vec3d vd_pos(vd->pos);
448 
 
 448 
 
449 
				vector<Vec3d> va_ring_bef;
 449 
				vector<Vec3d> va_ring_bef;
450 
				vector<Vec3d> va_ring_aft;
 450 
				vector<Vec3d> va_ring_aft;
451 
				vector<Vec3d> vb_ring_bef;
 451 
				vector<Vec3d> vb_ring_bef;
452 
				vector<Vec3d> vb_ring_aft;
 452 
				vector<Vec3d> vb_ring_aft;
Line 455... Line 455...
455 
				vector<Vec3d> vd_ring_bef;
 455 
				vector<Vec3d> vd_ring_bef;
456 
				vector<Vec3d> vd_ring_aft;
 456 
				vector<Vec3d> vd_ring_aft;
457 
 
 457 
 
458 
				for(VertexCirculator c(va); !c.end(); ++c)
 458 
				for(VertexCirculator c(va); !c.end(); ++c)
459 
					{
 459 
					{
460 
						va_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
 460 
						va_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
461 
						if(c.get_vertex() != vb)
 461 
						if(c.get_vertex() != vb)
462 
							va_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
 462 
							va_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
463 
					}
 463 
					}
464 
				for(VertexCirculator c(vb); !c.end(); ++c)
 464 
				for(VertexCirculator c(vb); !c.end(); ++c)
465 
					{
 465 
					{
466 
						vb_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
 466 
						vb_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
467 
						if(c.get_vertex() != va)
 467 
						if(c.get_vertex() != va)
468 
							vb_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
 468 
							vb_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
469 
					}
 469 
					}
470 
				for(VertexCirculator c(vc); !c.end(); ++c)
 470 
				for(VertexCirculator c(vc); !c.end(); ++c)
471 
					{
 471 
					{
472 
						vc_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
 472 
						vc_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
473 
						vc_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
 473 
						vc_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
474 
						if(c.get_vertex() == va)
 474 
						if(c.get_vertex() == va)
475 
							vc_ring_aft.push_back(Vec3d(vd->get_pos()));
 475 
							vc_ring_aft.push_back(Vec3d(vd->pos));
476 
					}
 476 
					}
477 
				for(VertexCirculator c(vd); !c.end(); ++c)
 477 
				for(VertexCirculator c(vd); !c.end(); ++c)
478 
					{
 478 
					{
479 
						vd_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
 479 
						vd_ring_bef.push_back(Vec3d(c.get_vertex()->pos));
480 
						vd_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
 480 
						vd_ring_aft.push_back(Vec3d(c.get_vertex()->pos));
481 
						if(c.get_vertex() == vb)
 481 
						if(c.get_vertex() == vb)
482 
							vd_ring_aft.push_back(Vec3d(vc->get_pos()));
 482 
							vd_ring_aft.push_back(Vec3d(vc->pos));
483 
					}
 483 
					}
484 
				double before =
 484 
				double before =
485 
					gauss_curv(va_pos,va_ring_bef) +
 485 
					gauss_curv(va_pos,va_ring_bef) +
486 
					gauss_curv(vb_pos,vb_ring_bef) +
 486 
					gauss_curv(vb_pos,vb_ring_bef) +
487 
					gauss_curv(vc_pos,vc_ring_bef) +
 487 
					gauss_curv(vc_pos,vc_ring_bef) +
Line 516... Line 516...
516 
 
 516 
 
517 
			void add_to_queue(priority_queue<PQElem>& Q,
 517 
			void add_to_queue(priority_queue<PQElem>& Q,
518 
												HalfEdgeIter he,
 518 
												HalfEdgeIter he,
519 
												const EnergyFun& efun)
 519 
												const EnergyFun& efun)
520 
			{
 520 
			{
521 
				if(!he->is_boundary())
 521 
				if(!is_boundary(he))
522 
					{
 522 
					{
523 
						double energy = efun.delta_energy(he);
 523 
						double energy = efun.delta_energy(he);
524 
						int t = he->touched + 1;
 524 
						int t = he->touched + 1;
525 
						he->touched = t;
 525 
						he->touched = t;
526 
						he->opp->touched = t;
 526 
						he->opp->touched = t;
Line 585... Line 585...
585 
		double tsum=0;
 585 
		double tsum=0;
586 
		{
 586 
		{
587 
			HalfEdgeIter he = m.halfedges_begin();
 587 
			HalfEdgeIter he = m.halfedges_begin();
588 
			for(;he != m.halfedges_end(); ++he)
 588 
			for(;he != m.halfedges_end(); ++he)
589 
				{
 589 
				{
590 
					if(!he->is_boundary())
 590 
					if(!is_boundary(he))
591 
						{
 591 
						{
592 
							double e = efun.delta_energy(he);
 592 
							double e = efun.delta_energy(he);
593 
							tmin = s_min(e, tmin);
 593 
							tmin = s_min(e, tmin);
594 
							tsum += efun.energy(he);
 594 
							tsum += efun.energy(he);
595 
						}
 595 
						}
Line 602... Line 602...
602 
			{
 602 
			{
603 
				vector<HalfEdgeIter>  halfedges;
 603 
				vector<HalfEdgeIter>  halfedges;
604 
				HalfEdgeIter he = m.halfedges_begin();
 604 
				HalfEdgeIter he = m.halfedges_begin();
605 
				for(;he != m.halfedges_end(); ++he)
 605 
				for(;he != m.halfedges_end(); ++he)
606 
					{
 606 
					{
607 
						if(!he->is_boundary())
 607 
						if(!is_boundary(he))
608 
							halfedges.push_back(he);
 608 
							halfedges.push_back(he);
609 
					}
 609 
					}
610 
				random_shuffle(halfedges.begin(), halfedges.end());
 610 
				random_shuffle(halfedges.begin(), halfedges.end());
611 
				swaps = 0;
 611 
				swaps = 0;
612 
				for(int i=0;i<halfedges.size();++i)
 612 
				for(int i=0;i<halfedges.size();++i)