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 module Main where 39 40 import System 41 42 epsilon, infinity :: Double 43 epsilon = 1.0e-6 44 infinity = 1.0e+20 45 46 47 48 49 50 type Vector = (Double, Double, Double) 51 52 vecadd, vecsub, vecmult :: Vector -> Vector -> Vector 53 vecadd (x1,y1,z1) (x2,y2,z2) = (x1+x2, y1+y2, z1+z2) 54 vecsub (x1,y1,z1) (x2,y2,z2) = (x1-x2, y1-y2, z1-z2) 55 vecmult (x1,y1,z1) (x2,y2,z2) = (x1*x2, y1*y2, z1*z2) 56 57 vecsum :: [Vector] -> Vector 58 vecsum = foldr vecadd (0,0,0) 59 60 vecnorm :: Vector -> (Vector, Double) 61 vecnorm (x,y,z) = ((x/len, y/len, z/len), len) 62 where len = sqrt (x^2 + y^2 + z^2) 63 64 vecscale :: Vector -> Double -> Vector 65 vecscale (x,y,z) a = (a*x, a*y, a*z) 66 67 vecdot :: Vector -> Vector -> Double 68 vecdot (x1,y1,z1) (x2,y2,z2) = x1*x2 + y1*y2 + z1*z2 69 70 veccross :: Vector -> Vector -> Vector 71 veccross (x1,y1,z1) (x2,y2,z2) = (y1*z2-y2*z1, z1*x2-z2*x1, x1*y2-x2*y1) 72 73 is_zerovector :: Vector -> Bool 74 is_zerovector (x,y,z) = x<epsilon && y<epsilon && z<epsilon 75 76 77 78 79 data Light = Directional Vector Vector -- direction, colour 80 | Point Vector Vector -- position, colour 81 82 lightpos :: Light -> Vector 83 lightpos (Point pos col) = pos 84 85 lightdir :: Light -> Vector 86 lightdir (Directional dir col) = fst (vecnorm dir) 87 88 lightcolour :: Light -> Vector 89 lightcolour (Directional dir col) = col 90 lightcolour (Point pos col) = col 91 92 93 data Surfspec = Ambient Vector -- all but specpow default to zero 94 | Diffuse Vector 95 | Specular Vector 96 | Specpow Double -- default 8 97 | Reflect Double 98 | Transmit Double 99 | Refract Double -- default 1, like air == no refraction 100 | Body Vector -- body colour, default (1, 1, 1) 101 102 103 ambientsurf :: [Surfspec] -> Vector 104 ambientsurf ss = head ([ s | Ambient s <- ss] ++ [(0,0,0)]) 105 106 diffusesurf :: [Surfspec] -> Vector 107 diffusesurf ss = head ([ s | Diffuse s <- ss] ++ [(0,0,0)]) 108 109 specularsurf :: [Surfspec] -> Vector 110 specularsurf ss = head ([ s | Specular s <- ss] ++ [(0,0,0)]) 111 112 specpowsurf :: [Surfspec] -> Double 113 specpowsurf ss = head ([ s | Specpow s <- ss] ++ [8]) 114 115 reflectsurf :: [Surfspec] -> Double 116 reflectsurf ss = head ([ s | Reflect s <- ss] ++ [0]) 117 118 transmitsurf :: [Surfspec] -> Double 119 transmitsurf ss = head ([ s | Transmit s <- ss] ++ [0]) 120 121 refractsurf :: [Surfspec] -> Double 122 refractsurf ss = head ([ s | Refract s <- ss] ++ [1]) 123 124 bodysurf :: [Surfspec] -> Vector 125 bodysurf ss = head ([ s | Body s <- ss] ++ [(1,1,1)]) 126 127 128 data Sphere = Sphere Vector Double [Surfspec] -- pos, radius, surface type 129 130 spheresurf :: Sphere -> [Surfspec] 131 spheresurf (Sphere pos rad surf) = surf 132 133 134 135 136 137 lookat, vup :: Vector 138 lookat = (0,0,0) 139 vup = (0,0,1) 140 141 fov :: Double 142 fov = 45 143 144 world :: [Sphere] 145 world = testspheres 146 147 redsurf :: [Surfspec] 148 redsurf = [Ambient (0.1, 0, 0), Diffuse (0.3, 0, 0), 149 Specular (0.8, 0.4, 0.4), Transmit 0.7] 150 151 greensurf :: [Surfspec] 152 greensurf = [Ambient (0, 0.1, 0), Diffuse (0, 0.3, 0), 153 Specular (0.4, 0.8, 0.4)] 154 155 bluesurf :: [Surfspec] 156 bluesurf = [Ambient (0, 0, 0.1), Diffuse (0, 0, 0.3), 157 Specular (0.4, 0.4, 0.8)] 158 159 160 161 s2 :: [Surfspec] 162 s2 = [Ambient (0.035, 0.0325, 0.025), Diffuse (0.5, 0.45, 0.35), 163 Specular (0.8, 0.8, 0.8), Specpow 3.0, Reflect 0.5] 164 165 testspheres :: [Sphere] 166 testspheres = [(Sphere (0,0,0) 0.5 s2), 167 (Sphere (0.272166,0.272166,0.544331) 0.166667 s2), 168 (Sphere (0.643951,0.172546,0) 0.166667 s2), 169 (Sphere (0.172546,0.643951,0) 0.166667 s2), 170 (Sphere (-0.371785,0.0996195,0.544331) 0.166667 s2), 171 (Sphere (-0.471405,0.471405,0) 0.166667 s2), 172 (Sphere (-0.643951,-0.172546,0) 0.166667 s2), 173 (Sphere (0.0996195,-0.371785,0.544331) 0.166667 s2), 174 (Sphere (-0.172546,-0.643951,0) 0.166667 s2), 175 (Sphere (0.471405,-0.471405,0) 0.166667 s2)] 176 177 178 testlights :: [Light] 179 testlights = [Point (4, 3, 2) (0.288675,0.288675,0.288675), 180 Point (1, -4, 4) (0.288675,0.288675,0.288675), 181 Point (-3,1,5) (0.288675,0.288675,0.288675)] 182 183 lookfrom, background :: Vector 184 lookfrom = (2.1, 1.3, 1.7) 185 background = (0.078, 0.361, 0.753) 186 187 188 189 190 191 ray :: Int -> [((Int, Int),Vector)] 192 ray winsize = [ ((i,j), f i j) | i<-[0..winsize-1], j<-[0..winsize-1]] 193 where 194 lights = testlights 195 (firstray, scrnx, scrny) = camparams lookfrom lookat vup fov (fromIntegral winsize) 196 f i j = tracepixel world lights (fromIntegral i) (fromIntegral j) firstray scrnx scrny 197 198 199 dtor :: Double -> Double 200 dtor x = x*pi / 180 201 202 203 camparams :: Vector -> Vector -> Vector -> Double -> Double 204 -> (Vector, Vector, Vector) 205 camparams lookfrom lookat vup fov winsize = (firstray, scrnx, scrny) 206 where 207 initfirstray = vecsub lookat lookfrom -- pre-normalized! 208 (lookdir, dist) = vecnorm initfirstray 209 (scrni, _) = vecnorm (veccross lookdir vup) 210 (scrnj, _) = vecnorm (veccross scrni lookdir) 211 xfov = fov 212 yfov = fov 213 xwinsize = winsize -- for now, square window 214 ywinsize = winsize 215 magx = 2 * dist * (tan (dtor (xfov/2))) / xwinsize 216 magy = 2 * dist * (tan (dtor (yfov/2))) / ywinsize 217 scrnx = vecscale scrni magx 218 scrny = vecscale scrnj magy 219 firstray = vecsub initfirstray 220 (vecadd (vecscale scrnx (0.5 * xwinsize)) 221 (vecscale scrny (0.5 * ywinsize))) 222 223 224 225 tracepixel :: [Sphere] -> [Light] -> Double -> Double -> Vector -> Vector 226 -> Vector -> Vector 227 tracepixel spheres lights x y firstray scrnx scrny 228 = if hit 229 then shade lights sp pos dir dist (1,1,1) 230 else background 231 where 232 pos = lookfrom 233 (dir, _) = vecnorm (vecadd (vecadd firstray (vecscale scrnx x)) 234 (vecscale scrny y)) 235 (hit, dist, sp) = trace spheres pos dir -- pick first intersection 236 237 238 239 240 241 trace :: [Sphere] -> Vector -> Vector -> (Bool,Double,Sphere) 242 trace spheres pos dir 243 = if (null dists) 244 then (False, infinity, head spheres) -- missed all 245 else (True, mindist, sp) -- pick the smallest one 246 where 247 (mindist, sp) = foldr f (head dists) (tail dists) 248 f (d1,s1) (d2,s2) | d1<d2 = (d1,s1) 249 | otherwise = (d2,s2) 250 -- make a list of the distances to intersection for each hit object 251 sphmap [] = [] 252 sphmap (x:xs) = if is_hit 253 then (where_hit, x):sphmap xs 254 else sphmap xs 255 where (is_hit, where_hit) = sphereintersect pos dir x 256 dists = sphmap spheres 257 258 259 260 261 262 263 264 265 266 267 268 shade :: [Light] -> Sphere -> Vector -> Vector -> Double -> Vector -> Vector 269 shade lights sp lookpos dir dist contrib = rcol 270 where 271 hitpos = vecadd lookpos (vecscale dir dist) 272 ambientlight = (1, 1, 1) -- full contribution as default 273 surf = spheresurf sp 274 amb = vecmult ambientlight (ambientsurf surf) 275 norm = spherenormal hitpos sp 276 refl = vecadd dir (vecscale norm (-2*(vecdot dir norm))) 277 -- diff is diffuse and specular contribution 278 diff = vecsum (map (\l->lightray l hitpos norm refl surf) lights) 279 transmitted = transmitsurf surf 280 simple = vecadd amb diff 281 -- calculate transmitted ray; it adds onto "simple" 282 trintensity = vecscale (bodysurf surf) transmitted 283 (is_tir, trcol) = if transmitted < epsilon 284 then (False, simple) 285 else transmitray lights simple hitpos dir 286 index trintensity contrib norm 287 where index = refractsurf surf 288 -- reflected ray; in case of TIR, add transmitted component 289 reflsurf = vecscale (specularsurf surf) (reflectsurf surf) 290 reflectiv = if (is_tir) 291 then (vecadd trintensity reflsurf) 292 else reflsurf 293 rcol = if is_zerovector reflectiv 294 then trcol 295 else reflectray hitpos refl lights reflectiv contrib trcol 296 297 298 299 300 301 transmitray :: [Light] -> Vector -> Vector -> Vector -> Double -> Vector 302 -> Vector -> Vector -> (Bool, Vector) 303 transmitray lights colour pos dir index intens contrib norm 304 = if is_zerovector newcontrib 305 then (False, colour) -- cutoff 306 else (False, vecadd (vecmult newcol intens) colour) 307 where 308 newcontrib = vecmult intens contrib 309 (is_tir, newdir) = refractray index dir norm 310 nearpos = vecadd pos (vecscale newdir epsilon) 311 (is_hit, dist, sp) = trace world nearpos newdir 312 newcol | is_hit = shade lights sp nearpos newdir dist newcontrib 313 | otherwise = background 314 315 316 317 318 reflectray :: Vector -> Vector -> [Light] -> Vector -> Vector -> Vector 319 -> Vector 320 reflectray pos newdir lights intens contrib colour 321 = if is_zerovector newcontrib 322 then colour 323 else vecadd colour (vecmult newcol intens) 324 where 325 newcontrib = vecmult intens contrib 326 nearpos = vecadd pos (vecscale newdir epsilon) 327 (is_hit, dist, sp) = trace world nearpos newdir 328 newcol = if is_hit 329 then shade lights sp nearpos newdir dist newcontrib 330 else background 331 332 333 334 335 336 337 refractray :: Double -> Vector -> Vector -> (Bool,Vector) 338 refractray newindex olddir innorm 339 = if disc < 0 340 then (True, (0,0,0)) -- total internal reflection 341 else (False, vecadd (vecscale norm t) (vecscale olddir nr)) 342 where 343 dotp = -(vecdot olddir innorm) 344 (norm, k, nr) = if dotp < 0 345 then (vecscale innorm (-1), -dotp, 1/newindex) 346 else (innorm, dotp, newindex) -- trans. only with air 347 disc = 1 - nr*nr*(1-k*k) 348 t = nr * k - (sqrt disc) 349 350 351 352 353 354 355 356 lightray :: Light -> Vector -> Vector -> Vector -> [Surfspec] -> Vector 357 lightray l pos norm refl surf = 358 let 359 (ldir, dist) = lightdirection l pos 360 cosangle = vecdot ldir norm -- lightray is this far off normal 361 (is_inshadow, lcolour) = shadowed pos ldir (lightcolour l) 362 in 363 if is_inshadow then (0,0,0) 364 else 365 let 366 diff = diffusesurf surf 367 spow = specpowsurf surf -- assumed trans is same as refl 368 in 369 if (cosangle <= 0) then -- opposite side 370 let bodycol = bodysurf surf 371 cosalpha = -(vecdot refl ldir) 372 diffcont = vecmult (vecscale diff (-cosangle)) lcolour 373 speccont = if cosalpha <= 0 then (0,0,0) 374 else vecmult (vecscale bodycol (cosalpha**spow)) lcolour 375 in vecadd diffcont speccont 376 else 377 let spec = specularsurf surf 378 cosalpha = vecdot refl ldir 379 -- this far off refl ray (for spec) 380 diffcont = vecmult (vecscale diff cosangle) lcolour 381 speccont | cosalpha <= 0 = (0,0,0) 382 | otherwise = vecmult (vecscale spec (cosalpha**spow)) lcolour 383 in vecadd diffcont speccont 384 385 386 387 388 lightdirection :: Light -> Vector -> (Vector, Double) 389 lightdirection (Directional dir col) pt = (fst (vecnorm dir), infinity) 390 lightdirection (Point pos col) pt = vecnorm (vecsub pos pt) 391 392 393 394 shadowed :: Vector -> Vector -> a -> (Bool,a) 395 shadowed pos dir lcolour = if not is_hit 396 then (False, lcolour) 397 else (True, lcolour) -- for now 398 where 399 (is_hit, dist, sp) = trace world (vecadd pos (vecscale dir epsilon)) dir 400 401 402 403 404 405 406 407 408 sphereintersect :: Vector -> Vector -> Sphere -> (Bool,Double) 409 sphereintersect pos dir sp = if disc < 0 410 then (False, 0) -- imaginary solns only 411 else if slo < 0 412 then if shi < 0 413 then (False, 0) 414 else (True, shi) 415 else (True, slo) 416 where 417 slo = -bm - (sqrt disc) 418 shi = -bm + (sqrt disc) 419 Sphere spos rad _ = sp 420 m = vecsub pos spos -- x-centre 421 m2 = vecdot m m -- (x-centre).(x-centre) 422 bm = vecdot m dir -- (x-centre).dir 423 disc = bm * bm - m2 + rad * rad -- discriminant 424 425 426 427 spherenormal :: Vector-> Sphere -> Vector 428 spherenormal pos sp = vecscale (vecsub pos spos) (1/rad) 429 where 430 Sphere spos rad _ = sp 431 432 433 434 435 436 437 438 439 440 441 main :: IO () 442 main = getArgs >>= \[winsize_string] -> 443 run (read winsize_string) 444 445 run :: Int -> IO () 446 run winsize = ppm winsize (ray winsize) 447 448 449 450 451 ppm :: Int -> [((Int, Int),Vector)] -> IO () 452 ppm winsize matrix = do putStrLn "P3" 453 putStr (show winsize); putStr " " 454 putStrLn (show winsize) 455 putStrLn "255" 456 pixels matrix 457 458 459 pixels :: [((Int, Int),Vector)] -> IO () 460 pixels [] = putStr "" 461 pixels ((point,colour):ps) = do rbg colour 462 putStrLn "" 463 pixels ps 464 465 466 rbg :: Vector -> IO () 467 rbg (r,g,b) = do putStr (eight_bit r); putStr " " 468 putStr (eight_bit g); putStr " " 469 putStr (eight_bit b) 470 where eight_bit = show . round . (255*) 471 472