1votos

Dibujar triangulo de sierpinski en Haskell

por josejuan hace 5 años

Hacen ya más de 16 años de mi primer raytracer, aquí uno muy simple para renderizar la alfombra de sierpinski. Paraleliza sobre N procesadores, pero necesita optimización (sólo paraleliza como el 40% de los recursos).

Gráfica el triangulo de sierpinski.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
{- 
 
  Es un pequeño raytracer, el código que ves es el módulo 
  que lo contiene, una función "main" para la línea de 
  comandos podría ser una como: 
 
    main = do 
      (configFile:outPutPngFile:_) <- getArgs 
      config <- readFile configFile >>= return . read 
      writePng outPutPngFile $ renderSierpinskiBox config 
 
  Así, con un archivo de configuración (de escena) como éste 
  (quitar líneas de comentario) 
 
  -- sierpinski.3d 
  Configuration { 
 
      -- resolución imagen de salida 
      resX = 700, 
      resY = 700, 
 
      -- cámara 
      cameraPos = Point3d {p3x = 2.6, p3y = 2.5, p3z = 3.0},  -- posición 
      cameraObj = Point3d {p3x = 0.0, p3y = 0.0, p3z = 0.0},  -- punto al que mirar 
      cameraUp  = Point3d {p3x = 0.0, p3y = -1.0, p3z = 0.0}, -- cenit de la cámara 
      cameraFov = 2.0,                                        -- apertura 
 
      -- se pueden definir todas las luces que se quieran 
      lights = [ 
          Light { 
              lightPos = Point3d {p3x = 2.0, p3y = 3.0, p3z = 1.5}, -- posición de la luz 
              lightCol = Point3d {p3x = 1.0, p3y = 0.5, p3z = 0.5}, -- color de la luz 
              lightAtK = 1,                                         -- atenuación 1/(k + c * d^2) 
              lightAtC = 0.25, 
              lightSpp = 2.25                                       -- factor especular 
          }, 
          Light { 
              lightPos = Point3d {p3x = 3.0, p3y = 2.0, p3z = 2.5}, 
              lightCol = Point3d {p3x = 0.5, p3y = 1.0, p3z = 0.5}, 
              lightAtK = 1, 
              lightAtC = 0.25, 
              lightSpp = 2.25 
          }, 
          Light { 
              lightPos = Point3d {p3x = 0.5, p3y = 0.5, p3z = 0.5}, 
              lightCol = Point3d {p3x = 0.5, p3y = 0.5, p3z = 1.0}, 
              lightAtK = 2, 
              lightAtC = 0.5, 
              lightSpp = 2.25 
      ], 
 
      minHeight = 0.01,                                             -- mínima altura para recursión sierpinski 
      backColor = Point3d{p3x = 0, p3y = 0, p3z = 0},               -- color de fondo de escena 
      matColor = Point3d{p3x = 0.5, p3y = 0.5, p3z = 1},            -- color del material (no se usa) 
      box = Box { 
          minCorner = Point3d {p3x = 0, p3y = 0, p3z = 0},          -- esquina inferior del cubo a renderizar 
          maxCorner = Point3d {p3x = 1, p3y = 1, p3z = 1}           -- la superior 
 
  Así un ejemplo de ejecución sería 
 
  $ sierpinski.exe sierpinski.3d sierpinski.png 
 
  Y el resultado es el que puedes ver (si sigue ahí) en los comentarios. 
 
-} 
 
{-# LANGUAGE StandaloneDeriving, TupleSections, BangPatterns #-} 
module SierpinskiBox where 
 
import Data.Maybe 
import Data.List 
import Control.Monad 
import Codec.Picture 
import qualified Data.Map as M 
import Control.Parallel 
import Control.Parallel.Strategies 
 
epsilon = 1e-6 
 
data Point3d = Point3d { p3x :: Double 
                       , p3y :: Double 
                       , p3z :: Double 
                       } deriving (Read, Show, Eq, Ord) 
 
data Ray = Ray { origin    :: Point3d 
               , direction :: Point3d 
               } deriving (Read, Show) 
 
data Box = Box { minCorner :: Point3d 
               , maxCorner :: Point3d 
               } deriving (Read, Show, Eq, Ord) 
 
data Collision = Collision { distance :: Double 
                           , position :: Point3d 
                           , normal   :: Normal3d 
                           } deriving (Read, Show, Eq, Ord) 
 
data Light = Light { lightPos :: Point3d 
                   , lightCol :: ColorRGB 
                   , lightAtK :: Double 
                   , lightAtC :: Double 
                   , lightSpp :: Double 
                   } deriving (Read, Show, Eq, Ord) 
 
data Configuration = Configuration { resX      :: Int 
                                   , resY      :: Int 
                                   , cameraPos :: Point3d 
                                   , cameraObj :: Point3d 
                                   , cameraUp  :: Point3d 
                                   , cameraFov :: Double 
                                   , lights    :: [Light] 
                                   , minHeight :: Double 
                                   , backColor :: ColorRGB 
                                   , matColor  :: ColorRGB 
                                   , box       :: Box 
                                   } deriving (Read, Show, Eq) 
 
type Normal3d = Point3d 
type ColorRGB = Point3d 
 
colorToPixel :: ColorRGB -> PixelRGB8 
colorToPixel (Point3d r g b) = PixelRGB8 (t r) (t g) (t b) where t c | c < 0      = 0 
                                                                     | c >= 254.5 = 255 
                                                                     | True       = round (255 * c) 
 
computeLightColor :: Configuration -> Ray -> Collision -> Light -> ColorRGB 
computeLightColor c (Ray o d) (Collision _ p n) (Light l i ak ac sp) = 
  if phi <= 0 || shadow sc then Point3d 0 0 0 else ((phi + phi2) / (ak + ac * ld * ld)) .* i 
  where lp' = l -. p 
        ld = norm lp' 
        lp = lp' ./ ld 
        phi = lp ^. n 
 
        -- diffuse 
        sc = sierpinskiCollision (minHeight c) (box c) (Ray l ((-1) .* lp)) 
        shadow Nothing = False 
        shadow (Just (Collision t _ _)) = t < ld - epsilon 
 
        -- specular 
        v = (2 * (n ^. d)) .* n 
        r = d -. v 
        phi2 = (r ^. n)**sp 
 
pixel3DColor :: Configuration -> Ray -> Maybe Collision -> PixelRGB8 
pixel3DColor c _ Nothing = colorToPixel $ backColor c 
pixel3DColor c r !(Just o) = colorToPixel $ foldl1 (+.) $ map (computeLightColor c r o) (lights c) 
 
renderSierpinskiBox :: Configuration -> Image PixelRGB8 
renderSierpinskiBox c@(Configuration rx ry cO cD cU cF _ minH _ _ box) = do 
  let rays = cameraRays rx ry cO cD cU cF 
      mpix = M.fromList $ (parMap rseq) (\(x, y, r) -> (key x y, pixel3DColor c r $ sierpinskiCollision minH box r)) rays 
      key x y = 10000 * y + x 
  rays `seq` mpix `seq` generateImage (\x y -> mpix M.! (key x y)) rx ry 
 
-- Una forma rápida sería generar recursivamente todos los subcubos y usar directamente 
-- `rayBoxesIntersection` pero es mucho más eficiente testear colisiones en jerarquía 
-- como si fuera un octree (en realidad es un `nonotree`) XD XD XD 
sierpinskiCollision :: Double -> Box -> Ray -> Maybe Collision 
sierpinskiCollision !minHeight !box !ray@(Ray o d) = 
  case rayBoxesIntersection ray $ splitSierpinskiBox box of 
    Nothing -> Nothing 
    Just (t, box'@(Box a b)) -> if p3y (b -. a) < minHeight 
                                  then Just $ Collision t p n 
                                  else sierpinskiCollision minHeight box' ray 
                                where p = o +. (t .* d) 
                                      n = boxNormal box' p 
 
splitSierpinskiBox :: Box -> [Box] 
splitSierpinskiBox !(Box a b) = catMaybes [cube x y z | x <- [1..3], y <- [1..3], z <- [1..3]] 
  where cube 2 _ 2 = Nothing 
        cube 2 2 _ = Nothing 
        cube _ 2 2 = Nothing 
        cube x y z = Just $ Box (a +. (p' *. s)) (a +. (p *. s)) 
                     where p = Point3d x y z 
                           p' = p -. Point3d 1 1 1 
                           s = (b -. a) ./ 3 
 
rayBoxesIntersection :: Ray -> [Box] -> Maybe (Double, Box) 
rayBoxesIntersection !ray = listToMaybe.sort.mapMaybe (\b -> rayBoxIntersection ray b >>= return.(,b)) 
 
rayBoxIntersection :: Ray -> Box -> Maybe Double 
rayBoxIntersection !(Ray (Point3d ox oy oz) (Point3d dx dy dz)) !(Box (Point3d ax ay az) (Point3d bx by bz)) = 
    if tmax < 0 || tmin > tmax then Nothing else Just tmin 
    where t1 = (ax - ox) / dx 
          t2 = (bx - ox) / dx 
          t3 = (ay - oy) / dy 
          t4 = (by - oy) / dy 
          t5 = (az - oz) / dz 
          t6 = (bz - oz) / dz 
          tmin = max (max (min t1 t2) (min t3 t4)) (min t5 t6) 
          tmax = min (min (max t1 t2) (max t3 t4)) (max t5 t6) 
 
boxNormal :: Box -> Point3d -> Normal3d 
boxNormal !(Box (Point3d ax ay az) (Point3d bx by bz)) !(Point3d x y z) = 
  snd $ head $ sortBy ((.fst).compare.fst) [ (abs (y - ay), Point3d    0 (-1)    0 ) 
                                           , (abs (y - by), Point3d    0    1    0 ) 
                                           , (abs (x - ax), Point3d (-1)    0    0 ) 
                                           , (abs (x - bx), Point3d    1    0    0 ) 
                                           , (abs (z - az), Point3d    0    0 (-1) ) 
                                           , (abs (z - bz), Point3d    0    0    1 ) 
 
 
cameraRays :: Int -> Int -> Point3d -> Point3d -> Point3d -> Double -> [(Int, Int, Ray)] 
cameraRays resx resy poscam dstcam upvec fov = [ray x y | x <- [0 .. resx - 1], y <- [0 .. resy - 1]] 
  where ray x y = (x, y, ray' (normRS x resx) (normRS y resy)) 
        ray' x y = Ray poscam $ normalized (dir +. (x .* hor) +. (y .* ver)) 
        tovec = dstcam -. poscam 
        hor = normalized $ upvec %. tovec 
        ver = normalized $ tovec %. hor 
        dir = fov .* normalized tovec 
        normRS c resc = (c' - resc'') / resc' where c' = fromIntegral c 
                                                    resc' = fromIntegral resc 
                                                    resc'' = 0.5 * resc' - 0.5 
 
-- Deberían prohibir crear un espacio vectorial sin las funciones básicas :P (arf!) 
(.*) :: Double -> Point3d -> Point3d 
s .* (Point3d x y z) = Point3d (s * x) (s * y) (s * z) 
 
(./) :: Point3d -> Double -> Point3d 
v ./ s = (1 / s) .* v 
 
(+.) :: Point3d -> Point3d -> Point3d 
(Point3d ax ay az) +. (Point3d bx by bz) = Point3d (ax + bx) (ay + by) (az + bz) 
 
(-.) :: Point3d -> Point3d -> Point3d 
(Point3d ax ay az) -. (Point3d bx by bz) = Point3d (ax - bx) (ay - by) (az - bz) 
 
(%.) :: Point3d -> Point3d -> Point3d 
(Point3d ax ay az) %. (Point3d bx by bz) = Point3d (ay * bz - az * by) (az * bx - ax * bz) (ax * by - ay * bx) 
 
(*.) :: Point3d -> Point3d -> Point3d 
(Point3d ax ay az) *. (Point3d bx by bz) = Point3d (ax * bx) (ay * by) (az * bz) 
 
(^.) :: Point3d -> Point3d -> Double 
a ^. b = x + y + z where Point3d x y z = a *. b 
 
norm2 :: Point3d -> Double 
norm2 p = x + y + z where Point3d x y z = p *. p 
 
norm :: Point3d -> Double 
norm = sqrt . norm2 
 
normalized :: Point3d -> Point3d 
normalized v = v ./ (norm v) 
1 comentario

Comenta la solución

Tienes que identificarte para poder publicar tu comentario.