WebSVN – gelsvn – Diff – /trunk/src/HMesh/mesh_optimization.cpp


/* ----------------------------------------------------------------------- *
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
 * Copyright (C) the authors and DTU Informatics
 * For license and list of authors, see ../../doc/intro.pdf
 * ----------------------------------------------------------------------- */
 
 
#include "mesh_optimization.h"
 
#include <cfloat>
#include <queue>
#include <vector>
 
#include "../CGLA/Vec3d.h"
 
//#include "Manifold.h"
#include "AttributeVector.h"
 
namespace HMesh
{
    using namespace std;
    using namespace CGLA;
	
    // Small utility functions
    namespace 
    {
        class LineSeg
        {
            const Vec3d p0;
            const Vec3d p1;
            const float l;
            const CGLA::Vec3d dir;
        public:
			
            LineSeg(const Vec3d& _p0, const Vec3d& _p1): 
			p0(_p0), p1(_p1), l(length(p1-p0)), dir((p1-p0)/l) {}
			
			bool inseg(const Vec3d& p) const
			{
				double t = dot(dir, p-p0);
				if(t<0)
					return false;
				if(t>l)
					return false;
				return true;
			}
        };
		
        Vec3d compute_normal(Vec3d* v)
        {
            Vec3d norm;
            for(int i = 0; i < 4; ++i)
            {
                norm[0] += (v[i][1]-v[(i+1)%4][1])*(v[i][2]+v[(i+1)%4][2]);
                norm[1] += (v[i][2]-v[(i+1)%4][2])*(v[i][0]+v[(i+1)%4][0]);
                norm[2] += (v[i][0]-v[(i+1)%4][0])*(v[i][1]+v[(i+1)%4][1]);
            }
            float l = norm.length();
            if(l>0.0f)
                norm /= l;
            return norm;
        }
		
