diff --git a/Bih.hs b/Bih.hs new file mode 100644 --- /dev/null +++ b/Bih.hs @@ -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 diff --git a/Bound.hs b/Bound.hs new file mode 100644 --- /dev/null +++ b/Bound.hs @@ -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 diff --git a/Box.hs b/Box.hs new file mode 100644 --- /dev/null +++ b/Box.hs @@ -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 diff --git a/Cone.hs b/Cone.hs new file mode 100644 --- /dev/null +++ b/Cone.hs @@ -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 diff --git a/Csg.hs b/Csg.hs new file mode 100644 --- /dev/null +++ b/Csg.hs @@ -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 diff --git a/Glome.hs b/Glome.hs --- a/Glome.hs +++ b/Glome.hs @@ -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" ] diff --git a/Plane.hs b/Plane.hs new file mode 100644 --- /dev/null +++ b/Plane.hs @@ -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 diff --git a/README b/README --- a/README +++ b/README @@ -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". diff --git a/Scene.hs b/Scene.hs new file mode 100644 --- /dev/null +++ b/Scene.hs @@ -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) +-} diff --git a/Solid.hs b/Solid.hs --- a/Solid.hs +++ b/Solid.hs @@ -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 diff --git a/SolidTexture.hs b/SolidTexture.hs --- a/SolidTexture.hs +++ b/SolidTexture.hs @@ -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 diff --git a/Spd.hs b/Spd.hs --- a/Spd.hs +++ b/Spd.hs @@ -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 diff --git a/Sphere.hs b/Sphere.hs new file mode 100644 --- /dev/null +++ b/Sphere.hs @@ -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 diff --git a/TestScene.hs b/TestScene.hs --- a/TestScene.hs +++ b/TestScene.hs @@ -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) diff --git a/Tex.hs b/Tex.hs new file mode 100644 --- /dev/null +++ b/Tex.hs @@ -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 diff --git a/Trace.hs b/Trace.hs --- a/Trace.hs +++ b/Trace.hs @@ -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) diff --git a/Triangle.hs b/Triangle.hs new file mode 100644 --- /dev/null +++ b/Triangle.hs @@ -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 diff --git a/Vec.hs b/Vec.hs --- a/Vec.hs +++ b/Vec.hs @@ -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 diff --git a/glome-hs.cabal b/glome-hs.cabal --- a/glome-hs.cabal +++ b/glome-hs.cabal @@ -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 diff --git a/make b/make --- a/make +++ b/make @@ -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 diff --git a/run b/run --- a/run +++ b/run @@ -1,4 +1,4 @@ #!/bin/bash #./Glome +RTS -N2 -sstderr -RTS #./Glome +RTS -./dist/build/glome/glome +./dist/build/glome/glome +RTS -N2