Subversion Repositories gelsvn

Rev

Rev 191 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 191 Rev 192
Line 29... Line 29...
29 
	using namespace std;
 29 
	using namespace std;
30 
 
 30 
 
31 
	namespace
31 
	namespace
32 
	{
 32 
	{
33 
 
 33 
 
34 
			template<class T>
 34 
		template<class T>
35 
			Vec3f push_vertex(const RGrid<T>& grid, const Vec3f& pos0, float iso)
 35 
		Vec3f push_vertex(const RGrid<T>& grid, const Vec3f& pos0, float iso)
36 
			{
 36 
		{
37 
					Geometry::TrilinFilter<Geometry::RGrid<T> > inter(&grid);
 37 
			Geometry::TrilinFilter<Geometry::RGrid<T> > inter(&grid);
- 
 
 38 
			Vec3f corner_min = pos0 - Vec3f(1);
- 
 
 39 
			Vec3f corner_max = pos0 + Vec3f(1);
38 
					Vec3f pos = pos0;
 40 
			Vec3f pos = pos0;
39 
					for(int i=0;i<100;++i)
 41 
			for(int i=0;i<1000;++i)
40 
							if(inter.in_domain(pos))
 42 
				if(inter.in_domain(pos))
41 
							{
 43 
					{
42 
									float val = inter(pos)-iso;
 44 
						float val = inter(pos)-iso;
43 
									Vec3f grad = inter.grad(pos);
 45 
						Vec3f grad = inter.grad(pos);
44 
									float sqlen_grad = sqr_length(grad);
 46 
						float sqlen_grad = sqr_length(grad);
- 
 
 47 
 
- 
 
 48 
						// Compute the projection in the gradient direction
45 
									Vec3f d1 = - val * grad/sqlen_grad;
 49 
						Vec3f d1 = - val * grad/sqlen_grad;
- 
 
 50 
						Vec3f p_new = pos + 0.5 * d1; // "0.5 *" to add damping
- 
 
 51 
 
- 
 
 52 
						// Don't move outside the box.
- 
 
 53 
						if(!(corner_max.all_g(p_new) &&
- 
 
 54 
								 corner_min.all_l(p_new)))
- 
 
 55 
							return pos;
- 
 
 56
- 
 
 57 
						// Move so that the gradient points in the direction
- 
 
 58 
						// of the original position
46 
									Vec3f p_half = pos + d1;
 59 
						Vec3f d2 = (pos0-pos) -
- 
 
 60 
							grad*dot((pos0-pos),grad)/sqlen_grad;
47
61
48 
									Vec3f d2 =
- 
 
49 
											(pos0-pos) -
- 
 
50 
											grad*dot((pos0-pos),grad)/sqlen_grad;
 62 
						// Compute new position and distance up to original position.
51
- 
 
52 
									pos = p_half + d2;
 63 
						p_new += 0.5 * d2; // "0.5 *" to add damping
53
64
- 
 
 65 
						// Don't move outside the box.
54 
									if(sqrt(sqr_length(d1)+sqr_length(d2))<.001)
 66 
						if(corner_max.all_g(p_new) &&
- 
 
 67 
							 corner_min.all_l(p_new))
55 
											return pos;
 68 
							pos = p_new;
56 
							}
 69 
						else
57 
					return pos;
 70 
							return pos;
58 
			}
 - 
 
59 
 
 71 
 
- 
 
 72 
						// Did we converge?
- 
 
 73 
						if(sqrt(sqr_length(d1)+sqr_length(d2))<.001)
- 
 
 74 
							return pos;
- 
 
 75 
					}
- 
 
 76 
			return pos;
- 
 
 77 
		}
- 
 
 78
60 
		void compute_dual(Manifold& m, Manifold& m2, const vector<Vec3f>& face_positions)
 79 
		void compute_dual(Manifold& m, Manifold& m2, const vector<Vec3f>& face_positions)
61 
		{
 80 
		{
62 
			// make sure every face knows its number
 81 
			// make sure every face knows its number
63 
			m.enumerate_faces();
 82 
			m.enumerate_faces();
64 
 
 83