comments
[nofib.git] / parallel / gray / Illumination.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 -- Modified to use stdout (for testing)
7
8 {-# LANGUAGE BangPatterns,CPP #-}
9 module Illumination
10 ( Object
11 , Light (..)
12 , light, pointlight, spotlight
13 , render
14 ) where
15
16 import Control.Parallel
17 import Control.Parallel.Strategies (withStrategy, parBuffer, rwhnf)
18
19 import Data.Array
20 import Data.Char(chr)
21 import Data.Maybe
22
23 import Geometry
24 import CSG
25 import Surface
26 import Misc
27
28 type Object = CSG (SurfaceFn Color Double)
29
30 data Cxt = Cxt {ambient::Color, lights::[Light], object::Object, depth::Int}
31 deriving Show
32
33 render :: (Matrix,Matrix) -> Color -> [Light] -> Object -> Int ->
34 Radian -> Int -> Int -> String -> IO ()
35 render (m,m') amb ls obj dep fov wid ht file
36 = do { debugging
37 ; putStrLn (showBitmap' wid ht (parallel pixels))
38 }
39 where
40 #ifdef STRATEGIES_2
41 parallel = parBuffer 100 rwhnf . map (\x -> seqList x `pseq` x)
42 #else
43 parallel = withStrategy (parBuffer 100 rwhnf) . map (\x -> seqList x `pseq` x)
44 #endif
45
46 debugging = return ()
47 {-
48 do { putStrLn (show cxt)
49 ; putStrLn (show (width, delta, aspect, left, top))
50 }
51 -}
52 obj' = transform (m',m) obj
53 ls' = [ transformLight m' l | l <- ls ]
54 pixelA = listArray ((1,1), (ht,wid))
55 [ illumination cxt (start,pixel i j)
56 | j <- take ht [0.5..]
57 , i <- take wid [0.5..] ]
58 antiA = pixelA //
59 [ (ix, superSample ix (pixelA ! ix))
60 | j <- [2 .. ht - 1], i <- [2 .. wid - 1]
61 , let ix = (j, i)
62 , contrast ix pixelA ]
63 pixels = [ [ illumination cxt (start,pixel i j) | i<- take wid [0.5..] ]
64 | j <- take ht [0.5..]
65 ]
66 cxt = Cxt {ambient=amb, lights=ls', object=obj', depth=dep}
67 start = point 0 0 (-1)
68 width = 2 * tan (fov/2)
69 delta = width / fromIntegral wid
70 aspect = fromIntegral ht / fromIntegral wid
71 left = - width / 2
72 top = - left * aspect
73 pixel i j = vector (left + i*delta) (top - j*delta) 1
74
75 superSample (y, x) col = avg $ col:
76 [ illumination cxt (start, pixel (fromIntegral x - 0.5 + xd) (fromIntegral y - 0.5 + yd))
77 | (xd, yd) <- [(-0.333, 0.0), (0.333, 0.0), (0.0, -0.333), (0.0, 0.333)]
78 ]
79
80 seqList :: [a] -> ()
81 seqList [] = ()
82 seqList (x:xs) = x `pseq` seqList xs
83
84 avg cs = divN (fromIntegral (length cs)) (uncolor (sumCC cs))
85 where divN n (r,g,b) = color (r / n) (g / n) (b / n)
86
87 contrast :: (Int, Int) -> Array (Int, Int) Color -> Bool
88 contrast (x, y) arr = any diffMax [ subCC cur (arr ! (x + xd, y + yd))
89 | xd <- [-1, 1], yd <- [-1, 1]
90 ]
91 where cur = arr ! (x, y)
92 diffMax col = (abs r) > 0.25 || (abs g) > 0.2 || (abs b) > 0.4
93 where
94 (r,g,b) = uncolor col
95
96
97 illumination :: Cxt -> Ray -> Color
98 illumination cxt (r,v)
99 | depth cxt <= 0 = black
100 | otherwise = case castRay (r,v) (object cxt) of
101 Nothing -> black
102 Just info -> illum (cxt{depth=(depth cxt)-1}) info v
103
104 illum :: Cxt -> (Point,Vector,Properties Color Double) -> Vector -> Color
105 illum cxt (pos,normV,(col,kd,ks,n)) v
106 = ambTerm `addCC` difTerm `addCC` spcTerm `addCC` recTerm
107 where
108 visibleLights = unobscured pos (object cxt) (lights cxt) normV
109 d = depth cxt
110 amb = ambient cxt
111 newV = subVV v (multSV (2 * dot normV v) normV)
112
113 ambTerm = multSC kd (multCC amb col)
114 difTerm = multSC kd (sumCC [multSC (dot normV lj) (multCC intensity col)
115 |(loc,intensity) <- visibleLights,
116 let lj = normalize ({- pos `subVV` -} loc)])
117 -- ZZ might want to avoid the phong, when you can...
118 spcTerm = multSC ks (sumCC [multSC ((dot normV hj) ** n ) (multCC intensity col)
119 |(loc,intensity) <- visibleLights,
120 -- ZZ note this is specific to the light at infinity
121 let lj = {- pos `subVV` -} normalize loc,
122 let hj = normalize (lj `subVV` normalize v)])
123 recTerm = if recCoeff `nearC` black then black else multCC recCoeff recRay
124 recCoeff = multSC ks col
125 recRay = illumination cxt (pos,newV)
126
127 showBitmapA :: Int -> Int -> Array (Int, Int) Color -> String
128 showBitmapA wid ht arr
129 = header ++ concatMap scaleColor (elems arr)
130 where
131 scaleColor col = [scalePixel r, scalePixel g, scalePixel b]
132 where (r,g,b) = uncolor col
133 header = "P6\n#Galois\n" ++ show wid ++ " " ++ show ht ++ "\n255\n"
134
135 showBitmap :: Int -> Int ->[[Color]] -> String
136 showBitmap wid ht pss
137 -- type of assert | length pss == ht && all (\ ps -> length ps == wid) pss
138 = header ++ concat [[scalePixel r,scalePixel g,scalePixel b]
139 | ps <- pss, (r,g,b) <- map uncolor ps]
140 where
141 header = "P6\n#Galois\n" ++ show wid ++ " " ++ show ht ++ "\n255\n"
142 showBitmap _ _ _ = error "incorrect length of bitmap string"
143
144 scalePixel :: Double -> Char
145 scalePixel p = chr (floor (clampf p * 255))
146
147 showBitmap' :: Int -> Int ->[[Color]] -> String
148 showBitmap' wid ht pss
149 -- type of assert | length pss == ht && all (\ ps -> length ps == wid) pss
150 = header
151 ++ unlines [ unwords [unwords [scalePixel' r,scalePixel' g,scalePixel' b]
152 | (r,g,b) <- map uncolor ps]
153 | ps <- pss ]
154 where
155 header = "P3\n#Galois\n" ++ show wid ++ " " ++ show ht ++ "\n255\n"
156 showBitmap' _ _ _ = error "incorrect length of bitmap string"
157
158 scalePixel' :: Double -> String
159 scalePixel' p = show (floor (clampf p * 255))
160
161 -- Lights
162
163 data Light = Light Vector Color
164 | PointLight Point Color
165 | SpotLight Point Point Color Radian Double
166 deriving Show
167
168 light :: Coords -> Color -> Light
169 light (x,y,z) color =
170 Light (normalize (vector (-x) (-y) (-z))) color
171 pointlight (x,y,z) color =
172 PointLight (point x y z) color
173 spotlight (x,y,z) (p,q,r) col cutoff exp =
174 SpotLight (point x y z) (point p q r) col cutoff exp
175
176 transformLight m (Light v c) = Light (multMV m v) c
177 transformLight m (PointLight p c) = PointLight (multMP m p) c
178 transformLight m (SpotLight p q c r d) = SpotLight (multMP m p) (multMP m q) c r d
179
180 unobscured :: Point -> Object -> [Light] -> Vector -> [(Vector,Color)]
181 unobscured pos obj lights normV = catMaybes (map (unobscure pos obj normV) lights)
182
183 unobscure :: Point -> Object -> Vector -> Light -> Maybe (Vector,Color)
184 unobscure pos obj normV (Light vec color)
185 -- ZZ probably want to make this faster
186 | vec `dot` normV < 0 = Nothing
187 | intersects (pos `addPV` (0.0001 `multSV` vec),vec) obj = Nothing
188 | otherwise = Just (vec,color)
189 unobscure pos obj normV (PointLight pp color)
190 | vec `dot` normV < 0 = Nothing
191 | intersectWithin (pos `addPV` (0.0001 `multSV` (normalize vec)), vec) obj = Nothing
192 | otherwise = Just (vec,is)
193 where vec = pp `subPP` pos
194 is = attenuate vec color
195 unobscure org obj normV (SpotLight pos at color cutoff exp)
196 | vec `dot` normV < 0 = Nothing
197 | intersectWithin (org `addPV` (0.0001 `multSV` (normalize vec)), vec) obj = Nothing
198 | angle > cutoff = Nothing
199 | otherwise = Just (vec, is)
200 where vec = pos `subPP` org
201 vec' = pos `subPP` at
202 angle = acos (normalize vec `dot` (normalize vec'))
203
204 asp = normalize (at `subPP` pos)
205 qsp = normalize (org `subPP` pos)
206 is = attenuate vec (((asp `dot` qsp) ** exp) `multSC` color)
207
208 attenuate :: Vector -> Color -> Color
209 attenuate vec color = (100 / (99 + sq (norm vec))) `multSC` color
210
211 --
212
213 castRay ray p
214 = case intersectRayWithObject ray p of
215 (True, _, _) -> Nothing -- eye is inside
216 (False, [], _) -> Nothing -- eye is inside
217 (False, (0, b, _) : _, _) -> Nothing -- eye is inside
218 (False, (i, False, _) : _, _) -> Nothing -- eye is inside
219 (False, (t, b, (s, p0)) : _, _) ->
220 let (v, prop) = surface s p0 in
221 Just (offsetToPoint ray t, v, prop)
222
223 intersects ray p
224 = case intersectRayWithObject ray p of
225 (True, _, _) -> False
226 (False, [], _) -> False
227 (False, (0, b, _) : _, _) -> False
228 (False, (i, False, _) : _, _) -> False
229 (False, (i, b, _) : _, _) -> True
230
231 intersectWithin :: Ray -> Object -> Bool
232 intersectWithin ray p
233 = case intersectRayWithObject ray p of
234 (True, _, _) -> False -- eye is inside
235 (False, [], _) -> False -- eye is inside
236 (False, (0, b, _) : _, _) -> False -- eye is inside
237 (False, (i, False, _) : _, _) -> False -- eye is inside
238 (False, (t, b, _) : _, _) -> t < 1.0