add Galois' Ray Tracer
[nofib.git] / parallel / gray / 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