Subversion Repositories gelsvn

Rev

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

Rev 39 Rev 59
Line 12... Line 12...
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;
Line 24... Line 24...
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 
		}
Line 74... Line 74...
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:
Line 306... Line 300...
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 
 
Line 518... Line 512...
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,