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 6module Intersections 7 ( intersectRayWithObject, 8 quadratic 9 ) where 10 11import Data.Maybe(isJust) 12 13import Construct 14import Geometry 15import Interval 16import 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 23intersectRayWithObject ray p 24 = clampIntervals is 25 where is = intersections ray p 26 27clampIntervals (True, [], True) = (False, [(0, True, undefined)], True) 28clampIntervals empty@(False, [], False) = empty 29clampIntervals (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) 34clampIntervals 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 41intersections ray (Union p q) 42 = unionIntervals is js 43 where is = intersections ray p 44 js = intersections ray q 45 46intersections ray (Intersect p q) 47 = intersectIntervals is js 48 where is = intersections ray p 49 js = intersections ray q 50 51intersections ray (Difference p q) 52 = differenceIntervals is (negateSurfaces js) 53 where is = intersections ray p 54 js = intersections ray q 55 56intersections 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 61intersections ray (Box box p) 62 | intersectWithBox ray box = intersections ray p 63 | otherwise = emptyIList 64 65intersections ray p@(Plane s) 66 = intersectPlane ray s 67 68intersections ray p@(Sphere s) 69 = intersectSphere ray s 70 71intersections ray p@(Cube s) 72 = intersectCube ray s 73 74intersections ray p@(Cylinder s) 75 = intersectCylinder ray s 76 77intersections ray p@(Cone s) 78 = intersectCone ray s 79 80negateSurfaces :: IList (Surface, Texture a) -> IList (Surface, Texture a) 81negateSurfaces = mapI negSurf 82 where negSurf (i, b, (s,t)) = (i, b, (negateSurface s, t)) 83 84negateSurface (Planar p0 v0 v1) 85 = Planar p0 v1 v0 86negateSurface (Spherical p0 v0 v1) 87 = Spherical p0 v1 v0 88negateSurface (Cylindrical p0 v0 v1) 89 = Cylindrical p0 v1 v0 90negateSurface (Conic p0 v0 v1) 91 = Conic p0 v1 v0 92 93transformSurface 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 99transformSurface 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 106transformSurface 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 112transformSurface 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 122intersectPlane :: Ray -> a -> IList (Surface, Texture a) 123intersectPlane ray texture = intersectXZPlane PlaneFace ray 0.0 texture 124 125intersectXZPlane :: Face -> Ray -> Double -> a -> IList (Surface, Texture a) 126intersectXZPlane 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 161intersectSphere :: Ray -> a -> IList (Surface, Texture a) 162intersectSphere 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 195intersectCube :: Ray -> a -> IList (Surface, Texture a) 196intersectCube 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 226intersectCubeSlab 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 232intersectSlabIval Nothing Nothing = Nothing 233intersectSlabIval Nothing (Just i) = Nothing 234intersectSlabIval (Just i) Nothing = Nothing 235intersectSlabIval (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 251intersectCylinder :: Ray -> a -> IList (Surface, Texture a) 252intersectCylinder 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 259intersectCylSide (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 299intersectCone :: Ray -> a -> IList (Surface, Texture a) 300intersectCone 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 307intersectConeSide (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 356quadratic :: Double -> Double -> Double -> Maybe (Double, Double) 357quadratic 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 373data MaybeInterval = Interval !Double !Double 374 | NoInterval 375 376isInterval (Interval _ _) = True 377isInterval _ = False 378 379intersectWithBox :: Ray -> Box -> Bool 380intersectWithBox (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 388intersectInterval :: MaybeInterval -> MaybeInterval -> MaybeInterval 389intersectInterval NoInterval _ = NoInterval 390intersectInterval _ NoInterval = NoInterval 391intersectInterval (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 #-} 396intersectRayWithSlab :: Double -> Double -> (Double,Double) -> MaybeInterval 397intersectRayWithSlab 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 404infInterval = Interval (-inf) inf 405