        bool would_flip(const Manifold& m, HalfEdgeID h)
        {
            Walker w = m.walker(h);
			
            VertexID hv = w.vertex();
            VertexID hov = w.opp().vertex();
            VertexID hnv = w.next().vertex();
            VertexID honv = w.opp().next().vertex();
			
            Vec3d v[4];
            v[0] = Vec3d(m.pos(hv));
            v[1] = Vec3d(m.pos(hov));
            v[2] = Vec3d(m.pos(hnv));
            v[3] = Vec3d(m.pos(honv));
			
            Vec3d dir = compute_normal(v);
			
            Vec3d n1a = cross(v[3]-v[0], v[2]-v[0]);
            Vec3d n2a = cross(v[2]-v[1], v[3]-v[1]);
			
            if(dot(normalize(n1a), dir) < 0)
                return true;
            if(dot(normalize(n2a), dir) < 0)
                return true;
            return false;
        }
    }
	
	
    double ValencyEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const
    {
        Walker w = m.walker(h);
        
        VertexID v1 = w.opp().vertex();
        VertexID v2 = w.vertex();
        VertexID vo1 = w.next().vertex();
        VertexID vo2 = w.opp().next().vertex();
        
        int val1  = valency(m, v1);
        int val2  = valency(m, v2);
        int valo1 = valency(m, vo1);
        int valo2 = valency(m, vo2);
        
        // The optimal valency is four for a boundary vertex
        // and six elsewhere.
        int t1 = boundary(m, v1) ? 4 : 6;
        int t2 = boundary(m, v2) ? 4 : 6;
        int to1 = boundary(m, vo1) ? 4 : 6;
        int to2 = boundary(m, vo2) ? 4 : 6;
        
        int before = 
        sqr(val1-t1)+sqr(val2-t2)+
        sqr(valo1-to1)+sqr(valo2-to2);
        int after = 
        sqr(valo1+1-to1)+sqr(val1-1-t1)+
        sqr(val2-1-t2)+sqr(valo2+1-to2);
        
        return static_cast<double>(after-before);
    }
	
	
	class RandomEnergy: public EnergyFun
	{
	public:
		double delta_energy(const Manifold& m, HalfEdgeID he) const
		{
			return static_cast<double>(gel_rand()/static_cast<float>(GEL_RAND_MAX));
		}
	};
	
	
	double MinAngleEnergy::min_angle(const Vec3d& v0, const Vec3d& v1, const Vec3d& v2) const
	{
		Vec3d a = normalize(v1-v0);
		Vec3d b = normalize(v2-v1);
		Vec3d c = normalize(v0-v2);
		
		return min(dot(a,-c), min(dot(b,-a), dot(c,-b)));
	}
	double MinAngleEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const
	{
		Walker w = m.walker(h);
		
		VertexID hv = w.vertex();
		VertexID hnv = w.next().vertex();
		VertexID hov= w.opp().vertex();
		VertexID honv = w.opp().next().vertex();
		
		Vec3d v0(m.pos(hv));
		Vec3d v1(m.pos(hnv));
		Vec3d v2(m.pos(hov));
		Vec3d v3(m.pos(honv));
		
		Vec3d n1a = normalize(cross(v1-v0,v3-v0));
		Vec3d n2a = normalize(cross(v3-v2,v1-v2));
		
		if(dot(n1a, n2a) > thresh){
			double before = min(min_angle(v0,v1,v2), min_angle(v0,v2,v3));
			double after = min(min_angle(v0,v1,v3), min_angle(v1,v2,v3));
			return -(after-before);
		}
		return DBL_MAX;
	}
	
	
	void DihedralEnergy::compute_angles(const Manifold & m, HalfEdgeID h) const
	{
		Walker w = m.walker(h);
		
		VertexID hv = w.vertex();
		VertexID hov = w.opp().vertex();
		VertexID hnv = w.next().vertex();
		VertexID honv = w.opp().next().vertex();
		
		Vec3d va(m.pos(hv));
		Vec3d vb(m.pos(hov));
		Vec3d vc(m.pos(hnv));
		Vec3d vd(m.pos(honv));
		
		FaceID fa = w.next().opp().face();
		FaceID fb = w.next().next().opp().face();
		FaceID fc = w.opp().next().opp().face();
		FaceID fd = w.opp().next().next().opp().face();
		
		Vec3d n1 = normalize(cross(vc-va, vb-va));
		Vec3d n2 = normalize(cross(vb-va, vd-va));
		
		Vec3d na = fa == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fa));
		Vec3d nb = fb == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fb));
		Vec3d nc = fc == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fc));
		Vec3d nd = fd == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fd));
		
		
		Vec3d fn1 = normalize(cross(vb-vc, vd-vc));
		Vec3d fn2 = normalize(cross(vd-vc, va-vc));
		
		ab_12 = cos_ang(n1,n2);
		ab_a1 = cos_ang(na,n1);
		ab_b1 = cos_ang(nb,n1);
		ab_2c = cos_ang(n2,nc);
		ab_2d = cos_ang(n2,nd);
		
		aa_12 = cos_ang(fn1,fn2);
		aa_b1 = cos_ang(nb,fn1);
		aa_c1 = cos_ang(nc, fn1);
		aa_2a = cos_ang(fn2, na);
		aa_2d = cos_ang(fn2,nd);
	}
	
	double DihedralEnergy::energy(const Manifold& m, HalfEdgeID h) const
	{
		Walker w = m.walker(h);
		
		FaceID hf = w.face();
		FaceID hof = w.opp().face();
		
		double a = cos_ang(Vec3d(normal(m, hf)), Vec3d(normal(m, hof)));
		
		VertexID hv = w.vertex();
		VertexID hov = w.opp().vertex();
		
		Vec3d va(m.pos(hv));
		Vec3d vb(m.pos(hov));
		
		if(use_alpha)
			return edge_alpha_energy(va,vb,a);
		
		return edge_c_energy(va,vb,a);
	}
	
	
	double DihedralEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const
	{
		compute_angles(m, h);
		
		Walker w = m.walker(h);
		
		VertexID hv = w.vertex();
		VertexID hov = w.opp().vertex();
		VertexID hnv = w.next().vertex();
		VertexID honv = w.opp().next().vertex();
		
		Vec3d va(m.pos(hv));
		Vec3d vb(m.pos(hov));
		Vec3d vc(m.pos(hnv));
		Vec3d vd(m.pos(honv));
		
		if(use_alpha){
			double before = 
			edge_alpha_energy(va,vb,ab_12)
			+edge_alpha_energy(va,vc,ab_a1)
			+edge_alpha_energy(vc,vb,ab_b1)
			+edge_alpha_energy(vd,vb,ab_2c)
			+edge_alpha_energy(vd,va,ab_2d);
			
			double after = 
			edge_alpha_energy(vd,vc,aa_12)
			+edge_alpha_energy(vb,vc,aa_b1)
			+edge_alpha_energy(vd,vb,aa_c1)
			+edge_alpha_energy(va,vc,aa_2a)
			+edge_alpha_energy(vd,va,aa_2d);
			
			return (after-before);
		}
		double before = 
		edge_c_energy(va,vb,ab_12)
		+edge_c_energy(va,vc,ab_a1)
		+edge_c_energy(vc,vb,ab_b1)
		+edge_c_energy(vd,vb,ab_2c)
		+edge_c_energy(vd,va,ab_2d);
		
		double after = 
		edge_c_energy(vd,vc,aa_12)
		+edge_c_energy(vb,vc,aa_b1)
		+edge_c_energy(vd,vb,aa_c1)
		+edge_c_energy(va,vc,aa_2a)
		+edge_c_energy(vd,va,aa_2d);
		
		return after-before;
	}
	
	
	
	double CurvatureEnergy::abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const
	{
		const size_t N = ring.size();
		
		double H = 0;
		for(int i = 0; i < N; ++i){
			Vec3d vnim1 = ring[(i+N-1)%N] - v;
			Vec3d vni   = ring[i] - v;
			Vec3d vnip1 = ring[(i+1)%N] - v;
			
			Vec3d Nm = normalize(cross(vni, vnim1));
			Vec3d Np = normalize(cross(vnip1, vni));
			
			double beta = acos(max(-1.0, min(1.0, dot(Nm, Np))));
			H += vni.length() * beta;
		}
		H /= 4;
		
		return H;
	}
	
	double CurvatureEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const
	{
		Walker w = m.walker(h);
		
		VertexID va = w.vertex();
		VertexID vb = w.opp().vertex();
		VertexID vc = w.next().vertex();
		VertexID vd = w.opp().next().vertex();
        
        
		
		Vec3d va_pos(m.pos(va));
		Vec3d vb_pos(m.pos(vb));
		Vec3d vc_pos(m.pos(vc));
		Vec3d vd_pos(m.pos(vd));
		
		vector<Vec3d> va_ring_bef;
		vector<Vec3d> va_ring_aft;
		vector<Vec3d> vb_ring_bef;
		vector<Vec3d> vb_ring_aft;
		vector<Vec3d> vc_ring_bef;
		vector<Vec3d> vc_ring_aft;
		vector<Vec3d> vd_ring_bef;
		vector<Vec3d> vd_ring_aft;
		
		for(Walker wv = m.walker(va); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
			VertexID v = wv.vertex();
			Vec3d pos(m.pos(v));
			
			va_ring_bef.push_back(pos);
			if(v != vb)
				va_ring_aft.push_back(pos);
		}
		for(Walker wv = m.walker(vb); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
			VertexID v = wv.vertex();
			Vec3d pos(m.pos(v));
			
			vb_ring_bef.push_back(pos);
			if(v != va)
				vb_ring_aft.push_back(pos);
		}
		for(Walker wv = m.walker(vc); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
			VertexID v = wv.vertex();
			Vec3d pos(m.pos(v));
			
			vc_ring_bef.push_back(pos);
			vc_ring_aft.push_back(pos);
			if(v == va)
				vc_ring_aft.push_back(vd_pos);
		}
		for(Walker wv = m.walker(vd); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
			VertexID v = wv.vertex();
			Vec3d pos(m.pos(v));
			
			vd_ring_bef.push_back(pos);
			vd_ring_aft.push_back(pos);
			if(v == vb)
				vd_ring_aft.push_back(vc_pos);
		}
		double before =
		abs_mean_curv(va_pos, va_ring_bef) +
		abs_mean_curv(vb_pos, vb_ring_bef) +
		abs_mean_curv(vc_pos, vc_ring_bef) +
		abs_mean_curv(vd_pos, vd_ring_bef);
		
		double after =
		abs_mean_curv(va_pos, va_ring_aft) +
		abs_mean_curv(vb_pos, vb_ring_aft) +
		abs_mean_curv(vc_pos, vc_ring_aft) +
		abs_mean_curv(vd_pos, vd_ring_aft);
		
		return after-before;
	}
	
	
	class GaussCurvatureEnergy: public EnergyFun
	{
		
		double gauss_curv(const Vec3d& v, const vector<Vec3d>& ring) const
		{
			const size_t N = ring.size();
			double asum=0.0f;
			double area_sum=0;
			for(int i = 0; i < N; ++i){
				const Vec3d& v1 = ring[i];
				const Vec3d& v2 = ring[(i+1)%N];
				Vec3d a = v1-v;
				Vec3d b = v2-v;
				asum += acos(max(-1.0, min(1.0, dot(a,b)/(length(a)*length(b)))));
				area_sum += 0.5 * length(cross(a,b));
			}
			return 3*abs(2 * M_PI - asum)/area_sum;
		}
		
	public:
		double delta_energy(const Manifold& m, HalfEdgeID h) const
		{
			Walker w = m.walker(h);
			
			VertexID va = w.vertex();
			VertexID vb = w.opp().vertex();
			VertexID vc = w.next().vertex();
			VertexID vd = w.opp().next().vertex();
			
			Vec3d va_pos(m.pos(va));
			Vec3d vb_pos(m.pos(vb));
			Vec3d vc_pos(m.pos(vc));
			Vec3d vd_pos(m.pos(vd));
			
			vector<Vec3d> va_ring_bef;
			vector<Vec3d> va_ring_aft;
			vector<Vec3d> vb_ring_bef;
			vector<Vec3d> vb_ring_aft;
			vector<Vec3d> vc_ring_bef;
			vector<Vec3d> vc_ring_aft;
			vector<Vec3d> vd_ring_bef;
			vector<Vec3d> vd_ring_aft;
			
			for(Walker wv = m.walker(va); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
				VertexID v = wv.vertex();
				Vec3d pos(m.pos(v));
				
				va_ring_bef.push_back(pos);
				if(v != vb)
					va_ring_aft.push_back(pos);
			}
			for(Walker wv = m.walker(vb); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
				VertexID v = wv.vertex();
				Vec3d pos(m.pos(v));
				
				vb_ring_bef.push_back(pos);
				if(v != va)
					vb_ring_aft.push_back(pos);
			}
			for(Walker wv = m.walker(vc); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
				VertexID v = wv.vertex();
				Vec3d pos(m.pos(v));
				
				vc_ring_bef.push_back(pos);
				vc_ring_aft.push_back(pos);
				if(v == va)
					vc_ring_aft.push_back(pos);
			}
			for(Walker wv = m.walker(vd); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
				VertexID v = wv.vertex();
				Vec3d pos(m.pos(v));
				
				vd_ring_bef.push_back(pos);
				vd_ring_aft.push_back(pos);
				if(v == vb)
					vd_ring_aft.push_back(pos);
			}
			double before =
			gauss_curv(va_pos,va_ring_bef) +
			gauss_curv(vb_pos,vb_ring_bef) +
			gauss_curv(vc_pos,vc_ring_bef) +
			gauss_curv(vd_pos,vd_ring_bef);
			
			double after =
			gauss_curv(va_pos,va_ring_aft) +
			gauss_curv(vb_pos,vb_ring_aft) +
			gauss_curv(vc_pos,vc_ring_aft) +
			gauss_curv(vd_pos,vd_ring_aft);
			
			return after-before;
		}
	};
	
	struct PQElem
	{
		double pri;
		HalfEdgeID h;
		int time;
		
		//PQElem() {}
		PQElem(double _pri, HalfEdgeID _h, int _time): 
		pri(_pri), h(_h), time(_time) {}
	};
	
	bool operator<(const PQElem& e0, const PQElem& e1)
	{
		return e0.pri > e1.pri;
	}
	
	
	void add_to_queue(const Manifold& m, HalfEdgeAttributeVector<int>& touched, priority_queue<PQElem>& Q, HalfEdgeID h, const EnergyFun& efun)
	{
		if(boundary(m, h))
			return;
		
		Walker w = m.walker(h);
		HalfEdgeID ho = w.opp().halfedge();
		
		double energy = efun.delta_energy(m, h);
		int t = touched[h] + 1;
		touched[h] = t;
		touched[ho] = t;
		if((energy<0) && (t < 10000)){
			Q.push(PQElem(energy, h, t));
		}
		
	}
	
	void add_one_ring_to_queue(const Manifold& m, HalfEdgeAttributeVector<int>& touched, priority_queue<PQElem>& Q, VertexID v, const EnergyFun& efun)
	{
		
		for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
			add_to_queue(m, touched, Q, w.halfedge(), efun);
		}
	}
	
	
	void priority_queue_optimization(Manifold& m, const EnergyFun& efun)
	{
		HalfEdgeAttributeVector<int> touched(m.allocated_halfedges(), 0);
		priority_queue<PQElem> Q;
		
		cout << "Building priority queue"<< endl;
		
		for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){
			if(!touched[*h])
				add_to_queue(m, touched, Q, *h, efun);
		}
		
		cout << "Emptying priority queue of size: " << Q.size() << " ";
		while(!Q.empty())
		{
			if(Q.size() % 1000 == 0)
				cout << ".";
			if(Q.size() % 10000 == 0)
				cout << Q.size();
			
			PQElem elem = Q.top();
			Q.pop();
			
			if(touched[elem.h] != elem.time)
				continue;
			if(!precond_flip_edge(m, elem.h))
				continue;
			
			m.flip_edge(elem.h);
			
			Walker w = m.walker(elem.h);
			add_one_ring_to_queue(m, touched, Q, w.vertex(), efun);
			add_one_ring_to_queue(m, touched, Q, w.next().vertex(), efun);
			add_one_ring_to_queue(m, touched, Q, w.opp().vertex(), efun);
			add_one_ring_to_queue(m, touched, Q, w.opp().next().vertex(), efun);
			
		}
		cout << endl;
	}
	
	void simulated_annealing_optimization(Manifold& m, const EnergyFun& efun, int max_iter)
	{
		gel_srand(0);
		int swaps;
		int iter = 0;
		double T = 1;
		
		double tmin=0;
		{
			for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){
				if(boundary(m, *h)) 
					continue;  
				double e = efun.delta_energy(m, *h);
				tmin = s_min(e, tmin);
			}
		}
		if (tmin < 0.0) 
			T = -2*tmin;
		
		if(max_iter>0){ 
			do{
				cout << "Temperature : " << T << endl;
				vector<HalfEdgeID>  halfedges;
				for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){
					if(boundary(m, *h))
						continue;
					halfedges.push_back(*h);
				}
				random_shuffle(halfedges.begin(), halfedges.end());
				swaps = 0;
				for(size_t i = 0; i < halfedges.size(); ++i){
					HalfEdgeID h = halfedges[i];
					DihedralEnergy dih_en;
					double dma = dih_en.min_angle(m, h);
					if(dma > -0.4){
						double delta = efun.delta_energy(m, h);
						if(delta < -1e-8){
							if(precond_flip_edge(m, h)){
								m.flip_edge(h);
								++swaps;
							}
						}
						else{
							delta = max(1e-8, delta);
							double prob = min(0.9999, exp(-delta/T));
							if(gel_rand()/double(GEL_RAND_MAX) < prob){
								if(precond_flip_edge(m, h)){
									m.flip_edge(h);
									++swaps;
								}
							}
						}
					}
				}
				cout << "Swaps = " << swaps << " T = " << T << endl;
				if(iter % 5 == 0 && iter > 0)
					T *= 0.9;
			}
			while(++iter < max_iter && swaps);
		}
		cout << "Iterations "  << iter << endl; 
		
	}
	
	
	void minimize_dihedral_angle(Manifold& m, 
								 int iter,
								 bool anneal, 
								 bool alpha,
								 double gamma)
	{
		DihedralEnergy energy_fun(gamma, alpha);
		if(anneal)
			simulated_annealing_optimization(m, energy_fun, iter);
		else
			priority_queue_optimization(m, energy_fun);
	}
	
	void randomize_mesh(Manifold& m, int max_iter)
	{
		RandomEnergy energy_fun;
		simulated_annealing_optimization(m, energy_fun, max_iter);
	}
	
	void minimize_curvature(Manifold& m, bool anneal)
	{
		CurvatureEnergy energy_fun;
		if(anneal)
			simulated_annealing_optimization(m, energy_fun);
		else
			priority_queue_optimization(m, energy_fun);
	}
	
	void minimize_gauss_curvature(Manifold& m, bool anneal)
	{
		GaussCurvatureEnergy energy_fun;
		if(anneal)
			simulated_annealing_optimization(m, energy_fun);
		else
			priority_queue_optimization(m, energy_fun);
	}
	
	void maximize_min_angle(Manifold& m, float thresh, bool anneal)
	{
		MinAngleEnergy energy_fun(thresh);
		if(anneal)
			simulated_annealing_optimization(m, energy_fun);
		else
			priority_queue_optimization(m, energy_fun);
	}
	
	void optimize_valency(Manifold& m, bool anneal)
	{
		ValencyEnergy energy_fun;
		if(anneal)
			simulated_annealing_optimization(m, energy_fun);
		else
			priority_queue_optimization(m, energy_fun);
	}
}
 

