add Galois' Ray Tracer
[nofib.git] / parallel / gray / Intersections.hs
1 -- Copyright (c) 2000 Galois Connections, Inc.
2 -- All rights reserved. This software is distributed as
3 -- free software under the license in the file "LICENSE",
4 -- which is included in the distribution.
5
6 module Intersections
7 ( intersectRayWithObject,
8 quadratic
9 ) where
10
11 import Maybe(isJust)
12
13 import Construct
14 import Geometry
15 import Interval
16 import Misc
17
18 -- This is factored into two bits. The main function `intersections'
19 -- intersects a line with an object.
20 -- The wrapper call `intersectRayWithObject' coerces this to an intersection
21 -- with a ray by clamping the result to start at 0.
22
23 intersectRayWithObject ray p
24 = clampIntervals is
25 where is = intersections ray p
26
27 clampIntervals (True, [], True) = (False, [(0, True, undefined)], True)
28 clampIntervals empty@(False, [], False) = empty
29 clampIntervals (True, is@((i, False, p) : is'), isOpen)
30 | i `near` 0 || i < 0
31 = clampIntervals (False, is', isOpen)
32 | otherwise
33 = (False, (0, True, undefined) : is, isOpen)
34 clampIntervals ivals@(False, is@((i, True, p) : is'), isOpen)
35 | i `near` 0 || i < 0
36 -- can unify this with first case...
37 = clampIntervals (True, is', isOpen)
38 | otherwise
39 = ivals
40
41 intersections ray (Union p q)
42 = unionIntervals is js
43 where is = intersections ray p
44 js = intersections ray q
45
46 intersections ray (Intersect p q)
47 = intersectIntervals is js
48 where is = intersections ray p
49 js = intersections ray q
50
51 intersections ray (Difference p q)
52 = differenceIntervals is (negateSurfaces js)
53 where is = intersections ray p
54 js = intersections ray q
55
56 intersections ray (Transform m m' p)
57 = mapI (xform m) is
58 where is = intersections (m' `multMR` ray) p
59 xform m (i, b, (s, p0)) = (i, b, (transformSurface m s, p0))
60
61 intersections ray (Box box p)
62 | intersectWithBox ray box = intersections ray p
63 | otherwise = emptyIList
64
65 intersections ray p@(Plane s)
66 = intersectPlane ray s
67
68 intersections ray p@(Sphere s)
69 = intersectSphere ray s
70
71 intersections ray p@(Cube s)
72 = intersectCube ray s
73
74 intersections ray p@(Cylinder s)
75 = intersectCylinder ray s
76
77 intersections ray p@(Cone s)
78 = intersectCone ray s
79
80 negateSurfaces :: IList (Surface, Texture a) -> IList (Surface, Texture a)
81 negateSurfaces = mapI negSurf
82 where negSurf (i, b, (s,t)) = (i, b, (negateSurface s, t))
83
84 negateSurface (Planar p0 v0 v1)
85 = Planar p0 v1 v0
86 negateSurface (Spherical p0 v0 v1)
87 = Spherical p0 v1 v0
88 negateSurface (Cylindrical p0 v0 v1)
89 = Cylindrical p0 v1 v0
90 negateSurface (Conic p0 v0 v1)
91 = Conic p0 v1 v0
92
93 transformSurface m (Planar p0 v0 v1)
94 = Planar p0' v0' v1'
95 where p0' = multMP m p0
96 v0' = multMV m v0
97 v1' = multMV m v1
98
99 transformSurface m (Spherical p0 v0 v1)
100 = Spherical p0' v0' v1'
101 where p0' = multMP m p0
102 v0' = multMV m v0
103 v1' = multMV m v1
104
105 -- ditto as above
106 transformSurface m (Cylindrical p0 v0 v1)
107 = Cylindrical p0' v0' v1'
108 where p0' = multMP m p0
109 v0' = multMV m v0
110 v1' = multMV m v1
111
112 transformSurface m (Conic p0 v0 v1)
113 = Conic p0' v0' v1'
114 where p0' = multMP m p0
115 v0' = multMV m v0
116 v1' = multMV m v1
117
118 --------------------------------
119 -- Plane
120 --------------------------------
121
122 intersectPlane :: Ray -> a -> IList (Surface, Texture a)
123 intersectPlane ray texture = intersectXZPlane PlaneFace ray 0.0 texture
124
125 intersectXZPlane :: Face -> Ray -> Double -> a -> IList (Surface, Texture a)
126 intersectXZPlane n (r,v) yoffset texture
127 | b `near` 0
128 = -- the ray is parallel to the plane - it's either all in, or all out
129 if y `near` yoffset || y < yoffset then openIList else emptyIList
130
131 -- The line intersects the plane. Find t such that
132 -- (x + at, y + bt, z + ct) intersects the X-Z plane.
133 -- t may be negative (the ray starts within the halfspace),
134 -- but we'll catch that later when we clamp the intervals
135
136 | b < 0 -- the ray is pointing downwards
137 = (False, [mkEntry (t0, (Planar p0 v0 v1, (n, p0, texture)))], True)
138
139 | otherwise -- the ray is pointing upwards
140 = (True, [mkExit (t0, (Planar p0 v0 v1, (n, p0, texture)))], False)
141
142 where t0 = (yoffset-y) / b
143 x0 = x + a * t0
144 z0 = z + c * t0
145 p0 = point x0 0 z0
146 v0 = vector 0 0 1
147 v1 = vector 1 0 0
148
149 x = xCoord r
150 y = yCoord r
151 z = zCoord r
152 a = xComponent v
153 b = yComponent v
154 c = zComponent v
155
156
157 --------------------------------
158 -- Sphere
159 --------------------------------
160
161 intersectSphere :: Ray -> a -> IList (Surface, Texture a)
162 intersectSphere ray@(r, v) texture
163 = -- Find t such that (x + ta, y + tb, z + tc) intersects the
164 -- unit sphere, that is, such that:
165 -- (x + ta)^2 + (y + tb)^2 + (z + tc)^2 = 1
166 -- This is a quadratic equation in t:
167 -- t^2(a^2 + b^2 + c^2) + 2t(xa + yb + zc) + (x^2 + y^2 + z^2 - 1) = 0
168 let c1 = sq a + sq b + sq c
169 c2 = 2 * (x * a + y * b + z * c)
170 c3 = sq x + sq y + sq z - 1
171 in
172 case quadratic c1 c2 c3 of
173 Nothing -> emptyIList
174 Just (t1, t2) -> entryexit (g t1) (g t2)
175 where x = xCoord r
176 y = yCoord r
177 z = zCoord r
178 a = xComponent v
179 b = yComponent v
180 c = zComponent v
181 g t = (t, (Spherical origin v1 v2, (SphereFace, p0, texture)))
182 where origin = point 0 0 0
183 x0 = x + t * a
184 y0 = y + t * b
185 z0 = z + t * c
186 p0 = point x0 y0 z0
187 v0 = vector x0 y0 z0
188 (v1, v2) = tangents v0
189
190
191 --------------------------------
192 -- Cube
193 --------------------------------
194
195 intersectCube :: Ray -> a -> IList (Surface, Texture a)
196 intersectCube ray@(r, v) texture
197 = -- The set of t such that (x + at, y + bt, z + ct) lies within
198 -- the unit cube satisfies:
199 -- 0 <= x + at <= 1, 0 <= y + bt <= 1, 0 <= z + ct <= 1
200 -- The minimum and maximum such values of t give us the two
201 -- intersection points.
202 case intersectSlabIval (intersectCubeSlab face2 face3 x a)
203 (intersectSlabIval (intersectCubeSlab face5 face4 y b)
204 (intersectCubeSlab face0 face1 z c)) of
205 Nothing -> emptyIList
206 Just (t1, t2) -> entryexit (g t1) (g t2)
207 where g ((n, v0, v1), t)
208 = (t, (Planar p0 v0 v1, (n, p0, texture)))
209 where p0 = offsetToPoint ray t
210 face0 = (CubeFront, vectorY, vectorX)
211 face1 = (CubeBack, vectorX, vectorY)
212 face2 = (CubeLeft, vectorZ, vectorY)
213 face3 = (CubeRight, vectorY, vectorZ)
214 face4 = (CubeTop, vectorZ, vectorX)
215 face5 = (CubeBottom, vectorX, vectorZ)
216 vectorX = vector 1 0 0
217 vectorY = vector 0 1 0
218 vectorZ = vector 0 0 1
219 x = xCoord r
220 y = yCoord r
221 z = zCoord r
222 a = xComponent v
223 b = yComponent v
224 c = zComponent v
225
226 intersectCubeSlab n m w d
227 | d `near` 0 = if (0 <= w) && (w <= 1)
228 then Just ((n, -inf), (m, inf)) else Nothing
229 | d > 0 = Just ((n, (-w)/d), (m, (1-w)/d))
230 | otherwise = Just ((m, (1-w)/d), (n, (-w)/d))
231
232 intersectSlabIval Nothing Nothing = Nothing
233 intersectSlabIval Nothing (Just i) = Nothing
234 intersectSlabIval (Just i) Nothing = Nothing
235 intersectSlabIval (Just (nu1@(n1, u1), mv1@(m1, v1)))
236 (Just (nu2@(n2, u2), mv2@(m2, v2)))
237 = checkInterval (nu, mv)
238 where nu = if u1 < u2 then nu2 else nu1
239 mv = if v1 < v2 then mv1 else mv2
240 checkInterval numv@(nu@(_, u), (m, v))
241 -- rounding error may force us to push v out a bit
242 | u `near` v = Just (nu, (m, u + epsilon))
243 | u < v = Just numv
244 | otherwise = Nothing
245
246
247 --------------------------------
248 -- Cylinder
249 --------------------------------
250
251 intersectCylinder :: Ray -> a -> IList (Surface, Texture a)
252 intersectCylinder ray texture
253 = isectSide `intersectIntervals` isectTop `intersectIntervals` isectBottom
254 where isectSide = intersectCylSide ray texture
255 isectTop = intersectXZPlane CylinderTop ray 1.0 texture
256 isectBottom = complementIntervals $ negateSurfaces $
257 intersectXZPlane CylinderBottom ray 0.0 texture
258
259 intersectCylSide (r, v) texture
260 = -- The ray (x + ta, y + tb, z + tc) intersects the sides of the
261 -- cylinder if:
262 -- (x + ta)^2 + (z + tc)^2 = 1 and 0 <= y + tb <= 1.
263 if (sq a + sq c) `near` 0
264 then -- The ray is parallel to the Y-axis, and does not intersect
265 -- the cylinder sides. It's either all in, or all out
266 if (sqxy `near` 1.0 || sqxy < 1.0) then openIList else emptyIList
267 else -- Find values of t that solve the quadratic equation
268 -- (a^2 + c^2)t^2 + 2(ax + cz)t + x^2 + z^2 - 1 = 0
269 let c1 = sq a + sq c
270 c2 = 2 * (x * a + z * c)
271 c3 = sq x + sq z - 1
272 in
273 case quadratic c1 c2 c3 of
274 Nothing -> emptyIList
275 Just (t1, t2) -> entryexit (g t1) (g t2)
276
277 where sqxy = sq x + sq y
278 g t = (t, (Cylindrical origin v1 v2, (CylinderSide, p0, texture)))
279 where origin = point 0 0 0
280 x0 = x + t * a
281 y0 = y + t * b
282 z0 = z + t * c
283 p0 = point x0 y0 z0
284 v0 = vector x0 0 z0
285 (v1, v2) = tangents v0
286
287 x = xCoord r
288 y = yCoord r
289 z = zCoord r
290 a = xComponent v
291 b = yComponent v
292 c = zComponent v
293
294
295 -------------------
296 -- Cone
297 -------------------
298
299 intersectCone :: Ray -> a -> IList (Surface, Texture a)
300 intersectCone ray texture
301 = isectSide `intersectIntervals` isectTop `intersectIntervals` isectBottom
302 where isectSide = intersectConeSide ray texture
303 isectTop = intersectXZPlane ConeBase ray 1.0 texture
304 isectBottom = complementIntervals $ negateSurfaces $
305 intersectXZPlane ConeBase ray 0.0 texture
306
307 intersectConeSide (r, v) texture
308 = -- Find the points where the ray intersects the cond side. At any points of
309 -- intersection, we must have:
310 -- (x + ta)^2 + (z + tc)^2 = (y + tb)^2
311 -- which is the following quadratic equation:
312 -- t^2(a^2-b^2+c^2) + 2t(xa-yb+cz) + (x^2-y^2+z^2) = 0
313 let c1 = sq a - sq b + sq c
314 c2 = 2 * (x * a - y * b + c * z)
315 c3 = sq x - sq y + sq z
316 in case quadratic c1 c2 c3 of
317 Nothing -> emptyIList
318 Just (t1, t2) ->
319 -- If either intersection strikes the middle, then the other
320 -- can only be off by rounding error, so we make a tangent
321 -- strike using the "good" value.
322 -- If the intersections straddle the origin, then it's
323 -- an exit/entry pair, otherwise it's an entry/exit pair.
324 let y1 = y + t1 * b
325 y2 = y + t2 * b
326 in if y1 `near` 0 then entryexit (g t1) (g t1)
327 else if y2 `near` 0 then entryexit (g t2) (g t2)
328 else if (y1 < 0) `xor` (y2 < 0) then exitentry (g t1) (g t2)
329 else entryexit (g t1) (g t2)
330
331 where g t = (t, (Conic origin v1 v2, (ConeSide, p0, texture)))
332 where origin = point 0 0 0
333 x0 = x + t * a
334 y0 = y + t * b
335 z0 = z + t * c
336 p0 = point x0 y0 z0
337 v0 = normalize $ vector x0 (-y0) z0
338 (v1, v2) = tangents v0
339
340 x = xCoord r
341 y = yCoord r
342 z = zCoord r
343 a = xComponent v
344 b = yComponent v
345 c = zComponent v
346
347 -- beyond me why this isn't defined in the prelude...
348 xor False b = b
349 xor True b = not b
350
351
352 -------------------
353 -- Solving quadratics
354 -------------------
355
356 quadratic :: Double -> Double -> Double -> Maybe (Double, Double)
357 quadratic a b c =
358 -- Solve the equation ax^2 + bx + c = 0 by using the quadratic formula.
359 let d = sq b - 4 * a * c
360 d' = if d `near` 0 then 0 else d
361 in if d' < 0
362 then Nothing -- There are no real roots.
363 else
364 if a > 0 then Just (((-b) - sqrt d') / (2 * a),
365 ((-b) + sqrt d') / (2 * a))
366 else Just (((-b) + sqrt d') / (2 * a),
367 ((-b) - sqrt d') / (2 * a))
368
369 -------------------
370 -- Bounding boxes
371 -------------------
372
373 data MaybeInterval = Interval !Double !Double
374 | NoInterval
375
376 isInterval (Interval _ _) = True
377 isInterval _ = False
378
379 intersectWithBox :: Ray -> Box -> Bool
380 intersectWithBox (r, v) (B x1 x2 y1 y2 z1 z2)
381 = isInterval interval
382 where x_interval = intersectRayWithSlab (xCoord r) (xComponent v) (x1, x2)
383 y_interval = intersectRayWithSlab (yCoord r) (yComponent v) (y1, y2)
384 z_interval = intersectRayWithSlab (zCoord r) (zComponent v) (z1, z2)
385 interval = intersectInterval x_interval
386 (intersectInterval y_interval z_interval)
387
388 intersectInterval :: MaybeInterval -> MaybeInterval -> MaybeInterval
389 intersectInterval NoInterval _ = NoInterval
390 intersectInterval _ NoInterval = NoInterval
391 intersectInterval (Interval a b) (Interval c d)
392 | b < c || d < a = NoInterval
393 | otherwise = Interval (a `max` c) (b `min` d)
394
395 {-# INLINE intersectRayWithSlab #-}
396 intersectRayWithSlab :: Double -> Double -> (Double,Double) -> MaybeInterval
397 intersectRayWithSlab xCoord alpha (x1, x2)
398 | alpha == 0 = if xCoord < x1 || xCoord > x2 then NoInterval else infInterval
399 | alpha > 0 = Interval a b
400 | otherwise = Interval b a
401 where a = (x1 - xCoord) / alpha
402 b = (x2 - xCoord) / alpha
403
404 infInterval = Interval (-inf) inf