Addding the Galois Connections ray tracer as an example for GHC.
--- /dev/null
+-- 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 CSG(module Construct,
+ module Geometry,
+ module Intersections,
+ module Interval,
+ module Misc) where
+
+import Construct
+import Geometry
+import Intersections
+import Interval
+import Misc
--- /dev/null
+-- 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 Construct
+ ( Surface (..)
+ , Face (..)
+ , CSG (..)
+ , Texture
+ , Transform
+ , union, intersect, difference
+ , plane, sphere, cube, cylinder, cone
+ , transform
+ , translate, translateX, translateY, translateZ
+ , scale, scaleX, scaleY, scaleZ, uscale
+ , rotateX, rotateY, rotateZ
+ , eye, translateEye
+ , rotateEyeX, rotateEyeY, rotateEyeZ
+ ) where
+
+import Geometry
+
+-- In each case, we model the surface by a point and a pair of tangent vectors.
+-- This gives us enough information to determine the surface
+-- normal at that point, which is all that is required by the current
+-- illumination model. We can't just save the surface normal because
+-- that isn't preserved by transformations.
+
+data Surface
+ = Planar Point Vector Vector
+ | Spherical Point Vector Vector
+ | Cylindrical Point Vector Vector
+ | Conic Point Vector Vector
+ deriving Show
+
+data Face
+ = PlaneFace
+ | SphereFace
+ | CubeFront
+ | CubeBack
+ | CubeLeft
+ | CubeRight
+ | CubeTop
+ | CubeBottom
+ | CylinderSide
+ | CylinderTop
+ | CylinderBottom
+ | ConeSide
+ | ConeBase
+ deriving Show
+
+data CSG a
+ = Plane a
+ | Sphere a
+ | Cylinder a
+ | Cube a
+ | Cone a
+ | Transform Matrix Matrix (CSG a)
+ | Union (CSG a) (CSG a)
+ | Intersect (CSG a) (CSG a)
+ | Difference (CSG a) (CSG a)
+ | Box Box (CSG a)
+ deriving (Show)
+
+-- the data returned for determining surface texture
+-- the Face tells which face of a primitive this is
+-- the Point is the point of intersection in object coordinates
+-- the a is application-specific texture information
+type Texture a = (Face, Point, a)
+
+union, intersect, difference :: CSG a -> CSG a -> CSG a
+
+union p@(Box b1 _) q@(Box b2 _) = Box (mergeBox b1 b2) (Union p q)
+union p q = Union p q
+
+-- rather pessimistic
+intersect p@(Box b1 _) q@(Box b2 _) = Box (mergeBox b1 b2) (Intersect p q)
+intersect p q = Intersect p q
+
+difference (Box b1 p) q = Box b1 (Difference p q)
+-- no need to box again inside
+-- difference p@(Box b1 _) q = Box b1 (Difference p q)
+difference p q = Difference p q
+
+mkBox b p = Box b p
+
+plane, sphere, cube, cylinder, cone :: a -> CSG a
+
+plane = Plane
+sphere s =
+ mkBox (B (-1 - epsilon) (1 + epsilon)
+ (-1 - epsilon) (1 + epsilon)
+ (-1 - epsilon) (1 + epsilon)) (Sphere s)
+cone s =
+ mkBox (B (-1 - epsilon) (1 + epsilon)
+ ( - epsilon) (1 + epsilon)
+ (-1 - epsilon) (1 + epsilon)) (Cone s)
+cube s =
+ mkBox (B (- epsilon) (1 + epsilon)
+ (- epsilon) (1 + epsilon)
+ (- epsilon) (1 + epsilon)) (Cube s)
+cylinder s =
+ mkBox (B (-1 - epsilon) (1 + epsilon)
+ ( - epsilon) (1 + epsilon)
+ (-1 - epsilon) (1 + epsilon)) (Cylinder s)
+
+----------------------------
+-- Object transformations
+----------------------------
+
+type Transform = (Matrix, Matrix)
+
+transform :: Transform -> CSG a -> CSG a
+
+transform (m, m') (Transform mp mp' p) = Transform (multMM m mp) (multMM mp' m') p
+transform mm' (Union p q) = Union (transform mm' p) (transform mm' q)
+transform mm' (Intersect p q) = Intersect (transform mm' p) (transform mm' q)
+transform mm' (Difference p q) = Difference (transform mm' p) (transform mm' q)
+transform mm'@(m,_) (Box box p) = Box (transformBox m box) (transform mm' p)
+transform (m, m') prim = Transform m m' prim
+
+translate :: Coords -> CSG a -> CSG a
+translateX, translateY, translateZ :: Double -> CSG a -> CSG a
+
+translate xyz = transform $ transM xyz
+translateX x = translate (x, 0, 0)
+translateY y = translate (0, y, 0)
+translateZ z = translate (0, 0, z)
+
+scale :: Coords -> CSG a -> CSG a
+scaleX, scaleY, scaleZ, uscale :: Double -> CSG a -> CSG a
+
+scale xyz = transform $ scaleM xyz
+scaleX x = scale (x, 1, 1)
+scaleY y = scale (1, y, 1)
+scaleZ z = scale (1, 1, z)
+uscale u = scale (u,u,u)
+
+rotateX, rotateY, rotateZ :: Radian -> CSG a -> CSG a
+
+rotateX a = transform $ rotxM a
+rotateY a = transform $ rotyM a
+rotateZ a = transform $ rotzM a
+
+unit = matrix
+ ( ( 1.0, 0.0, 0.0, 0.0 ),
+ ( 0.0, 1.0, 0.0, 0.0 ),
+ ( 0.0, 0.0, 1.0, 0.0 ),
+ ( 0.0, 0.0, 0.0, 1.0 ) )
+
+transM (x, y, z)
+ = ( matrix
+ ( ( 1, 0, 0, x ),
+ ( 0, 1, 0, y ),
+ ( 0, 0, 1, z ),
+ ( 0, 0, 0, 1 ) ),
+ matrix
+ ( ( 1, 0, 0, -x ),
+ ( 0, 1, 0, -y ),
+ ( 0, 0, 1, -z ),
+ ( 0, 0, 0, 1 ) ) )
+
+scaleM (x, y, z)
+ = ( matrix
+ ( ( x', 0, 0, 0 ),
+ ( 0, y', 0, 0 ),
+ ( 0, 0, z', 0 ),
+ ( 0, 0, 0, 1 ) ),
+ matrix
+ ( ( 1/x', 0, 0, 0 ),
+ ( 0, 1/y', 0, 0 ),
+ ( 0, 0, 1/z', 0 ),
+ ( 0, 0, 0, 1 ) ) )
+ where x' = nonZero x
+ y' = nonZero y
+ z' = nonZero z
+
+rotxM t
+ = ( matrix
+ ( ( 1, 0, 0, 0 ),
+ ( 0, cos t, -sin t, 0 ),
+ ( 0, sin t, cos t, 0 ),
+ ( 0, 0, 0, 1 ) ),
+ matrix
+ ( ( 1, 0, 0, 0 ),
+ ( 0, cos t, sin t, 0 ),
+ ( 0, -sin t, cos t, 0 ),
+ ( 0, 0, 0, 1 ) ) )
+
+rotyM t
+ = ( matrix
+ ( ( cos t, 0, sin t, 0 ),
+ ( 0, 1, 0, 0 ),
+ ( -sin t, 0, cos t, 0 ),
+ ( 0, 0, 0, 1 ) ),
+ matrix
+ ( ( cos t, 0, -sin t, 0 ),
+ ( 0, 1, 0, 0 ),
+ ( sin t, 0, cos t, 0 ),
+ ( 0, 0, 0, 1 ) ) )
+
+rotzM t
+ = ( matrix
+ ( ( cos t, -sin t, 0, 0 ),
+ ( sin t, cos t, 0, 0 ),
+ ( 0, 0, 1, 0 ),
+ ( 0, 0, 0, 1 ) ),
+ matrix
+ ( ( cos t, sin t, 0, 0 ),
+ ( -sin t, cos t, 0, 0 ),
+ ( 0, 0, 1, 0 ),
+ ( 0, 0, 0, 1 ) ) )
+
+-------------------
+-- Eye transformations
+
+-- These are used to specify placement of the eye.
+-- `eye' starts out at (0, 0, -1).
+-- These are implemented as inverse transforms of the model.
+-------------------
+
+eye :: Transform
+translateEye :: Coords -> Transform -> Transform
+rotateEyeX, rotateEyeY, rotateEyeZ :: Radian -> Transform -> Transform
+
+eye = (unit, unit)
+translateEye xyz (eye1, eye2)
+ = (multMM m1 eye1, multMM eye2 m2)
+ where (m1, m2) = transM xyz
+rotateEyeX t (eye1, eye2)
+ = (multMM m1 eye1, multMM eye2 m2)
+ where (m1, m2) = rotxM t
+rotateEyeY t (eye1, eye2)
+ = (multMM m1 eye1, multMM eye2 m2)
+ where (m1, m2) = rotyM t
+rotateEyeZ t (eye1, eye2)
+ = (multMM m1 eye1, multMM eye2 m2)
+ where (m1, m2) = rotzM t
+
+-------------------
+-- Bounding boxes
+-------------------
+
+mergeBox (B x11 x12 y11 y12 z11 z12) (B x21 x22 y21 y22 z21 z22) =
+ B (x11 `min` x21) (x12 `max` x22)
+ (y11 `min` y21) (y12 `max` y22)
+ (z11 `min` z21) (z12 `max` z22)
+
+transformBox t (B x1 x2 y1 y2 z1 z2)
+ = (B (foldr1 min (map xCoord pts'))
+ (foldr1 max (map xCoord pts'))
+ (foldr1 min (map yCoord pts'))
+ (foldr1 max (map yCoord pts'))
+ (foldr1 min (map zCoord pts'))
+ (foldr1 max (map zCoord pts')))
+ where pts' = map (multMP t) pts
+ pts = [point x1 y1 z1,
+ point x1 y1 z2,
+ point x1 y2 z1,
+ point x1 y2 z2,
+ point x2 y1 z1,
+ point x2 y1 z2,
+ point x2 y2 z1,
+ point x2 y2 z2]
--- /dev/null
+-- 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 Data where
+
+import Array
+import IOExts
+
+import CSG
+import Geometry
+import Illumination
+import Primitives
+import Surface
+
+-- Now the parsed (expresssion) language
+
+type Name = String
+
+type Code = [GMLToken]
+
+data GMLToken
+ -- All these can occur in parsed code
+ = TOp GMLOp
+ | TId Name
+ | TBind Name
+ | TBool Bool
+ | TInt Int
+ | TReal Double
+ | TString String
+ | TBody Code
+ | TArray Code
+ | TApply
+ | TIf
+ -- These can occur in optimized/transformed code
+ -- NONE (yet!)
+
+
+instance Show GMLToken where
+ showsPrec p (TOp op) = shows op
+ showsPrec p (TId id) = showString id
+ showsPrec p (TBind id) = showString ('/' : id)
+ showsPrec p (TBool bool) = shows bool
+ showsPrec p (TInt i) = shows i
+ showsPrec p (TReal d) = shows d
+ showsPrec p (TString s) = shows s
+ showsPrec p (TBody code) = shows code
+ showsPrec p (TArray code) = showString "[ "
+ . foldr (\ a b -> a . showChar ' ' . b) id (map shows code)
+ . showString "]"
+ showsPrec p (TApply) = showString "apply"
+ showsPrec p (TIf) = showString "if"
+
+ showList code = showString "{ "
+ . foldr (\ a b -> a . showChar ' ' . b) id (map shows code)
+ . showString "}"
+
+
+-- Now the value language, used inside the interpreter
+
+type Stack = [GMLValue]
+
+data GMLValue
+ = VBool !Bool
+ | VInt !Int
+ | VReal !Double
+ | VString String
+ | VClosure Env Code
+ | VArray (Array Int GMLValue) -- FIXME: Haskell array
+ -- uses the interpreter version of point
+ | VPoint { xPoint :: !Double
+ , yPoint :: !Double
+ , zPoint :: !Double
+ }
+ -- these are abstract to the interpreter
+ | VObject Object
+ | VLight Light
+ -- This is an abstract object, used by the abstract interpreter
+ | VAbsObj AbsObj
+
+
+-- There are only *3* basic abstract values,
+-- and the combinators also.
+
+data AbsObj
+ = AbsFACE
+ | AbsU
+ | AbsV
+ deriving (Show)
+
+instance Show GMLValue where
+ showsPrec p value = showString (showStkEle value)
+
+showStkEle :: GMLValue -> String
+showStkEle (VBool b) = show b ++ " :: Bool"
+showStkEle (VInt i) = show i ++ " :: Int"
+showStkEle (VReal r) = show r ++ " :: Real"
+showStkEle (VString s) = show s ++ " :: String"
+showStkEle (VClosure {}) = "<closure> :: Closure"
+showStkEle (VArray arr)
+ = "<array (" ++ show (succ (snd (bounds arr))) ++ " elements)> :: Array"
+showStkEle (VPoint x y z) = "(" ++ show x
+ ++ "," ++ show y
+ ++ "," ++ show z
+ ++ ") :: Point"
+showStkEle (VObject {}) = "<Object> :: Object"
+showStkEle (VLight {}) = "<Light> :: Object"
+showStkEle (VAbsObj vobs) = "{{ " ++ show vobs ++ "}} :: AbsObj"
+
+-- An abstract environment
+
+newtype Env = Env [(Name, GMLValue)] deriving Show
+
+emptyEnv :: Env
+emptyEnv = Env []
+
+extendEnv :: Env -> Name -> GMLValue -> Env
+extendEnv (Env e) n v = Env ((n, v):e)
+
+lookupEnv :: Env -> Name -> Maybe GMLValue
+lookupEnv (Env e) n = lookup n e
+
+-- All primitive operators
+--
+-- There is no Op_apply, Op_false, Op_true and Op_if
+-- (because they appear explcitly in the rules).
+
+data GMLOp
+ = Op_acos
+ | Op_addi
+ | Op_addf
+ | Op_asin
+ | Op_clampf
+ | Op_cone
+ | Op_cos
+ | Op_cube
+ | Op_cylinder
+ | Op_difference
+ | Op_divi
+ | Op_divf
+ | Op_eqi
+ | Op_eqf
+ | Op_floor
+ | Op_frac
+ | Op_get
+ | Op_getx
+ | Op_gety
+ | Op_getz
+ | Op_intersect
+ | Op_length
+ | Op_lessi
+ | Op_lessf
+ | Op_light
+ | Op_modi
+ | Op_muli
+ | Op_mulf
+ | Op_negi
+ | Op_negf
+ | Op_plane
+ | Op_point
+ | Op_pointlight
+ | Op_real
+ | Op_render
+ | Op_rotatex
+ | Op_rotatey
+ | Op_rotatez
+ | Op_scale
+ | Op_sin
+ | Op_sphere
+ | Op_spotlight
+ | Op_sqrt
+ | Op_subi
+ | Op_subf
+ | Op_trace -- non standard, for debugging GML programs
+ | Op_translate
+ | Op_union
+ | Op_uscale
+ deriving (Eq,Ord,Ix,Bounded)
+
+instance Show GMLOp where
+ showsPrec _ op = showString (opNameTable ! op)
+
+
+------------------------------------------------------------------------------
+
+-- And how we use the op codes (there names, there interface)
+
+-- These keywords include, "apply", "if", "true" and "false",
+-- they are not parsed as operators, but are
+-- captured by the parser as a special case.
+
+keyWords :: [String]
+keyWords = [ kwd | (kwd,_,_) <- opcodes ]
+
+-- Lookup has to look from the start (or else...)
+opTable :: [(Name,GMLToken)]
+opTable = [ (kwd,op) | (kwd,op,_) <- opcodes ]
+
+opNameTable :: Array GMLOp Name
+opNameTable = array (minBound,maxBound)
+ [ (op,name) | (name,TOp op,_) <- opcodes ]
+
+undef = error "undefined function"
+image = error "undefined function: talk to image group"
+
+-- typically, its best to have *one* opcode table,
+-- so that mis-alignments do not happen.
+
+opcodes :: [(String,GMLToken,PrimOp)]
+opcodes =
+ [ ("apply", TApply, error "incorrect use of apply")
+ , ("if", TIf, error "incorrect use of if")
+ , ("false", TBool False, error "incorrect use of false")
+ , ("true", TBool True, error "incorrect use of true")
+ ] ++ map (\ (a,b,c) -> (a,TOp b,c))
+ -- These are just invocation, any coersions need to occur between here
+ -- and before arriving at the application code (like deg -> rad).
+ [ ("acos", Op_acos, Real_Real (rad2deg . acos))
+ , ("addi", Op_addi, Int_Int_Int (+))
+ , ("addf", Op_addf, Real_Real_Real (+))
+ , ("asin", Op_asin, Real_Real (rad2deg . asin))
+ , ("clampf", Op_clampf, Real_Real clampf)
+ , ("cone", Op_cone, Surface_Obj cone)
+ , ("cos", Op_cos, Real_Real (cos . deg2rad))
+ , ("cube", Op_cube, Surface_Obj cube)
+ , ("cylinder", Op_cylinder, Surface_Obj cylinder)
+ , ("difference", Op_difference, Obj_Obj_Obj difference)
+ , ("divi", Op_divi, Int_Int_Int (ourQuot))
+ , ("divf", Op_divf, Real_Real_Real (/))
+ , ("eqi", Op_eqi, Int_Int_Bool (==))
+ , ("eqf", Op_eqf, Real_Real_Bool (==))
+ , ("floor", Op_floor, Real_Int floor)
+ , ("frac", Op_frac, Real_Real (snd . properFraction))
+ , ("get", Op_get, Arr_Int_Value ixGet)
+ , ("getx", Op_getx, Point_Real (\ x y z -> x))
+ , ("gety", Op_gety, Point_Real (\ x y z -> y))
+ , ("getz", Op_getz, Point_Real (\ x y z -> z))
+ , ("intersect", Op_intersect, Obj_Obj_Obj intersect)
+ , ("length", Op_length, Arr_Int (succ . snd . bounds))
+ , ("lessi", Op_lessi, Int_Int_Bool (<))
+ , ("lessf", Op_lessf, Real_Real_Bool (<))
+ , ("light", Op_light, Point_Color_Light light)
+ , ("modi", Op_modi, Int_Int_Int (ourRem))
+ , ("muli", Op_muli, Int_Int_Int (*))
+ , ("mulf", Op_mulf, Real_Real_Real (*))
+ , ("negi", Op_negi, Int_Int negate)
+ , ("negf", Op_negf, Real_Real negate)
+ , ("plane", Op_plane, Surface_Obj plane)
+ , ("point", Op_point, Real_Real_Real_Point VPoint)
+ , ("pointlight", Op_pointlight, Point_Color_Light pointlight)
+ , ("real", Op_real, Int_Real fromIntegral)
+ , ("render", Op_render, Render $ render eye)
+ , ("rotatex", Op_rotatex, Obj_Real_Obj (\ o d -> rotateX (deg2rad d) o))
+ , ("rotatey", Op_rotatey, Obj_Real_Obj (\ o d -> rotateY (deg2rad d) o))
+ , ("rotatez", Op_rotatez, Obj_Real_Obj (\ o d -> rotateZ (deg2rad d) o))
+ , ("scale", Op_scale, Obj_Real_Real_Real_Obj (\ o x y z -> scale (x,y,z) o))
+ , ("sin", Op_sin, Real_Real (sin . deg2rad))
+ , ("sphere", Op_sphere, Surface_Obj sphere') -- see comment at end of file
+ , ("spotlight", Op_spotlight, Point_Point_Color_Real_Real_Light mySpotlight)
+ , ("sqrt", Op_sqrt, Real_Real ourSqrt)
+ , ("subi", Op_subi, Int_Int_Int (-))
+ , ("subf", Op_subf, Real_Real_Real (-))
+ , ("trace", Op_trace, Value_String_Value mytrace)
+ , ("translate", Op_translate, Obj_Real_Real_Real_Obj (\ o x y z -> translate (x,y,z) o))
+ , ("union", Op_union, Obj_Obj_Obj union)
+ , ("uscale", Op_uscale, Obj_Real_Obj (\ o r -> uscale r o))
+ ]
+
+-- This enumerate all possible ways of calling the fixed primitives
+
+-- The datatype captures the type at the *interp* level,
+-- the type of the functional is mirrored on this (using Haskell types).
+
+data PrimOp
+
+ -- 1 argument
+ = Int_Int (Int -> Int)
+ | Real_Real (Double -> Double)
+ | Point_Real (Double -> Double -> Double -> Double)
+ | Surface_Obj (SurfaceFn Color Double -> Object)
+ | Real_Int (Double -> Int)
+ | Int_Real (Int -> Double)
+ | Arr_Int (Array Int GMLValue -> Int)
+
+ -- 2 arguments
+ | Int_Int_Int (Int -> Int -> Int)
+ | Int_Int_Bool (Int -> Int -> Bool)
+ | Real_Real_Real (Double -> Double -> Double)
+ | Real_Real_Bool (Double -> Double -> Bool)
+ | Arr_Int_Value (Array Int GMLValue -> Int -> GMLValue)
+
+ -- Many arguments, typically image mangling
+
+ | Obj_Obj_Obj (Object -> Object -> Object)
+ | Point_Color_Light (Coords -> Color -> Light)
+ | Real_Real_Real_Point (Double -> Double -> Double -> GMLValue)
+ | Obj_Real_Obj (Object -> Double -> Object)
+ | Obj_Real_Real_Real_Obj (Object -> Double -> Double -> Double -> Object)
+ | Value_String_Value (GMLValue -> String -> GMLValue)
+
+ | Point_Point_Color_Real_Real_Light
+ (Coords -> Coords -> Color -> Radian -> Radian -> Light)
+ -- And finally render
+ | Render (Color -> [Light] -> Object -> Int -> Double -> Int -> Int -> String -> IO ())
+
+data Type
+ = TyBool
+ | TyInt
+ | TyReal
+ | TyString
+ | TyCode
+ | TyArray
+ | TyPoint
+ | TyObject
+ | TyLight
+ | TyAlpha
+ | TyAbsObj
+ deriving (Eq,Ord,Ix,Bounded)
+
+typeTable =
+ [ ( TyBool, "Bool")
+ , ( TyInt, "Int")
+ , ( TyReal, "Real")
+ , ( TyString, "String")
+ , ( TyCode, "Code")
+ , ( TyArray, "Array")
+ , ( TyPoint, "Point")
+ , ( TyObject, "Object")
+ , ( TyLight, "Light")
+ , ( TyAlpha, "<anything>")
+ , ( TyAbsObj, "<abs>")
+ ]
+
+typeNames = array (minBound,maxBound) typeTable
+
+instance Show Type where
+ showsPrec _ op = showString (typeNames ! op)
+
+getPrimOpType :: PrimOp -> [Type]
+getPrimOpType (Int_Int _) = [TyInt]
+getPrimOpType (Real_Real _) = [TyReal]
+getPrimOpType (Point_Real _) = [TyPoint]
+getPrimOpType (Surface_Obj _) = [TyCode]
+getPrimOpType (Real_Int _) = [TyReal]
+getPrimOpType (Int_Real _) = [TyInt]
+getPrimOpType (Arr_Int _) = [TyArray]
+getPrimOpType (Int_Int_Int _) = [TyInt,TyInt]
+getPrimOpType (Int_Int_Bool _) = [TyInt,TyInt]
+getPrimOpType (Real_Real_Real _) = [TyReal,TyReal]
+getPrimOpType (Real_Real_Bool _) = [TyReal,TyReal]
+getPrimOpType (Arr_Int_Value _) = [TyArray,TyInt]
+getPrimOpType (Obj_Obj_Obj _) = [TyObject,TyObject]
+getPrimOpType (Point_Color_Light _) = [TyPoint,TyPoint]
+getPrimOpType (Real_Real_Real_Point _) = [TyReal,TyReal,TyReal]
+getPrimOpType (Obj_Real_Obj _) = [TyObject,TyReal]
+getPrimOpType (Obj_Real_Real_Real_Obj _) = [TyObject,TyReal,TyReal,TyReal]
+getPrimOpType (Value_String_Value _) = [TyAlpha,TyString]
+getPrimOpType (Point_Point_Color_Real_Real_Light _)
+ = [TyPoint,TyPoint,TyPoint,TyReal,TyReal]
+getPrimOpType (Render _) = [TyPoint,
+ TyLight,
+ TyObject,
+ TyInt,
+ TyReal,
+ TyReal,
+ TyReal,
+ TyString]
+
+
+-- Some primitives with better error message
+
+mytrace v s = trace (s ++" : "++ show v ++ "\n") v
+
+
+ixGet :: Array Int GMLValue -> Int -> GMLValue
+ixGet arr i
+ | inRange (bounds arr) i = arr ! i
+ | otherwise = error ("failed access with index value "
+ ++ show i
+ ++ " (should be between 0 and "
+ ++ show (snd (bounds arr)) ++ ")")
+
+ourQuot :: Int -> Int -> Int
+ourQuot _ 0 = error "attempt to use divi to divide by 0"
+ourQuot a b = a `quot` b
+
+ourRem :: Int -> Int -> Int
+ourRem _ 0 = error "attempt to use remi to divide by 0"
+ourRem a b = a `rem` b
+
+ourSqrt :: Double -> Double
+ourSqrt n | n < 0 = error "attempt to use sqrt on a negative number"
+ | otherwise = sqrt n
+
+
+mySpotlight p1 p2 col cutoff exp = spotlight p1 p2 col (deg2rad cutoff) exp
+
+-- The problem specification gets the mapping for spheres backwards
+-- (it maps the image from right to left).
+-- We've fixed that in the raytracing library so that it goes from left
+-- to right, but to keep the GML front compatible with the problem
+-- statement, we reverse it here.
+
+sphere' :: SurfaceFn Color Double -> CSG (SurfaceFn Color Double)
+sphere' (SFun f) = sphere (SFun (\i u v -> f i (1 - u) v))
+sphere' s = sphere s
--- /dev/null
+-- 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 Eval where
+
+import Array
+
+import IOExts
+
+import Geometry
+import CSG
+import Surface
+import Data
+import Parse (rayParse, rayParseF)
+
+class Monad m => MonadEval m where
+ doOp :: PrimOp -> GMLOp -> Stack -> m Stack
+ tick :: m ()
+ err :: String -> m a
+
+ tick = return ()
+
+newtype Pure a = Pure a deriving Show
+
+instance Monad Pure where
+ Pure x >>= k = k x
+ return = Pure
+ fail s = error s
+
+instance MonadEval Pure where
+ doOp = doPureOp
+ err s = error s
+
+instance MonadEval IO where
+ doOp prim op stk = do { -- putStrLn ("Calling " ++ show op
+ -- ++ " << " ++ show stk ++ " >>")
+ doAllOp prim op stk
+ }
+ err s = error s
+
+data State
+ = State { env :: Env
+ , stack :: Stack
+ , code :: Code
+ } deriving Show
+
+callback :: Env -> Code -> Stack -> Stack
+callback env code stk
+ = case eval (State { env = env, stack = stk, code = code}) of
+ Pure stk -> stk
+
+{-# SPECIALIZE eval :: State -> Pure Stack #-}
+{-# SPECIALIZE eval :: State -> IO Stack #-}
+
+eval :: MonadEval m => State -> m Stack
+eval st =
+ do { () <- return () -- $ unsafePerformIO (print st) -- Functional debugger
+ ; if moreCode st then
+ do { tick -- tick first, so as to catch loops on new eval.
+ ; st' <- step st
+ ; eval st'
+ }
+ else return (stack st)
+ }
+
+moreCode :: State -> Bool
+moreCode (State {code = []}) = False
+moreCode _ = True
+
+-- Step has a precondition that there *is* code to run
+{-# SPECIALIZE step :: State -> Pure State #-}
+{-# SPECIALIZE step :: State -> IO State #-}
+step :: MonadEval m => State -> m State
+
+-- Rule 1: Pushing BaseValues
+step st@(State{ stack = stack, code = (TBool b):cs })
+ = return (st { stack = (VBool b):stack, code = cs })
+step st@(State{ stack = stack, code = (TInt i):cs })
+ = return (st { stack = (VInt i):stack, code = cs })
+step st@(State{ stack = stack, code = (TReal r):cs })
+ = return (st { stack = (VReal r):stack, code = cs })
+step st@(State{ stack = stack, code = (TString s):cs })
+ = return (st { stack = (VString s):stack, code = cs })
+
+-- Rule 2: Name binding
+step st@(State{ env = env, stack = (v:stack), code = (TBind id):cs }) =
+ return (State { env = extendEnv env id v, stack = stack, code = cs })
+step st@(State{ env = env, stack = [], code = (TBind id):cs }) =
+ err "Attempt to bind the top of an empty stack"
+
+-- Rule 3: Name lookup
+step st@(State{ env = env, stack = stack, code = (TId id):cs }) =
+ case (lookupEnv env id) of
+ Just v -> return (st { stack = v:stack, code = cs })
+ Nothing -> err ("Cannot find value for identifier: " ++ id)
+
+-- Rule 4: Closure creation
+step st@(State{ env = env, stack = stack, code = (TBody body):cs }) =
+ return (st { stack = (VClosure env body):stack, code = cs })
+
+-- Rule 5: Application
+step st@(State{ env = env, stack = (VClosure env' code'):stack, code = TApply:cs }) =
+ do { stk <- eval (State {env = env', stack = stack, code = code'})
+ ; return (st { stack = stk, code = cs })
+ }
+step st@(State{ env = env, stack = [], code = TApply:cs }) =
+ err "Application with an empty stack"
+step st@(State{ env = env, stack = _:_, code = TApply:cs }) =
+ err "Application of a non-closure"
+
+-- Rule 6: Arrays
+step st@(State{ env = env, stack = stack, code = TArray code':cs }) =
+ do { stk <- eval (State {env = env, stack = [], code = code'})
+ ; let last = length stk-1
+ ; let arr = array (0,last) (zip [last,last-1..] stk)
+ ; return (st { stack = (VArray arr):stack, code = cs })
+ }
+
+-- Rule 7 & 8: If statement
+step st@(State{ env = env, stack = (VClosure e2 c2):(VClosure e1 c1):(VBool True):stack, code = TIf:cs }) =
+ do { stk <- eval (State {env = e1, stack = stack, code = c1})
+ ; return (st { stack = stk, code = cs })
+ }
+step st@(State{ env = env, stack = (VClosure e2 c2):(VClosure e1 c1):(VBool False):stack, code = TIf:cs }) =
+ do { stk <- eval (State {env = e2, stack = stack, code = c2})
+ ; return (st { stack = stk, code = cs })
+ }
+step st@(State{ env = env, stack = _, code = TIf:cs }) =
+ err "Incorrect use of if (bad and/or inappropriate values on the stack)"
+
+-- Rule 9: Operators
+step st@(State{ env = env, stack = stack, code = (TOp op):cs }) =
+ do { stk <- doOp (opFnTable ! op) op stack
+ ; return (st { stack = stk, code = cs })
+ }
+
+-- Rule Opps
+step _ = err "Tripped on sidewalk while stepping."
+
+
+--------------------------------------------------------------------------
+-- Operator code
+
+opFnTable :: Array GMLOp PrimOp
+opFnTable = array (minBound,maxBound)
+ [ (op,prim) | (_,TOp op,prim) <- opcodes ]
+
+
+
+
+doPureOp :: (MonadEval m) => PrimOp -> GMLOp -> Stack -> m Stack
+doPureOp _ Op_render _ =
+ err ("\nAttempting to call render from inside a purely functional callback.")
+doPureOp primOp op stk = doPrimOp primOp op stk -- call the purely functional operators
+
+{-# SPECIALIZE doPrimOp :: PrimOp -> GMLOp -> Stack -> Pure Stack #-}
+{-# SPECIALIZE doPrimOp :: PrimOp -> GMLOp -> Stack -> IO Stack #-}
+{-# SPECIALIZE doPrimOp :: PrimOp -> GMLOp -> Stack -> Abs Stack #-}
+
+doPrimOp :: (MonadEval m) => PrimOp -> GMLOp -> Stack -> m Stack
+
+-- 1 argument.
+
+doPrimOp (Int_Int fn) _ (VInt i1:stk)
+ = return ((VInt (fn i1)) : stk)
+doPrimOp (Real_Real fn) _ (VReal r1:stk)
+ = return ((VReal (fn r1)) : stk)
+doPrimOp (Point_Real fn) _ (VPoint x y z:stk)
+ = return ((VReal (fn x y z)) : stk)
+
+-- This is where the callbacks happen from...
+doPrimOp (Surface_Obj fn) _ (VClosure env code:stk)
+ = case absapply env code [VAbsObj AbsFACE,VAbsObj AbsU,VAbsObj AbsV] of
+ Just [VReal r3,VReal r2,VReal r1,VPoint c1 c2 c3] ->
+ let
+ res = prop (color c1 c2 c3) r1 r2 r3
+ in
+ return ((VObject (fn (SConst res))) : stk)
+ _ -> return ((VObject (fn (SFun call))) : stk)
+ where
+ -- The most general case
+ call i r1 r2 =
+ case callback env code [VReal r2,VReal r1,VInt i] of
+ [VReal r3,VReal r2,VReal r1,VPoint c1 c2 c3]
+ -> prop (color c1 c2 c3) r1 r2 r3
+ stk -> error ("callback failed: incorrectly typed return arguments"
+ ++ show stk)
+
+doPrimOp (Real_Int fn) _ (VReal r1:stk)
+ = return ((VInt (fn r1)) : stk)
+doPrimOp (Int_Real fn) _ (VInt r1:stk)
+ = return ((VReal (fn r1)) : stk)
+doPrimOp (Arr_Int fn) _ (VArray arr:stk)
+ = return ((VInt (fn arr)) : stk)
+
+-- 2 arguments.
+
+doPrimOp (Int_Int_Int fn) _ (VInt i2:VInt i1:stk)
+ = return ((VInt (fn i1 i2)) : stk)
+doPrimOp (Int_Int_Bool fn) _ (VInt i2:VInt i1:stk)
+ = return ((VBool (fn i1 i2)) : stk)
+doPrimOp (Real_Real_Real fn) _ (VReal r2:VReal r1:stk)
+ = return ((VReal (fn r1 r2)) : stk)
+doPrimOp (Real_Real_Bool fn) _ (VReal r2:VReal r1:stk)
+ = return ((VBool (fn r1 r2)) : stk)
+doPrimOp (Arr_Int_Value fn) _ (VInt i:VArray arr:stk)
+ = return ((fn arr i) : stk)
+
+
+ -- Many arguments, typically image mangling
+
+doPrimOp (Obj_Obj_Obj fn) _ (VObject o2:VObject o1:stk)
+ = return ((VObject (fn o1 o2)) : stk)
+doPrimOp (Point_Color_Light fn) _ (VPoint r g b:VPoint x y z : stk)
+ = return (VLight (fn (x,y,z) (color r g b)) : stk)
+doPrimOp (Point_Point_Color_Real_Real_Light fn) _
+ (VReal r2:VReal r1:VPoint r g b:VPoint x2 y2 z2:VPoint x1 y1 z1 : stk)
+ = return (VLight (fn (x1,y1,z1) (x2,y2,z2) (color r g b) r1 r2) : stk)
+doPrimOp (Real_Real_Real_Point fn) _ (VReal r3:VReal r2:VReal r1:stk)
+ = return ((fn r1 r2 r3) : stk)
+doPrimOp (Obj_Real_Obj fn) _ (VReal r:VObject o:stk)
+ = return (VObject (fn o r) : stk)
+doPrimOp (Obj_Real_Real_Real_Obj fn) _ (VReal r3:VReal r2:VReal r1:VObject o:stk)
+ = return (VObject (fn o r1 r2 r3) : stk)
+
+-- This one is our testing harness
+doPrimOp (Value_String_Value fn) _ (VString s:o:stk)
+ = res `seq` return (res : stk)
+ where
+ res = fn o s
+
+doPrimOp primOp op args
+ = err ("\n\ntype error when attempting to execute builtin primitive \"" ++
+ show op ++ "\"\n\n| " ++
+ show op ++ " takes " ++ show (length types) ++ " argument" ++ s
+ ++ " with" ++ the ++ " type" ++ s ++ "\n|\n|" ++
+ " " ++ unwords [ show ty | ty <- types ] ++ "\n|\n|" ++
+ " currently, the relevent argument" ++ s ++ " on the stack " ++
+ are ++ "\n|\n| " ++
+ unwords [ "(" ++ show arg ++ ")"
+ | arg <- reverse (take (length types) args) ] ++ "\n|\n| "
+ ++ " (top of stack is on the right hand side)\n\n")
+ where
+ len = length types
+ s = (if len /= 1 then "s" else "")
+ are = (if len /= 1 then "are" else "is")
+ the = (if len /= 1 then "" else " the")
+ types = getPrimOpType primOp
+
+
+-- Render is somewhat funny, becauase it can only get called at top level.
+-- All other operations are purely functional.
+
+doAllOp :: PrimOp -> GMLOp -> Stack -> IO Stack
+doAllOp (Render render) Op_render
+ (VString str:VInt ht:VInt wid:VReal fov
+ :VInt dep:VObject obj:VArray arr
+ :VPoint r g b : stk)
+ = do { render (color r g b) lights obj dep (fov * (pi / 180.0)) wid ht str
+ ; return stk
+ }
+ where
+ lights = [ light | (VLight light) <- elems arr ]
+
+doAllOp primOp op stk = doPrimOp primOp op stk -- call the purely functional operators
+
+------------------------------------------------------------------------------
+{-
+ - Abstract evaluation.
+ -
+ - The idea is you check for constant code that
+ - (1) does not look at its arguments
+ - (2) gives a fixed result
+ -
+ - We run for 100 steps.
+ -
+ -}
+
+absapply :: Env -> Code -> Stack -> Maybe Stack
+absapply env code stk =
+ case runAbs (eval (State env stk code)) 100 of
+ AbsState stk _ -> Just stk
+ AbsFail m -> Nothing
+
+newtype Abs a = Abs { runAbs :: Int -> AbsState a }
+data AbsState a = AbsState a !Int
+ | AbsFail String
+
+instance Monad Abs where
+ (Abs fn) >>= k = Abs (\ s -> case fn s of
+ AbsState r s' -> runAbs (k r) s'
+ AbsFail m -> AbsFail m)
+ return x = Abs (\ n -> AbsState x n)
+ fail s = Abs (\ n -> AbsFail s)
+
+instance MonadEval Abs where
+ doOp = doAbsOp
+ err = fail
+ tick = Abs (\ n -> if n <= 0
+ then AbsFail "run out of time"
+ else AbsState () (n-1))
+
+doAbsOp :: PrimOp -> GMLOp -> Stack -> Abs Stack
+doAbsOp _ Op_point (VReal r3:VReal r2:VReal r1:stk)
+ = return ((VPoint r1 r2 r3) : stk)
+ -- here, you could have an (AbsPoint :: AbsObj) which you put on the
+ -- stack, with any object in the three fields.
+doAbsOp _ op _ = err ("operator not understood (" ++ show op ++ ")")
+
+------------------------------------------------------------------------------
+-- Driver
+
+mainEval :: Code -> IO ()
+mainEval prog = do { stk <- eval (State emptyEnv [] prog)
+ ; return ()
+ }
+{-
+ * Oops, one of the example actually has something
+ * on the stack at the end.
+ * Oh well...
+ ; if null stk
+ then return ()
+ else do { putStrLn done
+ ; print stk
+ }
+-}
+
+done = "Items still on stack at (successfull) termination of program"
+
+------------------------------------------------------------------------------
+-- testing
+
+test :: String -> Pure Stack
+test is = eval (State emptyEnv [] (rayParse is))
+
+testF :: String -> IO Stack
+testF is = do prog <- rayParseF is
+ eval (State emptyEnv [] prog)
+
+testA :: String -> Either String (Stack,Int)
+testA is = case runAbs (eval (State emptyEnv
+ [VAbsObj AbsFACE,VAbsObj AbsU,VAbsObj AbsV]
+ (rayParse is))) 100 of
+ AbsState a n -> Right (a,n)
+ AbsFail m -> Left m
+
+abstest1 = "1.0 0.0 0.0 point /red { /v /u /face red 1.0 0.0 1.0 } apply"
+
+-- should be [3:: Int]
+et1 = test "1 /x { x } /f 2 /x f apply x addi"
+
+
+
+
+
--- /dev/null
+-- 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 Geometry
+ ( Coords
+ , Ray
+ , Point -- abstract
+ , Vector -- abstract
+ , Matrix -- abstract
+ , Color -- abstract
+ , Box(..)
+ , Radian
+ , matrix
+ , coord
+ , color
+ , uncolor
+ , xCoord , yCoord , zCoord
+ , xComponent , yComponent , zComponent
+ , point
+ , vector
+ , nearV
+ , point_to_vector
+ , vector_to_point
+ , dot
+ , cross
+ , tangents
+ , addVV
+ , addPV
+ , subVV
+ , negV
+ , subPP
+ , norm
+ , normalize
+ , dist2
+ , sq
+ , distFrom0Sq
+ , distFrom0
+ , multSV
+ , multMM
+ , transposeM
+ , multMV
+ , multMP
+ , multMQ
+ , multMR
+ , white
+ , black
+ , addCC
+ , subCC
+ , sumCC
+ , multCC
+ , multSC
+ , nearC
+ , offsetToPoint
+ , epsilon
+ , inf
+ , nonZero
+ , eqEps
+ , near
+ , clampf
+ ) where
+
+import List
+
+type Coords = (Double,Double,Double)
+
+type Ray = (Point,Vector) -- origin of ray, and unit vector giving direction
+
+data Point = P !Double !Double !Double -- implicit extra arg of 1
+ deriving (Show)
+data Vector = V !Double !Double !Double -- implicit extra arg of 0
+ deriving (Show, Eq)
+data Matrix = M !Quad !Quad !Quad !Quad
+ deriving (Show)
+
+data Color = C !Double !Double !Double
+ deriving (Show, Eq)
+
+data Box = B !Double !Double !Double !Double !Double !Double
+ deriving (Show)
+
+data Quad = Q !Double !Double !Double !Double
+ deriving (Show)
+
+type Radian = Double
+
+type Tup4 a = (a,a,a,a)
+
+--{-# INLINE matrix #-}
+matrix :: Tup4 (Tup4 Double) -> Matrix
+matrix ((m11, m12, m13, m14),
+ (m21, m22, m23, m24),
+ (m31, m32, m33, m34),
+ (m41, m42, m43, m44))
+ = M (Q m11 m12 m13 m14)
+ (Q m21 m22 m23 m24)
+ (Q m31 m32 m33 m34)
+ (Q m41 m42 m43 m44)
+
+coord x y z = (x, y, z)
+
+color r g b = C r g b
+
+uncolor (C r g b) = (r,g,b)
+
+{-# INLINE xCoord #-}
+xCoord (P x y z) = x
+{-# INLINE yCoord #-}
+yCoord (P x y z) = y
+{-# INLINE zCoord #-}
+zCoord (P x y z) = z
+
+{-# INLINE xComponent #-}
+xComponent (V x y z) = x
+{-# INLINE yComponent #-}
+yComponent (V x y z) = y
+{-# INLINE zComponent #-}
+zComponent (V x y z) = z
+
+point :: Double -> Double -> Double -> Point
+point x y z = P x y z
+
+vector :: Double -> Double -> Double -> Vector
+vector x y z = V x y z
+
+nearV :: Vector -> Vector -> Bool
+nearV (V a b c) (V d e f) = a `near` d && b `near` e && c `near` f
+
+point_to_vector :: Point -> Vector
+point_to_vector (P x y z) = V x y z
+
+vector_to_point :: Vector -> Point
+vector_to_point (V x y z) = P x y z
+
+{-# INLINE vector_to_quad #-}
+vector_to_quad :: Vector -> Quad
+vector_to_quad (V x y z) = Q x y z 0
+
+{-# INLINE point_to_quad #-}
+point_to_quad :: Point -> Quad
+point_to_quad (P x y z) = Q x y z 1
+
+{-# INLINE quad_to_point #-}
+quad_to_point :: Quad -> Point
+quad_to_point (Q x y z _) = P x y z
+
+{-# INLINE quad_to_vector #-}
+quad_to_vector :: Quad -> Vector
+quad_to_vector (Q x y z _) = V x y z
+
+--{-# INLINE dot #-}
+dot :: Vector -> Vector -> Double
+dot (V x1 y1 z1) (V x2 y2 z2) = x1 * x2 + y1 * y2 + z1 * z2
+
+cross :: Vector -> Vector -> Vector
+cross (V x1 y1 z1) (V x2 y2 z2)
+ = V (y1 * z2 - z1 * y2) (z1 * x2 - x1 * z2) (x1 * y2 - y1 * x2)
+
+-- assumption: the input vector is a normal
+tangents :: Vector -> (Vector, Vector)
+tangents v@(V x y z)
+ = (v1, v `cross` v1)
+ where v1 | x == 0 = normalize (vector 0 z (-y))
+ | otherwise = normalize (vector (-y) x 0)
+
+{-# INLINE dot4 #-}
+dot4 :: Quad -> Quad -> Double
+dot4 (Q x1 y1 z1 w1) (Q x2 y2 z2 w2) = x1 * x2 + y1 * y2 + z1 * z2 + w1 * w2
+
+addVV :: Vector -> Vector -> Vector
+addVV (V x1 y1 z1) (V x2 y2 z2)
+ = V (x1 + x2) (y1 + y2) (z1 + z2)
+
+addPV :: Point -> Vector -> Point
+addPV (P x1 y1 z1) (V x2 y2 z2)
+ = P (x1 + x2) (y1 + y2) (z1 + z2)
+
+subVV :: Vector -> Vector -> Vector
+subVV (V x1 y1 z1) (V x2 y2 z2)
+ = V (x1 - x2) (y1 - y2) (z1 - z2)
+
+negV :: Vector -> Vector
+negV (V x1 y1 z1)
+ = V (-x1) (-y1) (-z1)
+
+subPP :: Point -> Point -> Vector
+subPP (P x1 y1 z1) (P x2 y2 z2)
+ = V (x1 - x2) (y1 - y2) (z1 - z2)
+
+--{-# INLINE norm #-}
+norm :: Vector -> Double
+norm (V x y z) = sqrt (sq x + sq y + sq z)
+
+--{-# INLINE normalize #-}
+-- normalize a vector to a unit vector
+normalize :: Vector -> Vector
+normalize v@(V x y z)
+ | norm /= 0 = multSV (1/norm) v
+ | otherwise = error "normalize empty!"
+ where norm = sqrt (sq x + sq y + sq z)
+
+-- This does computes the distance *squared*
+dist2 :: Point -> Point -> Double
+dist2 us vs = sq x + sq y + sq z
+ where
+ (V x y z) = subPP us vs
+
+{-# INLINE sq #-}
+sq :: Double -> Double
+sq d = d * d
+
+{-# INLINE distFrom0Sq #-}
+distFrom0Sq :: Point -> Double -- Distance of point from origin.
+distFrom0Sq (P x y z) = sq x + sq y + sq z
+
+{-# INLINE distFrom0 #-}
+distFrom0 :: Point -> Double -- Distance of point from origin.
+distFrom0 p = sqrt (distFrom0Sq p)
+
+--{-# INLINE multSV #-}
+multSV :: Double -> Vector -> Vector
+multSV k (V x y z) = V (k*x) (k*y) (k*z)
+
+--{-# INLINE multMM #-}
+multMM :: Matrix -> Matrix -> Matrix
+multMM m1@(M q1 q2 q3 q4) m2
+ = M (multMQ m2' q1)
+ (multMQ m2' q2)
+ (multMQ m2' q3)
+ (multMQ m2' q4)
+ where
+ m2' = transposeM m2
+
+{-# INLINE transposeM #-}
+transposeM :: Matrix -> Matrix
+transposeM (M (Q e11 e12 e13 e14)
+ (Q e21 e22 e23 e24)
+ (Q e31 e32 e33 e34)
+ (Q e41 e42 e43 e44)) = (M (Q e11 e21 e31 e41)
+ (Q e12 e22 e32 e42)
+ (Q e13 e23 e33 e43)
+ (Q e14 e24 e34 e44))
+
+
+--multMM m1 m2 = [map (dot4 row) (transpose m2) | row <- m1]
+
+--{-# INLINE multMV #-}
+multMV :: Matrix -> Vector -> Vector
+multMV m v = quad_to_vector (multMQ m (vector_to_quad v))
+
+--{-# INLINE multMP #-}
+multMP :: Matrix -> Point -> Point
+multMP m p = quad_to_point (multMQ m (point_to_quad p))
+
+-- mat vec = map (dot4 vec) mat
+
+{-# INLINE multMQ #-}
+multMQ :: Matrix -> Quad -> Quad
+multMQ (M q1 q2 q3 q4) q
+ = Q (dot4 q q1)
+ (dot4 q q2)
+ (dot4 q q3)
+ (dot4 q q4)
+
+{-# INLINE multMR #-}
+multMR :: Matrix -> Ray -> Ray
+multMR m (r, v) = (multMP m r, multMV m v)
+
+white :: Color
+white = C 1 1 1
+black :: Color
+black = C 0 0 0
+
+addCC :: Color -> Color -> Color
+addCC (C a b c) (C d e f) = C (a+d) (b+e) (c+f)
+
+subCC :: Color -> Color -> Color
+subCC (C a b c) (C d e f) = C (a-d) (b-e) (c-f)
+
+sumCC :: [Color] -> Color
+sumCC = foldr addCC black
+
+multCC :: Color -> Color -> Color
+multCC (C a b c) (C d e f) = C (a*d) (b*e) (c*f)
+
+multSC :: Double -> Color -> Color
+multSC k (C a b c) = C (a*k) (b*k) (c*k)
+
+nearC :: Color -> Color -> Bool
+nearC (C a b c) (C d e f) = a `near` d && b `near` e && c `near` f
+
+offsetToPoint :: Ray -> Double -> Point
+offsetToPoint (r,v) i = r `addPV` (i `multSV` v)
+
+--
+
+epsilon, inf :: Double -- aproximate zero and infinity
+epsilon = 1.0e-10
+inf = 1.0e20
+
+nonZero :: Double -> Double -- Use before a division. It makes definitions
+nonZero x | x > epsilon = x -- more complete and I bet the errors that get
+ | x < -epsilon = x -- introduced will be undetectable if epsilon
+ | otherwise = epsilon -- is small enough
+
+
+eqEps x y = abs (x-y) < epsilon
+near = eqEps
+
+clampf :: Double -> Double
+clampf p | p < 0 = 0
+ | p > 1 = 1
+ | True = p
--- /dev/null
+-- 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.
+
+-- Modified to use stdout (for testing)
+
+module Illumination
+ ( Object
+ , Light (..)
+ , light, pointlight, spotlight
+ , render
+ ) where
+
+import Array
+import Char(chr)
+import IOExts
+import Maybe
+
+import Geometry
+import CSG
+import Surface
+import Misc
+
+type Object = CSG (SurfaceFn Color Double)
+
+data Cxt = Cxt {ambient::Color, lights::[Light], object::Object, depth::Int}
+ deriving Show
+
+render :: (Matrix,Matrix) -> Color -> [Light] -> Object -> Int ->
+ Radian -> Int -> Int -> String -> IO ()
+render (m,m') amb ls obj dep fov wid ht file
+ = do { debugging
+ ; putStrLn (showBitmap wid ht pixels)
+ }
+ where
+ debugging = return ()
+{-
+ do { putStrLn (show cxt)
+ ; putStrLn (show (width, delta, aspect, left, top))
+ }
+-}
+ obj' = transform (m',m) obj
+ ls' = [ transformLight m' l | l <- ls ]
+ pixelA = listArray ((1,1), (ht,wid))
+ [ illumination cxt (start,pixel i j)
+ | j <- take ht [0.5..]
+ , i <- take wid [0.5..] ]
+ antiA = pixelA //
+ [ (ix, superSample ix (pixelA ! ix))
+ | j <- [2 .. ht - 1], i <- [2 .. wid - 1]
+ , let ix = (j, i)
+ , contrast ix pixelA ]
+ pixels = [ [ illumination cxt (start,pixel i j) | i<- take wid [0.5..] ]
+ | j <- take ht [0.5..]
+ ]
+ cxt = Cxt {ambient=amb, lights=ls', object=obj', depth=dep}
+ start = point 0 0 (-1)
+ width = 2 * tan (fov/2)
+ delta = width / fromIntegral wid
+ aspect = fromIntegral ht / fromIntegral wid
+ left = - width / 2
+ top = - left * aspect
+ pixel i j = vector (left + i*delta) (top - j*delta) 1
+
+ superSample (y, x) col = avg $ col:
+ [ illumination cxt (start, pixel (fromIntegral x - 0.5 + xd) (fromIntegral y - 0.5 + yd))
+ | (xd, yd) <- [(-0.333, 0.0), (0.333, 0.0), (0.0, -0.333), (0.0, 0.333)]
+ ]
+
+avg cs = divN (fromIntegral (length cs)) (uncolor (sumCC cs))
+ where divN n (r,g,b) = color (r / n) (g / n) (b / n)
+
+contrast :: (Int, Int) -> Array (Int, Int) Color -> Bool
+contrast (x, y) arr = any diffMax [ subCC cur (arr ! (x + xd, y + yd))
+ | xd <- [-1, 1], yd <- [-1, 1]
+ ]
+ where cur = arr ! (x, y)
+ diffMax col = (abs r) > 0.25 || (abs g) > 0.2 || (abs b) > 0.4
+ where
+ (r,g,b) = uncolor col
+
+
+illumination :: Cxt -> Ray -> Color
+illumination cxt (r,v)
+ | depth cxt <= 0 = black
+ | otherwise = case castRay (r,v) (object cxt) of
+ Nothing -> black
+ Just info -> illum (cxt{depth=(depth cxt)-1}) info v
+
+illum :: Cxt -> (Point,Vector,Properties Color Double) -> Vector -> Color
+illum cxt (pos,normV,(col,kd,ks,n)) v
+ = ambTerm `addCC` difTerm `addCC` spcTerm `addCC` recTerm
+ where
+ visibleLights = unobscured pos (object cxt) (lights cxt) normV
+ d = depth cxt
+ amb = ambient cxt
+ newV = subVV v (multSV (2 * dot normV v) normV)
+
+ ambTerm = multSC kd (multCC amb col)
+ difTerm = multSC kd (sumCC [multSC (dot normV lj) (multCC intensity col)
+ |(loc,intensity) <- visibleLights,
+ let lj = normalize ({- pos `subVV` -} loc)])
+ -- ZZ might want to avoid the phong, when you can...
+ spcTerm = multSC ks (sumCC [multSC ((dot normV hj) ** n ) (multCC intensity col)
+ |(loc,intensity) <- visibleLights,
+ -- ZZ note this is specific to the light at infinity
+ let lj = {- pos `subVV` -} normalize loc,
+ let hj = normalize (lj `subVV` normalize v)])
+ recTerm = if recCoeff `nearC` black then black else multCC recCoeff recRay
+ recCoeff = multSC ks col
+ recRay = illumination cxt (pos,newV)
+
+showBitmapA :: Int -> Int -> Array (Int, Int) Color -> String
+showBitmapA wid ht arr
+ = header ++ concatMap scaleColor (elems arr)
+ where
+ scaleColor col = [scalePixel r, scalePixel g, scalePixel b]
+ where (r,g,b) = uncolor col
+ header = "P6\n#Galois\n" ++ show wid ++ " " ++ show ht ++ "\n255\n"
+
+showBitmap :: Int -> Int ->[[Color]] -> String
+showBitmap wid ht pss
+-- type of assert | length pss == ht && all (\ ps -> length ps == wid) pss
+ = header ++ concat [[scalePixel r,scalePixel g,scalePixel b]
+ | ps <- pss, (r,g,b) <- map uncolor ps]
+ where
+ header = "P6\n#Galois\n" ++ show wid ++ " " ++ show ht ++ "\n255\n"
+showBitmap _ _ _ = error "incorrect length of bitmap string"
+
+scalePixel :: Double -> Char
+scalePixel p = chr (floor (clampf p * 255))
+
+
+-- Lights
+
+data Light = Light Vector Color
+ | PointLight Point Color
+ | SpotLight Point Point Color Radian Double
+ deriving Show
+
+light :: Coords -> Color -> Light
+light (x,y,z) color =
+ Light (normalize (vector (-x) (-y) (-z))) color
+pointlight (x,y,z) color =
+ PointLight (point x y z) color
+spotlight (x,y,z) (p,q,r) col cutoff exp =
+ SpotLight (point x y z) (point p q r) col cutoff exp
+
+transformLight m (Light v c) = Light (multMV m v) c
+transformLight m (PointLight p c) = PointLight (multMP m p) c
+transformLight m (SpotLight p q c r d) = SpotLight (multMP m p) (multMP m q) c r d
+
+unobscured :: Point -> Object -> [Light] -> Vector -> [(Vector,Color)]
+unobscured pos obj lights normV = catMaybes (map (unobscure pos obj normV) lights)
+
+unobscure :: Point -> Object -> Vector -> Light -> Maybe (Vector,Color)
+unobscure pos obj normV (Light vec color)
+ -- ZZ probably want to make this faster
+ | vec `dot` normV < 0 = Nothing
+ | intersects (pos `addPV` (0.0001 `multSV` vec),vec) obj = Nothing
+ | otherwise = Just (vec,color)
+unobscure pos obj normV (PointLight pp color)
+ | vec `dot` normV < 0 = Nothing
+ | intersectWithin (pos `addPV` (0.0001 `multSV` (normalize vec)), vec) obj = Nothing
+ | otherwise = Just (vec,is)
+ where vec = pp `subPP` pos
+ is = attenuate vec color
+unobscure org obj normV (SpotLight pos at color cutoff exp)
+ | vec `dot` normV < 0 = Nothing
+ | intersectWithin (org `addPV` (0.0001 `multSV` (normalize vec)), vec) obj = Nothing
+ | angle > cutoff = Nothing
+ | otherwise = Just (vec, is)
+ where vec = pos `subPP` org
+ vec' = pos `subPP` at
+ angle = acos (normalize vec `dot` (normalize vec'))
+
+ asp = normalize (at `subPP` pos)
+ qsp = normalize (org `subPP` pos)
+ is = attenuate vec (((asp `dot` qsp) ** exp) `multSC` color)
+
+attenuate :: Vector -> Color -> Color
+attenuate vec color = (100 / (99 + sq (norm vec))) `multSC` color
+
+--
+
+castRay ray p
+ = case intersectRayWithObject ray p of
+ (True, _, _) -> Nothing -- eye is inside
+ (False, [], _) -> Nothing -- eye is inside
+ (False, (0, b, _) : _, _) -> Nothing -- eye is inside
+ (False, (i, False, _) : _, _) -> Nothing -- eye is inside
+ (False, (t, b, (s, p0)) : _, _) ->
+ let (v, prop) = surface s p0 in
+ Just (offsetToPoint ray t, v, prop)
+
+intersects ray p
+ = case intersectRayWithObject ray p of
+ (True, _, _) -> False
+ (False, [], _) -> False
+ (False, (0, b, _) : _, _) -> False
+ (False, (i, False, _) : _, _) -> False
+ (False, (i, b, _) : _, _) -> True
+
+intersectWithin :: Ray -> Object -> Bool
+intersectWithin ray p
+ = case intersectRayWithObject ray p of
+ (True, _, _) -> False -- eye is inside
+ (False, [], _) -> False -- eye is inside
+ (False, (0, b, _) : _, _) -> False -- eye is inside
+ (False, (i, False, _) : _, _) -> False -- eye is inside
+ (False, (t, b, _) : _, _) -> t < 1.0
--- /dev/null
+-- 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
--- /dev/null
+-- 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 Interval
+ ( IList
+ , Intersection
+ , emptyIList, openIList
+ , mkEntry, mkExit
+ , entryexit, exitentry
+ , mapI
+ , unionIntervals, intersectIntervals, differenceIntervals
+ , complementIntervals
+ ) where
+
+import Geometry
+
+-- The result of a ray trace is represented as a list of surface
+-- intersections. Each intersection is a point along the ray with
+-- a flag indicating whether this intersection is an entry or an
+-- exit from the solid. Each intersection also carries unspecified
+-- surface data for use by the illumination model.
+
+-- Just the list of intersections isn't enough, however. An empty
+-- list can denote either a trace that is always within the solid
+-- or never in the solid. To dissambiguate, an extra flag is kept
+-- that indicates whether we are starting inside or outside of the
+-- solid. As a convenience, we also keep an additional flag that
+-- indicates whether the last intersection ends inside or outside.
+
+type IList a = (Bool, [Intersection a], Bool)
+type Intersection a = (Double, Bool, a)
+
+emptyIList = (False, [], False)
+openIList = (True, [], True)
+
+mapI f (b1, is, b2) = (b1, map f is, b2)
+
+isEntry (_, entry, _) = entry
+isExit (_, entry, _) = not entry
+
+mkEntry (t, a) = (t, True, a)
+mkExit (t, a) = (t, False, a)
+
+entryexit w1 w2 = (False, [mkEntry w1, mkExit w2], False)
+exitentry w1 w2 = (True, [mkExit w1, mkEntry w2], True)
+arrange w1@(t1, _) w2@(t2, _) | t1 < t2 = entryexit w1 w2
+ | otherwise = entryexit w2 w1
+
+
+cmpI :: Intersection a -> Intersection a -> Ordering
+cmpI (i, _, _) (j, _, _)
+ | i `near` j = EQ
+ | i < j = LT
+ | otherwise = GT
+
+bad (b1, [], b2) = b1 /= b2
+bad (b1, is, b2) = bad' b1 is || b2 /= b3
+ where (_, b3, _) = last is
+
+bad' b [] = False
+bad' b ((_, c, _) : is) = b == c || bad' c is
+
+unionIntervals :: IList a -> IList a -> IList a
+unionIntervals (isStartOpen, is, isEndOpen) (jsStartOpen, js, jsEndOpen)
+ = (isStartOpen || jsStartOpen, uniIntervals is js, isEndOpen || jsEndOpen)
+ where uniIntervals is [] | jsEndOpen = []
+ | otherwise = is
+ uniIntervals [] js | isEndOpen = []
+ | otherwise = js
+ uniIntervals is@(i : is') js@(j : js')
+ = case cmpI i j of
+ EQ -> if isEntry i == isEntry j then i : uniIntervals is' js'
+ else uniIntervals is' js'
+ LT -> if isEntry j then i : uniIntervals is' js
+ else uniIntervals is' js
+ GT -> if isEntry i then j : uniIntervals is js'
+ else uniIntervals is js'
+
+intersectIntervals :: IList a -> IList a -> IList a
+intersectIntervals is js
+ = complementIntervals (unionIntervals is' js')
+ where is' = complementIntervals is
+ js' = complementIntervals js
+
+differenceIntervals :: IList a -> IList a -> IList a
+differenceIntervals is js
+ = complementIntervals (unionIntervals is' js)
+ where is' = complementIntervals is
+
+complementIntervals :: IList a -> IList a
+complementIntervals (o1, is, o2)
+ = (not o1, [ (i, not isentry, a) | (i, isentry, a) <- is ], not o2)
+
+-- tests...
+
+{-
+mkIn, mkOut :: Double -> Intersection a
+mkIn x = (x, True, undefined)
+mkOut x = (x, False, undefined)
+
+i1 = (False, [ mkIn 2, mkOut 7 ], False)
+i1' = (True, [ mkOut 2, mkIn 7 ], True)
+i2 = (False, [ mkIn 1, mkOut 3, mkIn 4, mkOut 5, mkIn 6, mkOut 8 ], False)
+
+t1 = unionIntervals i1 i2
+t2 = intersectIntervals i1 i2
+t3 = intersectIntervals i2 i1
+t4 = complementIntervals i1
+t5 = intersectIntervals i2 i1'
+t6 = differenceIntervals i2 i1
+t7 = differenceIntervals i2 i2
+
+sh (o1,is,o2) =
+ do if o1 then putStr "..." else return ()
+ putStr $ foldr1 (++) (map si is)
+ if o2 then putStr "..." else return ()
+si (i, True, _, _) = "<" ++ show i
+si (i, False, _, _) = " " ++ show i ++ ">"
+-}
--- /dev/null
+-- 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.
+
+-- Modified to read sample input directly from a file.
+
+module Main where
+
+import System
+
+import Parse
+import Eval
+
+main = do { str <- readFile "galois.gml"
+ ; mainEval (rayParse str)
+ }
--- /dev/null
+TOP = ..
+include $(TOP)/mk/boilerplate.mk
+
+SRC_HC_OPTS += -package text -package lang
+
+all :: runtest
+
+include $(TOP)/mk/target.mk
+
--- /dev/null
+-- 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 Misc where
+
+import IOExts
+
+debug s v = trace (s ++" : "++ show v ++ "\n") v
+-- debug s v = v
--- /dev/null
+-- 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 Parse where
+
+import Char
+import Parsec hiding (token)
+
+import Data
+
+
+program :: Parser Code
+program =
+ do { whiteSpace
+ ; ts <- tokenList
+ ; eof
+ ; return ts
+ }
+
+tokenList :: Parser Code
+tokenList = many token <?> "list of tokens"
+
+token :: Parser GMLToken
+token =
+ do { ts <- braces tokenList ; return (TBody ts) }
+ <|> do { ts <- brackets tokenList ; return (TArray ts) }
+ <|> (do { s <- gmlString ; return (TString s) } <?> "string")
+ <|> (do { t <- pident False ; return t } <?> "identifier")
+ <|> (do { char '/' -- No whitespace after slash
+ ; t <- pident True ; return t } <?> "binding identifier")
+ <|> (do { n <- number ; return n } <?> "number")
+
+pident :: Bool -> Parser GMLToken
+pident rebind =
+ do { id <- ident
+ ; case (lookup id opTable) of
+ Nothing -> if rebind then return (TBind id) else return (TId id)
+ Just t -> if rebind then error ("Attempted rebinding of identifier " ++ id) else return t
+ }
+
+ident :: Parser String
+ident = lexeme $
+ do { l <- letter
+ ; ls <- many (satisfy (\x -> isAlphaNum x || x == '-' || x == '_'))
+ ; return (l:ls)
+ }
+
+gmlString :: Parser String
+gmlString = lexeme $ between (char '"') (char '"') (many (satisfy (\x -> isPrint x && x /= '"')))
+
+-- Tests for numbers
+-- Hugs breaks on big exponents (> ~40)
+test_number = "1234 -1234 1 -0 0" ++
+ " 1234.5678 -1234.5678 1234.5678e12 1234.5678e-12 -1234.5678e-12" ++
+ " -1234.5678e12 -1234.5678E-12 -1234.5678E12" ++
+ " 1234e11 1234E33 -1234e33 1234e-33" ++
+ " 123e 123.4e 123ee 123.4ee 123E 123.4E 123EE 123.4EE"
+
+
+-- Always int or real
+number :: Parser GMLToken
+number = lexeme $
+ do { s <- optSign
+ ; n <- decimal
+ ; do { string "."
+ ; m <- decimal
+ ; e <- option "" exponent'
+ ; return (TReal (read (s ++ n ++ "." ++ m ++ e))) -- FIXME: Handle error conditions
+ }
+ <|> do { e <- exponent'
+ ; return (TReal (read (s ++ n ++ ".0" ++ e)))
+ }
+ <|> do { return (TInt (read (s ++ n))) }
+ }
+
+exponent' :: Parser String
+exponent' = try $
+ do { e <- oneOf "eE"
+ ; s <- optSign
+ ; n <- decimal
+ ; return (e:s ++ n)
+ }
+
+decimal = many1 digit
+
+optSign :: Parser String
+optSign = option "" (string "-")
+
+
+------------------------------------------------------
+-- Library for tokenizing.
+
+braces p = between (symbol "{") (symbol "}") p
+brackets p = between (symbol "[") (symbol "]") p
+
+symbol name = lexeme (string name)
+
+lexeme p = do{ x <- p; whiteSpace; return x }
+
+whiteSpace = skipMany (simpleSpace <|> oneLineComment <?> "")
+ where simpleSpace = skipMany1 (oneOf " \t\n\r\v")
+ oneLineComment =
+ do{ string "%"
+ ; skipMany (noneOf "\n\r\v")
+ ; return ()
+ }
+
+
+------------------------------------------------------------------------------
+
+rayParse :: String -> Code
+rayParse is = case (parse program "<stdin>" is) of
+ Left err -> error (show err)
+ Right x -> x
+
+rayParseF :: String -> IO Code
+rayParseF file =
+ do { r <- parseFromFile program file
+ ; case r of
+ Left err -> error (show err)
+ Right x -> return x
+ }
+
+run :: String -> IO ()
+run is = case (parse program "" is) of
+ Left err -> print err
+ Right x -> print x
+
+runF :: IO ()
+runF =
+ do { r <- parseFromFile program "simple.gml"
+ ; case r of
+ Left err -> print err
+ Right x -> print x
+ }
--- /dev/null
+-- 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 Pixmap where
+
+import Char
+import IO hiding (try)
+import Parsec
+
+readPPM f
+ = do h <- openFile f ReadMode
+ s <- hGetContents h
+ case (parse parsePPM f s) of
+ Left err -> error (show err)
+ Right x -> return x
+
+writePPM f ppm
+ = do h <- openFile f WriteMode
+ let s = showPPM (length (head ppm)) (length ppm) ppm
+ hPutStr h s
+
+-- parsing
+
+parsePPM
+ = do string "P6"
+ whiteSpace
+ width <- number
+ whiteSpace
+ height <- number
+ whiteSpace
+ colormax <- number
+ whiteSpace
+ cs <- getInput
+ return (chop width (chopColors cs))
+
+chopColors [] = []
+chopColors (a:b:c:ds) = (ord a, ord b, ord c) : chopColors ds
+
+chop n [] = []
+chop n xs = h : chop n t
+ where (h, t) = splitAt n xs
+
+number
+ = do ds <- many1 digit
+ return (read ds :: Int)
+
+whiteSpace
+ = skipMany (simpleSpace <|> oneLineComment <?> "")
+ where simpleSpace = skipMany1 (oneOf " \t\n\r\v")
+ oneLineComment =
+ do char '#'
+ skipMany (noneOf "\n\r\v")
+ return ()
+
+-- printing
+
+showPPM :: Int -> Int -> [[(Int,Int,Int)]] -> String
+showPPM wid ht pss
+ = header ++ concat [[chr r, chr g, chr b] | ps <- pss, (r, g, b) <-ps]
+ where
+ header = "P6\n#Galois\n" ++ show wid ++ " " ++ show ht ++ "\n255\n"
+showPPM _ _ _ = error "incorrect length of bitmap string"
--- /dev/null
+-- 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 Primitives where
+
+rad2deg :: Double -> Double
+rad2deg r = r * 180 / pi
+
+deg2rad :: Double -> Double
+deg2rad d = d * pi / 180
+
+addi :: Int -> Int -> Int
+addi = (+)
+
+addf :: Double -> Double -> Double
+addf = (+)
+
+acosD :: Double -> Double
+acosD x = acos x * 180 / pi
+
+asinD :: Double -> Double
+asinD x = asin x * 180 / pi
--- /dev/null
+-- 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 RayTrace(module Illumination, module Surface) where
+
+import Illumination
+import Surface
--- /dev/null
+-- 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 Surface
+ ( SurfaceFn (..)
+ , Properties
+ , sfun, sconst
+ , prop
+ , matte, shiny
+ , chgColor
+ , surface
+ ) where
+
+import Geometry
+import CSG
+import Misc
+
+-- the surface gets passed face then u then v.
+data SurfaceFn c v = SFun (Int -> Double -> Double -> Properties c v)
+ | SConst (Properties c v)
+
+sfun :: (Int -> Double -> Double -> Properties c v) -> SurfaceFn c v
+sfun = SFun
+sconst :: Properties c v -> SurfaceFn c v
+sconst = SConst
+
+type Properties c v = (c, v, v, v)
+
+prop c d s p = (c, d, s, p)
+
+matte = (white, 1.0, 0.0, 1.0)
+shiny = (white, 0.0, 1.0, 1.0)
+
+chgColor :: c -> Properties d v -> Properties c v
+chgColor c (_, d, s, p) = (c, d, s, p)
+
+instance (Show c, Show v) => Show (SurfaceFn c v) where
+ show (SFun _) = "Surface function"
+ -- show (SConst p) = "Surface constant: " ++ show p
+ show (SConst p) = "Surface constant"
+
+evalSurface :: SurfaceFn Color Double -> Int -> Double -> Double -> Properties Color Double
+evalSurface (SConst p) = \_ _ _ -> p
+evalSurface (SFun f) = f
+
+-- calculate surface properties, given the type of
+-- surface, and intersection point in object coordinates
+
+-- surface :: Surface SurfaceFn -> (Int, Point) -> (Vector, Properties)
+
+surface (Planar _ v0 v1) (n, p0, fn)
+ = (norm, evalSurface fn n' u v)
+ where norm = normalize $ cross v0 v1
+ (n', u, v) = planarUV n p0
+
+surface (Spherical _ v0 v1) (_, p0, fn)
+ = (norm, evalSurface fn 0 u v)
+ where x = xCoord p0
+ y = yCoord p0
+ z = zCoord p0
+ k = sqrt (1 - sq y)
+ theta = adjustRadian (atan2 (x / k) (z / k))
+ -- correct so that the image grows left-to-right
+ -- instead of right-to-left
+ u = 1.0 - clampf (theta / (2 * pi))
+ v = clampf ((y + 1) / 2)
+ norm = normalize $ cross v0 v1
+
+-- ZZ ignore the (incorrect) surface model, and estimate the normal
+-- from the intersection in object space
+surface (Cylindrical _ v0 v1) (_, p0, fn)
+ = (norm, evalSurface fn 0 u v)
+ where x = xCoord p0
+ y = yCoord p0
+ z = zCoord p0
+ u = clampf $ adjustRadian (atan2 x z) / (2 * pi)
+ v = y
+ norm = normalize $ cross v0 v1
+
+-- ZZ ignore the (incorrect) surface model, and estimate the normal
+-- from the intersection in object space
+surface (Conic _ v0 v1) (_, p0, fn)
+ = (norm, evalSurface fn 0 u v)
+ where x = xCoord p0
+ y = yCoord p0
+ z = zCoord p0
+ u = clampf $ adjustRadian (atan2 (x / y) (z / y)) / (2 * pi)
+ v = y
+ norm = normalize $ cross v0 v1
+
+planarUV face p0
+ = case face of
+ PlaneFace -> (0, x, z)
+
+ CubeFront -> (0, x, y)
+ CubeBack -> (1, x, y)
+ CubeLeft -> (2, z, y)
+ CubeRight -> (3, z, y)
+ CubeTop -> (4, x, z)
+ CubeBottom -> (5, x, z)
+
+ CylinderTop -> (1, (x + 1) / 2, (z + 1) / 2)
+ CylinderBottom -> (2, (x + 1) / 2, (z + 1) / 2)
+
+ ConeBase -> (1, (x + 1) / 2, (z + 1) / 2)
+ where x = xCoord p0
+ y = yCoord p0
+ z = zCoord p0
+
+-- misc
+
+adjustRadian :: Radian -> Radian
+adjustRadian r = if r > 0 then r else r + 2 * pi
--- /dev/null
+
+[ [97 95 73 50 89 97 99 99 99 99 99 99 99 99 99 98 99 99 99 99 99 99 99 99 99 99 99 99 99 99 99 98 99 99 99 99 99 98 99 99 98 97 99 99 99 99 99 99 99 99 99 99 97 97 96 96 96 96 96 96 99 99 99]
+ [88 96 66 53 52 86 99 99 99 99 99 99 99 99 99 99 99 99 99 98 99 99 99 99 99 99 99 99 99 99 98 96 98 99 99 99 99 99 99 99 99 97 98 99 99 99 99 99 99 99 99 99 96 96 96 98 97 96 96 96 97 97 96]
+ [89 92 79 50 54 45 91 98 99 99 99 99 99 99 99 99 99 99 99 99 99 99 99 99 99 99 99 98 98 99 98 96 98 98 98 98 97 98 99 99 99 99 99 99 99 99 99 99 99 99 99 97 96 96 97 99 99 96 96 98 98 97 97]
+ [88 91 96 81 40 35 39 91 95 99 99 99 99 99 99 99 99 99 99 99 99 99 99 99 99 98 97 97 98 99 98 99 99 96 95 95 95 96 97 99 99 98 97 99 99 97 98 99 99 98 96 96 96 94 96 98 99 96 96 96 97 95 96]
+ [83 92 91 48 54 33 62 64 98 99 99 99 99 99 99 99 97 98 99 99 99 99 99 99 99 98 98 99 99 99 99 99 99 98 97 97 98 99 97 99 99 99 99 99 99 99 99 99 99 99 96 97 98 96 96 97 99 98 96 95 96 96 96]
+ [91 93 64 78 94 75 57 50 81 97 99 99 99 99 99 98 94 96 99 99 99 99 99 99 99 97 98 99 99 99 98 98 99 98 99 98 99 99 98 99 99 99 99 99 99 99 99 99 99 99 97 99 98 96 91 96 97 98 98 96 96 96 99]
+ [95 63 85 94 84 95 72 61 44 84 96 98 99 99 99 99 98 98 98 99 99 99 99 99 99 97 96 98 99 99 97 96 98 99 99 97 98 99 98 97 99 99 99 99 99 99 99 99 99 99 97 99 99 96 93 98 96 97 96 96 96 96 95]
+ [63 80 88 96 96 88 90 52 64 52 95 98 99 99 99 99 99 99 97 98 99 99 99 97 97 99 98 97 97 97 99 96 98 99 99 97 98 99 99 98 98 98 99 99 99 99 99 99 99 99 98 98 99 96 96 96 94 97 98 99 96 92 95]
+ [92 84 90 92 91 88 89 75 50 58 64 96 99 99 99 99 99 99 99 99 99 99 99 98 99 99 98 99 99 98 99 98 99 99 99 99 99 99 99 99 98 98 99 99 99 96 99 98 99 99 99 97 96 97 96 96 92 96 99 98 95 94 95]
+ [91 80 85 85 92 96 93 87 81 49 66 88 99 99 99 99 99 99 99 99 99 99 99 98 99 99 96 98 99 99 99 99 99 98 99 99 99 98 99 99 99 99 99 99 99 98 99 98 99 99 98 97 98 96 96 96 93 96 99 98 96 97 97]
+ [70 90 96 96 95 95 97 93 60 73 64 67 93 97 99 99 99 99 99 99 99 99 99 98 97 97 98 99 99 99 99 99 98 94 97 98 99 98 99 99 99 99 99 99 99 99 98 98 98 97 98 99 99 96 96 96 96 96 99 97 96 96 95]
+ [93 93 97 97 94 88 85 89 90 57 72 43 82 97 99 99 99 99 99 99 99 99 99 99 98 96 97 99 99 99 99 99 99 96 96 96 98 99 99 98 99 97 98 99 99 98 95 81 88 84 98 98 95 96 96 95 96 96 98 95 94 94 92]
+ [87 96 91 94 96 97 98 94 75 66 76 60 67 83 99 99 99 99 99 99 99 99 99 99 99 98 98 99 99 99 99 99 99 97 97 97 95 96 97 95 96 76 70 66 73 83 92 60 88 58 88 95 95 95 96 94 95 96 98 96 97 97 97]
+ [90 96 86 84 89 85 93 92 96 96 94 84 56 85 98 98 99 99 99 99 99 99 99 99 99 99 99 99 99 99 99 98 99 97 56 40 57 71 69 66 78 73 84 55 34 39 39 41 44 46 31 61 90 98 97 92 94 98 99 97 98 98 98]
+ [93 91 94 89 66 81 86 94 89 87 97 82 84 65 84 82 87 89 97 96 97 99 99 99 99 99 99 99 99 99 99 97 96 40 4 15 9 9 6 2 7 14 23 9 8 8 10 3 3 8 13 12 17 42 90 93 86 93 96 97 96 96 96]
+ [85 82 82 76 91 90 86 83 84 86 54 37 26 49 31 20 13 43 40 55 80 95 98 98 99 99 99 99 99 98 99 97 66 10 3 4 12 19 4 4 7 18 7 4 11 11 6 11 24 11 12 8 18 9 30 85 93 90 98 99 96 93 96]
+ [77 68 78 66 72 66 79 92 73 57 57 73 69 39 73 68 70 57 22 15 17 21 50 83 93 97 98 99 99 97 94 66 12 0 0 0 3 4 6 21 7 3 5 2 3 14 17 7 4 6 28 39 39 24 20 53 95 80 98 99 98 96 97]
+ [69 75 86 82 89 86 64 51 66 93 93 80 90 92 89 87 92 94 68 34 24 13 17 6 20 55 91 98 98 95 66 10 18 4 3 3 1 3 12 10 9 11 17 3 15 14 17 23 15 26 25 25 21 43 21 10 65 62 97 98 99 98 96]
+ [95 94 94 68 45 48 41 69 84 72 80 94 78 78 94 96 88 93 95 94 83 67 35 13 28 38 15 40 87 80 17 7 4 7 6 9 37 20 19 19 22 9 5 7 11 13 26 14 29 27 41 35 44 19 20 16 43 70 96 96 97 99 99]
+ [93 96 88 54 57 57 75 74 70 79 94 96 94 97 97 88 88 96 83 97 97 72 65 75 46 21 17 9 28 29 19 22 4 9 8 5 24 34 17 19 7 8 25 16 9 5 8 12 15 32 54 43 30 14 18 11 45 62 96 96 96 97 98]
+ [97 91 56 59 44 70 63 78 88 82 73 82 96 97 78 94 94 96 99 99 98 91 54 21 4 4 3 11 14 35 14 17 7 26 36 21 11 35 42 24 3 9 17 13 9 20 9 21 20 45 28 35 42 19 38 14 36 67 98 98 97 97 96]
+ [80 42 56 39 41 71 74 62 54 57 62 68 96 87 92 95 93 97 98 98 67 21 7 4 5 19 22 25 10 8 24 7 6 26 42 43 38 42 32 13 36 60 81 88 67 62 37 63 28 15 26 44 32 30 25 38 62 70 73 78 87 95 97]
+ [33 60 39 44 80 73 40 59 58 55 70 71 90 91 86 97 98 98 92 43 8 9 2 6 42 65 50 45 27 27 24 9 13 37 66 66 88 80 60 72 92 90 96 98 94 84 43 43 18 31 34 35 19 53 28 46 31 23 49 56 40 78 96]
+ [74 32 63 76 67 52 85 93 81 74 85 83 65 80 88 96 98 68 30 4 2 2 2 24 82 56 40 25 44 44 30 52 70 83 87 91 90 73 69 51 86 97 99 99 97 92 44 39 28 26 26 22 49 23 36 67 37 17 23 27 39 54 91]
+ [44 62 77 88 81 94 95 94 55 78 57 42 77 79 90 92 39 7 3 9 1 4 19 68 57 39 49 54 62 53 74 92 88 97 99 99 90 64 82 53 84 99 99 99 97 99 87 76 52 39 48 37 35 26 51 32 8 6 19 6 42 78 84]
+ [61 82 82 81 96 90 94 66 64 77 65 81 96 97 81 20 2 5 9 10 21 23 61 67 55 48 53 61 48 87 97 96 98 98 98 96 69 50 85 60 65 98 99 99 98 97 94 62 45 50 57 54 21 11 16 11 6 11 21 6 26 71 75]
+ [93 77 71 97 90 93 51 85 65 77 94 96 94 77 15 1 0 0 8 9 8 53 73 61 64 59 59 53 87 98 98 98 96 94 98 75 37 34 84 51 37 97 98 99 99 97 57 5 5 4 5 4 7 3 18 4 6 15 9 4 8 64 76]
+ [96 82 93 79 95 82 81 79 85 96 91 74 44 11 12 9 2 0 5 5 24 69 55 63 84 82 75 78 92 98 99 97 95 91 89 33 38 41 70 55 30 94 98 99 98 90 11 8 4 6 3 12 41 2 3 1 6 6 6 6 16 30 58]
+ [89 66 92 84 96 98 97 92 96 75 45 29 9 38 6 10 3 10 1 8 69 67 64 70 85 87 93 94 97 99 99 99 94 94 70 44 23 60 78 55 28 98 99 99 99 78 6 9 16 15 5 2 8 4 7 1 18 13 18 16 4 23 66]
+ [48 79 85 82 90 97 96 95 64 20 25 26 27 12 5 8 6 7 6 42 74 64 49 80 94 96 98 99 99 99 99 98 94 96 71 50 7 55 92 31 49 98 99 99 98 62 7 10 13 8 12 11 9 3 5 4 16 9 11 9 13 44 68]
+ [93 72 87 92 96 98 95 62 38 9 23 64 5 6 12 32 7 9 15 70 69 54 75 97 96 97 99 99 99 99 99 98 98 96 71 79 62 66 69 20 41 98 99 99 98 63 9 29 38 33 28 17 17 40 10 9 18 10 22 10 8 39 74]
+ [93 69 84 88 95 93 72 41 17 31 69 54 35 6 6 26 18 36 39 68 77 82 95 98 98 97 99 98 99 99 99 99 99 97 91 89 92 63 31 11 44 98 99 99 98 77 19 8 13 22 30 28 35 14 18 21 15 7 37 46 23 53 91]
+ [77 85 92 96 88 65 24 23 58 63 44 35 24 17 9 24 18 21 43 54 87 95 97 98 98 97 99 97 96 98 99 99 99 99 97 89 64 47 42 16 68 99 99 99 99 94 31 17 37 25 20 51 13 7 8 7 6 34 40 39 65 63 94]
+ [58 83 97 92 90 42 34 66 76 38 42 35 22 12 10 51 8 33 54 80 98 96 94 90 56 88 94 95 70 68 98 98 98 97 96 94 89 66 30 48 97 99 99 99 98 98 61 13 24 10 20 19 18 6 15 10 11 25 43 27 63 53 92]
+ [87 71 66 73 79 65 69 44 74 47 44 38 18 15 24 29 9 38 67 93 98 96 79 53 46 96 95 88 23 61 87 95 95 95 87 79 56 35 46 89 99 99 99 99 99 98 67 23 13 10 21 49 31 21 26 16 17 36 46 40 46 67 96]
+ [96 80 47 69 60 77 59 22 29 74 61 50 56 41 34 27 12 48 74 96 99 97 56 40 23 95 96 63 39 77 85 82 77 75 69 68 71 77 90 97 99 99 99 99 99 98 66 29 44 17 7 23 23 18 24 27 16 33 67 18 30 87 96]
+ [97 96 92 70 45 52 21 22 30 18 38 74 82 25 24 21 12 35 74 96 99 93 31 57 33 95 94 35 17 53 50 48 77 89 93 96 96 98 95 97 99 99 99 99 99 99 72 7 16 28 22 31 14 10 20 23 16 13 39 25 54 95 92]
+ [90 64 52 70 92 83 73 49 25 8 36 40 44 25 39 28 42 26 72 96 98 96 28 41 43 95 88 18 8 75 92 89 91 95 98 93 75 61 94 97 99 99 99 99 99 98 82 13 12 20 34 25 9 23 73 55 45 19 21 15 93 95 96]
+ [22 38 82 96 93 95 93 81 89 51 45 46 26 36 28 24 36 15 31 94 98 95 55 44 30 80 93 34 19 87 88 89 83 94 94 79 30 16 33 90 97 99 99 99 98 98 74 26 9 32 38 21 14 54 74 41 46 46 18 29 91 92 96]
+ [79 92 98 98 98 95 86 82 85 91 94 81 47 28 13 16 47 7 6 69 93 95 92 70 40 75 95 83 44 42 43 83 94 96 81 61 26 26 7 32 85 98 99 99 98 95 62 38 8 13 20 30 23 62 42 31 3 11 11 59 94 97 97]
+ [96 96 97 98 98 85 82 94 93 88 89 96 94 92 59 40 19 12 3 10 37 84 90 94 79 60 93 67 75 65 68 94 97 96 66 40 9 11 48 4 35 93 98 99 97 62 50 36 4 3 15 53 33 63 28 6 2 15 56 72 97 97 95]
+ [98 94 97 95 86 94 91 89 90 95 95 96 95 93 98 88 75 58 34 29 7 21 54 83 84 89 75 79 85 94 94 95 98 82 43 25 6 21 60 24 15 89 97 99 70 66 47 48 2 0 3 22 21 7 4 1 6 55 83 96 98 95 93]
+ [93 96 96 91 60 65 73 87 93 96 95 91 91 88 94 97 97 97 90 76 74 25 2 19 37 59 72 71 83 91 87 90 92 95 67 47 69 37 51 45 3 76 98 98 46 83 71 21 3 4 2 19 13 2 0 7 60 89 97 99 98 92 90]
+ [90 88 92 95 97 95 95 93 87 88 82 86 95 93 92 91 97 96 49 73 58 58 24 6 10 58 77 49 56 84 82 89 87 94 87 46 71 37 47 67 9 54 98 98 56 50 58 13 14 18 0 10 20 4 2 31 96 96 99 99 99 96 95]
+ [78 55 84 95 97 89 73 71 89 92 91 89 90 90 83 82 93 93 68 54 41 27 24 49 85 96 97 89 72 70 48 70 76 91 95 92 43 25 40 87 27 69 97 96 75 33 13 12 10 5 2 14 6 3 15 86 99 97 99 99 97 96 97]
+ [43 61 66 70 93 79 86 94 89 84 80 84 90 96 90 92 94 92 81 37 34 70 94 96 99 98 99 98 96 95 85 77 57 62 63 59 55 50 52 77 33 63 91 71 18 6 4 5 4 4 4 12 4 15 77 98 99 99 99 98 95 95 96]
+ [77 82 75 56 62 67 86 85 84 89 92 90 80 77 52 61 56 59 70 89 98 98 99 99 99 99 99 99 99 98 98 98 92 92 91 69 61 47 28 22 4 9 14 7 4 3 7 11 4 9 4 4 22 79 97 96 97 99 98 96 93 95 96]
+ [97 96 93 97 94 86 69 64 73 59 71 71 71 64 87 94 95 96 99 99 99 99 99 99 99 97 99 99 99 99 99 99 98 96 96 96 96 95 96 93 77 35 11 4 10 2 7 6 2 1 18 59 97 97 99 97 97 99 99 97 95 95 95]
+ [96 97 98 94 97 95 89 75 65 85 89 93 97 98 99 98 98 99 99 99 99 99 99 97 97 99 98 99 99 97 96 96 95 94 95 97 96 96 96 97 98 93 75 45 19 4 2 4 18 49 90 98 99 99 99 95 96 98 96 96 96 95 96]
+ [96 98 99 97 96 95 78 59 78 97 98 98 99 99 99 99 99 99 99 99 99 99 99 98 98 99 97 97 98 96 92 91 92 95 96 98 98 95 96 97 99 98 96 97 93 82 80 87 98 97 98 99 99 98 97 96 96 98 96 96 96 97 96]
+ [95 84 92 97 95 94 50 80 95 98 99 99 99 99 99 97 96 99 98 97 98 99 99 97 97 99 97 97 96 96 90 82 89 92 97 94 94 94 93 96 99 99 98 99 97 99 98 99 99 97 97 97 98 98 96 97 96 96 96 96 96 98 96]
+ [97 90 92 98 92 52 39 93 98 99 99 99 99 99 98 95 95 98 98 96 97 99 99 97 97 99 99 98 98 95 93 89 87 92 94 88 90 90 88 94 96 97 98 98 97 99 99 99 99 97 97 96 95 96 95 96 96 97 98 98 99 98 96]
+ [96 94 90 95 69 43 75 97 99 99 99 99 99 99 97 96 94 94 95 96 99 99 97 94 95 98 99 97 95 95 98 95 87 93 89 89 87 88 93 94 90 94 96 96 98 98 99 98 99 98 96 97 97 95 93 95 96 98 99 99 99 98 99]
+ [97 81 85 68 51 53 96 97 99 99 99 99 99 99 99 97 94 95 95 96 96 98 98 97 95 94 95 95 84 90 96 95 93 95 94 88 79 88 92 97 95 95 95 96 98 99 97 97 99 99 99 98 95 94 96 98 97 99 98 98 98 99 99]
+ [81 86 70 56 71 72 97 99 99 99 99 99 99 99 99 99 97 98 97 96 96 96 98 98 95 94 96 96 92 94 92 94 95 94 93 91 91 93 95 97 92 87 93 95 98 99 95 95 99 98 98 98 94 88 95 95 96 98 96 97 97 96 98]
+ [67 62 30 43 90 73 95 98 98 96 98 98 97 98 95 98 95 95 96 95 98 94 95 96 94 95 96 95 95 96 90 94 96 96 96 92 91 92 94 96 97 91 91 95 94 97 95 93 98 96 96 96 97 95 96 93 95 96 96 96 97 95 96]
+ [53 35 46 37 71 72 94 97 99 99 98 97 98 98 95 94 96 97 96 96 98 87 92 93 96 95 95 92 93 96 95 93 93 89 94 90 87 89 92 90 95 96 97 90 90 96 95 94 97 97 96 96 96 97 97 94 95 96 96 97 96 97 96]
+ [74 88 56 35 91 90 88 94 97 97 95 95 97 98 95 94 96 97 96 96 98 96 95 95 92 88 91 92 92 94 93 94 90 86 88 95 94 93 93 92 97 96 96 94 95 96 96 96 98 99 98 99 98 98 98 97 95 97 98 98 96 94 93]
+ [92 93 63 65 92 97 93 95 99 96 89 98 99 98 96 98 96 93 93 95 97 94 91 94 97 94 95 96 93 88 85 90 92 86 90 97 96 96 96 93 94 97 95 96 94 96 99 96 96 96 96 98 98 96 95 95 94 97 96 93 92 95 93]
+ [82 94 95 81 92 95 92 97 99 98 94 93 92 97 94 97 97 98 94 90 97 95 92 94 94 90 91 96 92 88 94 94 89 83 90 97 96 96 95 96 97 95 94 97 96 96 97 94 93 93 96 98 94 96 96 92 89 90 94 94 94 92 88]
+ [78 85 76 94 97 95 96 97 99 99 98 96 97 97 96 97 97 95 95 96 98 96 96 96 88 86 93 96 94 93 89 88 88 90 90 94 94 97 97 96 96 97 98 98 98 96 92 87 96 96 96 96 96 94 92 93 88 88 93 95 94 90 87]
+ [83 88 91 94 97 97 99 98 98 98 98 99 99 98 99 98 99 97 98 97 98 96 96 95 96 94 96 95 95 91 85 90 90 93 94 94 92 94 95 96 98 97 98 98 97 97 97 96 96 96 95 96 94 95 95 94 93 93 94 95 93 88 91]
+ [95 90 94 94 98 96 98 99 96 98 97 97 98 98 99 99 96 95 97 97 99 95 96 98 93 96 96 96 93 95 89 93 93 95 96 96 97 97 97 97 99 98 97 96 98 99 99 99 98 93 93 96 96 96 95 94 91 92 90 93 94 96 96]
+] /galois
+
+{ /v /u /face % bind parameters
+ { % toIntCoord : float -> int
+ 63.0 mulf floor /i % i = floor(3.0*i)
+ i 63 eqi { 62 } { i } if % return max(2, i)
+ } /toIntCoord
+ galois u toIntCoord apply get % val = texture[u][v]
+ v toIntCoord apply get
+ real 100.0 divf /gal
+ gal gal gal point % b/w galois
+ 1.0 % kd = 1.0
+ 0.0 % ks = 0.0
+ 1.0 % n = 1.0
+} /galoisface
+
+
+galoisface cube
+-0.5 -0.5 -0.5 translate % center
+2.5 uscale % make it bigger
+-25.0 rotatex -25.0 rotatey % rotate
+0.0 -1.0 7.0 translate % move to final position
+
+%galoisface cylinder
+%-0.5 -0.5 -0.5 translate % center
+%1.5 uscale % make it bigger
+%0.0 rotatex 90.0 rotatey % rotate
+%0.0 0.0 5.0 translate % move to final position
+
+%galoisface sphere
+%-0.5 -0.5 -0.5 translate % center
+%1.5 uscale % make it bigger
+%-25.0 rotatex 25.0 rotatey % rotate
+%-3.0 0.0 5.0 translate % move to final position
+
+%union union % join the 3 together
+
+{ /v /u /face
+ v 5.0 divf /v
+ u 5.0 divf /u
+ v floor 2 modi 0 eqi
+ { 1.0 }
+ { 0.8 }
+ if /r
+ u floor 2 modi 0 eqi
+ { 1.0 }
+ { 0.8 }
+ if /g
+ v frac /v
+ u frac /u
+ v 0.0 lessf { v 1.0 addf } { v } if /v
+ u 0.0 lessf { u 1.0 addf } { u } if /u
+ { % toIntCoord : float -> int
+ 63.0 mulf floor /i % i = floor(3.0*i)
+ i 63 eqi { 62 } { i } if % return max(2, i)
+ } /toIntCoord
+ galois u toIntCoord apply get % val = texture[u][v]
+ v toIntCoord apply get
+ real 100.0 divf /gal
+ r gal mulf g gal mulf gal point % b/w galois
+ 0.0 % kd = 1.0
+ 1.0 % ks = 0.0
+ 1.0 % n = 1.0
+} plane /p
+
+p 0.0 -3.0 0.0 translate % plane at Y = -3
+
+union
+
+/scene
+ % directional light
+1.0 -1.0 0.0 point % direction
+1.0 1.0 1.0 point light /l % directional light
+
+1.0 0.5 0.0 point % ambient light
+[ l ] % lights
+scene % scene to render
+3 % tracing depth
+90.0 % field of view
+300 200
+"galois.ppm" % output file
+render
+
projects / ghc-hetmet.git / commitdiff