Rev 595	Rev 601
1	`/* ----------------------------------------------------------------------- *`	1	`/* ----------------------------------------------------------------------- *`
2	`* This file is part of GEL, http://www.imm.dtu.dk/GEL`	2	`* This file is part of GEL, http://www.imm.dtu.dk/GEL`
3	`* Copyright (C) the authors and DTU Informatics`	3	`* Copyright (C) the authors and DTU Informatics`
4	`* For license and list of authors, see ../../doc/intro.pdf`	4	`* For license and list of authors, see ../../doc/intro.pdf`
5	`* ----------------------------------------------------------------------- */`	5	`* ----------------------------------------------------------------------- */`
6		6
7		7
8	`#include "mesh_optimization.h"`	8	`#include "mesh_optimization.h"`
9		9
10	`#include <cfloat>`	10	`#include <cfloat>`
11	`#include <queue>`	11	`#include <queue>`
12	`#include <vector>`	12	`#include <vector>`
13		13
14	`#include <CGLA/Vec3d.h>`	14	`#include "../CGLA/Vec3d.h"`
15		15
16	`//#include "Manifold.h"`	16	`//#include "Manifold.h"`
17	`#include "AttributeVector.h"`	17	`#include "AttributeVector.h"`
18		18
19	`namespace HMesh`	19	`namespace HMesh`
20	`{`	20	`{`
21	`using namespace std;`	21	`using namespace std;`
22	`using namespace CGLA;`	22	`using namespace CGLA;`
23		23
24	`// Small utility functions`	24	`// Small utility functions`
25	`namespace`	25	`namespace`
26	`{`	26	`{`
27	`class LineSeg`	27	`class LineSeg`
28	`{`	28	`{`
29	`const Vec3d p0;`	29	`const Vec3d p0;`
30	`const Vec3d p1;`	30	`const Vec3d p1;`
31	`const float l;`	31	`const float l;`
32	`const CGLA::Vec3d dir;`	32	`const CGLA::Vec3d dir;`
33	`public:`	33	`public:`
34		34
35	`LineSeg(const Vec3d& _p0, const Vec3d& _p1):`	35	`LineSeg(const Vec3d& _p0, const Vec3d& _p1):`
36	`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) {}`
37		37
38	`bool inseg(const Vec3d& p) const`	38	`bool inseg(const Vec3d& p) const`
39	`{`	39	`{`
40	`double t = dot(dir, p-p0);`	40	`double t = dot(dir, p-p0);`
41	`if(t<0)`	41	`if(t<0)`
42	`return false;`	42	`return false;`
43	`if(t>l)`	43	`if(t>l)`
44	`return false;`	44	`return false;`
45	`return true;`	45	`return true;`
46	`}`	46	`}`
47	`};`	47	`};`
48		48
49	`Vec3d compute_normal(Vec3d* v)`	49	`Vec3d compute_normal(Vec3d* v)`
50	`{`	50	`{`
51	`Vec3d norm;`	51	`Vec3d norm;`
52	`for(int i = 0; i < 4; ++i)`	52	`for(int i = 0; i < 4; ++i)`
53	`{`	53	`{`
54	`norm[0] += (v[i][1]-v[(i+1)%4][1])*(v[i][2]+v[(i+1)%4][2]);`	54	`norm[0] += (v[i][1]-v[(i+1)%4][1])*(v[i][2]+v[(i+1)%4][2]);`
55	`norm[1] += (v[i][2]-v[(i+1)%4][2])*(v[i][0]+v[(i+1)%4][0]);`	55	`norm[1] += (v[i][2]-v[(i+1)%4][2])*(v[i][0]+v[(i+1)%4][0]);`
56	`norm[2] += (v[i][0]-v[(i+1)%4][0])*(v[i][1]+v[(i+1)%4][1]);`	56	`norm[2] += (v[i][0]-v[(i+1)%4][0])*(v[i][1]+v[(i+1)%4][1]);`
57	`}`	57	`}`
58	`float l = norm.length();`	58	`float l = norm.length();`
59	`if(l>0.0f)`	59	`if(l>0.0f)`
60	`norm /= l;`	60	`norm /= l;`
61	`return norm;`	61	`return norm;`
62	`}`	62	`}`
63		63
64	`bool would_flip(const Manifold& m, HalfEdgeID h)`	64	`bool would_flip(const Manifold& m, HalfEdgeID h)`
65	`{`	65	`{`
66	`Walker w = m.walker(h);`	66	`Walker w = m.walker(h);`
67		67
68	`VertexID hv = w.vertex();`	68	`VertexID hv = w.vertex();`
69	`VertexID hov = w.opp().vertex();`	69	`VertexID hov = w.opp().vertex();`
70	`VertexID hnv = w.next().vertex();`	70	`VertexID hnv = w.next().vertex();`
71	`VertexID honv = w.opp().next().vertex();`	71	`VertexID honv = w.opp().next().vertex();`
72		72
73	`Vec3d v[4];`	73	`Vec3d v[4];`
74	`v[0] = Vec3d(m.pos(hv));`	74	`v[0] = Vec3d(m.pos(hv));`
75	`v[1] = Vec3d(m.pos(hov));`	75	`v[1] = Vec3d(m.pos(hov));`
76	`v[2] = Vec3d(m.pos(hnv));`	76	`v[2] = Vec3d(m.pos(hnv));`
77	`v[3] = Vec3d(m.pos(honv));`	77	`v[3] = Vec3d(m.pos(honv));`
78		78
79	`Vec3d dir = compute_normal(v);`	79	`Vec3d dir = compute_normal(v);`
80		80
81	`Vec3d n1a = cross(v[3]-v[0], v[2]-v[0]);`	81	`Vec3d n1a = cross(v[3]-v[0], v[2]-v[0]);`
82	`Vec3d n2a = cross(v[2]-v[1], v[3]-v[1]);`	82	`Vec3d n2a = cross(v[2]-v[1], v[3]-v[1]);`
83		83
84	`if(dot(normalize(n1a), dir) < 0)`	84	`if(dot(normalize(n1a), dir) < 0)`
85	`return true;`	85	`return true;`
86	`if(dot(normalize(n2a), dir) < 0)`	86	`if(dot(normalize(n2a), dir) < 0)`
87	`return true;`	87	`return true;`
88	`return false;`	88	`return false;`
89	`}`	89	`}`
90	`}`	90	`}`
91		91
92		92
93	`double ValencyEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const`	93	`double ValencyEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const`
94	`{`	94	`{`
95	`Walker w = m.walker(h);`	95	`Walker w = m.walker(h);`
96		96
97	`VertexID v1 = w.opp().vertex();`	97	`VertexID v1 = w.opp().vertex();`
98	`VertexID v2 = w.vertex();`	98	`VertexID v2 = w.vertex();`
99	`VertexID vo1 = w.next().vertex();`	99	`VertexID vo1 = w.next().vertex();`
100	`VertexID vo2 = w.opp().next().vertex();`	100	`VertexID vo2 = w.opp().next().vertex();`
101		101
102	`int val1 = valency(m, v1);`	102	`int val1 = valency(m, v1);`
103	`int val2 = valency(m, v2);`	103	`int val2 = valency(m, v2);`
104	`int valo1 = valency(m, vo1);`	104	`int valo1 = valency(m, vo1);`
105	`int valo2 = valency(m, vo2);`	105	`int valo2 = valency(m, vo2);`
106		106
107	`// The optimal valency is four for a boundary vertex`	107	`// The optimal valency is four for a boundary vertex`
108	`// and six elsewhere.`	108	`// and six elsewhere.`
109	`int t1 = boundary(m, v1) ? 4 : 6;`	109	`int t1 = boundary(m, v1) ? 4 : 6;`
110	`int t2 = boundary(m, v2) ? 4 : 6;`	110	`int t2 = boundary(m, v2) ? 4 : 6;`
111	`int to1 = boundary(m, vo1) ? 4 : 6;`	111	`int to1 = boundary(m, vo1) ? 4 : 6;`
112	`int to2 = boundary(m, vo2) ? 4 : 6;`	112	`int to2 = boundary(m, vo2) ? 4 : 6;`
113		113
114	`int before =`	114	`int before =`
115	`sqr(val1-t1)+sqr(val2-t2)+`	115	`sqr(val1-t1)+sqr(val2-t2)+`
116	`sqr(valo1-to1)+sqr(valo2-to2);`	116	`sqr(valo1-to1)+sqr(valo2-to2);`
117	`int after =`	117	`int after =`
118	`sqr(valo1+1-to1)+sqr(val1-1-t1)+`	118	`sqr(valo1+1-to1)+sqr(val1-1-t1)+`
119	`sqr(val2-1-t2)+sqr(valo2+1-to2);`	119	`sqr(val2-1-t2)+sqr(valo2+1-to2);`
120		120
121	`return static_cast<double>(after-before);`	121	`return static_cast<double>(after-before);`
122	`}`	122	`}`
123		123
124		124
125	`class RandomEnergy: public EnergyFun`	125	`class RandomEnergy: public EnergyFun`
126	`{`	126	`{`
127	`public:`	127	`public:`
128	`double delta_energy(const Manifold& m, HalfEdgeID he) const`	128	`double delta_energy(const Manifold& m, HalfEdgeID he) const`
129	`{`	129	`{`
130	`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));`
131	`}`	131	`}`
132	`};`	132	`};`
133		133
134		134
135	`double MinAngleEnergy::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`
136	`{`	136	`{`
137	`Vec3d a = normalize(v1-v0);`	137	`Vec3d a = normalize(v1-v0);`
138	`Vec3d b = normalize(v2-v1);`	138	`Vec3d b = normalize(v2-v1);`
139	`Vec3d c = normalize(v0-v2);`	139	`Vec3d c = normalize(v0-v2);`
140		140
141	`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)));`
142	`}`	142	`}`
143	`double MinAngleEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const`	143	`double MinAngleEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const`
144	`{`	144	`{`
145	`Walker w = m.walker(h);`	145	`Walker w = m.walker(h);`
146		146
147	`VertexID hv = w.vertex();`	147	`VertexID hv = w.vertex();`
148	`VertexID hnv = w.next().vertex();`	148	`VertexID hnv = w.next().vertex();`
149	`VertexID hov= w.opp().vertex();`	149	`VertexID hov= w.opp().vertex();`
150	`VertexID honv = w.opp().next().vertex();`	150	`VertexID honv = w.opp().next().vertex();`
151		151
152	`Vec3d v0(m.pos(hv));`	152	`Vec3d v0(m.pos(hv));`
153	`Vec3d v1(m.pos(hnv));`	153	`Vec3d v1(m.pos(hnv));`
154	`Vec3d v2(m.pos(hov));`	154	`Vec3d v2(m.pos(hov));`
155	`Vec3d v3(m.pos(honv));`	155	`Vec3d v3(m.pos(honv));`
156		156
157	`Vec3d n1a = normalize(cross(v1-v0,v3-v0));`	157	`Vec3d n1a = normalize(cross(v1-v0,v3-v0));`
158	`Vec3d n2a = normalize(cross(v3-v2,v1-v2));`	158	`Vec3d n2a = normalize(cross(v3-v2,v1-v2));`
159		159
160	`if(dot(n1a, n2a) > thresh){`	160	`if(dot(n1a, n2a) > thresh){`
161	`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));`
162	`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));`
163	`return -(after-before);`	163	`return -(after-before);`
164	`}`	164	`}`
165	`return DBL_MAX;`	165	`return DBL_MAX;`
166	`}`	166	`}`
167		167
168		168
169	`void DihedralEnergy::compute_angles(const Manifold & m, HalfEdgeID h) const`	169	`void DihedralEnergy::compute_angles(const Manifold & m, HalfEdgeID h) const`
170	`{`	170	`{`
171	`Walker w = m.walker(h);`	171	`Walker w = m.walker(h);`
172		172
173	`VertexID hv = w.vertex();`	173	`VertexID hv = w.vertex();`
174	`VertexID hov = w.opp().vertex();`	174	`VertexID hov = w.opp().vertex();`
175	`VertexID hnv = w.next().vertex();`	175	`VertexID hnv = w.next().vertex();`
176	`VertexID honv = w.opp().next().vertex();`	176	`VertexID honv = w.opp().next().vertex();`
177		177
178	`Vec3d va(m.pos(hv));`	178	`Vec3d va(m.pos(hv));`
179	`Vec3d vb(m.pos(hov));`	179	`Vec3d vb(m.pos(hov));`
180	`Vec3d vc(m.pos(hnv));`	180	`Vec3d vc(m.pos(hnv));`
181	`Vec3d vd(m.pos(honv));`	181	`Vec3d vd(m.pos(honv));`
182		182
183	`FaceID fa = w.next().opp().face();`	183	`FaceID fa = w.next().opp().face();`
184	`FaceID fb = w.next().next().opp().face();`	184	`FaceID fb = w.next().next().opp().face();`
185	`FaceID fc = w.opp().next().opp().face();`	185	`FaceID fc = w.opp().next().opp().face();`
186	`FaceID fd = w.opp().next().next().opp().face();`	186	`FaceID fd = w.opp().next().next().opp().face();`
187		187
188	`Vec3d n1 = normalize(cross(vc-va, vb-va));`	188	`Vec3d n1 = normalize(cross(vc-va, vb-va));`
189	`Vec3d n2 = normalize(cross(vb-va, vd-va));`	189	`Vec3d n2 = normalize(cross(vb-va, vd-va));`
190		190
191	`Vec3d na = fa == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fa));`	191	`Vec3d na = fa == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fa));`
192	`Vec3d nb = fb == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fb));`	192	`Vec3d nb = fb == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fb));`
193	`Vec3d nc = fc == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fc));`	193	`Vec3d nc = fc == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fc));`
194	`Vec3d nd = fd == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fd));`	194	`Vec3d nd = fd == InvalidFaceID ? Vec3d(0) : Vec3d(normal(m, fd));`
195		195
196		196
197	`Vec3d fn1 = normalize(cross(vb-vc, vd-vc));`	197	`Vec3d fn1 = normalize(cross(vb-vc, vd-vc));`
198	`Vec3d fn2 = normalize(cross(vd-vc, va-vc));`	198	`Vec3d fn2 = normalize(cross(vd-vc, va-vc));`
199		199
200	`ab_12 = cos_ang(n1,n2);`	200	`ab_12 = cos_ang(n1,n2);`
201	`ab_a1 = cos_ang(na,n1);`	201	`ab_a1 = cos_ang(na,n1);`
202	`ab_b1 = cos_ang(nb,n1);`	202	`ab_b1 = cos_ang(nb,n1);`
203	`ab_2c = cos_ang(n2,nc);`	203	`ab_2c = cos_ang(n2,nc);`
204	`ab_2d = cos_ang(n2,nd);`	204	`ab_2d = cos_ang(n2,nd);`
205		205
206	`aa_12 = cos_ang(fn1,fn2);`	206	`aa_12 = cos_ang(fn1,fn2);`
207	`aa_b1 = cos_ang(nb,fn1);`	207	`aa_b1 = cos_ang(nb,fn1);`
208	`aa_c1 = cos_ang(nc, fn1);`	208	`aa_c1 = cos_ang(nc, fn1);`
209	`aa_2a = cos_ang(fn2, na);`	209	`aa_2a = cos_ang(fn2, na);`
210	`aa_2d = cos_ang(fn2,nd);`	210	`aa_2d = cos_ang(fn2,nd);`
211	`}`	211	`}`
212		212
213	`double DihedralEnergy::energy(const Manifold& m, HalfEdgeID h) const`	213	`double DihedralEnergy::energy(const Manifold& m, HalfEdgeID h) const`
214	`{`	214	`{`
215	`Walker w = m.walker(h);`	215	`Walker w = m.walker(h);`
216		216
217	`FaceID hf = w.face();`	217	`FaceID hf = w.face();`
218	`FaceID hof = w.opp().face();`	218	`FaceID hof = w.opp().face();`
219		219
220	`double a = cos_ang(Vec3d(normal(m, hf)), Vec3d(normal(m, hof)));`	220	`double a = cos_ang(Vec3d(normal(m, hf)), Vec3d(normal(m, hof)));`
221		221
222	`VertexID hv = w.vertex();`	222	`VertexID hv = w.vertex();`
223	`VertexID hov = w.opp().vertex();`	223	`VertexID hov = w.opp().vertex();`
224		224
225	`Vec3d va(m.pos(hv));`	225	`Vec3d va(m.pos(hv));`
226	`Vec3d vb(m.pos(hov));`	226	`Vec3d vb(m.pos(hov));`
227		227
228	`if(use_alpha)`	228	`if(use_alpha)`
229	`return edge_alpha_energy(va,vb,a);`	229	`return edge_alpha_energy(va,vb,a);`
230		230
231	`return edge_c_energy(va,vb,a);`	231	`return edge_c_energy(va,vb,a);`
232	`}`	232	`}`
233		233
234		234
235	`double DihedralEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const`	235	`double DihedralEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const`
236	`{`	236	`{`
237	`compute_angles(m, h);`	237	`compute_angles(m, h);`
238		238
239	`Walker w = m.walker(h);`	239	`Walker w = m.walker(h);`
240		240
241	`VertexID hv = w.vertex();`	241	`VertexID hv = w.vertex();`
242	`VertexID hov = w.opp().vertex();`	242	`VertexID hov = w.opp().vertex();`
243	`VertexID hnv = w.next().vertex();`	243	`VertexID hnv = w.next().vertex();`
244	`VertexID honv = w.opp().next().vertex();`	244	`VertexID honv = w.opp().next().vertex();`
245		245
246	`Vec3d va(m.pos(hv));`	246	`Vec3d va(m.pos(hv));`
247	`Vec3d vb(m.pos(hov));`	247	`Vec3d vb(m.pos(hov));`
248	`Vec3d vc(m.pos(hnv));`	248	`Vec3d vc(m.pos(hnv));`
249	`Vec3d vd(m.pos(honv));`	249	`Vec3d vd(m.pos(honv));`
250		250
251	`if(use_alpha){`	251	`if(use_alpha){`
252	`double before =`	252	`double before =`
253	`edge_alpha_energy(va,vb,ab_12)`	253	`edge_alpha_energy(va,vb,ab_12)`
254	`+edge_alpha_energy(va,vc,ab_a1)`	254	`+edge_alpha_energy(va,vc,ab_a1)`
255	`+edge_alpha_energy(vc,vb,ab_b1)`	255	`+edge_alpha_energy(vc,vb,ab_b1)`
256	`+edge_alpha_energy(vd,vb,ab_2c)`	256	`+edge_alpha_energy(vd,vb,ab_2c)`
257	`+edge_alpha_energy(vd,va,ab_2d);`	257	`+edge_alpha_energy(vd,va,ab_2d);`
258		258
259	`double after =`	259	`double after =`
260	`edge_alpha_energy(vd,vc,aa_12)`	260	`edge_alpha_energy(vd,vc,aa_12)`
261	`+edge_alpha_energy(vb,vc,aa_b1)`	261	`+edge_alpha_energy(vb,vc,aa_b1)`
262	`+edge_alpha_energy(vd,vb,aa_c1)`	262	`+edge_alpha_energy(vd,vb,aa_c1)`
263	`+edge_alpha_energy(va,vc,aa_2a)`	263	`+edge_alpha_energy(va,vc,aa_2a)`
264	`+edge_alpha_energy(vd,va,aa_2d);`	264	`+edge_alpha_energy(vd,va,aa_2d);`
265		265
266	`return (after-before);`	266	`return (after-before);`
267	`}`	267	`}`
268	`double before =`	268	`double before =`
269	`edge_c_energy(va,vb,ab_12)`	269	`edge_c_energy(va,vb,ab_12)`
270	`+edge_c_energy(va,vc,ab_a1)`	270	`+edge_c_energy(va,vc,ab_a1)`
271	`+edge_c_energy(vc,vb,ab_b1)`	271	`+edge_c_energy(vc,vb,ab_b1)`
272	`+edge_c_energy(vd,vb,ab_2c)`	272	`+edge_c_energy(vd,vb,ab_2c)`
273	`+edge_c_energy(vd,va,ab_2d);`	273	`+edge_c_energy(vd,va,ab_2d);`
274		274
275	`double after =`	275	`double after =`
276	`edge_c_energy(vd,vc,aa_12)`	276	`edge_c_energy(vd,vc,aa_12)`
277	`+edge_c_energy(vb,vc,aa_b1)`	277	`+edge_c_energy(vb,vc,aa_b1)`
278	`+edge_c_energy(vd,vb,aa_c1)`	278	`+edge_c_energy(vd,vb,aa_c1)`
279	`+edge_c_energy(va,vc,aa_2a)`	279	`+edge_c_energy(va,vc,aa_2a)`
280	`+edge_c_energy(vd,va,aa_2d);`	280	`+edge_c_energy(vd,va,aa_2d);`
281		281
282	`return after-before;`	282	`return after-before;`
283	`}`	283	`}`
284		284
285		285
286		286
287	`double CurvatureEnergy::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`
288	`{`	288	`{`
289	`const size_t N = ring.size();`	289	`const size_t N = ring.size();`
290		290
291	`double H = 0;`	291	`double H = 0;`
292	`for(int i = 0; i < N; ++i){`	292	`for(int i = 0; i < N; ++i){`
293	`Vec3d vnim1 = ring[(i+N-1)%N] - v;`	293	`Vec3d vnim1 = ring[(i+N-1)%N] - v;`
294	`Vec3d vni = ring[i] - v;`	294	`Vec3d vni = ring[i] - v;`
295	`Vec3d vnip1 = ring[(i+1)%N] - v;`	295	`Vec3d vnip1 = ring[(i+1)%N] - v;`
296		296
297	`Vec3d Nm = normalize(cross(vni, vnim1));`	297	`Vec3d Nm = normalize(cross(vni, vnim1));`
298	`Vec3d Np = normalize(cross(vnip1, vni));`	298	`Vec3d Np = normalize(cross(vnip1, vni));`
299		299
300	`double beta = acos(max(-1.0, min(1.0, dot(Nm, Np))));`	300	`double beta = acos(max(-1.0, min(1.0, dot(Nm, Np))));`
301	`H += vni.length() * beta;`	301	`H += vni.length() * beta;`
302	`}`	302	`}`
303	`H /= 4;`	303	`H /= 4;`
304		304
305	`return H;`	305	`return H;`
306	`}`	306	`}`
307		307
308	`double CurvatureEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const`	308	`double CurvatureEnergy::delta_energy(const Manifold& m, HalfEdgeID h) const`
309	`{`	309	`{`
310	`Walker w = m.walker(h);`	310	`Walker w = m.walker(h);`
311		311
312	`VertexID va = w.vertex();`	312	`VertexID va = w.vertex();`
313	`VertexID vb = w.opp().vertex();`	313	`VertexID vb = w.opp().vertex();`
314	`VertexID vc = w.next().vertex();`	314	`VertexID vc = w.next().vertex();`
315	`VertexID vd = w.opp().next().vertex();`	315	`VertexID vd = w.opp().next().vertex();`
316		316
317		317
318		318
319	`Vec3d va_pos(m.pos(va));`	319	`Vec3d va_pos(m.pos(va));`
320	`Vec3d vb_pos(m.pos(vb));`	320	`Vec3d vb_pos(m.pos(vb));`
321	`Vec3d vc_pos(m.pos(vc));`	321	`Vec3d vc_pos(m.pos(vc));`
322	`Vec3d vd_pos(m.pos(vd));`	322	`Vec3d vd_pos(m.pos(vd));`
323		323
324	`vector<Vec3d> va_ring_bef;`	324	`vector<Vec3d> va_ring_bef;`
325	`vector<Vec3d> va_ring_aft;`	325	`vector<Vec3d> va_ring_aft;`
326	`vector<Vec3d> vb_ring_bef;`	326	`vector<Vec3d> vb_ring_bef;`
327	`vector<Vec3d> vb_ring_aft;`	327	`vector<Vec3d> vb_ring_aft;`
328	`vector<Vec3d> vc_ring_bef;`	328	`vector<Vec3d> vc_ring_bef;`
329	`vector<Vec3d> vc_ring_aft;`	329	`vector<Vec3d> vc_ring_aft;`
330	`vector<Vec3d> vd_ring_bef;`	330	`vector<Vec3d> vd_ring_bef;`
331	`vector<Vec3d> vd_ring_aft;`	331	`vector<Vec3d> vd_ring_aft;`
332		332
333	`for(Walker wv = m.walker(va); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`	333	`for(Walker wv = m.walker(va); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`
334	`VertexID v = wv.vertex();`	334	`VertexID v = wv.vertex();`
335	`Vec3d pos(m.pos(v));`	335	`Vec3d pos(m.pos(v));`
336		336
337	`va_ring_bef.push_back(pos);`	337	`va_ring_bef.push_back(pos);`
338	`if(v != vb)`	338	`if(v != vb)`
339	`va_ring_aft.push_back(pos);`	339	`va_ring_aft.push_back(pos);`
340	`}`	340	`}`
341	`for(Walker wv = m.walker(vb); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`	341	`for(Walker wv = m.walker(vb); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`
342	`VertexID v = wv.vertex();`	342	`VertexID v = wv.vertex();`
343	`Vec3d pos(m.pos(v));`	343	`Vec3d pos(m.pos(v));`
344		344
345	`vb_ring_bef.push_back(pos);`	345	`vb_ring_bef.push_back(pos);`
346	`if(v != va)`	346	`if(v != va)`
347	`vb_ring_aft.push_back(pos);`	347	`vb_ring_aft.push_back(pos);`
348	`}`	348	`}`
349	`for(Walker wv = m.walker(vc); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`	349	`for(Walker wv = m.walker(vc); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`
350	`VertexID v = wv.vertex();`	350	`VertexID v = wv.vertex();`
351	`Vec3d pos(m.pos(v));`	351	`Vec3d pos(m.pos(v));`
352		352
353	`vc_ring_bef.push_back(pos);`	353	`vc_ring_bef.push_back(pos);`
354	`vc_ring_aft.push_back(pos);`	354	`vc_ring_aft.push_back(pos);`
355	`if(v == va)`	355	`if(v == va)`
356	`vc_ring_aft.push_back(vd_pos);`	356	`vc_ring_aft.push_back(vd_pos);`
357	`}`	357	`}`
358	`for(Walker wv = m.walker(vd); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`	358	`for(Walker wv = m.walker(vd); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`
359	`VertexID v = wv.vertex();`	359	`VertexID v = wv.vertex();`
360	`Vec3d pos(m.pos(v));`	360	`Vec3d pos(m.pos(v));`
361		361
362	`vd_ring_bef.push_back(pos);`	362	`vd_ring_bef.push_back(pos);`
363	`vd_ring_aft.push_back(pos);`	363	`vd_ring_aft.push_back(pos);`
364	`if(v == vb)`	364	`if(v == vb)`
365	`vd_ring_aft.push_back(vc_pos);`	365	`vd_ring_aft.push_back(vc_pos);`
366	`}`	366	`}`
367	`double before =`	367	`double before =`
368	`abs_mean_curv(va_pos, va_ring_bef) +`	368	`abs_mean_curv(va_pos, va_ring_bef) +`
369	`abs_mean_curv(vb_pos, vb_ring_bef) +`	369	`abs_mean_curv(vb_pos, vb_ring_bef) +`
370	`abs_mean_curv(vc_pos, vc_ring_bef) +`	370	`abs_mean_curv(vc_pos, vc_ring_bef) +`
371	`abs_mean_curv(vd_pos, vd_ring_bef);`	371	`abs_mean_curv(vd_pos, vd_ring_bef);`
372		372
373	`double after =`	373	`double after =`
374	`abs_mean_curv(va_pos, va_ring_aft) +`	374	`abs_mean_curv(va_pos, va_ring_aft) +`
375	`abs_mean_curv(vb_pos, vb_ring_aft) +`	375	`abs_mean_curv(vb_pos, vb_ring_aft) +`
376	`abs_mean_curv(vc_pos, vc_ring_aft) +`	376	`abs_mean_curv(vc_pos, vc_ring_aft) +`
377	`abs_mean_curv(vd_pos, vd_ring_aft);`	377	`abs_mean_curv(vd_pos, vd_ring_aft);`
378		378
379	`return after-before;`	379	`return after-before;`
380	`}`	380	`}`
381		381
382		382
383	`class GaussCurvatureEnergy: public EnergyFun`	383	`class GaussCurvatureEnergy: public EnergyFun`
384	`{`	384	`{`
385		385
386	`double gauss_curv(const Vec3d& v, const vector<Vec3d>& ring) const`	386	`double gauss_curv(const Vec3d& v, const vector<Vec3d>& ring) const`
387	`{`	387	`{`
388	`const size_t N = ring.size();`	388	`const size_t N = ring.size();`
389	`double asum=0.0f;`	389	`double asum=0.0f;`
390	`double area_sum=0;`	390	`double area_sum=0;`
391	`for(int i = 0; i < N; ++i){`	391	`for(int i = 0; i < N; ++i){`
392	`const Vec3d& v1 = ring[i];`	392	`const Vec3d& v1 = ring[i];`
393	`const Vec3d& v2 = ring[(i+1)%N];`	393	`const Vec3d& v2 = ring[(i+1)%N];`
394	`Vec3d a = v1-v;`	394	`Vec3d a = v1-v;`
395	`Vec3d b = v2-v;`	395	`Vec3d b = v2-v;`
396	`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)))));`
397	`area_sum += 0.5 * length(cross(a,b));`	397	`area_sum += 0.5 * length(cross(a,b));`
398	`}`	398	`}`
399	`return 3abs(2 M_PI - asum)/area_sum;`	399	`return 3abs(2 M_PI - asum)/area_sum;`
400	`}`	400	`}`
401		401
402	`public:`	402	`public:`
403	`double delta_energy(const Manifold& m, HalfEdgeID h) const`	403	`double delta_energy(const Manifold& m, HalfEdgeID h) const`
404	`{`	404	`{`
405	`Walker w = m.walker(h);`	405	`Walker w = m.walker(h);`
406		406
407	`VertexID va = w.vertex();`	407	`VertexID va = w.vertex();`
408	`VertexID vb = w.opp().vertex();`	408	`VertexID vb = w.opp().vertex();`
409	`VertexID vc = w.next().vertex();`	409	`VertexID vc = w.next().vertex();`
410	`VertexID vd = w.opp().next().vertex();`	410	`VertexID vd = w.opp().next().vertex();`
411		411
412	`Vec3d va_pos(m.pos(va));`	412	`Vec3d va_pos(m.pos(va));`
413	`Vec3d vb_pos(m.pos(vb));`	413	`Vec3d vb_pos(m.pos(vb));`
414	`Vec3d vc_pos(m.pos(vc));`	414	`Vec3d vc_pos(m.pos(vc));`
415	`Vec3d vd_pos(m.pos(vd));`	415	`Vec3d vd_pos(m.pos(vd));`
416		416
417	`vector<Vec3d> va_ring_bef;`	417	`vector<Vec3d> va_ring_bef;`
418	`vector<Vec3d> va_ring_aft;`	418	`vector<Vec3d> va_ring_aft;`
419	`vector<Vec3d> vb_ring_bef;`	419	`vector<Vec3d> vb_ring_bef;`
420	`vector<Vec3d> vb_ring_aft;`	420	`vector<Vec3d> vb_ring_aft;`
421	`vector<Vec3d> vc_ring_bef;`	421	`vector<Vec3d> vc_ring_bef;`
422	`vector<Vec3d> vc_ring_aft;`	422	`vector<Vec3d> vc_ring_aft;`
423	`vector<Vec3d> vd_ring_bef;`	423	`vector<Vec3d> vd_ring_bef;`
424	`vector<Vec3d> vd_ring_aft;`	424	`vector<Vec3d> vd_ring_aft;`
425		425
426	`for(Walker wv = m.walker(va); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`	426	`for(Walker wv = m.walker(va); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`
427	`VertexID v = wv.vertex();`	427	`VertexID v = wv.vertex();`
428	`Vec3d pos(m.pos(v));`	428	`Vec3d pos(m.pos(v));`
429		429
430	`va_ring_bef.push_back(pos);`	430	`va_ring_bef.push_back(pos);`
431	`if(v != vb)`	431	`if(v != vb)`
432	`va_ring_aft.push_back(pos);`	432	`va_ring_aft.push_back(pos);`
433	`}`	433	`}`
434	`for(Walker wv = m.walker(vb); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`	434	`for(Walker wv = m.walker(vb); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`
435	`VertexID v = wv.vertex();`	435	`VertexID v = wv.vertex();`
436	`Vec3d pos(m.pos(v));`	436	`Vec3d pos(m.pos(v));`
437		437
438	`vb_ring_bef.push_back(pos);`	438	`vb_ring_bef.push_back(pos);`
439	`if(v != va)`	439	`if(v != va)`
440	`vb_ring_aft.push_back(pos);`	440	`vb_ring_aft.push_back(pos);`
441	`}`	441	`}`
442	`for(Walker wv = m.walker(vc); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`	442	`for(Walker wv = m.walker(vc); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`
443	`VertexID v = wv.vertex();`	443	`VertexID v = wv.vertex();`
444	`Vec3d pos(m.pos(v));`	444	`Vec3d pos(m.pos(v));`
445		445
446	`vc_ring_bef.push_back(pos);`	446	`vc_ring_bef.push_back(pos);`
447	`vc_ring_aft.push_back(pos);`	447	`vc_ring_aft.push_back(pos);`
448	`if(v == va)`	448	`if(v == va)`
449	`vc_ring_aft.push_back(pos);`	449	`vc_ring_aft.push_back(pos);`
450	`}`	450	`}`
451	`for(Walker wv = m.walker(vd); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`	451	`for(Walker wv = m.walker(vd); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`
452	`VertexID v = wv.vertex();`	452	`VertexID v = wv.vertex();`
453	`Vec3d pos(m.pos(v));`	453	`Vec3d pos(m.pos(v));`
454		454
455	`vd_ring_bef.push_back(pos);`	455	`vd_ring_bef.push_back(pos);`
456	`vd_ring_aft.push_back(pos);`	456	`vd_ring_aft.push_back(pos);`
457	`if(v == vb)`	457	`if(v == vb)`
458	`vd_ring_aft.push_back(pos);`	458	`vd_ring_aft.push_back(pos);`
459	`}`	459	`}`
460	`double before =`	460	`double before =`
461	`gauss_curv(va_pos,va_ring_bef) +`	461	`gauss_curv(va_pos,va_ring_bef) +`
462	`gauss_curv(vb_pos,vb_ring_bef) +`	462	`gauss_curv(vb_pos,vb_ring_bef) +`
463	`gauss_curv(vc_pos,vc_ring_bef) +`	463	`gauss_curv(vc_pos,vc_ring_bef) +`
464	`gauss_curv(vd_pos,vd_ring_bef);`	464	`gauss_curv(vd_pos,vd_ring_bef);`
465		465
466	`double after =`	466	`double after =`
467	`gauss_curv(va_pos,va_ring_aft) +`	467	`gauss_curv(va_pos,va_ring_aft) +`
468	`gauss_curv(vb_pos,vb_ring_aft) +`	468	`gauss_curv(vb_pos,vb_ring_aft) +`
469	`gauss_curv(vc_pos,vc_ring_aft) +`	469	`gauss_curv(vc_pos,vc_ring_aft) +`
470	`gauss_curv(vd_pos,vd_ring_aft);`	470	`gauss_curv(vd_pos,vd_ring_aft);`
471		471
472	`return after-before;`	472	`return after-before;`
473	`}`	473	`}`
474	`};`	474	`};`
475		475
476	`struct PQElem`	476	`struct PQElem`
477	`{`	477	`{`
478	`double pri;`	478	`double pri;`
479	`HalfEdgeID h;`	479	`HalfEdgeID h;`
480	`int time;`	480	`int time;`
481		481
482	`//PQElem() {}`	482	`//PQElem() {}`
483	`PQElem(double _pri, HalfEdgeID _h, int _time):`	483	`PQElem(double _pri, HalfEdgeID _h, int _time):`
484	`pri(_pri), h(_h), time(_time) {}`	484	`pri(_pri), h(_h), time(_time) {}`
485	`};`	485	`};`
486		486
487	`bool operator<(const PQElem& e0, const PQElem& e1)`	487	`bool operator<(const PQElem& e0, const PQElem& e1)`
488	`{`	488	`{`
489	`return e0.pri > e1.pri;`	489	`return e0.pri > e1.pri;`
490	`}`	490	`}`
491		491
492		492
493	`void add_to_queue(const Manifold& m, HalfEdgeAttributeVector<int>& touched, priority_queue<PQElem>& Q, HalfEdgeID h, const EnergyFun& efun)`	493	`void add_to_queue(const Manifold& m, HalfEdgeAttributeVector<int>& touched, priority_queue<PQElem>& Q, HalfEdgeID h, const EnergyFun& efun)`
494	`{`	494	`{`
495	`if(boundary(m, h))`	495	`if(boundary(m, h))`
496	`return;`	496	`return;`
497		497
498	`Walker w = m.walker(h);`	498	`Walker w = m.walker(h);`
499	`HalfEdgeID ho = w.opp().halfedge();`	499	`HalfEdgeID ho = w.opp().halfedge();`
500		500

Subversion Repositories gelsvn

(root)/trunk/src/HMesh/mesh_optimization.cpp – Rev 595 → 601