ghc/tests/programs/galois_raytrace/Geometry.hs

   1 -- Copyright (c) 2000 Galois Connections, Inc.
   2 -- All rights reserved.  This software is distributed as
   3 -- free software under the license in the file "LICENSE",
   4 -- which is included in the distribution.
   5
   6 module Geometry
   7     ( Coords
   8     , Ray
   9     , Point  -- abstract
  10     , Vector -- abstract
  11     , Matrix -- abstract
  12     , Color  -- abstract
  13     , Box(..)
  14     , Radian
  15     , matrix
  16     , coord
  17     , color
  18     , uncolor
  19     , xCoord , yCoord , zCoord
  20     , xComponent , yComponent , zComponent
  21     , point
  22     , vector
  23     , nearV
  24     , point_to_vector
  25     , vector_to_point
  26     , dot
  27     , cross
  28     , tangents
  29     , addVV
  30     , addPV
  31     , subVV
  32     , negV
  33     , subPP
  34     , norm
  35     , normalize
  36     , dist2
  37     , sq
  38     , distFrom0Sq
  39     , distFrom0
  40     , multSV
  41     , multMM
  42     , transposeM
  43     , multMV
  44     , multMP
  45     , multMQ
  46     , multMR
  47     , white
  48     , black
  49     , addCC
  50     , subCC
  51     , sumCC
  52     , multCC
  53     , multSC
  54     , nearC
  55     , offsetToPoint
  56     , epsilon
  57     , inf
  58     , nonZero
  59     , eqEps
  60     , near
  61     , clampf
  62     ) where
  63
  64 import List
  65
  66 type Coords = (Double,Double,Double)
  67
  68 type Ray = (Point,Vector)    -- origin of ray, and unit vector giving direction
  69
  70 data Point  = P !Double !Double !Double -- implicit extra arg of 1
  71     deriving (Show)
  72 data Vector = V !Double !Double !Double -- implicit extra arg of 0
  73     deriving (Show, Eq)
  74 data Matrix = M !Quad   !Quad   !Quad   !Quad
  75     deriving (Show)
  76
  77 data Color  = C !Double !Double !Double
  78     deriving (Show, Eq)
  79
  80 data Box = B !Double !Double !Double !Double !Double !Double
  81     deriving (Show)
  82
  83 data Quad   = Q !Double !Double !Double !Double
  84     deriving (Show)
  85
  86 type Radian = Double
  87
  88 type Tup4 a = (a,a,a,a)
  89
  90 --{-# INLINE matrix #-}
  91 matrix :: Tup4 (Tup4 Double) -> Matrix
  92 matrix ((m11, m12, m13, m14),
  93           (m21, m22, m23, m24),
  94           (m31, m32, m33, m34),
  95           (m41, m42, m43, m44))
  96   = M (Q m11 m12 m13 m14)
  97       (Q m21 m22 m23 m24)
  98       (Q m31 m32 m33 m34)
  99       (Q m41 m42 m43 m44)
 100
 101 coord x y z = (x, y, z)
 102
 103 color r g b = C r g b
 104
 105 uncolor (C r g b) = (r,g,b)
 106
 107 {-# INLINE xCoord #-}
 108 xCoord (P x y z) = x
 109 {-# INLINE yCoord #-}
 110 yCoord (P x y z) = y
 111 {-# INLINE zCoord #-}
 112 zCoord (P x y z) = z
 113
 114 {-# INLINE xComponent #-}
 115 xComponent (V x y z) = x
 116 {-# INLINE yComponent #-}
 117 yComponent (V x y z) = y
 118 {-# INLINE zComponent #-}
 119 zComponent (V x y z) = z
 120
 121 point :: Double -> Double -> Double -> Point
 122 point x y z = P x y z
 123
 124 vector :: Double -> Double -> Double -> Vector
 125 vector x y z = V x y z
 126
 127 nearV :: Vector -> Vector -> Bool
 128 nearV (V a b c) (V d e f) = a `near` d && b `near` e && c `near` f
 129
 130 point_to_vector :: Point -> Vector
 131 point_to_vector (P x y z) = V x y z
 132
 133 vector_to_point :: Vector -> Point
 134 vector_to_point (V x y z)  = P x y z
 135
 136 {-# INLINE vector_to_quad #-}
 137 vector_to_quad :: Vector -> Quad
 138 vector_to_quad (V x y z) = Q x y z 0
 139
 140 {-# INLINE point_to_quad #-}
 141 point_to_quad :: Point -> Quad
 142 point_to_quad (P x y z) = Q x y z 1
 143
 144 {-# INLINE quad_to_point #-}
 145 quad_to_point :: Quad -> Point
 146 quad_to_point (Q x y z _) = P x y z
 147
 148 {-# INLINE quad_to_vector #-}
 149 quad_to_vector :: Quad -> Vector
 150 quad_to_vector (Q x y z _) = V x y z
 151
 152 --{-# INLINE dot #-}
 153 dot :: Vector -> Vector -> Double
 154 dot (V x1 y1 z1) (V x2 y2 z2) = x1 * x2 + y1 * y2 + z1 * z2
 155
 156 cross :: Vector -> Vector -> Vector
 157 cross (V x1 y1 z1) (V x2 y2 z2)
 158   = V (y1 * z2 - z1 * y2) (z1 * x2 - x1 * z2) (x1 * y2 - y1 * x2)
 159
 160 -- assumption: the input vector is a normal
 161 tangents :: Vector -> (Vector, Vector)
 162 tangents v@(V x y z)
 163   = (v1, v `cross` v1)
 164   where v1 | x == 0    = normalize (vector 0 z (-y))
 165            | otherwise = normalize (vector (-y) x 0)
 166
 167 {-# INLINE dot4 #-}
 168 dot4 :: Quad -> Quad -> Double
 169 dot4 (Q x1 y1 z1 w1) (Q x2 y2 z2 w2) = x1 * x2 + y1 * y2 + z1 * z2 + w1 * w2
 170
 171 addVV :: Vector -> Vector -> Vector
 172 addVV (V x1 y1 z1) (V x2 y2 z2)
 173     = V (x1 + x2) (y1 + y2) (z1 + z2)
 174
 175 addPV :: Point -> Vector -> Point
 176 addPV (P x1 y1 z1) (V x2 y2 z2)
 177     = P (x1 + x2) (y1 + y2) (z1 + z2)
 178
 179 subVV :: Vector -> Vector -> Vector
 180 subVV (V x1 y1 z1) (V x2 y2 z2)
 181     = V (x1 - x2) (y1 - y2) (z1 - z2)
 182
 183 negV :: Vector -> Vector
 184 negV (V x1 y1 z1)
 185     = V (-x1) (-y1) (-z1)
 186
 187 subPP :: Point -> Point -> Vector
 188 subPP (P x1 y1 z1) (P x2 y2 z2)
 189     = V (x1 - x2) (y1 - y2) (z1 - z2)
 190
 191 --{-# INLINE norm #-}
 192 norm :: Vector -> Double
 193 norm (V x y z) = sqrt (sq x + sq y + sq z)
 194
 195 --{-# INLINE normalize #-}
 196 -- normalize a vector to a unit vector
 197 normalize :: Vector -> Vector
 198 normalize v@(V x y z)
 199              | norm /= 0 = multSV (1/norm) v
 200              | otherwise = error "normalize empty!"
 201     where norm = sqrt (sq x + sq y + sq z)
 202
 203 -- This does computes the distance *squared*
 204 dist2 :: Point -> Point -> Double
 205 dist2 us vs = sq x + sq y + sq z
 206     where
 207        (V x y z) = subPP us vs
 208
 209 {-# INLINE sq #-}
 210 sq :: Double -> Double
 211 sq d = d * d
 212
 213 {-# INLINE distFrom0Sq #-}
 214 distFrom0Sq :: Point -> Double  -- Distance of point from origin.
 215 distFrom0Sq (P x y z) = sq x + sq y + sq z
 216
 217 {-# INLINE distFrom0 #-}
 218 distFrom0 :: Point -> Double  -- Distance of point from origin.
 219 distFrom0 p = sqrt (distFrom0Sq p)
 220
 221 --{-# INLINE multSV #-}
 222 multSV :: Double -> Vector -> Vector
 223 multSV k (V x y z) = V (k*x) (k*y) (k*z)
 224
 225 --{-# INLINE multMM #-}
 226 multMM :: Matrix -> Matrix -> Matrix
 227 multMM m1@(M q1 q2 q3 q4) m2
 228      = M (multMQ m2' q1)
 229          (multMQ m2' q2)
 230          (multMQ m2' q3)
 231          (multMQ m2' q4)
 232   where
 233      m2' = transposeM m2
 234
 235 {-# INLINE transposeM #-}
 236 transposeM :: Matrix -> Matrix
 237 transposeM (M (Q e11  e12  e13  e14)
 238               (Q e21  e22  e23  e24)
 239               (Q e31  e32  e33  e34)
 240               (Q e41  e42  e43  e44)) = (M (Q e11  e21  e31  e41)
 241                                            (Q e12  e22  e32  e42)
 242                                            (Q e13  e23  e33  e43)
 243                                            (Q e14  e24  e34  e44))
 244
 245
 246 --multMM m1 m2 = [map (dot4 row) (transpose m2) | row <- m1]
 247
 248 --{-# INLINE multMV #-}
 249 multMV :: Matrix -> Vector -> Vector
 250 multMV m v = quad_to_vector (multMQ m (vector_to_quad v))
 251
 252 --{-# INLINE multMP #-}
 253 multMP :: Matrix -> Point -> Point
 254 multMP m p = quad_to_point (multMQ m (point_to_quad p))
 255
 256 -- mat vec = map (dot4 vec) mat
 257
 258 {-# INLINE multMQ #-}
 259 multMQ :: Matrix -> Quad -> Quad
 260 multMQ (M q1 q2 q3 q4) q
 261        = Q (dot4 q q1)
 262            (dot4 q q2)
 263            (dot4 q q3)
 264            (dot4 q q4)
 265
 266 {-# INLINE multMR #-}
 267 multMR :: Matrix -> Ray -> Ray
 268 multMR m (r, v) = (multMP m r, multMV m v)
 269
 270 white :: Color
 271 white = C 1 1 1
 272 black :: Color
 273 black = C 0 0 0
 274
 275 addCC :: Color -> Color -> Color
 276 addCC (C a b c) (C d e f) = C (a+d) (b+e) (c+f)
 277
 278 subCC :: Color -> Color -> Color
 279 subCC (C a b c) (C d e f) = C (a-d) (b-e) (c-f)
 280
 281 sumCC :: [Color] -> Color
 282 sumCC = foldr addCC black
 283
 284 multCC :: Color -> Color -> Color
 285 multCC (C a b c) (C d e f) = C (a*d) (b*e) (c*f)
 286
 287 multSC :: Double -> Color -> Color
 288 multSC k       (C a b c) = C (a*k) (b*k) (c*k)
 289
 290 nearC :: Color -> Color -> Bool
 291 nearC (C a b c) (C d e f) = a `near` d && b `near` e && c `near` f
 292
 293 offsetToPoint :: Ray -> Double -> Point
 294 offsetToPoint (r,v) i = r `addPV` (i `multSV` v)
 295
 296 --
 297
 298 epsilon, inf :: Double      -- aproximate zero and infinity
 299 epsilon = 1.0e-10
 300 inf = 1.0e20
 301
 302 nonZero :: Double -> Double         -- Use before a division. It makes definitions
 303 nonZero x | x > epsilon  = x        -- more complete and I bet the errors that get
 304           | x < -epsilon = x        -- introduced will be undetectable if epsilon
 305           | otherwise    = epsilon  -- is small enough
 306
 307
 308 eqEps x y = abs (x-y) < epsilon
 309 near = eqEps
 310
 311 clampf :: Double -> Double
 312 clampf p | p < 0 = 0
 313          | p > 1 = 1
 314          | True  = p