glome-hs 0.4.1 → 0.5
raw patch · 21 files changed
+1874/−1430 lines, 21 filesdep +binary
Dependencies added: binary
Files
- Bih.hs +284/−0
- Bound.hs +56/−0
- Box.hs +68/−0
- Cone.hs +262/−0
- Csg.hs +103/−0
- Glome.hs +113/−93
- Plane.hs +45/−0
- README +5/−1
- Scene.hs +73/−0
- Solid.hs +380/−1178
- SolidTexture.hs +8/−8
- Spd.hs +13/−15
- Sphere.hs +75/−0
- TestScene.hs +76/−70
- Tex.hs +47/−0
- Trace.hs +63/−36
- Triangle.hs +115/−0
- Vec.hs +83/−24
- glome-hs.cabal +3/−3
- make +1/−1
- run +1/−1
+ Bih.hs view
@@ -0,0 +1,284 @@+module Bih (bih) where+import Vec+import Solid+import Data.List hiding (group) -- for "partition"++-- Bounding Interval Heirarchy+-- http://en.wikipedia.org/wiki/Bounding_interval_hierarchy++data Bih = Bih {bihbb :: Bbox, bihroot :: BihNode} deriving Show+data BihNode = BihLeaf !SolidItem + | BihBranch {lmax :: !Flt, rmin :: !Flt, ax :: !Int, + l :: BihNode, r :: BihNode} deriving Show++-- bih construction+build_leaf :: [(Bbox, SolidItem)] -> BihNode+build_leaf objs =+ BihLeaf (group (map snd objs))++-- tuning parameter that controls threshold for separating+-- large objects from small objects instead of usual left/right+-- sorting +max_bih_sa = 0.3 :: Flt++build_rec :: [(Bbox,SolidItem)] -> Bbox -> Bbox -> Int -> BihNode+build_rec objs nodebox splitbox depth = + -- if (null objs) || (null $ tail objs) || + -- (null $ tail $ tail objs)+ if length objs < 2+ then build_leaf objs+ else+ let (Bbox nodeboxp1 nodeboxp2) = nodebox+ (Bbox splitboxp1 splitboxp2) = splitbox+ axis = vmaxaxis (vsub splitboxp2 splitboxp1)+ bbmin = va splitboxp1 axis+ bbmax = va splitboxp2 axis+ candidate = (bbmin + bbmax) * 0.5+ in+ if candidate > (va nodeboxp2 axis) then+ build_rec objs nodebox + (Bbox splitboxp1 (vset splitboxp2 axis candidate)) + depth+ else+ if candidate < (va nodeboxp1 axis) then+ build_rec objs nodebox (+ Bbox (vset splitboxp1 axis candidate) splitboxp2) + depth+ else+ -- not sure if this is a big win+ let nbsa = bbsa nodebox+ (big,small) = partition (\ (bb,_) -> + (bbsa bb) > (nbsa * max_bih_sa)) objs+ in + if (not $ null big) && ((length big) < ((length small)*2))+ then (BihBranch (va nodeboxp2 0) (va nodeboxp1 0) 0+ (build_rec big nodebox splitbox (depth+1))+ (build_rec small nodebox splitbox (depth+1)) )+ else+ let (l,r) = partition (\((Bbox bbp1 bbp2),_)-> + (((va bbp1 axis)+(va bbp2 axis))*0.5) + < candidate ) objs+ lmax = foldl fmax (-infinity) (map (\((Bbox _ p2),_) -> va p2 axis) l)+ rmin = foldl fmin infinity (map (\((Bbox p1 _),_) -> va p1 axis) r)+ (lsplit,rsplit) = bbsplit splitbox axis candidate+ lnb = (Bbox nodeboxp1 (vset nodeboxp2 axis lmax))+ rnb = (Bbox (vset nodeboxp1 axis rmin) nodeboxp2)+ in+ -- stop if there's no progress being made+ if ((null l) && (rmin <= bbmin)) ||+ ((null r) && (lmax >= bbmax))+ then build_leaf objs+ else+ (BihBranch (lmax+delta) (rmin-delta) axis+ (build_rec l lnb lsplit (depth+1))+ (build_rec r rnb rsplit (depth+1)) )++bih :: [SolidItem] -> SolidItem+bih [] = SolidItem Void+-- bih (sld:[]) = sld -- sometimes we'd like to be able to use a+ -- single object bih just for its aabb+bih slds =+ let objs = map (\x -> ((bound x),x)) (flatten_group slds)+ bb = foldl bbjoin empty_bbox (map (\(b,_)->b) objs)+ root = build_rec objs bb bb 0+ (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bb+ in+ if p1x == (-infinity) || p1y == (-infinity) || p1z == (-infinity) ||+ p2x == infinity || p2y == infinity || p2z == infinity+ then+ error $ "bih: infinite bounding box " ++ (show objs)+ else+ SolidItem (Bih bb root)++rayint_bih :: Bih -> Ray -> Flt -> Texture -> Rayint +rayint_bih (Bih bb root) r d t =+ let Ray orig dir = r+ dir_rcp = vrcp dir+ Interval near far = bbclip r bb+ traverse (BihLeaf s) near far = rayint s r (fmin d far) t+ traverse (BihBranch lsplit rsplit axis l r) near far =+ let dirr = va dir_rcp axis+ o = va orig axis+ dl = (lsplit - o) * dirr+ dr = (rsplit - o) * dirr+ in + if near > far + then RayMiss+ else+ if dirr > 0+ then + (nearest+ (if near < dl+ then traverse l near (fmin dl far)+ else RayMiss)+ (if dr < far+ then traverse r (fmax dr near) far+ else RayMiss))+ else+ (nearest+ (if near < dr+ then traverse r near (fmin dr far)+ else RayMiss)+ (if dl < far+ then traverse l (fmax dl near) far+ else RayMiss))++ in+ traverse root near far++-- This is unwieldy, but the performance gains+-- make it worthwhile. By testing 4 rays against +-- each cell, we do 1/4 the memory accesses. ++-- One simplifying assumption we make that adds a +-- little bit of overhead: If one ray hits a cell, +-- we act as though they all do. For that reason,+-- this only works well with coherent rays.++packetint_bih :: Bih -> Ray -> Ray -> Ray -> Ray -> Flt -> Texture -> PacketResult+packetint_bih (Bih bb root) r1 r2 r3 r4 d t =+ let bih = Bih bb root+ Ray orig1 dir1 = r1+ Ray orig2 dir2 = r2+ Ray orig3 dir3 = r3+ Ray orig4 dir4 = r4++ dir_rcp1 = vrcp dir1+ dir_rcp2 = vrcp dir2+ dir_rcp3 = vrcp dir3+ dir_rcp4 = vrcp dir4+ in+ -- We want all the ray components to have+ -- at least the same sign.+ if not $ veqsign dir_rcp1 dir_rcp2 &&+ veqsign dir_rcp1 dir_rcp3 &&+ veqsign dir_rcp1 dir_rcp4+ then+ PacketResult (rayint bih r1 d t)+ (rayint bih r2 d t)+ (rayint bih r3 d t)+ (rayint bih r4 d t)+ else + let Interval near1 far1 = bbclip r1 bb+ Interval near2 far2 = bbclip r2 bb+ Interval near3 far3 = bbclip r3 bb+ Interval near4 far4 = bbclip r4 bb++ near = fmin (fmin near1 near2) (fmin near3 near4)+ far = fmax (fmax far1 far2) (fmax far3 far4)++ traverse (BihLeaf s) near far = packetint s r1 r2 r3 r4 (fmin d far) t+ traverse (BihBranch lsplit rsplit axis l r) near far =+ if near > far + then packetmiss+ else+ let dirr1 = va dir_rcp1 axis+ dirr2 = va dir_rcp2 axis+ dirr3 = va dir_rcp3 axis+ dirr4 = va dir_rcp4 axis+ + o1 = va orig1 axis+ o2 = va orig2 axis+ o3 = va orig3 axis+ o4 = va orig4 axis++ dl1 = (lsplit - o1) * dirr1+ dl2 = (lsplit - o2) * dirr2+ dl3 = (lsplit - o3) * dirr3+ dl4 = (lsplit - o4) * dirr4++ dr1 = (rsplit - o1) * dirr1+ dr2 = (rsplit - o2) * dirr2+ dr3 = (rsplit - o3) * dirr3+ dr4 = (rsplit - o4) * dirr4++ in + if dirr1 > 0 -- true for all, since signs match+ then + let dl = fmax4 dl1 dl2 dl3 dl4+ dr = fmin4 dr1 dr2 dr3 dr4+ in+ (nearest_packetresult+ (if near < dl+ then traverse l near (fmin dl far)+ else packetmiss)+ (if dr < far+ then traverse r (fmax dr near) far+ else packetmiss))+ else+ let dl = fmin4 dl1 dl2 dl3 dl4+ dr = fmax4 dr1 dr2 dr3 dr4+ in+ (nearest_packetresult+ (if near < dr+ then traverse r near (fmin dr far)+ else packetmiss)+ (if dl < far+ then traverse l (fmax dl near) far+ else packetmiss))++ in+ traverse root near far++shadow_bih :: Bih -> Ray -> Flt -> Bool+shadow_bih (Bih bb root) r d =+ let (Ray orig dir) = r+ dir_rcp = vrcp dir+ Interval near far = bbclip r bb+ traverse (BihLeaf s) near far = shadow s r (fmin d far)+ traverse (BihBranch lsplit rsplit axis l r) near far =+ let dirr = va dir_rcp axis+ o = va orig axis+ dl = (lsplit - o) * dirr+ dr = (rsplit - o) * dirr+ in + if near > far + then False+ else+ if dirr > 0+ then+ ((if near < dl+ then traverse l near (fmin dl far)+ else False) + ||+ (if dr < far+ then traverse r (fmax dr near) far+ else False))+ else+ ((if near < dr+ then traverse r near (fmin dr far)+ else False)+ ||+ (if dl < far+ then traverse l (fmax dl near) far+ else False))++ in traverse root near far++inside_bih :: Bih -> Vec -> Bool+inside_bih (Bih (Bbox (Vec x1 y1 z1) (Vec x2 y2 z2)) root) pt =+ let (Vec x y z) = pt+ traverse (BihLeaf s) = inside s pt+ traverse (BihBranch lsplit rsplit axis l r) =+ let o = va pt axis+ in (if o < lsplit+ then (traverse l)+ else False) + ||+ (if o > rsplit + then (traverse r)+ else False)+ in+ (x > x1) && (x < x2) && + (y > y1) && (y < y2) && + (z > z1) && (z < z2) && (traverse root)++bound_bih :: Bih -> Bbox+bound_bih (Bih bb root) = bb++instance Solid Bih where+ rayint = rayint_bih+ packetint = packetint_bih+ shadow = shadow_bih+ inside = inside_bih+ bound = bound_bih
+ Bound.hs view
@@ -0,0 +1,56 @@+module Bound (bound_object) where+import Vec+import Solid++-- Bounding objects: we can use any object as a bounding+-- object for any other object; if a ray misses the+-- bounding object, we can assume it missed the bounded+-- object as well. Unlike bih, setting up bounds is a manual+-- process. It is important that the bounded object is+-- completely inside the bounding object.++-- The bounding object should have a cheaper intersection test than+-- the bounded object for this to be useful.++-- The first SolidItem is the bounding object, the second+-- is the bounded object.+data Bound = Bound SolidItem SolidItem deriving Show++bound_object :: SolidItem -> SolidItem -> SolidItem+bound_object a b = SolidItem $ Bound a b++rayint_bound :: Bound -> Ray -> Flt -> Texture -> Rayint+rayint_bound (Bound sa sb) r d t =+ let (Ray orig _) = r+ in if inside sa orig || shadow sa r d+ then rayint sb r d t+ else RayMiss++shadow_bound :: Bound -> Ray -> Flt -> Bool+shadow_bound (Bound sa sb) r d =+ let (Ray orig _ ) = r+ in if inside sa orig || shadow sa r d+ then shadow sb r d+ else False++inside_bound :: Bound -> Vec -> Bool+inside_bound (Bound sa sb) pt = inside sa pt && inside sb pt++-- if this is too slow, we could just take the bounding box for sa+bound_bound :: Bound -> Bbox+bound_bound (Bound sa sb) = bboverlap (bound sa) (bound sb)++-- remove bounding objects when we flatten transformations+-- (this is so that the accelleration structure can +-- build an automatic bounding hierarchy rather than+-- a manual one)++flatten_transform_bound :: Bound -> SolidItem+flatten_transform_bound (Bound sa sb) = flatten_transform sb++instance Solid Bound where+ rayint = rayint_bound+ shadow = shadow_bound+ inside = inside_bound+ bound = bound_bound+ flatten_transform = flatten_transform_bound
+ Box.hs view
@@ -0,0 +1,68 @@+module Box (box) where+import Vec+import Solid++-- Simple, axis-aligned bounding box defined with two points at opposing corners.++data Box = Box !Bbox deriving Show++box :: Vec -> Vec -> SolidItem+box (Vec x1 y1 z1) (Vec x2 y2 z2) =+ SolidItem (Box (Bbox (Vec (fmin x1 x2) (fmin y1 y2) (fmin z1 z2))+ (Vec (fmax x1 x2) (fmax y1 y2) (fmax z1 z2))))++-- this could be optimized a bit more+rayint_box :: Box -> Ray -> Flt -> Texture -> Rayint+rayint_box (Box (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z))) r d t =+ let (Ray orig dir) = r+ (Vec ox oy oz) = orig+ (Vec dx dy dz) = dir+ dxrcp = 1/dx+ dyrcp = 1/dy+ dzrcp = 1/dz+ Interval inx outx = if dx > 0 + then Interval ((p1x-ox)*dxrcp) ((p2x-ox)*dxrcp)+ else Interval ((p2x-ox)*dxrcp) ((p1x-ox)*dxrcp)+ Interval iny outy = if dy > 0+ then Interval ((p1y-oy)*dyrcp) ((p2y-oy)*dyrcp)+ else Interval ((p2y-oy)*dyrcp) ((p1y-oy)*dyrcp)+ Interval inz outz = if dz > 0+ then Interval ((p1z-oz)*dzrcp) ((p2z-oz)*dzrcp)+ else Interval ((p2z-oz)*dzrcp) ((p1z-oz)*dzrcp)+ lastin = (fmax3 inx iny inz)+ firstout = (fmin3 outx outy outz)+ in if lastin > firstout || firstout < 0 || lastin > d+ then RayMiss+ else + let n = if inx == lastin + then if dx > 0 then nvx else vx+ else if iny == lastin+ then if dy > 0 then nvy else vy+ else if dz > 0 then nvz else vz+ norm = if lastin > 0 then n else vinvert n+ hitdepth = fmax 0 lastin+ in+ RayHit hitdepth (vscaleadd orig dir hitdepth) norm t ++shadow_box :: Box -> Ray -> Flt -> Bool+shadow_box (Box box) r d =+ let Interval near far = bbclip r box + in+ if (near > far) || far <= 0 || far > d+ then False+ else True++inside_box :: Box -> Vec -> Bool+inside_box (Box (Bbox (Vec x1 y1 z1) (Vec x2 y2 z2))) (Vec x y z) =+ x > x1 && x < x2 && + y > y1 && y < y2 && + z > z1 && z < z2++bound_box :: Box -> Bbox+bound_box (Box box) = box++instance Solid Box where+ rayint = rayint_box+ shadow = shadow_box+ inside = inside_box+ bound = bound_box
+ Cone.hs view
@@ -0,0 +1,262 @@+module Cone (disc, cone, cylinder) where+import Vec+import Solid+import Sphere -- for disc bounding box++-- We define "Cone", "Cylinder", and "Disc" in this module.+-- A Cone is really a tapered cylinder with a different radius+-- at each end, though the underlying representation is a+-- clipped cone.++-- We represent Cylinders and Cones as transformations of axis-aligned+-- primitives.++-- Todo: cylinder shadow test++data Disc = Disc !Vec !Vec !Flt deriving Show -- position, normal, r*r+data Cylinder = Cylinder !Flt !Flt !Flt deriving Show -- radius height1 height2+data Cone = Cone !Flt !Flt !Flt !Flt deriving Show -- r clip1 clip2 height++-- CONSTRUCTORS --++disc :: Vec -> Vec -> Flt -> SolidItem+disc pos norm r =+ SolidItem $ Disc pos norm (r*r)++cylinder_z :: Flt -> Flt -> Flt -> SolidItem+cylinder_z r h1 h2 = SolidItem (Cylinder r h1 h2)++cone_z :: Flt -> Flt -> Flt -> Flt -> SolidItem+cone_z r h1 h2 height = SolidItem (Cone r h1 h2 height)++-- construct a general cylinder from p1 to p2 with radius r+cylinder :: Vec -> Vec -> Flt -> SolidItem+cylinder p1 p2 r =+ let axis = vsub p2 p1+ len = vlen axis+ ax1 = vscale axis (1/len)+ (ax2,ax3) = orth ax1 + in transform (cylinder_z r 0 len)+ [ (xyz_to_uvw ax2 ax3 ax1),+ (translate p1) ]+ +-- similar for cone+cone :: Vec -> Flt -> Vec -> Flt -> SolidItem+cone p1 r1 p2 r2 =+ if r1 < r2 + then cone p2 r2 p1 r1+ else if r1-r2 < delta+ then cylinder p1 p2 r2+ else+ let axis = vsub p2 p1+ len = vlen axis+ ax1 = vscale axis (1/len)+ (ax2,ax3) = orth ax1 + height = (r1*len)/(r1-r2) -- distance to end point+ in+ transform (cone_z r1 0 len height)+ [ (xyz_to_uvw ax2 ax3 ax1),+ (translate p1) ] ++rayint_disc :: Disc -> Ray -> Flt -> Texture -> Rayint+rayint_disc (Disc point norm radius_sqr) r d t =+ let (Ray orig dir) = r+ dist = plane_int_dist r point norm + in if dist < 0 || dist > d + then RayMiss+ else let pos = vscaleadd orig dir dist+ offset = vsub pos point+ in + if (vdot offset offset) > radius_sqr+ then RayMiss+ else RayHit dist pos norm t++shadow_disc :: Disc -> Ray -> Flt -> Bool+shadow_disc (Disc point norm radius_sqr) r d =+ let (Ray orig dir) = r+ dist = plane_int_dist r point norm + in if dist < 0 || dist > d + then False+ else let pos = vscaleadd orig dir dist+ offset = vsub pos point+ in + if (vdot offset offset) > radius_sqr+ then False+ else True++bound_disc :: Disc -> Bbox+bound_disc (Disc pos norm rsqr) =+ bound (sphere pos (sqrt rsqr))++instance Solid Disc where+ rayint = rayint_disc+ shadow = shadow_disc+ inside (Disc _ _ _) _ = False+ bound = bound_disc+++rayint_cylinder :: Cylinder -> Ray -> Flt -> Texture -> Rayint+rayint_cylinder (Cylinder r h1 h2) (Ray orig dir) d t =+ let Vec ox oy oz = orig+ Vec dx dy dz = dir+ a = dx*dx + dy*dy+ b = 2*(dx*ox + dy*oy)+ c = ox*ox + oy*oy - r*r+ disc = b*b - 4*a*c+ in if disc < 0 + then RayMiss+ else + let discsqrt = sqrt disc + q = if b < 0 + then (b-discsqrt)*(-0.5)+ else (b+discsqrt)*(-0.5)+ t0' = q/a+ t1' = c/q+ t0 = fmin t0' t1'+ t1 = fmax t0' t1'+ in if t1 < 0 || t0 > d + then RayMiss+ else let dist = if t0 < 0+ then t1+ else t0+ in if dist < 0 || dist > d+ then RayMiss+ else let pos = vscaleadd orig dir dist+ Vec posx posy posz = pos+ in if posz > h1 && posz < h2+ then RayHit dist pos (Vec (posx/r) (posy/r) 0) t+ else if dz > 0 -- ray pointing up from bottom+ then if oz < h1+ then rayint_disc (Disc (Vec 0 0 h1) nvz (r*r)) (Ray orig dir) d t+ --then rayint_aadisc h1 r (Ray orig dir) d t+ else RayMiss+ else if oz > h2+ then rayint_disc (Disc (Vec 0 0 h2) vz (r*r)) (Ray orig dir) d t+ --rayint_aadisc h2 r (Ray orig dir) d t -- todo: fix normal+ else RayMiss++inside_cylinder :: Cylinder -> Vec -> Bool+inside_cylinder (Cylinder r h1 h2) (Vec x y z) =+ z > h1 && z < h2 && x*x + y*y < r*r+ +bound_cylinder :: Cylinder -> Bbox+bound_cylinder (Cylinder r h1 h2) =+ Bbox (Vec (-r) (-r) h1) (Vec r r h2)++instance Solid Cylinder where+ rayint = rayint_cylinder+ inside = inside_cylinder+ bound = bound_cylinder+++rayint_cone :: Cone -> Ray -> Flt -> Texture -> Rayint+rayint_cone (Cone r clip1 clip2 height) (Ray orig dir) d t =+ let Vec ox oy oz = orig+ Vec dx dy dz = dir+ k' = (r/height)+ k = k'*k'+ a = dx*dx + dy*dy - k*dz*dz+ b = 2*(dx*ox + dy*oy - k*dz*(oz-height))+ c = ox*ox + oy*oy - k*(oz-height)*(oz-height)+ disc = b*b - 4*a*c+ in if disc < 0+ then RayMiss+ else+ let discsqrt = sqrt disc+ q = if b < 0+ then (b-discsqrt)*(-0.5)+ else (b+discsqrt)*(-0.5)+ t0' = q/a+ t1' = c/q+ t0 = fmin t0' t1'+ t1 = fmax t0' t1'+ in if t1 < 0 || t0 > d + then RayMiss+ else let dist = if t0 < 0+ then t1+ else t0+ in if dist < 0 || dist > d+ then RayMiss+ else+ let pos = vscaleadd orig dir dist+ Vec posx posy posz = pos+ in if posz > clip1 && posz < clip2+ then let invhyp = 1 / (sqrt (height*height + r*r))+ up = r * invhyp+ out = height * invhyp+ r_ = sqrt (posx*posx + posy*posy)+ correction = (out)/(r_)+ in RayHit dist pos (Vec (posx*correction) (posy*correction) up) t+ else + if dz > 0 -- ray pointing up from bottom+ then if oz < clip1+ then rayint_disc (Disc (Vec 0 0 clip1) nvz (r*r)) (Ray orig dir) d t+ else RayMiss+ else if oz > clip2+ then let r2 = r*(1-((clip2-clip1)/(height)))+ in rayint_disc (Disc (Vec 0 0 clip2) vz (r2*r2)) (Ray orig dir) d t+ --rayint_aadisc clip2 r2 (Ray orig dir) d t+ else RayMiss+ -- then rayint_aadisc clip1 r (Ray orig dir) d t+ -- else RayMiss -- rayint_aadisc clip2 + -- (r*((clip2-clip1)/height)) + -- (Ray orig dir) d t -- todo: fix normal++shadow_cone :: Cone -> Ray -> Flt -> Bool+shadow_cone (Cone r clip1 clip2 height) (Ray orig dir) d =+ let Vec ox oy oz = orig+ Vec dx dy dz = dir+ k' = (r/height)+ k = k'*k'+ a = dx*dx + dy*dy - k*dz*dz+ b = 2*(dx*ox + dy*oy - k*dz*(oz-height))+ c = ox*ox + oy*oy - k*(oz-height)*(oz-height)+ disc = b*b - 4*a*c+ in if disc < 0+ then False+ else+ let discsqrt = sqrt disc+ q = if b < 0+ then (b-discsqrt)*(-0.5)+ else (b+discsqrt)*(-0.5)+ t0' = q/a+ t1' = c/q+ t0 = fmin t0' t1'+ t1 = fmax t0' t1'+ in if t1 < 0 || t0 > d + then False+ else let dist = if t0 < 0+ then t1+ else t0+ in if dist < 0 || dist > d+ then False+ else+ let pos = vscaleadd orig dir dist+ Vec posx posy posz = pos+ in if posz > clip1 && posz < clip2+ then True+ else + if dz > 0 -- ray pointing up from bottom+ then if oz < clip1+ then shadow (Disc (Vec 0 0 clip1) nvz (r*r)) (Ray orig dir) d+ else False+ else if oz > clip2+ then let r2 = r*(1-((clip2-clip1)/(height)))+ in shadow (Disc (Vec 0 0 clip2) vz (r2*r2)) (Ray orig dir) d+ else False+++inside_cone :: Cone -> Vec -> Bool+inside_cone (Cone rbase h1 h2 height) (Vec x y z) =+ let r = rbase*(1-(((z-h1)/height)))+ in z > h1 && z < h2 && x*x + y*y < r*r++bound_cone :: Cone -> Bbox+bound_cone (Cone r h1 h2 height) =+ Bbox (Vec (-r) (-r) h1) (Vec r r h2)++instance Solid Cone where+ rayint = rayint_cone+ shadow = shadow_cone+ inside = inside_cone+ bound = bound_cone
+ Csg.hs view
@@ -0,0 +1,103 @@+module Csg (difference, intersection) where+import Vec+import Solid+import Data.List++-- Constructive Solid Geometry+-- (boolean operations for solids)++-- todo: implement shadow tests++data Difference = Difference SolidItem SolidItem deriving Show+data Intersection = Intersection [SolidItem] deriving Show++--Difference--+-- csg of object b subtracted from object a --+difference :: SolidItem -> SolidItem -> SolidItem+difference a b = SolidItem $ Difference a b++rayint_difference :: Difference -> Ray -> Flt -> Texture -> Rayint+rayint_difference dif r d t =+ let Difference sa sb = dif+ Ray orig dir = r+ ria = rayint sa r d t+ in+ case ria of+ RayMiss -> RayMiss+ RayHit ad ap an at ->+ if inside sb orig + then+ case rayint sb r d t of+ RayMiss -> RayMiss + RayHit bd bp bn bt ->+ if bd < ad + then if inside sa bp + then RayHit bd bp (vinvert bn) bt+ else rayint_advance (SolidItem dif) r d t bd+ else rayint_advance (SolidItem dif) r d t bd+ else + if inside sb ap+ then rayint_advance (SolidItem dif) r d t ad+ else RayHit ad ap an at+++--Intersection--+intersection :: [SolidItem] -> SolidItem+intersection slds = SolidItem $ Intersection slds++-- fixme: there's some numerical instability near edges+rayint_intersection :: Intersection -> Ray -> Flt -> Texture -> Rayint+rayint_intersection (Intersection slds) r d t =+ let (Ray orig dir) = r + in+ if null slds || d < 0+ then RayMiss+ else + let s = head slds + in case tail slds of+ [] -> rayint s r d t+ ss -> if inside s orig+ then case rayint s r d t of + RayMiss -> rayint (Intersection ss) r d t+ RayHit sd sp sn st -> + case rayint (Intersection ss) r sd t of+ RayMiss -> rayint_advance (SolidItem (Intersection slds)) + r d t sd + hit -> hit+ else case rayint s r d t of+ RayMiss -> RayMiss+ RayHit sd sp sn st ->+ if inside (Intersection ss) sp+ then RayHit sd sp sn st+ else rayint_advance (SolidItem (Intersection slds))+ r d t sd++inside_difference :: Difference -> Vec -> Bool+inside_difference (Difference sa sb) pt =+ (inside sa pt) && (not $ inside sb pt)++-- note: inside is True for an empty intersection.+-- this is actually the preferred semantics in +-- some cases, strange as it may seem.+inside_intersection :: Intersection -> Vec -> Bool+inside_intersection (Intersection slds) pt =+ foldl' (&&) True (map (\x -> inside x pt) slds) ++bound_difference :: Difference -> Bbox+bound_difference (Difference sa sb) = bound sa++bound_intersection :: Intersection -> Bbox+bound_intersection (Intersection slds) =+ if null slds + then empty_bbox+ else foldl' bboverlap everything_bbox (map bound slds)++instance Solid Difference where+ rayint = rayint_difference+ inside = inside_difference+ bound = bound_difference++instance Solid Intersection where+ rayint = rayint_intersection+ inside = inside_intersection+ bound = bound_intersection
Glome.hs view
@@ -1,6 +1,4 @@-import Vec-import Clr-import Solid+import Scene import Trace import Spd import TestScene@@ -14,122 +12,108 @@ import System import System.Console.GetOpt import Data.Maybe( fromMaybe )+-- import OpenEXR -- work in progress -- import Debug.Trace -- import Data.ByteString -get_color :: Flt -> Flt -> Scene -> Clr.Color---+maxdepth = 2 -- recursion depth for reflection/refraction++-- compute ray, invoke trace function, return color+get_color :: Flt -> Flt -> Scene -> (Scene.Color,Flt) get_color x y scn = let (Scene sld lights (Camera pos fwd up right) dtex bgcolor) = scn dir = vnorm $ vadd3 fwd (vscale right (-x)) (vscale up y) ray = (Ray pos dir) in- Trace.trace scn ray infinity 2---}+ Trace.trace_depth scn ray infinity maxdepth -{----- for testing screen-draw overhead-get_color x y scn =- Clr.Color (1*x) (1*y) 0 ---}+-- compute a packet of four rays from corners of box+get_packet :: Flt -> Flt -> Flt -> Flt -> Scene -> PacketColor+get_packet x1 y1 x2 y2 scn =+ let (Scene sld lights (Camera pos fwd up right) dtex bgcolor) = scn+ dir1 = vnorm $ vadd3 fwd (vscale right (-x1)) (vscale up y1)+ dir2 = vnorm $ vadd3 fwd (vscale right (-x2)) (vscale up y1)+ dir3 = vnorm $ vadd3 fwd (vscale right (-x1)) (vscale up y2)+ dir4 = vnorm $ vadd3 fwd (vscale right (-x2)) (vscale up y2)+ ray1 = Ray pos dir1+ ray2 = Ray pos dir2+ ray3 = Ray pos dir3+ ray4 = Ray pos dir4+ in trace_packet scn ray1 ray2 ray3 ray4 infinity maxdepth +-- convert trace result to +-- appropriate float type for OpenGL fc :: Flt -> Float---fc x = if x == x then clamp 0 x 0.5 else error "nan" fc x = realToFrac x---fc x = x+delta --- double nested for loop--- for rendering a sub-box of the screen-gen_pixels curx cury stopx stopy maxx maxy scene =+-- given a block of screen coordinates, return list of pixels+gen_pixel_list :: Flt -> Flt -> Flt -> Flt -> Flt -> Flt -> Scene -> [(Flt,Flt,Flt,Flt,Flt,Flt)]+gen_pixel_list curx cury stopx stopy maxx maxy scene = let midx = maxx/2 midy = maxy/2 gp x y =- if y>=stopy then- return ()+ if y >= stopy then+ [] else- if x>=stopx + if x >= stopx then gp curx (y+1) else - do- let scx = (x-midx) / midx- let scy = (y-midy) / midy- -- let (Clr.Color r g b) = get_color scx (scy*(midy/midx)) scene- let (Clr.Color r g b) = get_color (scx*(midx/midy)) scy scene- currentColor $= Color4 (fc r) (fc g) (fc b) 1- vertex$Vertex3 scx scy 0- gp (x+1) y+ let scx = (x-midx) / midx+ scy = (y-midy) / midy+ --(Clr.Color r g b) = get_color scx (scy*(midy/midx)) scene+ --(Clr.Color r g b) = get_color (scx*(midx/midy)) scy scene+ ((Scene.Color r g b),d) = get_color (scx*(midx/midy)) scy scene+ in+ (scx,scy,r,g,b,0) : (gp (x+1) y) in gp curx cury --- same as above, but non-monadic, instead return list of pixels-gen_pixel_list :: Flt -> Flt -> Flt -> Flt -> Flt -> Flt -> Scene -> [(Flt,Flt,Flt,Flt,Flt)]-gen_pixel_list curx cury stopx stopy maxx maxy scene =+-- same, but trace packets instead of mono-rays+gen_pixel_list_packet :: Flt -> Flt -> Flt -> Flt -> Flt -> Flt -> Scene -> [(Flt,Flt,Flt,Flt,Flt,Flt)]+gen_pixel_list_packet curx cury stopx stopy maxx maxy scene = let midx = maxx/2 midy = maxy/2 gp x y =- if y>=stopy then+ if y >= stopy then [] else- if x>=stopx + if x >= stopx then- gp curx (y+1)+ gp curx (y+2) else - let scx = (x-midx) / midx- scy = (y-midy) / midy- (Clr.Color r g b) = get_color scx (scy*(midy/midx)) scene- --(Clr.Color r g b) = get_color (scx*(midx/midy)) scy scene+ let scx1 = (x-midx) / midx+ scy1 = (y-midy) / midy+ scx2 = ((x+1)-midx) / midx+ scy2 = ((y+1)-midy) / midy++ PacketColor (Scene.Color r1 g1 b1)+ (Scene.Color r2 g2 b2)+ (Scene.Color r3 g3 b3)+ (Scene.Color r4 g4 b4) = get_packet (scx1*(midx/midy)) scy1+ (scx2*(midx/midy)) scy2 scene in- (scx,scy,r,g,b) : (gp (x+1) y)+ [(scx1,scy1,r1,g1,b1,0),+ (scx2,scy1,r2,g2,b2,0),+ (scx1,scy2,r3,g3,b3,0),+ (scx2,scy2,r4,g4,b4,0)] ++ (gp (x+2) y) in gp curx cury --- split screen into little blocks for parallel rendering-{--gen_blocks maxx maxy block_size scene =- let foo = 1- bar = 2- gb x y = - if y>=maxy then- return ()- else- if x>=maxx then- gb 0 (y+block_size)- else- do- gen_pixels x y (x+block_size-1) (y+block_size-1) maxx maxy scene- gb (x+block_size) y - in gb 0 0--}--{---- this doesn't seem to parallelize-gen_blocks maxx maxy block_size scene =- let xblocks = maxx/block_size- yblocks = maxy/block_size- blocks = Prelude.map (\x -> Prelude.map (\y -> (x*block_size,y*block_size) ) [0..yblocks-1] ) [0..xblocks-1]- in- do - -- mapM_ (\(x,y) -> gen_pixels x y (x+block_size) (y+block_size) maxx maxy scene) (concat blocks)- -- sequence_ $ Prelude.map (\(x,y) -> gen_pixels x y (x+block_size) (y+block_size) maxx maxy scene) (concat blocks)- sequence_ $ (parMap rwhnf) - (\(x,y) -> gen_pixels x y (x+block_size) (y+block_size) maxx maxy scene) - (concat blocks)--}---- as above, but parallel code is pure--- odd, this seems to be faster, even without multicore--- rnf is faster than rwhnf gen_blocks_list maxx maxy block_size scene = let xblocks = maxx/block_size yblocks = maxy/block_size- blocks = Prelude.concat $ Prelude.map (\x -> Prelude.map (\y -> (x*block_size,y*block_size) ) [0..yblocks-1] ) [0..xblocks-1]+ blocks = Prelude.concat $ Prelude.map + (\x -> Prelude.map + (\y -> (x*block_size,y*block_size) ) + [0..yblocks-1] ) + [0..xblocks-1] pixels = map -- (parMap rnf) - (\(x,y) -> gen_pixel_list x y (x+block_size) (y+block_size) maxx maxy scene)+ (\(x,y) -> gen_pixel_list_packet x y (x+block_size) (y+block_size) maxx maxy scene) (blocks) in do- mapM_ (\pix -> mapM_ (\(x,y,r,g,b) -> do currentColor $= Color4 (fc r) (fc g) (fc b) 1- vertex$Vertex3 (fc x) (fc y) 0 + mapM_ (\pix -> mapM_ (\(x,y,r,g,b,d) -> do currentColor $= Color4 (fc r) (fc g) (fc b) 1+ vertex$Vertex3 (fc x) (fc y) (fc d) ) pix) pixels @@ -149,44 +133,79 @@ args <- getArgs let (flags, nonOpts, msgs) = getOpt RequireOrder options args print $ "recognized options: " ++ (show (length flags))- scene <- getscene flags+ t1 <- getPOSIXTime+ scene <- getscene flags + -- print $ "(primitives,transforms,bounding objects): " ++ (show (primcount_scene scene))+ t2 <- getPOSIXTime+ print $ "scene setup: " ++ (show (t2-t1)) let sx = 720 :: GLsizei let sy = 480 :: GLsizei+ let sizex = fromIntegral sx+ let sizey = fromIntegral sy (name, _) <- getArgsAndInitialize- initialDisplayMode $= []--- initialDisplayMode $= [DoubleBuffered]+ --initialDisplayMode $= []+ initialDisplayMode $= [DoubleBuffered]+ pointSmooth $= Enabled++ -- create window createWindow name windowSize $= Size sx sy- displayCallback $= display scene sx sy++ -- set up camera+ -- why is the z-value (-100 < z < 0)? + -- I don't know, it just works this way for some reason++ {-+ matrixMode $= Projection+ loadIdentity+ ortho (-1) 1 (-1) 1 (-10000) 0+ matrixMode $= Modelview 0+ -}++ -- create display list+ t1 <- getPOSIXTime+ dlist <- defineNewList CompileAndExecute $ do + renderPrimitive Points $ gen_blocks_list sizex sizey 32 scene+ t2 <- getPOSIXTime+ print $ "render: " ++ (show (t2-t1))+ displayCallback $= display dlist keyboardMouseCallback $= Just (keyboard scene) mainLoop -display scene sx sy = do+display dlist = do+ clear [ColorBuffer]+ callList dlist+ swapBuffers++-- dodo: make this do some kind of antialiasing+display_aa scene sx sy = do t1 <- getPOSIXTime- clearColor $= Color4 0 0 0 1+ -- clearColor $= Color4 0 0 0 1 clear [ColorBuffer]- --(Size sx sy) <- GLUT.get windowSize let sizex = fromIntegral sx let sizey = fromIntegral sy -- renderPrimitive Points $ gen_blocks_list 512 512 128 scene- -- renderPrimitive Points $ gen_blocks_list 720 480 80 scene- renderPrimitive Points $ gen_pixels 0 0 sizex sizey sizex sizey scene- swapBuffers+ -- renderPrimitive Points $ gen_blocks_list sizex sizey 80 scene+ -- renderPrimitive Points $ gen_pixels 0 0 sizex sizey sizex sizey scene+ -- swapBuffers+ -- GLUT.rotate (fc (deg 1)) $Vector3 0 (1::GLfloat) 0 t2 <- getPOSIXTime- -- print ((fromInteger (t2-t1))/1000000000) print (t2-t1) + keyboard _ (Char 'q') Down _ _ = do exitWith ExitSuccess +-- for debugging, print a full scene dump keyboard s (Char 's') Down _ _ = do print (show s) keyboard _ _ _ _ _ = return () -+-- if a scene has been specified on the command line, render that;+-- otherwise, render whatever we find in TestScene.hs getscene :: [Flag] -> IO Scene getscene flags = case flags of @@ -196,10 +215,11 @@ (scene,s) <- return $ Prelude.head $ reads filestring return scene --data Flag = Filename String | Res Int Int deriving Show+-- todo: add argument for screen resolution+data Flag = Filename String | Res Int Int | Time Flt deriving Show options :: [OptDescr Flag] options = [ Option ['n'] ["filename"] (ReqArg Filename "FILE") "input NFF scene"+ --, Option ['t'] ["time"] (ReqArg Time 0) "time value for scene generation" ]
+ Plane.hs view
@@ -0,0 +1,45 @@+module Plane (plane, plane_offset) where+import Vec+import Solid++-- A plane is effectively a half-space; everything below the plane is+-- "inside", everything above is "outside".++data Plane = Plane Vec Flt deriving Show -- normal, perpendicular offset from origin++-- Usually, the most convenient way to define a plane is +-- by specifying a point on the plane and a normal++plane :: Vec -> Vec -> SolidItem+plane orig norm_ = SolidItem $ Plane norm d+ where norm = vnorm norm_+ d = vdot orig norm++-- we can also specify a point and a perpindicular offset:++plane_offset :: Vec -> Flt -> SolidItem+plane_offset pt off = SolidItem $ Plane pt off++rayint_plane :: Plane -> Ray -> Flt -> Texture -> Rayint+rayint_plane (Plane norm offset) (Ray orig dir) d t =+ let hit = -(((vdot norm orig)-offset) / (vdot norm dir))+ in if hit < 0 || hit > d + then RayMiss+ else RayHit hit (vscaleadd orig dir hit) norm t++inside_plane :: Plane -> Vec -> Bool+inside_plane (Plane norm offset) pt =+ let onplane = (vscale norm offset)+ newvec = vsub onplane pt+ in vdot newvec norm > 0++-- Note: attempting to use an infinite object (such as+-- a plane) inside a bih will cause an exception.++bound_plane :: Plane -> Bbox+bound_plane (Plane norm offset) = everything_bbox++instance Solid Plane where+ rayint = rayint_plane+ inside = inside_plane+ bound = bound_plane
README view
@@ -26,7 +26,11 @@ - perlin noise and a few basic textures are implemented - uses a bounding interval heirarchy acceleration structure or you can construct a BVH manually with the "Bound" primitive-- multiprocessor support is currently disabled+- multiprocessor support is available, but it leaks memory+ and scaling is sub-linear+- packet tracing (2x2) of primary rays is enabled by default+- as of 0.5, glome uses type classes, and new primitives can+ be defined in their own module Using: to load an NFF scene, run "./Glome -n [filename]". Otherwise, a default scene is rendered, defined in "TestScene.hs".
+ Scene.hs view
@@ -0,0 +1,73 @@+module Scene (Scene(Scene), Light(Light), Camera(Camera),+ scene, camera, light, + sld, lits, cam, dtex, bground,+ module Clr,+ module Vec,+ module Solid,+ module Sphere,+ module Triangle,+ module Bih,+ module Csg,+ module Plane,+ module Box,+ module Bound,+ module Cone,+ module Tex) where+import Clr+import Vec+import Solid+import Sphere+import Triangle+import Bih+import Csg+import Plane+import Box+import Bound+import Cone+import Tex++-- This is the module to import if you want to have+-- access to all the Solid constructors and scene+-- defininition code.++--LIGHTS--+data Light = Light {litpos :: !Vec,+ litcol :: !Color} deriving Show++light :: Vec -> Color -> Light+light pos clr = Light pos clr++-- CAMERA --+data Camera = Camera {campos, fwd, up, right :: !Vec} + deriving Show++default_cam = (Camera (vec 0.0 0.0 (-3.0)) + (vec 0.0 0.0 1.0) + (vec 0.0 1.0 0.0) + (vec 1.0 0.0 0.0) )++camera :: Vec -> Vec -> Vec -> Flt -> Camera+camera pos at up angle =+ let fwd = vnorm $ vsub at pos+ right = vnorm $ vcross up fwd+ up_ = vnorm $ vcross fwd right+ cam_scale = tan ((pi/180)*(angle/2))+ in+ Camera pos fwd+ (vscale up_ cam_scale) + (vscale right cam_scale)++--SCENE--+data Scene = Scene {sld :: SolidItem,+ lits :: [Light], + cam :: Camera, + dtex :: Texture, + bground :: Color} deriving Show++scene :: SolidItem -> [Light] -> Camera -> Texture -> Color -> Scene+scene s l cam t clr = Scene s l cam t clr++{-+default_scene = (Scene (sphere (vec 0.0 0.0 0.0) 1.0) + [] default_cam t_white c_white)+-}
Solid.hs view
@@ -3,1181 +3,383 @@ import Clr import Data.List hiding (group) ---COMMON DATATYPES AND FUNCTIONS---data Bbox = Bbox {p1 :: !Vec, p2 :: !Vec} deriving Show-data Interval = Interval !Flt !Flt deriving Show -- used instead of a tuple----union of two bounding boxes-bbjoin :: Bbox -> Bbox -> Bbox-bbjoin (Bbox p1a p2a) (Bbox p1b p2b) =- (Bbox (vmin p1a p1b) (vmax p2a p2b))----overlap of two bounding boxes-bboverlap :: Bbox -> Bbox -> Bbox-bboverlap (Bbox p1a p2a) (Bbox p1b p2b) =- (Bbox (vmax p1a p1b) (vmin p2a p2b))----split a bounding box into two-bbsplit :: Bbox -> Int -> Flt -> (Bbox,Bbox)-bbsplit (Bbox p1 p2) axis offset =- if (offset < (va p1 axis)) || (offset > (va p2 axis))- then error "degenerate bounding box split"- else ((Bbox p1 (vset p2 axis offset)),- (Bbox (vset p1 axis offset) p2))---- generate a bounding box from a list of points-bbpts :: [Vec] -> Bbox-bbpts [] = empty_bbox-bbpts ((Vec x y z):[]) =- Bbox (Vec (x-delta) (y-delta) (z-delta)) - (Vec (x+delta) (y+delta) (z+delta))--bbpts ((Vec x y z):pts) =- let (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bbpts pts- minx = fmin (x-delta) p1x- miny = fmin (y-delta) p1y- minz = fmin (z-delta) p1z- maxx = fmax (x+delta) p2x- maxy = fmax (y+delta) p2y- maxz = fmax (z+delta) p2z in- Bbox (Vec minx miny minz) (Vec maxx maxy maxz)---- surface area, volume-bbsa :: Bbox -> Flt-bbsa (Bbox p1 p2) =- let Vec dx dy dz = vsub p2 p1 - in dx*dy + dx*dz + dy*dz--bbvol :: Bbox -> Flt-bbvol (Bbox p1 p2) =- let (Vec dx dy dz) = vsub p2 p1- in dx*dy*dz--empty_bbox = - Bbox (Vec infinity infinity infinity) - (Vec (-infinity) (-infinity) (-infinity))--everything_bbox =- Bbox (Vec (-infinity) (-infinity) (-infinity))- (Vec infinity infinity infinity)---bbclip :: Ray -> Bbox -> Interval-bbclip (Ray (Vec ox oy oz) (Vec dx dy dz)) - (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) =- let dxrcp = 1/dx- dyrcp = 1/dy- dzrcp = 1/dz- Interval inx outx = if dx > 0 - then Interval ((p1x-ox)*dxrcp) ((p2x-ox)*dxrcp)- else Interval ((p2x-ox)*dxrcp) ((p1x-ox)*dxrcp)- Interval iny outy = if dy > 0- then Interval ((p1y-oy)*dyrcp) ((p2y-oy)*dyrcp)- else Interval ((p2y-oy)*dyrcp) ((p1y-oy)*dyrcp)- Interval inz outz = if dz > 0- then Interval ((p1z-oz)*dzrcp) ((p2z-oz)*dzrcp)- else Interval ((p2z-oz)*dzrcp) ((p1z-oz)*dzrcp)- in- Interval (fmax3 inx iny inz) (fmin3 outx outy outz)--data Rayint = RayHit {- depth :: !Flt,- pos :: !Vec,- norm :: !Vec,- texture :: !Texture-} | RayMiss deriving Show--nearest :: Rayint -> Rayint -> Rayint-nearest a RayMiss = a-nearest RayMiss b = b-nearest (RayHit da pa na ta) (RayHit db pb nb tb) =- if da < db- then RayHit da pa na ta- else RayHit db pb nb tb--furthest :: Rayint -> Rayint -> Rayint-furthest a RayMiss = RayMiss-furthest RayMiss b = RayMiss-furthest (RayHit da pa na ta) (RayHit db pb nb tb) =- if da > db- then RayHit da pa na ta- else RayHit db pb nb tb--hit :: Rayint -> Bool-hit (RayHit _ _ _ _) = True-hit RayMiss = False--dist :: Rayint -> Flt-dist RayMiss = infinity-dist (RayHit d _ _ _) = d-----LIGHTS---data Light = Light {litpos :: !Vec,- litcol :: !Color} deriving Show----MATERIALS---data Material = Material {clr :: Color, - reflect, refract, ior, - kd, shine :: !Flt} deriving Show-type Texture = Rayint -> Material---- this is sort of a no-op; we don't have a--- good way to show an arbitrary function-showTexture :: Texture -> String-showTexture t = show $ t RayMiss--instance Show Texture where- show = showTexture--m_white = (Material c_white 0 0 0 1 2)-t_white ri = m_white--t_uniform :: Material -> Texture-t_uniform m = \x -> m--interp :: Flt -> Flt -> Flt -> Flt-interp scale a b =- scale*a + (1-scale)*b----not really correct, but we'll go with it for now-m_interp :: Material -> Material -> Flt -> Material-m_interp m1 m2 scale =- let (Material m1c m1refl m1refr m1ior m1kd m1shine) = m1- (Material m2c m2refl m2refr m2ior m2kd m2shine) = m2- intp = interp scale- c = cadd (cscale m1c scale) (cscale m2c (1-scale))- refl = intp m1refl m2refl- refr = intp m1refr m2refr- ior = intp m1ior m2ior- kd = intp m1kd m2kd- shine = intp m1shine m2shine- in (Material c refl refr ior kd shine)----SOLID TYPES----data Solid = Sphere {center :: !Vec, - radius, invradius :: !Flt}- | Triangle {v1, v2, v3 :: Vec}- | TriangleNorm {v1, v2, v3, n1, n2, n3 :: Vec}- | Disc !Vec !Vec !Flt -- position, normal, r*r- | Cylinder !Flt !Flt !Flt -- radius height1 height2- | Cone !Flt !Flt !Flt !Flt -- r clip1 clip2 height- | Plane Vec Flt -- normal, offset from origin- | Box !Bbox- | Group ![Solid]- | Intersection ![Solid]- | Bevel !Solid !Flt- | Bound Solid Solid- | Difference !Solid !Solid- | Bih {bihbb :: !Bbox, bihroot :: !BihNode}- | Instance !Solid !Xfm- | Tex !Solid Texture- | Portal Solid Solid -- if we hit a, intersect with b- | Photon Vec Vec Color Flt -- pos, incident ray, color, radius- | Nothing deriving Show -- conflicts with - -- Nothing :: Maybe a from prelude--data BihNode = BihLeaf !Solid - | BihBranch {lmax :: !Flt, rmin :: !Flt, ax :: !Int, - l :: BihNode, r :: BihNode} deriving Show----CONSTRUCTORS---sphere :: Vec -> Flt -> Solid-sphere c r =- Sphere c r (1.0/r)--triangle :: Vec -> Vec -> Vec -> Solid-triangle v1 v2 v3 =- Triangle v1 v2 v3----simple tesselation-triangles :: [Vec] -> [Solid]-triangles (v1:vs) =- zipWith (\v2 v3 -> triangle v1 v2 v3) vs (tail vs) --trianglenorm v1 v2 v3 n1 n2 n3 =- -- Triangle v1 v2 v3- TriangleNorm v1 v2 v3 n1 n2 n3--trianglesnorms :: [(Vec,Vec)] -> [Solid]-trianglesnorms (vn1:vns) =- zipWith (\vn2 vn3 -> trianglenorm (fst vn1) (fst vn2) (fst vn3) - (snd vn1) (snd vn2) (snd vn3))- vns (tail vns)--box :: Vec -> Vec -> Solid-box p1 p2 =- Box (Bbox p1 p2)--disc :: Vec -> Vec -> Flt -> Solid-disc pos norm r =- Disc pos norm (r*r)--cylinder_z :: Flt -> Flt -> Flt -> Solid-cylinder_z r h1 h2 = Cylinder r h1 h2--cone_z :: Flt -> Flt -> Flt -> Flt -> Solid-cone_z r h1 h2 height = Cone r h1 h2 height---- construct a general cylinder from p1 to p2 with radius r-cylinder :: Vec -> Vec -> Flt -> Solid-cylinder p1 p2 r =- let axis = vsub p2 p1- len = vlen axis- ax1 = vscale axis (1/len)- (ax2,ax3) = orth ax1 - in Instance (cylinder_z r 0 len)- (compose [ (xyz_to_uvw ax2 ax3 ax1),- (translate p1) ])- --- similar for cone-cone :: Vec -> Flt -> Vec -> Flt -> Solid-cone p1 r1 p2 r2 =- if r1 < r2 - then cone p2 r2 p1 r1- else if r1-r2 < delta- then cylinder p1 p2 r2- else- let axis = vsub p2 p1- len = vlen axis- ax1 = vscale axis (1/len)- (ax2,ax3) = orth ax1 - height = (r1*len)/(r1-r2) -- distance to end point- in- Instance (cone_z r1 0 len height)- (compose [ (xyz_to_uvw ax2 ax3 ax1),- (translate p1) ]) --plane :: Vec -> Vec -> Solid-plane orig norm_ = Plane norm d- where norm = vnorm norm_- d = vdot orig norm---- flatten tree of groups into a single group-flatten_group :: [Solid] -> [Solid]-flatten_group ((Group slds):xs) =- (flatten_group slds) ++ xs-flatten_group ((Solid.Nothing):xs) = xs-flatten_group x = x--group :: [Solid] -> Solid-group [] = Solid.Nothing-group (sld:[]) = sld-group slds =- Group (flatten_group slds)--transform :: Solid -> [Xfm] -> Solid-transform (Instance s xfm2) xfm1 = - transform s [compose ([xfm2] ++ xfm1)]--transform s xfm =- Instance s (compose xfm)---- push all the transforms out to the leaves--- and throw away pre-existing bounding volumes--- so we can run the bih constructor on the--- resulting group-flatten_transform :: Solid -> [Solid]-flatten_transform (Group slds) =- flatten_group $ concat (map flatten_transform slds)--flatten_transform (Instance s xfm) =- case s of - Group slds -> flatten_transform $ group (map (\x -> transform x [xfm]) slds)- Bound sa sb -> flatten_transform (transform sb [xfm])- Instance sa xfm2 -> flatten_transform (transform s [xfm])- _ -> [transform s [xfm]]--flatten_transform (Bound sa sb) = flatten_transform sb---- bih construction-build_leaf objs =- BihLeaf (group (map snd objs))--max_bih_sa = 0.3 :: Flt--build_rec :: [(Bbox,Solid)] -> Bbox -> Bbox -> Int -> BihNode-build_rec objs nodebox splitbox depth = - -- if (null objs) || (null $ tail objs) || - -- (null $ tail $ tail objs)- if length objs < 2- then build_leaf objs- else- let (Bbox nodeboxp1 nodeboxp2) = nodebox- (Bbox splitboxp1 splitboxp2) = splitbox- axis = vmaxaxis (vsub splitboxp2 splitboxp1)- bbmin = va splitboxp1 axis- bbmax = va splitboxp2 axis- candidate = (bbmin + bbmax) * 0.5- in- if candidate > (va nodeboxp2 axis) then- build_rec objs nodebox - (Bbox splitboxp1 (vset splitboxp2 axis candidate)) - depth- else- if candidate < (va nodeboxp1 axis) then- build_rec objs nodebox (- Bbox (vset splitboxp1 axis candidate) splitboxp2) - depth- else- -- not sure if this is a big win- let nbsa = bbsa nodebox- (big,small) = partition (\ (bb,_) -> - (bbsa bb) > (nbsa * max_bih_sa)) objs- in - if (not $ null big) && ((length big) < ((length small)*2))- then (BihBranch (va nodeboxp2 0) (va nodeboxp1 0) 0- (build_rec big nodebox splitbox (depth+1))- (build_rec small nodebox splitbox (depth+1)) )- else- let (l,r) = partition (\((Bbox bbp1 bbp2),_)-> - (((va bbp1 axis)+(va bbp2 axis))*0.5) - < candidate ) objs- lmax = foldl fmax (-infinity) (map (\((Bbox _ p2),_) -> va p2 axis) l)- rmin = foldl fmin infinity (map (\((Bbox p1 _),_) -> va p1 axis) r)- (lsplit,rsplit) = bbsplit splitbox axis candidate- lnb = (Bbox nodeboxp1 (vset nodeboxp2 axis lmax))- rnb = (Bbox (vset nodeboxp1 axis rmin) nodeboxp2)- in- -- stop if there's no progress being made- if ((null l) && (rmin <= bbmin)) ||- ((null r) && (lmax >= bbmax))- then build_leaf objs- else- (BihBranch (lmax+delta) (rmin-delta) axis- (build_rec l lnb lsplit (depth+1))- (build_rec r rnb rsplit (depth+1)) )--bih :: [Solid] -> Solid-bih [] = Solid.Nothing--- bih (sld:[]) = sld -- sometimes we'd like to be able to use a- -- single object bih just for its aabb-bih slds =- let objs = map (\x -> ((bound x),x)) (flatten_group slds)- bb = foldl bbjoin empty_bbox (map (\(b,_)->b) objs)- root = build_rec objs bb bb 0- (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bb- in- if p1x == (-infinity) || p1y == (-infinity) || p1z == (-infinity) ||- p2x == infinity || p2y == infinity || p2z == infinity- then- error $ "bih: infinite bounding box " ++ (show objs)- else- (Bih bb root)----INTERSECTION TESTS---rayint :: Solid -> Ray -> Flt -> Texture -> Rayint----Basic Primitives-----Triangle----- adaptation of Moller and Trumbore from pbrt page 127-rayint (Triangle p1 p2 p3) (Ray o dir) dist tex =- let e1 = vsub p2 p1- e2 = vsub p3 p1- s1 = vcross dir e2- divisor = vdot s1 e1- in - if (divisor == 0)- then RayMiss- else- let invdivisor = 1.0 / divisor- d = vsub o p1 - b1 = (vdot d s1) * invdivisor- in- if (b1 < 0) || (b1 > 1) - then RayMiss - else- let s2 = vcross d e1- b2 = (vdot dir s2) * invdivisor- in- if (b2 < 0) || (b1 + b2 > 1) - then RayMiss- else- let t = (vdot e2 s2) * invdivisor- in- if (t < 0) || (t > dist)- then RayMiss- else- RayHit t (vscaleadd o dir t) (vnorm (vcross e1 e2)) tex--rayint (TriangleNorm p1 p2 p3 n1 n2 n3) (Ray o dir) dist tex =- let e1 = vsub p2 p1- e2 = vsub p3 p1- s1 = vcross dir e2- divisor = vdot s1 e1- in - if (divisor == 0)- then RayMiss- else- let invdivisor = 1.0 / divisor- d = vsub o p1 - b1 = (vdot d s1) * invdivisor- in- if (b1 < 0) || (b1 > 1) - then RayMiss - else- let s2 = vcross d e1- b2 = (vdot dir s2) * invdivisor- in- if (b2 < 0) || (b1 + b2 > 1) - then RayMiss- else- let t = (vdot e2 s2) * invdivisor- in- if (t < 0) || (t > dist)- then RayMiss- else- let n1scaled = (vscale n1 (1-(b1+b2))) - n2scaled = (vscale n2 b1)- n3scaled = (vscale n3 b2)- norm = vnorm $ vadd3 n1scaled n2scaled n3scaled- in RayHit t (vscaleadd o dir t) norm tex----Sphere----- adapted from graphics gems volume 1-rayint (Sphere center r invr) (Ray e dir) dist t = - let eo = vsub center e- v = vdot eo dir- in- if (dist >= (v - r)) && (v > 0.0)- then- let vsqr = v*v- csqr = vdot eo eo- rsqr = r*r- disc = rsqr - (csqr - vsqr) in- if disc < 0.0 then- RayMiss- else- let d = sqrt disc- hitdist = if (v-d) > 0 then (v-d) else (v+d)- in if (hitdist < 0) || (hitdist > dist)- then RayMiss- else- let p = vscaleadd e dir hitdist- -- n = vscale (vsub p center) invr in- -- n = vsub (vscale p invr) (vscale center invr) in- n = vnorm (vsub p center) - in RayHit hitdist p n t- else- RayMiss---- nice and simple-rayint (Disc point norm radius_sqr) r d t =- let (Ray orig dir) = r- dist = plane_int_dist r point norm - in if dist < 0 || dist > d - then RayMiss- else let pos = vscaleadd orig dir dist- offset = vsub pos point- in - if (vdot offset offset) > radius_sqr- then RayMiss- else RayHit dist pos norm t---- cylinder aligned to z axis--- no end caps--- adapted from pbrt-rayint (Cylinder r h1 h2) (Ray orig dir) d t =- let Vec ox oy oz = orig- Vec dx dy dz = dir- a = dx*dx + dy*dy- b = 2*(dx*ox + dy*oy)- c = ox*ox + oy*oy - r*r- disc = b*b - 4*a*c- in if disc < 0 - then RayMiss- else - let discsqrt = sqrt disc - q = if b < 0 - then (b-discsqrt)*(-0.5)- else (b+discsqrt)*(-0.5)- t0' = q/a- t1' = c/q- t0 = fmin t0' t1'- t1 = fmax t0' t1'- in if t1 < 0 || t0 > d - then RayMiss- else let dist = if t0 < 0- then t1- else t0- in if dist < 0 || dist > d- then RayMiss- else let pos = vscaleadd orig dir dist- Vec posx posy posz = pos- in if posz > h1 && posz < h2- then RayHit dist pos (Vec (posx/r) (posy/r) 0) t- else if dz > 0 -- ray pointing up from bottom- then if oz < h1- then rayint (Disc (Vec 0 0 h1) nvz (r*r)) (Ray orig dir) d t- --then rayint_aadisc h1 r (Ray orig dir) d t- else RayMiss- else if oz > h2- then rayint (Disc (Vec 0 0 h2) vz (r*r)) (Ray orig dir) d t- --rayint_aadisc h2 r (Ray orig dir) d t -- todo: fix normal- else RayMiss---- cone centered on z axiz, height of hp, clipped at h1 and h2-rayint (Cone r clip1 clip2 height) (Ray orig dir) d t =- let Vec ox oy oz = orig- Vec dx dy dz = dir- k' = (r/height)- k = k'*k'- a = dx*dx + dy*dy - k*dz*dz- b = 2*(dx*ox + dy*oy - k*dz*(oz-height))- c = ox*ox + oy*oy - k*(oz-height)*(oz-height)- disc = b*b - 4*a*c- in if disc < 0- then RayMiss- else- let discsqrt = sqrt disc- q = if b < 0- then (b-discsqrt)*(-0.5)- else (b+discsqrt)*(-0.5)- t0' = q/a- t1' = c/q- t0 = fmin t0' t1'- t1 = fmax t0' t1'- in if t1 < 0 || t0 > d - then RayMiss- else let dist = if t0 < 0- then t1- else t0- in if dist < 0 || dist > d- then RayMiss- else- let pos = vscaleadd orig dir dist- Vec posx posy posz = pos- in if posz > clip1 && posz < clip2- then let invhyp = 1 / (sqrt (height*height + r*r))- up = r * invhyp- out = height * invhyp- r_ = sqrt (posx*posx + posy*posy)- correction = (out)/(r_)- in RayHit dist pos (Vec (posx*correction) (posy*correction) up) t- else - if dz > 0 -- ray pointing up from bottom- then if oz < clip1- then rayint (Disc (Vec 0 0 clip1) nvz (r*r)) (Ray orig dir) d t- else RayMiss- else if oz > clip2- then let r2 = r*(1-((clip2-clip1)/(height)))- in rayint (Disc (Vec 0 0 clip2) vz (r2*r2)) (Ray orig dir) d t- --rayint_aadisc clip2 r2 (Ray orig dir) d t- else RayMiss- -- then rayint_aadisc clip1 r (Ray orig dir) d t- -- else RayMiss -- rayint_aadisc clip2 - -- (r*((clip2-clip1)/height)) - -- (Ray orig dir) d t -- todo: fix normal----Plane---rayint (Plane norm offset) (Ray orig dir) d t =- let hit = -(((vdot norm orig)-offset) / (vdot norm dir))- in if hit < 0 || hit > d - then RayMiss- else RayHit hit (vscaleadd orig dir hit) norm t----Box----- this could be optimized a bit more-rayint (Box (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z))) r d t =- let (Ray orig dir) = r- (Vec ox oy oz) = orig- (Vec dx dy dz) = dir- dxrcp = 1/dx- dyrcp = 1/dy- dzrcp = 1/dz- Interval inx outx = if dx > 0 - then Interval ((p1x-ox)*dxrcp) ((p2x-ox)*dxrcp)- else Interval ((p2x-ox)*dxrcp) ((p1x-ox)*dxrcp)- Interval iny outy = if dy > 0- then Interval ((p1y-oy)*dyrcp) ((p2y-oy)*dyrcp)- else Interval ((p2y-oy)*dyrcp) ((p1y-oy)*dyrcp)- Interval inz outz = if dz > 0- then Interval ((p1z-oz)*dzrcp) ((p2z-oz)*dzrcp)- else Interval ((p2z-oz)*dzrcp) ((p1z-oz)*dzrcp)- lastin = (fmax3 inx iny inz)- firstout = (fmin3 outx outy outz)- in if lastin > firstout || firstout < 0 || lastin > d- then RayMiss- else - let n = if inx == lastin - then if dx > 0 then nvx else vx- else if iny == lastin- then if dy > 0 then nvy else vy- else if dz > 0 then nvz else vz- norm = if lastin > 0 then n else vinvert n- hitdepth = fmax 0 lastin- in- RayHit hitdepth (vscaleadd orig dir hitdepth) norm t ----Composite objects-----Instance----- transforming the distance is a little awkward--- the normal shouldn't have to be re-normalized, should it?-rayint (Instance sld xfm) (Ray orig dir) d t =- let newdir = invxfm_vec xfm dir- neworig = invxfm_point xfm orig- lenscale = vlen newdir- invlenscale = 1/lenscale- in- case (rayint_check sld (Ray neworig (vscale newdir invlenscale)) (d*lenscale) t) of- RayMiss -> RayMiss- RayHit depth pos n tex -> RayHit (depth*invlenscale) (xfm_point xfm pos) (vnorm (invxfm_norm xfm n)) tex- ---Group----- rayint (Group lst) r d t = foldl nearest RayMiss (map (\x -> rayint x r d t) lst)-rayint (Group xs) r d t =- let rig [] = RayMiss- rig (x:xs) = nearest (rayint_check x r d t) (rig xs)- in rig xs----Difference----- csg of object a - object b ---rayint (Difference sa sb) r d t =- let dif = Difference sa sb- Ray orig dir = r- ria = rayint sa r d t- in- case ria of- RayMiss -> RayMiss- RayHit ad ap an at ->- if inside sb orig - then- case rayint sb r d t of- RayMiss -> RayMiss - RayHit bd bp bn bt ->- if bd < ad - then if inside sa bp - then RayHit bd bp (vinvert bn) bt- else rayint_advance dif r d t bd- else rayint_advance dif r d t bd- else - if inside sb ap- then rayint_advance dif r d t ad- else RayHit ad ap an at-----Intersection----- fixme: there's some numerical instability near edges-rayint (Intersection slds) r d t =- let (Ray orig dir) = r - in- if null slds || d < 0- then RayMiss- else - let s = head slds - in case tail slds of- [] -> rayint_check s r d t- ss -> if inside s orig- then case rayint s r d t of - RayMiss -> rayint (Intersection ss) r d t- RayHit sd sp sn st -> - case rayint_check (Intersection ss) r sd t of- RayMiss -> rayint_advance (Intersection slds) - r d t sd - hit -> hit- else case rayint s r d t of- RayMiss -> RayMiss- RayHit sd sp sn st ->- if inside (Intersection ss) sp- then RayHit sd sp sn st- else rayint_advance (Intersection slds)- r d t sd---rayint (Bound sa sb) r d t =- let (Ray orig _) = r- in if inside sa orig || shadow sa r d- then rayint sb r d t- else RayMiss----Bih---rayint (Bih bb root) r d t =- let Ray orig dir = r- dir_rcp = vrcp dir- Interval near far = bbclip r bb- traverse (BihLeaf s) near far = rayint_check s r (fmin d far) t- traverse (BihBranch lsplit rsplit axis l r) near far =- let dirr = va dir_rcp axis- o = va orig axis- dl = (lsplit - o) * dirr- dr = (rsplit - o) * dirr- in - if near > far - then RayMiss- else- -- this is ugly and verbose, - -- but it does what it needs to- if dirr > 0- then - (nearest- (if near < dl- then traverse l near (fmin dl far)- else RayMiss)- (if dr < far- then traverse r (fmax dr near) far- else RayMiss))- else- (nearest- (if near < dr- then traverse r near (fmin dr far)- else RayMiss)- (if dl < far- then traverse l (fmax dl near) far- else RayMiss)) --}-- in- traverse root near far---- anything that hits it gets teleported--- one-way door to another world--- this is dangerous and doesn't really work right-rayint (Portal port world) r d t =- case rayint port r d t of- RayMiss -> RayMiss- RayHit depth _ _ _ -> - rayint_advance world r d t depth----Tex----- this is a little odd; rather than associate--- a texture with each primitive, we use a--- container object; everything inside has that--- texture, unless it's overridden by a nested--- texture-rayint (Tex s tex) r d t = rayint_check s r d tex----Nothing---rayint (Solid.Nothing) _ _ _ = RayMiss---- default case: miss--- rayint _ _ _ _ = RayMiss---- various specialized ray intersections, used as helper functions---- used by cylinder / cone code--- broken, do not use-rayint_aadisc :: Flt -> Flt-> Ray -> Flt -> Texture -> Rayint-rayint_aadisc height radius r d t =- let Ray orig dir = r- -- Vec _ _ oz = orig- -- Vec _ _ dz = dir- -- dist = (height-oz)/dz- dist = plane_int_dist r (Vec 0 0 height) vx- in if dist <= 0 || dist >= d || isNaN dist - then RayMiss- else let pos = vscaleadd orig dir dist- (Vec px py _) = pos- in - if (px*px + py+py) > radius*radius- then RayMiss- else RayHit dist pos vz t---- move ray forward, intersect, fix result--- useful in csg-rayint_advance :: Solid -> Ray -> Flt -> Texture -> Flt -> Rayint-rayint_advance s r d t adv =- let a = adv+delta- in- case (rayint s (ray_move r a) (d-a) t) of- RayMiss -> RayMiss- RayHit depth pos norm tex -> RayHit (depth+a) pos norm tex---- check results of a ray-intersection test-rayint_check s r d t =- let Ray orig dir = r in- case rayint s r d t of- RayMiss -> RayMiss- RayHit depth pos norm tex ->- if depth < 0 - then error $ "rayint depth < 0 " ++ (show depth) ++ " " ++ (show s)- else- if depth > d- then error $ "rayint depth (" ++ (show depth) ++ ") > d (" ++ (show d) ++ ") " ++ (show s)- else - if not $ veq pos (vscaleadd orig dir depth)- then error $ "rayint position doesn't match depth" ++ - (show pos) ++ (show $ vscaleadd orig dir depth) ++ (show depth) ++ (show s)- else - if (vlen norm) < 1-delta - then error "normal too short"- else - if (vlen norm) > 1+delta- then error $ "normal too long " ++ (show norm) ++ " " ++ (show s)- else RayHit depth pos norm tex-----SHADOW---shadow :: Solid -> Ray -> Flt -> Bool----Sphere---shadow (Sphere center r invr) (Ray e dir) dist = - let eo = vsub center e- v = vdot eo dir- in- if (dist >= (v - r)) && (v > 0.0)- then- let vsqr = v*v- csqr = vdot eo eo- rsqr = r*r- disc = rsqr - (csqr - vsqr) in- if disc < 0.0 then- False- else- let d = sqrt disc- hitdist = if (v-d) > 0 then (v-d) else (v+d)- in if (hitdist < 0) || (hitdist > dist)- then False- else True- else- False---shadow (Triangle p1 p2 p3) (Ray o dir) dist =- let e1 = vsub p2 p1- e2 = vsub p3 p1- s1 = vcross dir e2- divisor = vdot s1 e1- in - if (divisor == 0)- then False- else- let invdivisor = 1.0 / divisor- d = vsub o p1 - b1 = (vdot d s1) * invdivisor- in- if (b1 < 0) || (b1 > 1) - then False- else- let s2 = vcross d e1- b2 = (vdot dir s2) * invdivisor- in- if (b2 < 0) || (b1 + b2 > 1) - then False- else- let t = (vdot e2 s2) * invdivisor- in- if (t < 0) || (t > dist)- then False- else True--shadow (TriangleNorm p1 p2 p3 n1 n2 n3) r d =- shadow (Triangle p1 p2 p3) r d--shadow (Disc point norm radius_sqr) r d =- let (Ray orig dir) = r- dist = plane_int_dist r point norm - in if dist < 0 || dist > d - then False- else let pos = vscaleadd orig dir dist- offset = vsub pos point- in - if (vdot offset offset) > radius_sqr- then False- else True--shadow (Box box) r d =- let Interval near far = bbclip r box - in- if (near > far) || far <= 0 || far > d- then False- else True---- cone centered on z axiz, height of hp, clipped at h1 and h2-shadow (Cone r clip1 clip2 height) (Ray orig dir) d =- let Vec ox oy oz = orig- Vec dx dy dz = dir- k' = (r/height)- k = k'*k'- a = dx*dx + dy*dy - k*dz*dz- b = 2*(dx*ox + dy*oy - k*dz*(oz-height))- c = ox*ox + oy*oy - k*(oz-height)*(oz-height)- disc = b*b - 4*a*c- in if disc < 0- then False- else- let discsqrt = sqrt disc- q = if b < 0- then (b-discsqrt)*(-0.5)- else (b+discsqrt)*(-0.5)- t0' = q/a- t1' = c/q- t0 = fmin t0' t1'- t1 = fmax t0' t1'- in if t1 < 0 || t0 > d - then False- else let dist = if t0 < 0- then t1- else t0- in if dist < 0 || dist > d- then False- else- let pos = vscaleadd orig dir dist- Vec posx posy posz = pos- in if posz > clip1 && posz < clip2- then True- else - if dz > 0 -- ray pointing up from bottom- then if oz < clip1- then shadow (Disc (Vec 0 0 clip1) nvz (r*r)) (Ray orig dir) d- else False- else if oz > clip2- then let r2 = r*(1-((clip2-clip1)/(height)))- in shadow (Disc (Vec 0 0 clip2) vz (r2*r2)) (Ray orig dir) d- else False--shadow (Instance sld xfm) (Ray orig dir) d =- let newdir = invxfm_vec xfm dir- neworig = invxfm_point xfm orig- lenscale = vlen newdir- invlenscale = 1/lenscale- in- shadow sld (Ray neworig (vscale newdir invlenscale)) (d*lenscale)--shadow (Tex s t) r d = shadow s r d--shadow (Group xs) r d =- let sg [] = False- sg (x:xs) = (shadow x r d) || (sg xs)- in sg xs--shadow (Bound sa sb) r d =- let (Ray orig _ ) = r- in if inside sa orig || shadow sa r d- then shadow sb r d- else False--shadow (Bih bb root) r d =- let (Ray orig dir) = r- dir_rcp = vrcp dir- Interval near far = bbclip r bb- traverse (BihLeaf s) near far = shadow s r (fmin d far)- traverse (BihBranch lsplit rsplit axis l r) near far =- let dirr = va dir_rcp axis- o = va orig axis- dl = (lsplit - o) * dirr- dr = (rsplit - o) * dirr- in - if near > far - then False- else- if dirr > 0- then- ((if near < dl- then traverse l near (fmin dl far)- else False) - ||- (if dr < far- then traverse r (fmax dr near) far- else False))- else- ((if near < dr- then traverse r near (fmin dr far)- else False)- ||- (if dl < far- then traverse l (fmax dl near) far- else False))-- in traverse root near far---- default shadow test, in case an optimized test --- isn't implemented yet, we use the regular --- trace function; we have to be careful with --- container objects, though; using the slow --- test on the container means doing the slow--- test against everything it contains-shadow s r d =- case (rayint s r d t_white) of- RayHit _ _ _ _ -> True- RayMiss -> False----INSIDE---inside :: Solid -> Vec -> Bool-inside (Sphere center r invr) pt =- let offset = vsub center pt - in (vdot offset offset) < r*r--inside (Group slds) pt =- foldl' (||) False (map (\x -> inside x pt) slds)--inside (Difference sa sb) pt =- (inside sa pt) && (not $ inside sb pt)---- note: inside is True for an empty intersection.--- this is actually the preferred semantics in --- some cases, strange as it may seem.-inside (Intersection slds) pt =- foldl' (&&) True (map (\x -> inside x pt) slds)--inside (Instance s xfm) pt =- inside s (xfm_point xfm pt)--inside (Plane norm offset) pt =- let onplane = (vscale norm offset)- newvec = vsub onplane pt- in vdot newvec norm > 0--inside (Box (Bbox (Vec x1 y1 z1) (Vec x2 y2 z2))) (Vec x y z) =- x > x1 && x < x2 && y > y1 && y < y2 && z > z1 && z < z2--inside (Tex s t) pt = inside s pt --inside (Bih (Bbox (Vec x1 y1 z1) (Vec x2 y2 z2)) root) pt =- let (Vec x y z) = pt- traverse (BihLeaf s) = inside s pt- traverse (BihBranch lsplit rsplit axis l r) =- let o = va pt axis- in (if o < lsplit- then (traverse l)- else False) - ||- (if o > rsplit - then (traverse r)- else False)- in- (x > x1) && (x < x2) && (y > y1) && (y < y2) && (z > z1) && (z < z2) && (traverse root)--inside (Cylinder r h1 h2) (Vec x y z) =- z > h1 && z < h2 && x*x + y*y < r*r- -inside (Cone rbase h1 h2 height) (Vec x y z) =- let r = rbase*(1-(((z-h1)/height)))- in z > h1 && z < h2 && x*x + y*y < r*r--inside (Bound sa sb) pt = inside sa pt && inside sb pt--inside _ _ = False---- return distance to surface, positive if inside, negative if outside--- this isn't used for anything in particular yet-power :: Solid -> Vec -> Flt--power (Sphere center r invr) pt =- let offset = vsub center pt- in r - (vlen offset)--power (Group slds) pt =- foldl' (max) (-infinity) (map (\x -> power x pt) slds)---- not accurate-power (Instance s xfm) pt =- power s (xfm_point xfm pt)--power (Plane norm offset) pt =- let onplane = (vscale norm offset)- newvec = vsub onplane pt- in vdot newvec norm-----BOUND---bound :: Solid -> Bbox-bound (Sphere center r invr) =- let offset = (vec r r r) in- (Bbox (vsub center offset) (vadd center offset))--bound (Triangle (Vec v1x v1y v1z) - (Vec v2x v2y v2z) - (Vec v3x v3y v3z)) =- Bbox- (Vec ((fmin (fmin v1x v2x) v3x) - delta)- ((fmin (fmin v1y v2y) v3y) - delta)- ((fmin (fmin v1z v2z) v3z) - delta) )-- (Vec ((fmax (fmax v1x v2x) v3x) + delta)- ((fmax (fmax v1y v2y) v3y) + delta)- ((fmax (fmax v1z v2z) v3z) + delta) )--bound (TriangleNorm v1 v2 v3 n1 n2 n3) =- bound (Triangle v1 v2 v3)----dangerous to use inside a bih-bound (Plane norm offset) = everything_bbox--bound (Box box) = box---- this could be a tighter fit-bound (Disc pos norm rsqr) =- bound (sphere pos (sqrt rsqr))--bound (Cylinder r h1 h2) =- Bbox (Vec (-r) (-r) h1) (Vec r r h2)--bound (Cone r h1 h2 height) =- Bbox (Vec (-r) (-r) h1) (Vec r r h2)--bound (Group slds) = - foldl' bbjoin empty_bbox (map bound slds)--bound (Difference sa sb) = bound sa---bound (Intersection slds) =- if null slds - then empty_bbox- else foldl' bboverlap everything_bbox (map bound slds)------ the reason the following doesn't work:--- bound (Bound sa sb) = bboverlap (bound sa) (bound sb)--- is that the ray has to hit solid a, and it might not --- if the bound is smaller than sa.---- update: this is probably no longer true, as I added an "inside"--- test to the bounds check as well--bound (Bound sa sb) = bound sa---- not optimal, but it does the job-bound (Instance sld xfm) =- let (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bound sld- pxfm = xfm_point xfm- in- bbpts [(pxfm (Vec x y z)) | x <- [p1x,p2x],- y <- [p1y,p2y],- z <- [p1z,p2z]]---bound (Bih bb root) = bb--bound (Portal port world) = bound port--bound (Tex s t) = bound s --bound (Photon p indir clr r) = bound (sphere p r)--bound Solid.Nothing = empty_bbox---- no default case: we want an exception --- if there is no match-----CAMERA---data Camera = Camera {campos, fwd, up, right :: !Vec} - deriving Show--default_cam = (Camera (vec 0.0 0.0 (-3.0)) - (vec 0.0 0.0 1.0) - (vec 0.0 1.0 0.0) - (vec 1.0 0.0 0.0) )--camera :: Vec -> Vec -> Vec -> Flt -> Camera-camera pos at up angle =- let fwd = vnorm $ vsub at pos- right = vnorm $ vcross up fwd- up_ = vnorm $ vcross fwd right- cam_scale = tan ((pi/180)*(angle/2))- in- Camera pos fwd- (vscale up_ cam_scale) - (vscale right cam_scale)----SCENE---data Scene = Scene {sld :: !Solid, - lights :: ![Light], - cam :: !Camera, - dtex :: !Texture, - bground :: !Color} deriving Show--default_scene = (Scene (sphere (vec 0.0 0.0 0.0) 1.0) - [] default_cam t_white c_white)+--COMMON DATATYPES AND UTILITY FUNCTIONS--+data Bbox = Bbox {p1 :: !Vec, p2 :: !Vec} deriving Show+data Interval = Interval !Flt !Flt deriving Show -- used instead of a tuple++--union of two bounding boxes+bbjoin :: Bbox -> Bbox -> Bbox+bbjoin (Bbox p1a p2a) (Bbox p1b p2b) =+ (Bbox (vmin p1a p1b) (vmax p2a p2b))++--overlap of two bounding boxes+bboverlap :: Bbox -> Bbox -> Bbox+bboverlap (Bbox p1a p2a) (Bbox p1b p2b) =+ (Bbox (vmax p1a p1b) (vmin p2a p2b))++--split a bounding box into two+bbsplit :: Bbox -> Int -> Flt -> (Bbox,Bbox)+bbsplit (Bbox p1 p2) axis offset =+ if (offset < (va p1 axis)) || (offset > (va p2 axis))+ then error "degenerate bounding box split"+ else ((Bbox p1 (vset p2 axis offset)),+ (Bbox (vset p1 axis offset) p2))++-- generate a bounding box from a list of points+bbpts :: [Vec] -> Bbox+bbpts [] = empty_bbox+bbpts ((Vec x y z):[]) =+ Bbox (Vec (x-delta) (y-delta) (z-delta)) + (Vec (x+delta) (y+delta) (z+delta))++bbpts ((Vec x y z):pts) =+ let (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bbpts pts+ minx = fmin (x-delta) p1x+ miny = fmin (y-delta) p1y+ minz = fmin (z-delta) p1z+ maxx = fmax (x+delta) p2x+ maxy = fmax (y+delta) p2y+ maxz = fmax (z+delta) p2z in+ Bbox (Vec minx miny minz) (Vec maxx maxy maxz)++-- surface area, volume of bounding boxes+bbsa :: Bbox -> Flt+bbsa (Bbox p1 p2) =+ let Vec dx dy dz = vsub p2 p1 + in dx*dy + dx*dz + dy*dz++bbvol :: Bbox -> Flt+bbvol (Bbox p1 p2) =+ let (Vec dx dy dz) = vsub p2 p1+ in dx*dy*dz++empty_bbox = + Bbox (Vec infinity infinity infinity) + (Vec (-infinity) (-infinity) (-infinity))++everything_bbox =+ Bbox (Vec (-infinity) (-infinity) (-infinity))+ (Vec infinity infinity infinity)++-- Find a ray's entrance and exit from a bounding +-- box. If last entrance is before the first exit,+-- we hit. Otherwise, we miss. (It's up to the +-- caller to figure that out.)++bbclip :: Ray -> Bbox -> Interval+bbclip (Ray (Vec ox oy oz) (Vec dx dy dz)) + (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) =+ let dxrcp = 1/dx+ dyrcp = 1/dy+ dzrcp = 1/dz+ Interval inx outx = if dx > 0 + then Interval ((p1x-ox)*dxrcp) ((p2x-ox)*dxrcp)+ else Interval ((p2x-ox)*dxrcp) ((p1x-ox)*dxrcp)+ Interval iny outy = if dy > 0+ then Interval ((p1y-oy)*dyrcp) ((p2y-oy)*dyrcp)+ else Interval ((p2y-oy)*dyrcp) ((p1y-oy)*dyrcp)+ Interval inz outz = if dz > 0+ then Interval ((p1z-oz)*dzrcp) ((p2z-oz)*dzrcp)+ else Interval ((p2z-oz)*dzrcp) ((p1z-oz)*dzrcp)+ in+ Interval (fmax3 inx iny inz) (fmin3 outx outy outz)++data Rayint = RayHit {+ depth :: !Flt,+ pos :: !Vec,+ norm :: !Vec,+ texture :: !Texture+} | RayMiss deriving Show++nearest :: Rayint -> Rayint -> Rayint+nearest a RayMiss = a+nearest RayMiss b = b+nearest (RayHit da pa na ta) (RayHit db pb nb tb) =+ if da < db+ then RayHit da pa na ta+ else RayHit db pb nb tb++furthest :: Rayint -> Rayint -> Rayint+furthest a RayMiss = RayMiss+furthest RayMiss b = RayMiss+furthest (RayHit da pa na ta) (RayHit db pb nb tb) =+ if da > db+ then RayHit da pa na ta+ else RayHit db pb nb tb++hit :: Rayint -> Bool+hit (RayHit _ _ _ _) = True+hit RayMiss = False++dist :: Rayint -> Flt+dist RayMiss = infinity+dist (RayHit d _ _ _) = d++--Packet Types--+data PacketResult = PacketResult Rayint Rayint Rayint Rayint+packetmiss = PacketResult RayMiss RayMiss RayMiss RayMiss+++nearest_packetresult :: PacketResult -> PacketResult -> PacketResult+nearest_packetresult (PacketResult a1 a2 a3 a4) (PacketResult b1 b2 b3 b4) =+ PacketResult (nearest a1 b1)+ (nearest a2 b2)+ (nearest a3 b3)+ (nearest a4 b4)++-- move ray forward, intersect, fix result+-- useful in csg+rayint_advance :: SolidItem -> Ray -> Flt -> Texture -> Flt -> Rayint+rayint_advance s r d t adv =+ let a = adv+delta+ in+ case (rayint s (ray_move r a) (d-a) t) of+ RayMiss -> RayMiss+ RayHit depth pos norm tex -> RayHit (depth+a) pos norm tex+++--MATERIALS--+data Material = Material {clr :: Color, + refl, refr, ior, + kd, shine :: !Flt} deriving Show+type Texture = Rayint -> Material++-- this is sort of a no-op; we don't have a+-- good way to show an arbitrary function+showTexture :: Texture -> String+showTexture t = show $ t RayMiss++instance Show Texture where+ show = showTexture++m_white = (Material c_white 0 0 0 1 2)+t_white ri = m_white++t_uniform :: Material -> Texture+t_uniform m = \x -> m++interp :: Flt -> Flt -> Flt -> Flt+interp scale a b =+ scale*a + (1-scale)*b++--not really correct, but we'll go with it for now+m_interp :: Material -> Material -> Flt -> Material+m_interp m1 m2 scale =+ let (Material m1c m1refl m1refr m1ior m1kd m1shine) = m1+ (Material m2c m2refl m2refr m2ior m2kd m2shine) = m2+ intp = interp scale+ c = cadd (cscale m1c scale) (cscale m2c (1-scale))+ refl = intp m1refl m2refl+ refr = intp m1refr m2refr+ ior = intp m1ior m2ior+ kd = intp m1kd m2kd+ shine = intp m1shine m2shine+ in (Material c refl refr ior kd shine)++--SOLID CLASS--++class (Show a) => Solid a where+ rayint :: a -> Ray -> Flt -> Texture -> Rayint+ packetint :: a -> Ray -> Ray -> Ray -> Ray -> Flt -> Texture -> PacketResult + shadow :: a -> Ray -> Flt -> Bool+ inside :: a -> Vec -> Bool+ bound :: a -> Bbox+ tolist :: a -> [SolidItem]+ transform :: a -> [Xfm] -> SolidItem+ flatten_transform :: a -> SolidItem++ -- Sometimes, we can improve performance by + -- intersecting 4 rays at once. This is + -- especially true of acceleration structures.+ -- By default, we fall back on mono-rays.+ + packetint s r1 r2 r3 r4 d t = + PacketResult (rayint s r1 d t)+ (rayint s r2 d t)+ (rayint s r3 d t)+ (rayint s r4 d t)++ -- if there is no shadow function, we fall back on rayint+ shadow s r d =+ case (rayint s r d t_white) of+ RayHit _ _ _ _ -> True+ RayMiss -> False++ -- This is here so we can flatten a group of groups+ -- into a single group; the default is fine for everything+ -- but groups and Void and SolidItem+ tolist a = [SolidItem (a)]+ + -- Method to transform an object; the default works fine+ -- except for instances themselves, which will want to+ -- collapse the two transformations into a sigle transform.+ transform a xfm = SolidItem $ Instance (SolidItem a) (compose xfm)++ -- This prepares a composite primitive to be fed into the bih constructor+ -- by pushing all the transformations out to the leaves and + -- throwing away manual bounding structures.+ flatten_transform a = SolidItem a++--Existential type so we can make a heterogeneous list of solids+--http://notes-on-haskell.blogspot.com/2007/01/proxies-and-delegation-vs-existential.html++data SolidItem = forall a. Solid a => SolidItem a++instance Solid SolidItem where+ rayint (SolidItem s) r d t = rayint s r d t+ shadow (SolidItem s) r d = shadow s r d+ inside (SolidItem s) pt = inside s pt+ bound (SolidItem s) = bound s+ tolist s = [s] -- don't wrap in a redundant SolidItem like everything else+ transform (SolidItem s) xfm = transform s xfm -- same here+ flatten_transform (SolidItem s) = (SolidItem (flatten_transform s)) -- and here++instance Show SolidItem where+ show (SolidItem s) = "SI " ++ show s++-- we implement "group", "void", and "instance" here because they're+-- used by some of the other primitives++-- GROUP --++group :: [SolidItem] -> SolidItem+group [] = SolidItem Void+group (sld:[]) = sld+group slds = SolidItem (flatten_group slds)++-- smash a group of groups into a single group,+-- so we can build an efficient bounding heirarchy+flatten_group :: [SolidItem] -> [SolidItem]+flatten_group slds = concat (map tolist slds)++-- this lets us treat lists of SolidItems as regular Solids+rayint_group :: [SolidItem] -> Ray -> Flt -> Texture -> Rayint+rayint_group [] _ _ _ = RayMiss+rayint_group (x:xs) r d t = nearest (rayint x r d t) (rayint_group xs r d t)++shadow_group :: [SolidItem] -> Ray -> Flt -> Bool+shadow_group [] r d = False+shadow_group (x:xs) r d = (shadow x r d) || (shadow_group xs r d)++inside_group :: [SolidItem] -> Vec -> Bool+inside_group slds pt =+ foldl' (||) False (map (\x -> inside x pt) slds)++bound_group :: [SolidItem] -> Bbox+bound_group slds = + foldl' bbjoin empty_bbox (map bound slds)++flatten_transform_group :: [SolidItem] -> SolidItem+flatten_transform_group slds =+ SolidItem $ map flatten_transform slds++instance Solid [SolidItem] where+ rayint = rayint_group+ shadow = shadow_group+ inside = inside_group+ bound = bound_group+ tolist a = concat $ map tolist a++-- VOID --+-- non-object (originally called "Nothing", but that+-- conflicted with the prelude maybe type, so we call+-- it "Void" instead) +data Void = Void deriving Show++nothing = SolidItem Void++instance Solid Void where+ rayint Void r d t = RayMiss+ shadow Void r d = False+ inside Void pt = False+ bound Void = empty_bbox+ tolist Void = [] +++-- INSTANCE --+-- this would be better in its own module, but we need+-- "Instance" to be defined here for the default implementation+-- of "transform". (I tried mutually recursive modules, it+-- didn't work. http://www.haskell.org/ghc/docs/latest/html/+-- users_guide/separate-compilation.html#mutual-recursion ) ++-- Another good reason to include Instance in Solid.hs+-- is that it's referenced from Cone.hs++-- An instance is a primitive that has been modified+-- by a transformation (i.e. some combination of+-- translation, rotation, and scaling). This is a+-- reasonably space-efficient way of making multiple copies+-- of a complex object.++-- It's unfortunate that "instance" is also a reserved word. +-- "instance Solid Instance where..." is a little confusing.++data Instance = Instance SolidItem Xfm deriving Show++rayint_instance :: Instance -> Ray -> Flt -> Texture -> Rayint+rayint_instance (Instance sld xfm) (Ray orig dir) d t =+ let newdir = invxfm_vec xfm dir+ neworig = invxfm_point xfm orig+ lenscale = vlen newdir+ invlenscale = 1/lenscale+ in+ case (rayint sld (Ray neworig (vscale newdir invlenscale)) (d*lenscale) t) of+ RayMiss -> RayMiss+ RayHit depth pos n tex -> RayHit (depth*invlenscale) + (xfm_point xfm pos) + (vnorm (invxfm_norm xfm n)) + tex++shadow_instance :: Instance -> Ray -> Flt -> Bool+shadow_instance (Instance sld xfm) (Ray orig dir) d =+ let newdir = invxfm_vec xfm dir+ neworig = invxfm_point xfm orig+ lenscale = vlen newdir+ invlenscale = 1/lenscale+ in+ shadow sld (Ray neworig (vscale newdir invlenscale)) (d*lenscale)++inside_instance :: Instance -> Vec -> Bool+inside_instance (Instance s xfm) pt =+ inside s (xfm_point xfm pt)++bound_instance :: Instance -> Bbox+bound_instance (Instance sld xfm) =+ let (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bound sld+ pxfm = xfm_point xfm+ in+ bbpts [(pxfm (Vec x y z)) | x <- [p1x,p2x],+ y <- [p1y,p2y],+ z <- [p1z,p2z]]++-- If we try to create a transformation of+-- a transformation, we can merge those+-- into a single transformation.++-- This ought to be tested to verify this+-- is really applying transforms in the+-- correct order...++transform_instance :: Instance -> [Xfm] -> SolidItem+transform_instance (Instance s xfm2) xfm1 =+ transform s [compose ([xfm2]++xfm1) ]++-- Flatten_transform attempts to push all transformations +-- in a heirarchy out to the leaf nodes. The case we're+-- interested in here is an instance of a group, and we +-- want to replace that with a group of individually +-- transformed instances. This could be construed as a+-- waste of memory, but in some cases it's necessary.++flatten_transform_instance :: Instance -> SolidItem+flatten_transform_instance (Instance s xfm) = + group $ map (\x -> transform (flatten_transform x) [xfm]) (tolist s)++instance Solid Instance where+ rayint = rayint_instance+ shadow = shadow_instance+ inside = inside_instance+ bound = bound_instance+ transform = transform_instance+ flatten_transform = flatten_transform_instance
SolidTexture.hs view
@@ -49,7 +49,7 @@ -- but t^5 works and t^6 doesn't. omega :: Flt -> Flt omega t_ = - let t = abs t_+ let t = fabs t_ tsqr = t*t tcube = tsqr*t in (-6)*tcube*tsqr + 15*tcube*t - 10*tcube + 1@@ -67,9 +67,9 @@ gamma :: Int -> Int -> Int -> Vec gamma i j k =- let a = phi!(mod (abs k) 12)- b = phi!(mod (abs (j+a)) 12)- c = phi!(mod (abs (i+b)) 12)+ let a = phi!(mod (iabs k) 12)+ b = phi!(mod (iabs (j+a)) 12)+ c = phi!(mod (iabs (i+b)) 12) in grad!c knot :: Int -> Int -> Int -> Vec -> Flt@@ -79,15 +79,15 @@ intGamma :: Int -> Int -> Int intGamma i j =- let a = phi!(mod (abs j) 16)- b = phi!(mod (abs (i+a)) 16)+ let a = phi!(mod (iabs j) 16)+ b = phi!(mod (iabs (i+a)) 16) in b turbulence :: Vec -> Int -> Flt-turbulence p 1 = abs(noise(p))+turbulence p 1 = fabs(noise(p)) turbulence p n = let newp = vscale p 0.5- t = abs (noise p)+ t = fabs (noise p) in t + (0.5 * (turbulence newp (n-1))) noise :: Vec -> Flt
Spd.hs view
@@ -1,7 +1,5 @@ module Spd where-import Vec-import Clr-import Solid+import Scene -- NFF file format description: -- http://tog.acm.org/resources/SPD/NFF.TXT@@ -155,7 +153,7 @@ -- vert1.x vert1.y vert1.z -- [etc. for total_vertices vertices] -readsSpdSolid :: ReadS Solid+readsSpdSolid :: ReadS SolidItem readsSpdSolid s = [((sphere center radius),s3) | ("s", s1) <- lexcr s, (center,s2) <- reads s1 :: [(Vec,String)], (radius,s3) <- reads s2 :: [(Flt,String)] ]@@ -174,7 +172,7 @@ (n,s2) <- reads s1 :: [(Int,String)], (vns,s3) <- readsSpdVecsNorms s2 :: [([(Vec,Vec)],String)] ] {- ++- [(Tex(Solid.Nothing,t),s1) | (t,s1) <- reads s :: [(Texture,String)]] -}+ [(tex(Void,t),s1) | (t,s1) <- reads s :: [(Texture,String)]] -} -- instance Read Solid where@@ -182,9 +180,9 @@ -- same as readSpdVecs, just different types-readsSpdPrims :: ReadS [Solid]+readsSpdPrims :: ReadS [SolidItem] readsSpdPrims s =- let parses = readsSpdSolid s :: [(Solid,String)]+ let parses = readsSpdSolid s :: [(SolidItem,String)] in if null parses then [([],s)]@@ -193,25 +191,25 @@ (vs,returns) = head (readsSpdPrims rest) in [((v:vs),returns)] -instance Read [Solid] where+instance Read [SolidItem] where readsPrec _ = readsSpdPrims -readsSpdTextureGroup :: ReadS Solid+readsSpdTextureGroup :: ReadS SolidItem readsSpdTextureGroup s =- [((Tex (bih prims) t),s2) | (t,s1) <- reads s :: [(Texture,String)],- (prims,s2) <- readsSpdPrims s1 :: [([Solid],String)] ]+ [((tex (bih prims) t),s2) | (t,s1) <- reads s :: [(Texture,String)],+ (prims,s2) <- readsSpdPrims s1 :: [([SolidItem],String)] ] -instance Read Solid where+instance Read SolidItem where readsPrec _ = readsSpdTextureGroup -accum_rss :: [Camera] -> [Light] -> [Solid] -> [BgColor] -> String -> ([Camera],[Light],[Solid],[BgColor],String)+accum_rss :: [Camera] -> [Light] -> [SolidItem] -> [BgColor] -> String -> ([Camera],[Light],[SolidItem],[BgColor],String) accum_rss cams lights prims background s = if null s then (cams,lights,prims,background,s) else let cam = reads s :: [(Camera,String)]- sld = reads s :: [(Solid,String)]+ sld = reads s :: [(SolidItem,String)] lit = reads s :: [(Light,String)] bgc = reads s :: [(BgColor,String)] in@@ -240,7 +238,7 @@ readsSpdScene :: ReadS Scene readsSpdScene s = let ((cam:cams),lights,prims,(BgColor(bgc):bgcs),s1) = accum_rss [] [] [] [] s- in [((Scene (bih prims) lights cam t_white bgc),s1)]+ in [((scene (bih prims) lights cam t_white bgc),s1)] instance Read Scene where readsPrec _ = readsSpdScene
+ Sphere.hs view
@@ -0,0 +1,75 @@+module Sphere (sphere) where+import Vec+import Solid++data Sphere = Sphere Vec Flt Flt deriving Show++sphere :: Vec -> Flt -> SolidItem+sphere c r =+ SolidItem (Sphere c r (1.0/r))++-- adapted from graphics gems volume 1+rayint_sphere :: Sphere -> Ray -> Flt -> Texture -> Rayint+rayint_sphere (Sphere center r invr) (Ray e dir) dist t = + let eo = vsub center e+ v = vdot eo dir+ in+ if (dist >= (v - r)) && (v > 0.0)+ then+ let vsqr = v*v+ csqr = vdot eo eo+ rsqr = r*r+ disc = rsqr - (csqr - vsqr) in+ if disc < 0.0 then+ RayMiss+ else+ let d = sqrt disc+ hitdist = if (v-d) > 0 then (v-d) else (v+d)+ in if (hitdist < 0) || (hitdist > dist)+ then RayMiss+ else+ let p = vscaleadd e dir hitdist+ -- n = vscale (vsub p center) invr in+ -- n = vsub (vscale p invr) (vscale center invr) in+ n = vnorm (vsub p center) + in RayHit hitdist p n t+ else+ RayMiss++shadow_sphere :: Sphere -> Ray -> Flt -> Bool+shadow_sphere (Sphere center r invr) (Ray e dir) dist = + let eo = vsub center e+ v = vdot eo dir+ in+ if (dist >= (v - r)) && (v > 0.0)+ then+ let vsqr = v*v+ csqr = vdot eo eo+ rsqr = r*r+ disc = rsqr - (csqr - vsqr) in+ if disc < 0.0 then+ False+ else+ let d = sqrt disc+ hitdist = if (v-d) > 0 then (v-d) else (v+d)+ in if (hitdist < 0) || (hitdist > dist)+ then False+ else True+ else+ False++inside_sphere :: Sphere -> Vec -> Bool+inside_sphere (Sphere center r invr) pt =+ let offset = vsub center pt + in (vdot offset offset) < r*r++bound_sphere :: Sphere -> Bbox+bound_sphere (Sphere center r invr) =+ let offset = (vec r r r) in+ (Bbox (vsub center offset) (vadd center offset))++instance Solid Sphere where + rayint = rayint_sphere+ shadow = shadow_sphere+ inside = inside_sphere+ bound = bound_sphere
TestScene.hs view
@@ -1,24 +1,19 @@-module TestScene where-import Vec-import Solid-import Clr+module TestScene (scn) where+import Scene import Data.List hiding (group) import SolidTexture import System.Random -{--lits = [Light {litpos = Vec {x = -3.0, y = 1.7, z = 5.0}, - litcol = Color 0.9 1 1 },- Light {litpos = Vec {x = 1.1, y = -4.2, z = 4.0}, - litcol = Color 1 0.9 1 },- Light {litpos = Vec {x = 4.3, y = 3.0, z = 2.3}, - litcol = Color 1 1 0.9}]--}--lits = [ Light (Vec (-100) 70 (140)) (cscale (Color 1 0.8 0.8) 1500)+lights = [ Light (Vec (-100) 70 (140)) (cscale (Color 1 0.8 0.8) 2500) , Light (Vec (-3) 5 8) (Color 1.5 2 2) ] +lattice = + let n = 15 :: Flt+ in bih [sphere (vec x y z) 0.2 | x <- [(-n)..n],+ y <- [(-n)..n],+ z <- [(-n)..n]]+ icosahedron pos r = let gr = ((1+(sqrt 5))/2) -- golden ratio, 1.618033988749895 n11 = [(-r),r]@@ -33,9 +28,9 @@ y <- ngrgr] ++ [Vec x 0 z | x <- ngrgr, z <- grrcp]- pln x = (Plane (vnorm x) (r+(vdot (vnorm x) pos)))+ pln x = (plane_offset (vnorm x) (r+(vdot (vnorm x) pos))) in- Intersection ((sphere pos (1.26*r)):(map pln points))+ intersection ((sphere pos (1.26*r)):(map pln points)) dodecahedron pos r = let gr = (1+(sqrt 5))/2 -- golden ratio, 1.618033988749895@@ -44,90 +39,92 @@ points = [Vec 0 y z | y <- n11, z <- ngrgr] ++ [Vec x 0 z | z <- n11, x <- ngrgr] ++ [Vec x y 0 | x <- n11, y <- ngrgr]- pln x = (Plane (vnorm x) (r+(vdot (vnorm x) pos)))+ pln x = (plane_offset (vnorm x) (r+(vdot (vnorm x) pos))) in- Intersection ((sphere pos (1.26*r)):(map pln points))--cust_cam = camera (vec (-2) (4.3) (15)) (vec 0 2 0) (vec 0 1 0) 45-+ intersection ((sphere pos (1.26*r)):(map pln points)) spiral = [ ((Vec ((sin (rot n))*n) ((cos (rot n))*n) (n-3)), (n/15)) | n <- [0, 0.01..6]] -coil = bih (zipWith (\ (p1,r1) (p2,r2) -> (Solid.group [(cone p1 r1 p2 r2), - (sphere p1 r1)] )) +coil = bih (zipWith (\ (p1,r1) (p2,r2) -> (group [(cone p1 r1 p2 r2), + (sphere p1 r1)] )) spiral (tail spiral)) + -- we branch once per year -- not really a plausible oak, but it's getting there-oak :: Flt -> StdGen -> Solid+oak :: Flt -> StdGen -> SolidItem oak age rng = if age < 0 - then Solid.Nothing+ then nothing else let year :: Int = floor age season = age-(fromIntegral year)- thickness = 0.03- minbranch = deg 10- maxbranch = deg 25- tree 0 r = Solid.Nothing+ thickness = 0.025+ minbranch = deg 12+ maxbranch = deg 18+ tree 0 r = nothing tree 1 r = -- cone (Vec 0 0 0) thickness (Vec 0 season 0) 0- Tex (sphere (Vec 0 0 0) season) (t_matte (Color 0.2 1 0.4))- tree n r = let nf = fromIntegral n - height = nf- (rng1,rng2) = split r- (rng3,rng4) = split rng1- (r1,rng5) = randomR (0,0.5) rng4- (r2,rng6) = randomR (minbranch,maxbranch) rng5- (r3,rng7) = randomR (0.8,0.95) rng6- seglen = 0.5 + r1- branchang = r2- scaling = r3+ tex (sphere (Vec 0 0 0) season) (t_matte (Color 0.2 1 0.4))+ tree n_ r = let nf = fromIntegral n_ + height_ = nf+ (rng1,rng2) = split r+ (rng3,rng4) = split rng1+ (r1,rng5) = randomR (0,0.5) rng4+ (r2,rng6) = randomR (minbranch,maxbranch) rng5+ (r3,rng7) = randomR (0.75,0.95) rng6+ (r4,rng8) = randomR (0.0,1.0) rng7+ seglen = 0.5 + r1+ branchang = r2+ scaling = r3+ (height,n) = if r4 > (1 :: Float)+ then ((height_/2),(ceiling (nf/2)))+ else (height_, n_)+ -- we make our own manual bounding heirarchy -- (bih doesn't know what to do with heirachies -- of transformed objects)- in Bound (sphere (Vec 0 (height/2) 0) (height/2))- (group [ cone (Vec 0 0 0) (thickness*nf) (Vec 0 seglen 0) (thickness*(nf-1)*scaling)- , transform (tree (n-1) rng2) [(scale (Vec scaling scaling scaling)),- (rotate (Vec 0 0 1) branchang),- (rotate (Vec 0 1 0) (deg 30)),- (translate (Vec 0 seglen 0))]- , transform (tree (n-1) rng3) [(scale (Vec scaling scaling scaling)),- (rotate (Vec 0 0 1) (-branchang)),- (rotate (Vec 0 1 0) (deg 30)),- (translate (Vec 0 seglen 0))]- ])- in Tex (bih (flatten_transform (tree year rng))) (t_matte (Color 0.8 0.5 0.4))+ in bound_object (sphere (Vec 0 (height/2) 0) (height/2))+ (group [ cone (Vec 0 0 0) (thickness*height) (Vec 0 seglen 0) (thickness*(height-1)*scaling)+ , transform (tree (n-1) rng2) [(scale (Vec scaling scaling scaling)),+ (rotate (Vec 0 0 1) branchang),+ (rotate (Vec 0 1 0) (deg 30)),+ (translate (Vec 0 seglen 0))]+ , transform (tree (n-1) rng3) [(scale (Vec scaling scaling scaling)),+ (rotate (Vec 0 0 1) (-branchang)),+ (rotate (Vec 0 1 0) (deg 30)),+ (translate (Vec 0 seglen 0))]+ ])+ in tex (bih (tolist (flatten_transform (tree year rng)))) (t_matte (Color 0.8 0.5 0.4)) -sphereint = Intersection [ (sphere (Vec (-1) 0 0) 2), +sphereint = intersection [ (sphere (Vec (-1) 0 0) 2), (sphere (Vec 1 0 0) 2), (sphere (Vec 0 (-1) 0) 2), (sphere (Vec 0 1 0) 2) ] --- scene with everything--- can't have the plane inside a bih structure because it's infinite-geom = group [ Tex (plane (Vec 0 0 0) (Vec 0 1 0)) (t_matte (Color 0 0.8 0.3))- , bih [ Tex (dodecahedron (Vec (-6) 3 0) 1) t_stripe- , Tex (icosahedron (Vec 4 1.5 3) 1.5) t_mottled- --,(oak 11.6 (mkStdGen 42))- , transform (oak 11.6 (mkStdGen 42)) [ scale (Vec 1.5 1.5 1.5)]- , Tex (transform (coil) [ scale (Vec (1/3) (1/3) (1/3))+geom = group [ tex (plane (Vec 0 0 0) (Vec 0 1 0)) (t_matte (Color 0 0.8 0.3))+ , bih [ tex (dodecahedron (Vec (-6) 3 0) 1) t_stripe+ , tex (transform (icosahedron (Vec 4 1.5 3) 1.5) [rotate vz (deg 11)+ ,rotate vx (deg 7) ] ) t_mottled+ + , transform (oak 4.6 (mkStdGen 42)) [ scale (Vec 1.5 1.5 1.5)]+ , tex (transform (coil) [ scale (Vec (1/3) (1/3) (1/3)) , rotate (Vec 0 1 0) (deg 65) , translate (Vec (-3.5) 1 (5)) ]) t_mirror --} , cone (Vec (-6) 0 0) 1 (Vec (-6) 3 0) 0- , Tex (Difference (sphere (Vec 0 (-4) 5) 4.7) (sphere (Vec 1.5 (1.5) 5.2) 1.6)) t_mirror- , transform (Tex sphereint (t_matte (Color 0.5 0 1))) [ scale (Vec 0.6 0.6 0.6),+ , tex (difference (sphere (Vec 0 (-4) 5) 4.7) (sphere (Vec 1.5 (1.5) 5.2) 1.6)) t_mirror+ , transform (tex sphereint (t_matte (Color 0.5 0 1))) [ scale (Vec 0.6 0.6 0.6), translate (Vec (-5.2) 1 5)] ] ] -geom1 = transform (oak 7.2 (mkStdGen 42)) [ scale (Vec 1.5 1.5 1.5)]---- color reflect refract ior kd shine +cust_cam = camera (vec (-2) (5.3) (20)) (vec 0 5 0) (vec 0 1 0) 45+ +-- some textures m_shiny_white :: Material m_shiny_white = (Material c_white 0.3 0 0 0.7 10) @@ -143,13 +140,13 @@ in if scale < 0 then error "foo" else if scale > 1 then error "bar"- else m_interp m_shiny_white m_dull_gray scale+ else m_interp m_mirror (m_matte (Color 0 0 1)) scale --shouldn't happen t_mottled RayMiss = m_shiny_white t_stripe (RayHit _ pos norm _) =- let scale = (stripe (Vec 2 4 3) sine_wave) pos+ let scale = (stripe (Vec 4 8 5) triangle_wave) pos in if scale < 0 then error "foo" else if scale > 1 then error "bar"@@ -158,11 +155,20 @@ --shouldn't happen t_stripe RayMiss = m_shiny_white ++m_matte c = (Material c 0 0 0 1 2)+ t_matte c = (\ri -> (Material c 0 0 0 1 2)) +m_mirror = (Material (Color 0.8 0.8 1) 1 0 0 0.2 1000) t_mirror = - (\ri -> (Material (Color 0.8 0.8 1) 1 0 0 0.2 100))+ (\ri -> m_mirror) +c_sky = (Color 0.4 0.5 0.8)+ scn :: IO Scene-scn = return (Scene geom lits cust_cam (t_matte (Color 0.8 0.5 0.4)) (Color 0.4 0.5 0.8))+scn = return (Scene geom+ lights cust_cam + (t_matte (Color 0.8 0.5 0.4)) + c_sky)
+ Tex.hs view
@@ -0,0 +1,47 @@+module Tex (tex) where+import Vec+import Solid++-- Textured objects --+-- The type "Texture" is used elsewhere, so+-- we just call a textured object a "Tex".++-- How textured objects work is a little strange:+-- instead of having a texture applied to every object,+-- which seems rather wastefull, we use the container +-- object "Tex"; anything contained in that Tex has +-- that texture.++-- In the case of nested Tex objects, the innermost +-- texture has precedence. Textures are implemented+-- by passing a Texture in to the rayint function.+-- Most objects just return the Texture unchanged (as+-- part of the RayHit record) but Tex overwrites the +-- texture with its own.++data Tex = Tex SolidItem Texture deriving Show++tex :: SolidItem -> Texture -> SolidItem+tex s t = SolidItem $ Tex s t++rayint_tex :: Tex -> Ray -> Flt -> Texture -> Rayint+rayint_tex (Tex s tex) r d t = rayint s r d tex++packetint_tex :: Tex -> Ray -> Ray -> Ray -> Ray -> Flt -> Texture -> PacketResult+packetint_tex (Tex s tx) r1 r2 r3 r4 d t = packetint s r1 r2 r3 r4 d tx++shadow_tex :: Tex -> Ray -> Flt -> Bool+shadow_tex (Tex s _) r d = shadow s r d++inside_tex :: Tex -> Vec -> Bool+inside_tex (Tex s _) pt = inside s pt++bound_tex :: Tex -> Bbox +bound_tex (Tex s _) = bound s++instance Solid Tex where+ rayint = rayint_tex+ packetint = packetint_tex+ shadow = shadow_tex+ inside = inside_tex+ bound = bound_tex
Trace.hs view
@@ -1,17 +1,19 @@ module Trace where-import Vec-import Clr-import Solid+import Scene import Data.List import Control.Concurrent.MVar import System.IO.Unsafe+--import Packet {-- We put lighting code in this file because it needs to be - mutually recursive with the trace function, for refraction- and reflection.+We put lighting code in this file because it needs to be +mutually recursive with the trace function, for refraction+and reflection. -} +data PacketColor = PacketColor !Color !Color !Color !Color++ {- simple_shade :: Rayint -> [Light] -> Solid -> Color -> Color simple_shade ri lights s bg =@@ -22,21 +24,32 @@ (RayMiss) -> bg -} +-- set rgb to normal's xyz coordinates+-- as a debugging aid debug_norm_shade :: Rayint -> Ray -> Scene -> Int -> Int -> Color debug_norm_shade ri (Ray o indir) scn recurs debug = case ri of- RayHit d p (Vec nx ny nz) t -> (Color (abs $ nx/2) (abs $ ny/2) (abs $ nz/2))+ RayHit d p (Vec nx ny nz) t -> (Color (fabs $ nx/2) (fabs $ ny/2) (fabs $ nz/2)) RayMiss -> bground scn --- handles diffuse light, shadows, and reflection--- todo: specular highlights, refraction+-- no shadows, reflection, or lighting+flat_shade :: Rayint -> Ray -> Scene -> Int -> Int -> Color+flat_shade ri (Ray o indir) scn recurs debug =+ case ri of+ RayMiss -> bground scn+ RayHit d p n t -> + let (Material clr refl refr ior kd shine) = t ri+ in clr++-- handles diffuse light, shadows, specular highlights and reflection+-- todo: refraction shade :: Rayint -> Ray -> Scene -> Int -> Int -> Color shade ri (Ray o indir) scn recurs !debug = case ri of (RayHit d p n t) ->- let (Material clr refl refr ior kd shine) = t ri+ let (Material clr refl_ refr ior kd shine) = t ri s = sld scn- lits = lights scn+ lights = lits scn direct = foldl' cadd c_black (map (\ (Light lp lc) -> let eyedir = vinvert indir@@ -52,10 +65,9 @@ intensity = 5.0 / (llen*llen) --intensity = 0.2 in- --if blinn /= blinn - --then error $ "nan " ++ (show (vdot halfangle n)) ++ " " ++ - -- (show $ shine*3) ++ " " ++ (show blinn)- --else+ if vdot n lvec < 0 + then c_black+ else if not $ shadow s (Ray (vscaleadd p n delta) ldir) (llen-(2*delta)) then cadd @@ -66,41 +78,56 @@ (cscale lc $ blinn_correct * intensity) -- c_black else - c_black) lits)- reflect = - if (refl > delta) && (recurs > 0)- then let outdir = Vec.reflect indir n + c_black) lights)+ reflect_ = + if (refl_ > delta) && (recurs > 0)+ then let outdir = reflect indir n in cscale (trace scn (Ray (vscaleadd p outdir delta) outdir) - infinity (recurs-1) ) refl+ infinity (recurs-1) ) refl_ else c_black refract = if (refr > delta) && (recurs > 0) then c_black else c_black in- cadd direct $ cadd reflect refract+ cadd direct $ cadd reflect_ refract (RayMiss) -> bground scn - trace :: Scene -> Ray -> Flt -> Int -> Color- trace scn ray depth recurs = let (Scene sld lights cam dtex bgcolor) = scn - in shade (rayint_check sld ray depth dtex) ray scn recurs 0+ in shade (rayint sld ray depth dtex) ray scn recurs 0 +-- return depth as well as color, for post-processing effects+trace_depth :: Scene -> Ray -> Flt -> Int -> (Color,Flt)+trace_depth scn ray depth recurs =+ let (Scene sld lights cam dtex bgcolor) = scn + ri = rayint sld ray depth dtex + d = case ri of+ RayHit d_ _ _ _ -> d_+ RayMiss -> infinity+ clr = shade ri ray scn recurs 0+ in (clr,d) -{- experiments in MVar usage-trace scn ray depth recurs =+-- return hit position as well as color+trace_pos :: Scene -> Ray -> Flt -> Int -> (Color,Vec)+trace_pos scn ray depth recurs = let (Scene sld lights cam dtex bgcolor) = scn - in- unsafePerformIO $- do- debug_start <- (readMVar bihctr)- let ri = rayint sld ray depth dtex- debug_end <- (readMVar bihctr)- let bihhits = debug_end - debug_start- print bihhits- return $ shade ri ray scn recurs bihhits--}+ ri = rayint sld ray depth dtex + p = case ri of+ RayHit _ p _ _ -> p+ RayMiss -> (Vec 0 0 0) -- fixme+ clr = shade ri ray scn recurs 0+ in (clr,p)+++trace_packet :: Scene -> Ray -> Ray -> Ray -> Ray -> Flt -> Int -> PacketColor+trace_packet scn ray1 ray2 ray3 ray4 depth recurs =+ let (Scene sld lights cam dtex bgcolor) = scn+ PacketResult ri1 ri2 ri3 ri4 = packetint sld ray1 ray2 ray3 ray4 depth dtex+ in PacketColor (shade ri1 ray1 scn recurs 0)+ (shade ri2 ray2 scn recurs 0)+ (shade ri3 ray3 scn recurs 0)+ (shade ri4 ray4 scn recurs 0)
+ Triangle.hs view
@@ -0,0 +1,115 @@+module Triangle where+import Vec+import Solid++data Triangle = Triangle Vec Vec Vec deriving Show+data TriangleNorm = TriangleNorm Vec Vec Vec Vec Vec Vec deriving Show++triangle :: Vec -> Vec -> Vec -> SolidItem+triangle v1 v2 v3 =+ SolidItem (Triangle v1 v2 v3)++triangles :: [Vec] -> [SolidItem]+triangles (v1:vs) =+ zipWith (\v2 v3 -> triangle v1 v2 v3) vs (tail vs) ++trianglenorm v1 v2 v3 n1 n2 n3 =+ SolidItem (TriangleNorm v1 v2 v3 n1 n2 n3)++trianglesnorms :: [(Vec,Vec)] -> [SolidItem]+trianglesnorms (vn1:vns) =+ zipWith (\vn2 vn3 -> trianglenorm (fst vn1) (fst vn2) (fst vn3)+ (snd vn1) (snd vn2) (snd vn3))+ vns (tail vns)++-- adaptation of Moller and Trumbore from pbrt page 127+rayint_triangle :: Triangle -> Ray -> Flt -> Texture -> Rayint+rayint_triangle (Triangle p1 p2 p3) (Ray o dir) dist tex =+ let e1 = vsub p2 p1+ e2 = vsub p3 p1+ s1 = vcross dir e2+ divisor = vdot s1 e1+ in + if (divisor == 0)+ then RayMiss+ else+ let invdivisor = 1.0 / divisor+ d = vsub o p1 + b1 = (vdot d s1) * invdivisor+ in+ if (b1 < 0) || (b1 > 1) + then RayMiss + else+ let s2 = vcross d e1+ b2 = (vdot dir s2) * invdivisor+ in+ if (b2 < 0) || (b1 + b2 > 1) + then RayMiss+ else+ let t = (vdot e2 s2) * invdivisor+ in+ if (t < 0) || (t > dist)+ then RayMiss+ else+ RayHit t (vscaleadd o dir t) (vnorm (vcross e1 e2)) tex++rayint_trianglenorm :: TriangleNorm -> Ray -> Flt -> Texture -> Rayint+rayint_trianglenorm (TriangleNorm p1 p2 p3 n1 n2 n3) (Ray o dir) dist tex =+ let e1 = vsub p2 p1+ e2 = vsub p3 p1+ s1 = vcross dir e2+ divisor = vdot s1 e1+ in + if (divisor == 0)+ then RayMiss+ else+ let invdivisor = 1.0 / divisor+ d = vsub o p1 + b1 = (vdot d s1) * invdivisor+ in+ if (b1 < 0) || (b1 > 1) + then RayMiss + else+ let s2 = vcross d e1+ b2 = (vdot dir s2) * invdivisor+ in+ if (b2 < 0) || (b1 + b2 > 1) + then RayMiss+ else+ let t = (vdot e2 s2) * invdivisor+ in+ if (t < 0) || (t > dist)+ then RayMiss+ else+ let n1scaled = (vscale n1 (1-(b1+b2))) + n2scaled = (vscale n2 b1)+ n3scaled = (vscale n3 b2)+ norm = vnorm $ vadd3 n1scaled n2scaled n3scaled+ in RayHit t (vscaleadd o dir t) norm tex++bound_triangle :: Triangle -> Bbox+bound_triangle (Triangle (Vec v1x v1y v1z) + (Vec v2x v2y v2z) + (Vec v3x v3y v3z)) =+ Bbox+ (Vec ((fmin (fmin v1x v2x) v3x) - delta)+ ((fmin (fmin v1y v2y) v3y) - delta)+ ((fmin (fmin v1z v2z) v3z) - delta) )++ (Vec ((fmax (fmax v1x v2x) v3x) + delta)+ ((fmax (fmax v1y v2y) v3y) + delta)+ ((fmax (fmax v1z v2z) v3z) + delta) )++bound_trianglenorm :: TriangleNorm -> Bbox+bound_trianglenorm (TriangleNorm v1 v2 v3 n1 n2 n3) =+ bound (Triangle v1 v2 v3)++instance Solid Triangle where+ rayint = rayint_triangle+ inside _ _ = False+ bound = bound_triangle++instance Solid TriangleNorm where+ rayint = rayint_trianglenorm+ inside _ _ = False+ bound = bound_trianglenorm
Vec.hs view
@@ -82,11 +82,34 @@ then b else c +fmin4 :: Flt -> Flt -> Flt -> Flt -> Flt+fmin4 a b c d = fmin (fmin a b) (fmin c d)++fmax4 :: Flt -> Flt -> Flt -> Flt -> Flt+fmax4 a b c d = fmax (fmax a b) (fmax c d)++fabs :: Flt -> Flt+fabs !a = + if a < 0 then (-a) else a++iabs :: Int -> Int+iabs !a =+ if a < 0 then (-a) else a++abs a = error "use non-polymorphic version, fabs"+ -- true if a and b are "almost" equal+-- the (abs $ a-b) test doesn't work if+-- a and b are large about_equal :: Flt -> Flt -> Bool-about_equal !a !b =- (abs $ a - b) < (delta*10)+about_equal a b =+ if a > 1 + then+ fabs (1 - (a/b)) < (delta*10) + else+ (fabs $ a - b) < (delta*10) + data Vec = Vec {x, y, z :: !Flt} deriving Show data Ray = Ray {origin, dir :: !Vec} deriving Show --data Plane = Plane {norm :: !Vec, offset :: !Flt} deriving Show@@ -169,6 +192,12 @@ (y1 - y2) (z1 - z2) +vmul :: Vec -> Vec -> Vec+vmul !(Vec !x1 !y1 !z1) !(Vec !x2 !y2 !z2) =+ Vec (x1 * x2)+ (y1 * y2)+ (z1 * z2)+ vinc :: Vec -> Flt -> Vec vinc v1 n = Vec ((x v1) + n)@@ -217,16 +246,16 @@ vnorm :: Vec -> Vec vnorm (Vec x1 y1 z1) = - let len = 1.0 / (sqrt ((x1*x1)+(y1*y1)+(z1*z1))) in- Vec (x1*len) (y1*len) (z1*len)+ let invlen = 1.0 / (sqrt ((x1*x1)+(y1*y1)+(z1*z1))) in+ Vec (x1*invlen) (y1*invlen) (z1*invlen) assert_norm :: Vec -> Vec assert_norm v = let l = vdot v v in if l > (1+delta) - then error "vector too long"+ then error $ "vector too long" ++ (show v) else if l < (1-delta)- then error "vector too short"+ then error $ "vector too short: " ++ (show v) else v bisect :: Vec -> Vec -> Vec@@ -250,6 +279,11 @@ veq (Vec ax ay az) (Vec bx by bz) = (about_equal ax bx) && (about_equal ay by) && (about_equal az bz) +--returns false on zero value+veqsign :: Vec -> Vec -> Bool+veqsign (Vec ax ay az) (Vec bx by bz) =+ ax*bx > 0 && ay*by > 0 && az*bz > 0+ -- translate a ray's origin in ray's direction by d amount ray_move :: Ray -> Flt -> Ray ray_move (Ray orig dir) d =@@ -325,13 +359,13 @@ xfm_mult (Xfm a inva) (Xfm b invb) = Xfm (mat_mult a b) (mat_mult invb inva) --- UTILITY FUNCTIONS --+-- TRANSFORM UTILITY FUNCTIONS -- -- If we multiply two transformation matricies, we get -- a transformation matrix equivalent to applying the --- second then the first to a vector.+-- second then the first. --- By reversing the list, we apply transforms in expected order.+-- By reversing the list, the transforms are applied in the expected order. compose :: [Xfm] -> Xfm compose xfms = check_xfm $ foldr xfm_mult ident_xfm (reverse xfms) @@ -346,8 +380,27 @@ ae m10 0 && ae m11 1 && ae m12 0 && ae m13 0 && ae m20 0 && ae m21 0 && ae m22 1 && ae m23 0 then (Xfm m i)- else error $ "corrupt matrix " ++ (show (Xfm m i)) + else error $ "corrupt matrix " ++ (show (Xfm m i)) ++ "\n" ++ (show (mat_mult m i)) +-- rotate point (or vector) a about ray b by angle c+vrotate :: Vec -> Ray -> Flt -> Vec+vrotate pt (Ray orig axis_) angle =+ let axis = assert_norm axis_+ transform = compose [ translate (vinvert orig)+ , rotate axis angle+ , translate orig+ ]+ new_pt = xfm_point transform pt+ in if about_equal (vlen (vsub orig pt)) (vlen (vsub orig new_pt))+ then new_pt+ else error $ "something is wrong with vrotate" ++ + (show $ vlen (vsub orig pt)) ++ " " ++ + (show $ vlen (vsub orig new_pt))+++-- TRANSFORM APPLICATION --+-- these need to be fast+ -- point is treated as (x y z 1) xfm_point :: Xfm -> Vec -> Vec xfm_point (Xfm (Matrix m00 m01 m02 m03 @@ -405,35 +458,41 @@ Ray (invxfm_point xfm orig) (vnorm (invxfm_vec xfm dir)) -- BASIC TRANSFORMS --+-- move translate :: Vec -> Xfm translate (Vec x y z) = check_xfm $ Xfm (Matrix 1 0 0 x 0 1 0 y 0 0 1 z) (Matrix 1 0 0 (-x) 0 1 0 (-y) 0 0 1 (-z)) +-- strectch along three axes (if x==y==z, then it's uniform scaling) scale :: Vec -> Xfm scale (Vec x y z) = check_xfm $ Xfm (Matrix x 0 0 0 0 y 0 0 0 0 z 0) (Matrix (1/x) 0 0 0 0 (1/y) 0 0 0 0 (1/z) 0) +-- rotate about an arbitrary axis and angle rotate :: Vec -> Flt -> Xfm rotate (Vec x y z) angle =- let s = sin angle- c = cos angle + if not $ (vlen (Vec x y z)) `about_equal` 1+ then error $ "please use a normalized vector for rotation: " ++ (show (vlen (Vec x y z)))+ else + let s = sin angle+ c = cos angle - m00 = ((x*x)+((1-(x*x))*c)) - m01 = (((x*y)*(1-c))-(z*s)) - m02 = ((x*z*(1-c))+(y*s))+ m00 = ((x*x)+((1-(x*x))*c)) + m01 = (((x*y)*(1-c))-(z*s)) + m02 = ((x*z*(1-c))+(y*s)) - m10 = (((x*y)*(1-c))+(z*s))- m11 = ((y*y)+((1-(y*y))*c))- m12 = ((y*z*(1-c))-(x*s))+ m10 = (((x*y)*(1-c))+(z*s))+ m11 = ((y*y)+((1-(y*y))*c))+ m12 = ((y*z*(1-c))-(x*s)) - m20 = ((x*z*(1-c))-(y*s))- m21 = ((y*z*(1-c))+(x*s))- m22 = ((z*z)+((1-(z*z))*c))- in- check_xfm $ Xfm (Matrix m00 m01 m02 0 m10 m11 m12 0 m20 m21 m22 0)- (Matrix m00 m10 m20 0 m01 m11 m21 0 m02 m12 m22 0)+ m20 = ((x*z*(1-c))-(y*s))+ m21 = ((y*z*(1-c))+(x*s))+ m22 = ((z*z)+((1-(z*z))*c))+ in+ check_xfm $ Xfm (Matrix m00 m01 m02 0 m10 m11 m12 0 m20 m21 m22 0)+ (Matrix m00 m10 m20 0 m01 m11 m21 0 m02 m12 m22 0) -- convert canonical coordinates to uvw coordinates xyz_to_uvw :: Vec -> Vec -> Vec -> Xfm
glome-hs.cabal view
@@ -1,5 +1,5 @@ Name: glome-hs-Version: 0.4.1+Version: 0.5 Synopsis: ray tracer Description: Ray Tracer capable of rendering a variety of primitives, with support for CSG (difference and intersection of solids),@@ -13,10 +13,10 @@ Homepage: http://syn.cs.pdx.edu/~jsnow/glome Stability: experimental Category: graphics-Build-Depends: base,haskell98,time,parallel,GLUT,OpenGL,random,array+Build-Depends: base,haskell98,time,parallel,GLUT,OpenGL,random,array,binary build-type: Simple Executable: glome-ghc-options: -fglasgow-exts -funbox-strict-fields+ghc-options: -fglasgow-exts -funbox-strict-fields -threaded extensions: BangPatterns Main-is: Glome.hs
make view
@@ -3,7 +3,7 @@ #ghc -O3 --make Glome.hs #ghc Glome.hs --make -O2 -threaded -fasm -optc-march=athlon64 -XFlexibleInstances -XTypeSynonymInstances #ghc Glome.hs --make -O2 -fasm -fglasgow-exts -funbox-strict-fields -fbang-patterns -fexcess-precision -optc-ffast-math -optc-O2 -optc-mfpmath=sse -optc-msse2-#ghc Glome.hs --make -fasm -fglasgow-exts -funbox-strict-fields -fbang-patterns -fexcess-precision -prof -auto+#ghc Glome.hs --make -fasm -fglasgow-exts -funbox-strict-fields -fbang-patterns -fexcess-precision -prof -auto-all runhaskell Setup.lhs configure --prefix=$HOME --user runhaskell Setup.lhs build
run view
@@ -1,4 +1,4 @@ #!/bin/bash #./Glome +RTS -N2 -sstderr -RTS #./Glome +RTS-./dist/build/glome/glome+./dist/build/glome/glome +RTS -N2