|
|
1.1 root 1: To: research!dbm
2: Status: R
3:
4: (* Find the intersection of a ray with a cylinder in three-space. If
5: there is more than one intersection, return the parameter of
6: intersection closest to the start of the ray. If there is no
7: intersection, return NONE.
8:
9: The ray is specified by its starting point (rayStart) and direction
10: vector (rayVec). The cylinder is specified by a starting point (cStart)
11: a length vector (cVec), and a radius (cRad). The cylinder is finite;
12: it does not extend beyond the direction vector at either end.
13:
14: *)
15:
16: local open Real
17: val EP = 0.0000001
18:
19: fun newton f x =
20: let val y = f x
21: in if abs(y)>EP then newton f (x+EP*y/(y-f(x+EP))) else x
22: end
23:
24: infix 7 @ fun (a,b,c) @ (d,e,f) = a*d+b*e+c*f (* dot product *)
25: infix 7 *! fun (a,b,c) *! t = (a*t,b*t,c*t) (* vector * scalar *)
26: infix 6 +! fun (a,b,c) +! (d,e,f) = (a+d,b+e,c+f) (* vector addition *)
27: infix 6 -! fun (a,b,c) -! (d,e,f) = (a-d,b-e,c-f) (* vector subtraction *)
28:
29: fun along (C,D) t = C +! (D *! t)
30:
31: in
32:
33: fun cylhit (ray as (rayStart, rayVec), cyl as (cStart, cVec), cRad) =
34: let val rayStart = rayStart -! cStart
35: fun cylpos p = (cVec@p) / (cVec@cVec)
36: (* parameter of closest point on cVec to the point p *)
37:
38: fun cyldist p = let val v = p -! (along cyl (cylpos p)) in v@v end
39: (* distance of p from closest point on cVec *)
40:
41: fun cyldist' p = p*!(cVec@cVec) -! cVec*!(p@cVec)
42: (* gradient of the cyldist function *)
43:
44: val tf = newton(fn t => cyldist' (along ray t) @ rayVec) 0.0
45: (* parameter along the ray of the closest point to cVec *)
46:
47: in if cyldist (along ray tf) > cRad*cRad then NONE else
48: let val t0 = newton(fn t => cyldist(along ray t)-(cRad*cRad)) 0.0
49: (* one of the two intersection points *)
50: val t1 = newton(fn t => cyldist(along ray t)-(cRad*cRad)) (tf+tf-t0)
51: (* the other intersection point *)
52: fun inRange v = v >= 0.0 andalso v <= 1.0
53: in case (t0 > 0.0 andalso inRange (cylpos (along ray t0)),
54: t1 > 0.0 andalso inRange (cylpos (along ray t1)))
55: of (false,false) => NONE
56: | (true,false) => SOME t0
57: | (false,true) => SOME t1
58: | (true,true) => SOME (if (t0<t1) then t0 else t1)
59: end
60: end
61:
62: end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.