Subversion Repositories gelsvn

Rev

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

Rev 191 Rev 192
Line 33... Line 33...
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) &&
46 
									Vec3f p_half = pos + d1;
 54 
								 corner_min.all_l(p_new)))
- 
 
 55 
							return pos;
47
56
- 
 
 57 
						// Move so that the gradient points in the direction
48 
									Vec3f d2 =
58 
						// of the original position
49 
											(pos0-pos) -
59 
						Vec3f d2 = (pos0-pos) -
50 
											grad*dot((pos0-pos),grad)/sqlen_grad;
 60 
							grad*dot((pos0-pos),grad)/sqlen_grad;
51
61
- 
 
 62 
						// Compute new position and distance up to original position.
- 
 
 63 
						p_new += 0.5 * d2; // "0.5 *" to add damping
- 
 
 64
- 
 
 65 
						// Don't move outside the box.
- 
 
 66 
						if(corner_max.all_g(p_new) &&
- 
 
 67 
							 corner_min.all_l(p_new))
52 
									pos = p_half + d2;
 68 
							pos = p_new;
- 
 
 69 
						else
- 
 
 70 
							return pos;
53
71 
 
- 
 
 72 
						// Did we converge?
54 
									if(sqrt(sqr_length(d1)+sqr_length(d2))<.001)
 73 
						if(sqrt(sqr_length(d1)+sqr_length(d2))<.001)
55 
											return pos;
 74 
							return pos;
56 
							}
 75 
					}
57 
					return pos;
 76 
			return pos;
58 
			}
 77 
		}