ghc/tests/programs/galois_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