1 files changed, 404 insertions, 0 deletions
diff --git a/testsuite/tests/programs/galois_raytrace/Intersections.hs b/testsuite/tests/programs/galois_raytrace/Intersections.hs
new file mode 100644
index 0000000000..c7fe003eb3
--- /dev/null
+++ b/testsuite/tests/programs/galois_raytrace/Intersections.hs
@@ -0,0 +1,404 @@
+-- Copyright (c) 2000 Galois Connections, Inc.
+-- All rights reserved.  This software is distributed as
+-- free software under the license in the file "LICENSE",
+-- which is included in the distribution.
+
+module Intersections 
+    ( intersectRayWithObject,
+      quadratic
+    ) where
+
+import Maybe(isJust)
+
+import Construct
+import Geometry
+import Interval
+import Misc
+
+-- This is factored into two bits.  The main function `intersections'
+-- intersects a line with an object.
+-- The wrapper call `intersectRayWithObject' coerces this to an intersection
+-- with a ray by clamping the result to start at 0.
+
+intersectRayWithObject ray p
+  = clampIntervals is
+  where is = intersections ray p
+
+clampIntervals (True, [], True) = (False, [(0, True, undefined)], True)
+clampIntervals empty@(False, [], False) = empty
+clampIntervals (True, is@((i, False, p) : is'), isOpen)
+  | i `near` 0 || i < 0
+  = clampIntervals (False, is', isOpen)
+  | otherwise
+  = (False, (0, True, undefined) : is, isOpen)
+clampIntervals ivals@(False, is@((i, True, p) : is'), isOpen)
+  | i `near` 0 || i < 0
+  -- can unify this with first case...
+  = clampIntervals (True, is', isOpen)
+  | otherwise
+  = ivals
+
+intersections ray (Union p q)
+  = unionIntervals is js
+  where is = intersections ray p
+	js = intersections ray q
+
+intersections ray (Intersect p q)
+  = intersectIntervals is js
+  where is = intersections ray p
+	js = intersections ray q
+
+intersections ray (Difference p q)
+  = differenceIntervals is (negateSurfaces js)
+  where is = intersections ray p
+	js = intersections ray q
+
+intersections ray (Transform m m' p)
+  = mapI (xform m) is
+  where is = intersections (m' `multMR` ray) p
+	xform m (i, b, (s, p0)) = (i, b, (transformSurface m s, p0))
+
+intersections ray (Box box p)
+  | intersectWithBox ray box = intersections ray p
+  | otherwise = emptyIList
+
+intersections ray p@(Plane s)
+  = intersectPlane ray s
+
+intersections ray p@(Sphere s)
+  = intersectSphere ray s
+
+intersections ray p@(Cube s)
+  = intersectCube ray s
+
+intersections ray p@(Cylinder s)
+  = intersectCylinder ray s
+
+intersections ray p@(Cone s)
+  = intersectCone ray s
+
+negateSurfaces :: IList (Surface, Texture a) -> IList (Surface, Texture a)
+negateSurfaces = mapI negSurf
+  where negSurf (i, b, (s,t)) = (i, b, (negateSurface s, t))
+
+negateSurface (Planar p0 v0 v1)
+  = Planar p0 v1 v0
+negateSurface (Spherical p0 v0 v1)
+  = Spherical p0 v1 v0
+negateSurface (Cylindrical p0 v0 v1)
+  = Cylindrical p0 v1 v0
+negateSurface (Conic p0 v0 v1)
+  = Conic p0 v1 v0
+
+transformSurface m (Planar p0 v0 v1)
+  = Planar p0' v0' v1'
+  where p0' = multMP m p0
+	v0' = multMV m v0
+	v1' = multMV m v1
+
+transformSurface m (Spherical p0 v0 v1)
+  = Spherical p0' v0' v1'
+  where p0' = multMP m p0
+	v0' = multMV m v0
+	v1' = multMV m v1
+
+-- ditto as above
+transformSurface m (Cylindrical p0 v0 v1)
+  = Cylindrical p0' v0' v1'
+  where p0' = multMP m p0
+	v0' = multMV m v0
+	v1' = multMV m v1
+
+transformSurface m (Conic p0 v0 v1)
+  = Conic p0' v0' v1'
+  where p0' = multMP m p0
+	v0' = multMV m v0
+	v1' = multMV m v1
+
+--------------------------------
+-- Plane
+--------------------------------
+
+intersectPlane :: Ray -> a -> IList (Surface, Texture a)
+intersectPlane ray texture = intersectXZPlane PlaneFace ray 0.0 texture
+
+intersectXZPlane :: Face -> Ray -> Double -> a -> IList (Surface, Texture a)
+intersectXZPlane n (r,v) yoffset texture
+  | b `near` 0
+  = -- the ray is parallel to the plane - it's either all in, or all out
+    if y `near` yoffset || y < yoffset then openIList else emptyIList
+
+    -- The line intersects the plane. Find t such that
+    -- (x + at, y + bt, z + ct) intersects the X-Z plane.
+    -- t may be negative (the ray starts within the halfspace),
+    -- but we'll catch that later when we clamp the intervals
+
+  | b < 0	-- the ray is pointing downwards
+  = (False, [mkEntry (t0, (Planar p0 v0 v1, (n, p0, texture)))], True)
+
+  | otherwise	-- the ray is pointing upwards
+  = (True,  [mkExit (t0, (Planar p0 v0 v1, (n, p0, texture)))],  False)
+
+  where t0 = (yoffset-y) / b
+	x0 = x + a * t0
+	z0 = z + c * t0
+	p0 = point x0 0 z0
+	v0 = vector 0 0 1
+	v1 = vector 1 0 0
+
+	x = xCoord r
+	y = yCoord r
+	z = zCoord r
+	a = xComponent v
+	b = yComponent v
+	c = zComponent v
+
+
+--------------------------------
+-- Sphere
+--------------------------------
+
+intersectSphere :: Ray -> a -> IList (Surface, Texture a)
+intersectSphere ray@(r, v) texture
+  = -- Find t such that (x + ta, y + tb, z + tc) intersects the
+    -- unit sphere, that is, such that:
+    --   (x + ta)^2 + (y + tb)^2 + (z + tc)^2 = 1
+    -- This is a quadratic equation in t:
+    --   t^2(a^2 + b^2 + c^2) + 2t(xa + yb + zc) + (x^2 + y^2 + z^2 - 1) = 0
+    let c1 = sq a + sq b + sq c
+	c2 = 2 * (x * a + y * b + z * c)
+	c3 = sq x + sq y + sq z - 1
+    in
+	case quadratic c1 c2 c3 of
+        Nothing -> emptyIList
+        Just (t1, t2) -> entryexit (g t1) (g t2)
+    where x = xCoord r
+	  y = yCoord r
+	  z = zCoord r
+	  a = xComponent v
+	  b = yComponent v
+	  c = zComponent v
+	  g t = (t, (Spherical origin v1 v2, (SphereFace, p0, texture)))
+	      where origin = point 0 0 0
+		    x0 = x + t * a
+		    y0 = y + t * b
+		    z0 = z + t * c
+		    p0 = point  x0 y0 z0
+		    v0 = vector x0 y0 z0
+		    (v1, v2) = tangents v0
+
+
+--------------------------------
+-- Cube
+--------------------------------
+
+intersectCube :: Ray -> a -> IList (Surface, Texture a)
+intersectCube ray@(r, v) texture
+  = -- The set of t such that (x + at, y + bt, z + ct) lies within
+    -- the unit cube satisfies:
+    --    0 <= x + at <= 1,  0 <= y + bt <= 1,  0 <= z + ct <= 1
+    -- The minimum and maximum such values of t give us the two
+    -- intersection points.
+    case intersectSlabIval (intersectCubeSlab face2 face3 x a)
+	(intersectSlabIval (intersectCubeSlab face5 face4 y b)
+			   (intersectCubeSlab face0 face1 z c)) of
+    Nothing -> emptyIList
+    Just (t1, t2) -> entryexit (g t1) (g t2)
+  where g ((n, v0, v1), t)
+	  = (t, (Planar p0 v0 v1, (n, p0, texture)))
+	  where p0 = offsetToPoint ray t
+	face0 = (CubeFront,  vectorY, vectorX)
+	face1 = (CubeBack,   vectorX, vectorY)
+	face2 = (CubeLeft,   vectorZ, vectorY)
+	face3 = (CubeRight,  vectorY, vectorZ)
+	face4 = (CubeTop,    vectorZ, vectorX)
+	face5 = (CubeBottom, vectorX, vectorZ)
+	vectorX = vector 1 0 0
+	vectorY = vector 0 1 0
+	vectorZ = vector 0 0 1
+	x = xCoord r
+	y = yCoord r
+	z = zCoord r
+	a = xComponent v
+	b = yComponent v
+	c = zComponent v
+
+intersectCubeSlab n m w d
+  | d `near` 0 = if (0 <= w) && (w <= 1)
+		 then Just ((n, -inf), (m, inf)) else Nothing
+  | d > 0      = Just ((n,  (-w)/d), (m, (1-w)/d))
+  | otherwise  = Just ((m, (1-w)/d), (n,  (-w)/d))
+
+intersectSlabIval Nothing Nothing  = Nothing
+intersectSlabIval Nothing (Just i) = Nothing
+intersectSlabIval (Just i) Nothing = Nothing
+intersectSlabIval (Just (nu1@(n1, u1), mv1@(m1, v1)))
+		  (Just (nu2@(n2, u2), mv2@(m2, v2)))
+  = checkInterval (nu, mv)
+  where nu = if u1 < u2 then nu2 else nu1
+	mv = if v1 < v2 then mv1 else mv2
+	checkInterval numv@(nu@(_, u), (m, v))
+	  -- rounding error may force us to push v out a bit
+	  | u `near` v = Just (nu, (m, u + epsilon))
+	  | u    <   v = Just numv
+	  | otherwise  = Nothing
+
+
+--------------------------------
+-- Cylinder
+--------------------------------
+
+intersectCylinder :: Ray -> a -> IList (Surface, Texture a)
+intersectCylinder ray texture
+  = isectSide `intersectIntervals` isectTop `intersectIntervals` isectBottom
+  where isectSide   = intersectCylSide ray texture
+	isectTop    = intersectXZPlane CylinderTop ray 1.0 texture
+	isectBottom = complementIntervals $ negateSurfaces $
+		      intersectXZPlane CylinderBottom ray 0.0 texture
+
+intersectCylSide (r, v) texture
+  = -- The ray (x + ta, y + tb, z + tc) intersects the sides of the
+    -- cylinder if:
+    --    (x + ta)^2 + (z + tc)^2 = 1  and 0 <= y + tb <= 1.
+    if (sq a + sq c) `near` 0
+    then -- The ray is parallel to the Y-axis, and does not intersect
+	 -- the cylinder sides.  It's either all in, or all out
+	if (sqxy `near` 1.0 || sqxy < 1.0) then openIList else emptyIList
+   else -- Find values of t that solve the quadratic equation
+	--   (a^2 + c^2)t^2 + 2(ax + cz)t + x^2 + z^2 - 1 = 0
+        let c1 = sq a + sq c
+            c2 = 2 * (x * a + z * c)
+            c3 = sq x + sq z - 1
+	in
+	case quadratic c1 c2 c3 of
+        Nothing -> emptyIList
+        Just (t1, t2) -> entryexit (g t1) (g t2)
+
+  where sqxy = sq x + sq y
+	g t = (t, (Cylindrical origin v1 v2, (CylinderSide, p0, texture)))
+	    where origin = point 0 0 0
+		  x0 = x + t * a
+		  y0 = y + t * b
+		  z0 = z + t * c
+		  p0 = point  x0 y0 z0
+		  v0 = vector x0 0 z0
+		  (v1, v2) = tangents v0
+
+	x = xCoord r
+	y = yCoord r
+	z = zCoord r
+	a = xComponent v
+	b = yComponent v
+	c = zComponent v
+
+
+-------------------
+-- Cone
+-------------------
+
+intersectCone :: Ray -> a -> IList (Surface, Texture a)
+intersectCone ray texture
+  = isectSide `intersectIntervals` isectTop `intersectIntervals` isectBottom
+  where isectSide   = intersectConeSide ray texture
+	isectTop    = intersectXZPlane ConeBase ray 1.0 texture
+	isectBottom = complementIntervals $ negateSurfaces $
+		      intersectXZPlane ConeBase ray 0.0 texture
+
+intersectConeSide (r, v) texture
+  = -- Find the points where the ray intersects the cond side.  At any points of
+    -- intersection, we must have:
+    --    (x + ta)^2 + (z + tc)^2 = (y + tb)^2
+    -- which is the following quadratic equation:
+    --    t^2(a^2-b^2+c^2) + 2t(xa-yb+cz) + (x^2-y^2+z^2) = 0
+    let c1 = sq a - sq b + sq c
+	c2 = 2 * (x * a - y * b + c * z)
+	c3 = sq x - sq y + sq z
+    in  case quadratic c1 c2 c3 of
+	Nothing -> emptyIList
+	Just (t1, t2) ->
+	    -- If either intersection strikes the middle, then the other
+	    -- can only be off by rounding error, so we make a tangent
+	    -- strike using the "good" value.
+	    -- If the intersections straddle the origin, then it's
+	    -- an exit/entry pair, otherwise it's an entry/exit pair.
+	    let y1 = y + t1 * b
+		y2 = y + t2 * b
+	    in  if y1 `near` 0                  then entryexit (g t1) (g t1)
+	        else if y2 `near` 0             then entryexit (g t2) (g t2)
+		else if (y1 < 0) `xor` (y2 < 0) then exitentry (g t1) (g t2)
+		else                                 entryexit (g t1) (g t2)
+
+  where g t = (t, (Conic origin v1 v2, (ConeSide, p0, texture)))
+	    where origin = point 0 0 0
+		  x0 = x + t * a
+		  y0 = y + t * b
+		  z0 = z + t * c
+		  p0 = point  x0 y0 z0
+		  v0 = normalize $ vector x0 (-y0) z0
+		  (v1, v2) = tangents v0
+
+	x = xCoord r
+	y = yCoord r
+	z = zCoord r
+	a = xComponent v
+	b = yComponent v
+	c = zComponent v
+
+	-- beyond me why this isn't defined in the prelude...
+	xor False b = b
+	xor True  b = not b
+
+
+-------------------
+-- Solving quadratics
+-------------------
+
+quadratic :: Double -> Double -> Double -> Maybe (Double, Double)
+quadratic a b c =
+  -- Solve the equation ax^2 + bx + c = 0 by using the quadratic formula.
+  let d = sq b - 4 * a * c
+      d' = if d `near` 0 then 0 else d
+  in if d' < 0
+     then Nothing -- There are no real roots.
+     else
+	if a > 0 then Just (((-b) - sqrt d') / (2 * a),
+			    ((-b) + sqrt d') / (2 * a))
+		 else Just (((-b) + sqrt d') / (2 * a),
+			    ((-b) - sqrt d') / (2 * a))
+
+-------------------
+-- Bounding boxes
+-------------------
+
+data MaybeInterval = Interval !Double !Double 
+		   | NoInterval
+
+isInterval (Interval _ _) = True
+isInterval _              = False
+
+intersectWithBox :: Ray -> Box -> Bool
+intersectWithBox (r, v) (B x1 x2 y1 y2 z1 z2)
+  = isInterval interval
+  where x_interval = intersectRayWithSlab (xCoord r) (xComponent v) (x1, x2)
+	y_interval = intersectRayWithSlab (yCoord r) (yComponent v) (y1, y2)
+	z_interval = intersectRayWithSlab (zCoord r) (zComponent v) (z1, z2)
+	interval = intersectInterval x_interval
+		   (intersectInterval y_interval z_interval)
+
+intersectInterval :: MaybeInterval -> MaybeInterval -> MaybeInterval
+intersectInterval NoInterval _ = NoInterval
+intersectInterval _ NoInterval = NoInterval
+intersectInterval (Interval a b) (Interval c d)
+  | b < c || d < a = NoInterval
+  | otherwise = Interval (a `max` c) (b `min` d)
+
+{-# INLINE intersectRayWithSlab #-}
+intersectRayWithSlab :: Double -> Double -> (Double,Double) -> MaybeInterval
+intersectRayWithSlab xCoord alpha (x1, x2)
+  | alpha == 0 = if xCoord < x1 || xCoord > x2 then NoInterval else infInterval
+  | alpha >  0 = Interval a b
+  | otherwise  = Interval b a 
+  where a = (x1 - xCoord) / alpha
+	b = (x2 - xCoord) / alpha
+
+infInterval = Interval (-inf) inf