← Back to branch summary

~ubuntu-branches/ubuntu/precise/ghc/precise

« back to all changes in this revision

Viewing changes to libraries/hpc/tests/raytrace/Intersections.hs

  • Committer: Bazaar Package Importer
  • Author(s): Joachim Breitner
  • Date: 2011-01-17 12:49:24 UTC
  • Revision ID: james.westby@ubuntu.com-20110117124924-do1pym1jlf5o636m
Tags: upstream-7.0.1
Import upstream version 7.0.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
libraries/hpc/tests/raytrace/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