GlomeTrace (empty) → 0.1.1
raw patch · 18 files changed
+2429/−0 lines, 18 filesdep +GlomeVecdep +arraydep +basesetup-changed
Dependencies added: GlomeVec, array, base
Files
- Data/Glome/Bih.hs +402/−0
- Data/Glome/Bound.hs +79/−0
- Data/Glome/Box.hs +65/−0
- Data/Glome/Clr.hs +41/−0
- Data/Glome/Cone.hs +259/−0
- Data/Glome/Csg.hs +116/−0
- Data/Glome/Plane.hs +46/−0
- Data/Glome/Scene.hs +90/−0
- Data/Glome/Solid.hs +552/−0
- Data/Glome/Spd.hs +250/−0
- Data/Glome/Sphere.hs +88/−0
- Data/Glome/Tex.hs +57/−0
- Data/Glome/Trace.hs +174/−0
- Data/Glome/Triangle.hs +141/−0
- GlomeTrace.cabal +34/−0
- LICENSE +14/−0
- README.txt +15/−0
- Setup.hs +6/−0
+ Data/Glome/Bih.hs view
@@ -0,0 +1,402 @@+{-+Copyright (c) 2008 Jim Snow+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.+2. Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.+3. The name of the author may not be used to endorse or promote products+ derived from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR+IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.+IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT+NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF+THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.++-}++{-# OPTIONS_GHC -funbox-strict-fields #-}+{-# LANGUAGE BangPatterns #-}++module Data.Glome.Bih (bih) where+import Data.Glome.Vec+import Data.Glome.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+-- create a leaf node from a list of objects+-- we use "group" so we can treat a bunch of objects as a single object+build_leaf :: [(Bbox, SolidItem)] -> BihNode+build_leaf objs =+ BihLeaf (group (map snd objs))++-- return surface area of a bounding box that encloses bounding boxes+-- divided by the surface area of the nodebox++-- this doesn't seem to be much of a win++optimality :: [(Bbox, SolidItem)] -> Bbox -> Flt+optimality objs bb =+ let bbsurf = bbsa bb+ go [] accbb = accbb+ go ((obb,_):xs) accbb = go xs (bbjoin obb accbb)+ obbsurf = bbsa $! bboverlap (go objs empty_bbox) bb+ in+ obbsurf / bbsurf++-- tuning parameter that controls threshold for separating+-- large objects from small objects instead of usual left/right+-- sorting ++-- was 0.3++max_bih_sa = 0.4 :: Flt++-- Recursive constructor, it looks like quicksort if you squint hard enough.+-- We split along the splitbox's axis of greatest extent, then sort objects+-- to one side or the other (they can overlap the center), then construct the+-- branch node and recurse.++-- I added a nonstandard heuristic: if there's a few very large objects and a lot+-- of small ones, we create one branch with big objects and the other with small+-- objects, instead of sorting by location.++build_rec :: [(Bbox,SolidItem)] -> Bbox -> Bbox -> Int -> BihNode+build_rec objs nodebox@(Bbox nodeboxp1 nodeboxp2) splitbox@(Bbox splitboxp1 splitboxp2) depth = ++ if (length (take 3 objs) < 2) -- && (optimality objs nodebox) > 0.2+ then build_leaf objs+ else+ let 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)) )++-- | The bih constructor creates a Bounding Interval Heirarchy+-- from a list of primitives. BIH is a type of data structure+-- that groups primitives into a heirarchy of bounding objects,+-- so that a ray need not be tested against every single+-- primitive. This can make the difference betweeen a rendering+-- job that takes days or seconds. BIH usually performs a little+-- worse than a SAH-based KD-tree, but construction time is much+-- better.+--+-- See http://en.wikipedia.org/wiki/Bounding_interval_hierarchy++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)++-- Standard ray traversal.+rayint_bih :: Bih -> Ray -> Flt -> Texture -> Rayint +rayint_bih (Bih bb root) !r@(Ray orig dir) !d t =+ let 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++-- Ray traversal with debug counter. The counter gets incremented+-- when we hit a box.+rayint_debug_bih :: Bih -> Ray -> Flt -> Texture -> (Rayint,Int) +rayint_debug_bih (Bih bb root) r@(Ray orig dir) d t =+ let dir_rcp = vrcp dir+ Interval near far = bbclip r bb+ traverse (BihLeaf s) near far = rayint_debug 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 + debug_wrap + (if near > far + then (RayMiss,0)+ else+ if dirr > 0+ then + (nearest_debug+ (if near < dl+ then traverse l near (fmin dl far)+ else (RayMiss,0))+ (if dr < far+ then traverse r (fmax dr near) far+ else (RayMiss,0)))+ else+ (nearest_debug+ (if near < dr+ then traverse r near (fmin dr far)+ else (RayMiss,0))+ (if dl < far+ then traverse l (fmax dl near) far+ else (RayMiss,0))))+ 1 -- increment the debug value for every box we hit+ in+ traverse root near far++-- This is unwieldy, but the performance gains+-- sometimes make it worthwhile. By testing 4 rays against +-- each cell, we (theoretically) do ~1/4 the +-- memory accesses. ++-- This originally made a big difference, but after switching+-- everything to typeclasses, it doesn't perform any better+-- than regular traversal.++-- 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@(Bih bb root) + !r1@(Ray orig1 dir1) + !r2@(Ray orig2 dir2) + !r3@(Ray orig3 dir3) + !r4@(Ray orig4 dir4) !d t =+ let 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 = fmin4 near1 near2 near3 near4+ far = fmax4 far1 far2 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@(Ray orig dir) d =+ let 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/outside test; essentially a point traversal.+-- We test if the point is inside any of the objects contained in+-- the bih.++inside_bih :: Bih -> Vec -> Bool+inside_bih (Bih (Bbox (Vec x1 y1 z1) (Vec x2 y2 z2)) root) pt@(Vec x y z) =+ let 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)++-- We already have a bounding box computed.+bound_bih :: Bih -> Bbox+bound_bih (Bih bb root) = bb++primcount_bih :: Bih -> Pcount+primcount_bih (Bih bb root) = pcadd (bihcount root) pcsinglebound+ where bihcount (BihLeaf s) = primcount s+ bihcount (BihBranch _ _ _ l r) = + pcadd (pcadd (bihcount l) (bihcount r)) pcsinglebound++instance Solid Bih where+ rayint = rayint_bih+ rayint_debug = rayint_debug_bih+ packetint = packetint_bih+ shadow = shadow_bih+ inside = inside_bih+ bound = bound_bih+ primcount = primcount_bih
+ Data/Glome/Bound.hs view
@@ -0,0 +1,79 @@+module Data.Glome.Bound (bound_object) where+import Data.Glome.Vec+import Data.Glome.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++-- | Use the first object as a bounding volume for the second+-- object. If a ray misses the first object, it is assumed to+-- miss the second object. Used primarily to improve performance.+-- In general, bih will usually perform better than +-- manually-constructed bounds, though.++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++rayint_debug_bound :: Bound -> Ray -> Flt -> Texture -> (Rayint,Int)+rayint_debug_bound (Bound sa sb) r d t =+ let (Ray orig _) = r+ in if inside sa orig || shadow sa r d+ then (debug_wrap (rayint_debug sb r d t) 1)+ else (RayMiss,0)++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)++transform_leaf_bound :: Bound -> [Xfm] -> SolidItem+transform_leaf_bound (Bound sa sb) xfms =+ transform_leaf sb xfms++flatten_transform_bound :: Bound -> [SolidItem]+flatten_transform_bound (Bound sa sb) = flatten_transform sb++primcount_bound :: Bound -> Pcount+primcount_bound (Bound sa sb) = pcadd (asbound (primcount sa)) (primcount sb)++instance Solid Bound where+ rayint = rayint_bound+ rayint_debug = rayint_debug_bound+ shadow = shadow_bound+ inside = inside_bound+ bound = bound_bound+ flatten_transform = flatten_transform_bound+ transform_leaf = transform_leaf_bound+ primcount = primcount_bound
+ Data/Glome/Box.hs view
@@ -0,0 +1,65 @@+module Data.Glome.Box (box) where+import Data.Glome.Vec+import Data.Glome.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@(Ray orig@(Vec ox oy oz) dir@(Vec dx dy dz)) d t =+ 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)+ 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
+ Data/Glome/Clr.hs view
@@ -0,0 +1,41 @@+module Data.Glome.Clr where++type CFlt = Double+data Color = Color {r,g,b :: !CFlt} deriving Show++c_black = Color 0 0 0+c_white = Color 1 1 1+c_red = Color 1 0 0+c_green = Color 0 1 0+c_blue = Color 0 0 1++cadd :: Color -> Color -> Color+cadd (Color r1 g1 b1) (Color r2 g2 b2) =+ Color (r1+r2) (g1+g2) (b1+b2)++cdiv :: Color -> CFlt -> Color+cdiv c1 div =+ cscale c1 (1/div)++cscale :: Color -> CFlt -> Color+cscale (Color r g b) mul =+ Color (r * mul)+ (g * mul)+ (b * mul)++cmul :: Color -> Color -> Color+cmul (Color r1 g1 b1) (Color r2 g2 b2) =+ Color (r1*r2) (g1*g2) (b1*b2)++cavg :: Color -> Color -> Color+cavg c1 c2 = cscale (cadd c1 c2) 0.5++cscaleadd :: Color -> Color -> CFlt -> Color+cscaleadd (Color r1 g1 b1) (Color r2 g2 b2) mul =+ Color (r1+(r2*mul)) (g1+(g2*mul)) (b1+(b2*mul))++cclamp :: Color -> Color+cclamp (Color r g b) = + Color (if r > 0.0 then r else 0.0)+ (if g > 0.0 then g else 0.0)+ (if b > 0.0 then b else 0.0)
+ Data/Glome/Cone.hs view
@@ -0,0 +1,259 @@+{-# OPTIONS_GHC -funbox-strict-fields #-}+{-# LANGUAGE BangPatterns #-}++module Data.Glome.Cone (disc, cone, cylinder) where+import Data.Glome.Vec+import Data.Glome.Solid+import Data.Glome.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 --++-- | Create a disc. These are used as the end-caps on cones and cylinders,+-- but they can be constructed by themselves as well.+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) ]+ +-- | Construct a cone from p1 to p2. R1 and r2 are the radii at each+-- end. A cone need not come to a point at either end.+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@(Ray orig dir) d t =+ let 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@(Ray orig dir) !d =+ let 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@(Ray orig@(Vec ox oy oz) dir@(Vec dx dy dz)) d t =+ let 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@(Vec posx posy posz) = vscaleadd orig dir dist+ 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 d t+ --then rayint_aadisc h1 r ray d t+ else RayMiss+ else if oz > h2+ then rayint_disc (Disc (Vec 0 0 h2) vz (r*r)) ray d t+ --rayint_aadisc h2 r ray 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@(Ray orig@(Vec ox oy oz) dir@(Vec dx dy dz)) d t =+ let 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 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 d t+ --rayint_aadisc clip2 r2 ray d t+ else RayMiss+ -- then rayint_aadisc clip1 r ray 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@(Ray orig@(Vec ox oy oz) dir@(Vec dx dy dz)) d =+ let 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 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 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
+ Data/Glome/Csg.hs view
@@ -0,0 +1,116 @@+module Data.Glome.Csg (difference, intersection) where+import Data.Glome.Vec+import Data.Glome.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--+-- | Create a new object based on the subtraction of the second item+-- from the first. This only works if the items have a well-defined+-- inside and outside. Triangles and discs, for instance, have no +-- volume, so subtracting them from anything won't do anything.+difference :: SolidItem -> SolidItem -> SolidItem+difference a b = SolidItem $ Difference a b++rayint_difference :: Difference -> Ray -> Flt -> Texture -> Rayint+rayint_difference dif@(Difference sa sb) r@(Ray orig dir) d t =+ let 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--++-- | Create a new item from the boolean intersection of a+-- list of solids. A point is inside the object iff it is+-- inside every primitive. We can construct polyhedra from+-- intersections of planes, but this isn't the most efficient+-- way to do that.+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@(Ray orig dir) d t =+ 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)++primcount_difference :: Difference -> Pcount+primcount_difference (Difference sa sb) = pcadd (primcount sa) (primcount sb)++primcount_intersection :: Intersection -> Pcount+primcount_intersection (Intersection slds) = foldl (pcadd) pcnone (map primcount slds)++instance Solid Difference where+ rayint = rayint_difference+ inside = inside_difference+ bound = bound_difference+ primcount = primcount_difference++instance Solid Intersection where+ rayint = rayint_intersection+ inside = inside_intersection+ bound = bound_intersection+ primcount = primcount_intersection
+ Data/Glome/Plane.hs view
@@ -0,0 +1,46 @@+module Data.Glome.Plane (plane, plane_offset) where+import Data.Glome.Vec+import Data.Glome.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++-- | Construct a plane (or, more accurately, a half-space)+-- by specifying a point on the plane and a normal.+-- The normal points towards the outside of the plane.+-- Planes are often useful within CSG objects.+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
+ Data/Glome/Scene.hs view
@@ -0,0 +1,90 @@+module Data.Glome.Scene (+ Scene(Scene), Light(Light), Camera(Camera),+ scene, camera, light, + sld, lits, cam, dtex, bground,+ primcount_scene,+ module Data.Glome.Clr,+ module Data.Glome.Vec,+ module Data.Glome.Solid,+ module Data.Glome.Sphere,+ module Data.Glome.Triangle,+ module Data.Glome.Bih,+ module Data.Glome.Csg,+ module Data.Glome.Plane,+ module Data.Glome.Box,+ module Data.Glome.Bound,+ module Data.Glome.Cone,+ module Data.Glome.Tex) where+import Data.Glome.Clr+import Data.Glome.Vec+import Data.Glome.Solid+import Data.Glome.Sphere+import Data.Glome.Triangle+import Data.Glome.Bih+import Data.Glome.Csg+import Data.Glome.Plane+import Data.Glome.Box+import Data.Glome.Bound+import Data.Glome.Cone+import Data.Glome.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+++-- | Construct a light given a center location and a color.+light :: Vec -> Color -> Light+light pos clr = Light pos clr++-- CAMERA --+data Camera = Camera {campos, fwd, up, right :: !Vec} + deriving Show++-- | Construct a camera pointing in some default direction.+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) )++-- | Construct a camera, given a position, a forward vector, +-- a point that the camera should be pointed towards, an up vector,+-- and a right vector. The up and right vectors don't have to be+-- normalized or perfectly orthogonal.+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++-- | Create a scene from an item (which can be a composite item, such +-- as a bih or group), a list of lights, a camera, a default texture,+-- and a default background color.+scene :: SolidItem -> [Light] -> Camera -> Texture -> Color -> Scene+scene s l cam t clr = Scene s l cam t clr++-- | Count the primitives in the scene. See docs for primcount +-- in Solid.hs.+primcount_scene :: Scene -> Pcount+primcount_scene (Scene sld _ _ _ _) = primcount sld++{-+default_scene = (Scene (sphere (vec 0.0 0.0 0.0) 1.0) + [] default_cam t_white c_white)+-}
+ Data/Glome/Solid.hs view
@@ -0,0 +1,552 @@+{-# OPTIONS_GHC -fexcess-precision #-}+{-# OPTIONS_GHC -funbox-strict-fields #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE ExistentialQuantification #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}++module Data.Glome.Solid where+import Data.Glome.Vec+import Data.Glome.Clr+import Data.List hiding (group)++-- | Ray intersection type. If we hit, we store the distance from the ray+-- origin, the position, the normal, and the texture attached to the object.+-- We could just as easily have created a hit type and wrapped it in a Maybe.++data Rayint = RayHit {+ depth :: !Flt,+ pos :: !Vec,+ norm :: !Vec,+ texture :: Texture+} | RayMiss deriving Show++-- | Pick the closest of two Rayints+nearest :: Rayint -> Rayint -> Rayint+nearest a RayMiss = a+nearest RayMiss b = b+nearest !rha@(RayHit da _ _ _) !rhb@(RayHit db _ _ _) =+ if da < db+ then rha+ else rhb++-- | Pick the furthest of two Rayints+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++-- | Test if a Rayint is a hit or a miss+hit :: Rayint -> Bool+hit (RayHit _ _ _ _) = True+hit RayMiss = False++-- | Extract a distance from a Rayint, with infinity for a miss+dist :: Rayint -> Flt+dist RayMiss = infinity+dist (RayHit d _ _ _) = d++--Packet Types--++-- | Sometimes, it's more efficient to trace multiple rays against an +-- acceleration structure at the same time, provided the rays are almost+-- identical. A PacketResult is the result of tracing 4 rays at once.++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 a ray forward and test the new ray against an object.+-- Fix the depth of the 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--++-- | Surface properties at a point on an object's surface. We have color, +-- reflection amount, refraction amount index of refraction, kd, ks, and shine.+-- These are parameters to a Whitted - style illumination model.++data Material = Material {clr :: !Color, + refl, refr, ior, + kd, ks, shine :: !Flt} deriving Show++-- | A texture is a function that takes a Rayint and returns a Material.+-- In other words, textures can vary based on location, normal, etc...+-- in arbitrary ways.+type Texture = Rayint -> Material++-- | This is sort of a no-op; textures are functions, and 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++-- | Uniform white material+m_white = (Material c_white 0 0 0 1 0 2)+t_white ri = m_white++-- | Uniform texture+t_uniform :: Material -> Texture+t_uniform m = \x -> m++interp :: Flt -> Flt -> Flt -> Flt+interp scale a b =+ scale*a + (1-scale)*b++-- | Interpolate between textures. +-- 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 m1ks m1shine) = m1+ (Material m2c m2refl m2refr m2ior m2kd m2ks 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+ ks = intp m1ks m2ks+ shine = intp m1shine m2shine+ in (Material c refl refr ior kd ks shine)++--utility functions for "primcount"+newtype Pcount = Pcount (Int,Int,Int) deriving Show++pcadd :: Pcount -> Pcount -> Pcount+pcadd (Pcount (a1,a2,a3)) (Pcount (b1,b2,b3)) = Pcount (a1+b1, a2+b2, a3+b3)++asbound :: Pcount -> Pcount+asbound (Pcount (a,b,c)) = Pcount (0,b,a+c)++pcsinglexfm :: Pcount+pcsinglexfm = Pcount (0,1,0)++pcsingleprim :: Pcount+pcsingleprim = Pcount (1,0,0)++pcsinglebound :: Pcount+pcsinglebound = Pcount (0,0,1)++pcnone :: Pcount+pcnone = Pcount (0,0,0)++-- utility functions for rayint_debug+debug_wrap :: (Rayint,Int) -> Int -> (Rayint,Int)+debug_wrap (ri,a) b = (ri,(a+b))++nearest_debug :: (Rayint,Int) -> (Rayint,Int) -> (Rayint,Int)+nearest_debug (ari, an) (bri, bn) = ((nearest ari bri),(an+bn))++--SOLID CLASS--++-- | A solid is something we can test a ray against or do inside/outside tests.+-- Some of these are simple solids like Sphere or Triangle, but others+-- are composite solids than have other solids as children.++class (Show a) => Solid a where++ -- | Test a ray against a solid, returning a ray intersection.+ -- The distance parameter is used to specify a max distance.+ -- If it's further away, we aren't interested in the intersection.+ -- The texture parameter is a default texture we use, if it's not+ -- overridden by a more specific texture.+ rayint :: a -- ^ object to test against+ -> Ray -- ^ ray+ -> Flt -- ^ maximum distance we care about+ -> Texture -- ^ default texture+ -> Rayint -- ^ we return a Rayint describing the hit location++ -- | Same as rayint, but return a count of the number of+ -- primitives checked. Useful for optimizing acceleration structures.+ rayint_debug :: a -> Ray -> Flt -> Texture -> (Rayint, Int)++ -- | Trace four rays at once against a solid.+ packetint :: a -> Ray -> Ray -> Ray -> Ray -> Flt -> Texture -> PacketResult ++ -- | Shadow test - we just return a Bool rather than return a + -- a full Rayint.+ shadow :: a -> Ray -> Flt -> Bool++ -- | Test if a point is inside an object. Useful for CSG.+ -- Objects with no volume just return False.+ inside :: a -> Vec -> Bool++ -- | Generate an axis-aligned bounding box than completely encloses+ -- the object. For performance, it is important that this fits as + -- tight as possible.+ bound :: a -> Bbox++ -- | Most simple objects just return themselves as a singleton list,+ -- but for composite objects, we flatten the structure out and + -- return a list. We usually do this prior to re-building a + -- composite object in a (hopefully) more efficient fashion.+ tolist :: a -> [SolidItem]++ -- | Create a new object transformed by some transformation. The+ -- reason this method exists is so we can override it for the+ -- Instance type - if we transform a transformation, we should+ -- combine the two matricies into one.+ -- Most objects can use the default implementation.+ transform :: a -> [Xfm] -> SolidItem++ -- | Used by flatten_transform. I don't really remember how it works. + transform_leaf :: a -> [Xfm] -> SolidItem++ -- | Take a composite object inside a transform, and turn it into+ -- a group of individually-transformed objects. Most objects + -- can use the defaut implementation.+ flatten_transform :: a -> [SolidItem]++ -- | Count the number of primitives, transforms, and bounding+ -- objects in a scene. Simple objects can just use the default,+ -- which is to return a single primitive.+ primcount :: a -> Pcount++ -- | This is for counting bih split planes and the like, for+ -- performance tuning and debugging. Most objects can use+ -- the default implementation.+ rayint_debug s !r !d t = ((rayint s r d t),0)++ -- | Sometimes, we can improve performance by + -- intersecting 4 rays at once. This is + -- especially true of acceleration structures.+ -- The default implementation is to 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++ -- There's a name for what a bunch of these functions+ -- are trying to do (but poorly): what I really want is+ -- a "catamorphism".++ -- 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 is used by flatten_transform below. For simple objects, it + -- works the same as transform, but for groups it transforms all the+ -- objects individually.+ transform_leaf = transform++ -- 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. For simple primitives, this+ -- is a no-op.+ flatten_transform = tolist++ -- Figure out how complicated the scene really is.+ -- Returns (primitives, matricies, bounding objects/planes).+ -- Also, it forces the full construction of acceleration structures.+ primcount s = pcsingleprim+++-- | We create an existential type for solids so we can emded them+-- in composite types without know what kind of solid it is.+-- 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+ packetint (SolidItem s) !r1 !r2 !r3 !r4 !d t = packetint s r1 r2 r3 r4 d t+ rayint_debug (SolidItem s) r d t = rayint_debug 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 (SolidItem s) = tolist s -- don't wrap in a redundant SolidItem like everything else+ transform (SolidItem s) xfm = transform s xfm -- same here+ transform_leaf (SolidItem s) xfm = transform_leaf s xfm -- and here+ flatten_transform (SolidItem s) = [SolidItem (flatten_transform s)] -- and here+ primcount (SolidItem s) = primcount s++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 --+--+-- | A group is just a list of objects. Sometimes its convenient to be +-- able to treat a group as if it were a single object, and that is +-- exactly what we do here. The ray intersection routine tests the ray +-- against each object in turn. Not very efficient+-- for large groups, but this is a useful building block for+-- constructing the leaves of acceleration structures. (See the bih+-- module.)++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)++{-- this is not measurably faster+rayint_group slds r d t = go slds RayMiss+ where go [] res = res+ go (x:xs) res = go xs $ nearest (rayint x r d t) res+--}++packetint_group :: [SolidItem] -> Ray -> Ray -> Ray -> Ray -> Flt -> Texture -> PacketResult+packetint_group [] !r1 !r2 !r3 !r4 !d t = packetmiss+packetint_group (x:xs) !r1 !r2 !r3 !r4 !d t = + nearest_packetresult (packetint x r1 r2 r3 r4 d t) + (packetint_group xs r1 r2 r3 r4 d t)++rayint_debug_group :: [SolidItem] -> Ray -> Flt -> Texture -> (Rayint,Int)+rayint_debug_group [] _ _ _ = (RayMiss,0)+rayint_debug_group (x:xs) !r !d t = + nearest_debug (rayint_debug x r d t) + (rayint_debug_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)++transform_leaf_group :: [SolidItem] -> [Xfm] -> SolidItem+transform_leaf_group slds xfms =+ SolidItem $ map (\x -> transform_leaf x xfms) (tolist slds)++primcount_group :: [SolidItem] -> Pcount+primcount_group slds = foldl (pcadd) (Pcount (0,0,0)) (map primcount slds)++instance Solid [SolidItem] where+ rayint = rayint_group+ packetint = packetint_group+ rayint_debug = rayint_debug_group+ shadow = shadow_group+ inside = inside_group+ bound = bound_group+ tolist a = concat $ map tolist a+ transform_leaf = transform_leaf_group+ flatten_transform a = concat $ map flatten_transform a+ primcount = primcount_group++-- VOID --++-- | A Void is a non-object, that we treat as if it were+-- one. This is functionally equivalent to an empty Group.+-- (Originally I called this "Nothing", but that+-- conflicted with the prelude maybe type, so I call+-- it "Void" instead) +data Void = Void deriving Show++nothing = SolidItem Void++instance Solid Void where+ rayint Void _ _ _ = RayMiss+ packetint Void _ _ _ _ _ _ = packetmiss+ shadow Void _ _ = False+ inside Void _ = False+ bound Void = empty_bbox+ tolist Void = []+ transform Void xfms = SolidItem Void ++-- INSTANCE --+--+-- | 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.+--+-- Usually, the application doesn't need to create an +-- instance directly, but should use "transform" on an+-- existing object.+-- +-- It's unfortunate that "instance" is also a reserved word. +-- "instance Solid Instance where..." is a little confusing.+-- +-- This would be better in its own module, but we need+-- "Instance" to be defined here so we can define the default +-- implementation of "transform" in terms on Instance.+-- (Mutually recursive modules would be useful, if I could+-- get them to work.)+--+-- Another good reason to include Instance in Solid.hs+-- is that it's referenced from Cone.hs++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++packetint_instance :: Instance -> Ray -> Ray -> Ray -> Ray -> Flt -> Texture -> PacketResult+packetint_instance !(Instance sld xfm) !(Ray orig1 dir1) !(Ray orig2 dir2) + !(Ray orig3 dir3) !(Ray orig4 dir4) d t =+ let newdir1 = invxfm_vec xfm dir1+ newdir2 = invxfm_vec xfm dir2+ newdir3 = invxfm_vec xfm dir3+ newdir4 = invxfm_vec xfm dir4+ neworig1 = invxfm_point xfm orig1+ neworig2 = invxfm_point xfm orig2+ neworig3 = invxfm_point xfm orig3+ neworig4 = invxfm_point xfm orig4+ lenscale1 = vlen newdir1+ lenscale2 = vlen newdir2+ lenscale3 = vlen newdir3+ lenscale4 = vlen newdir4+ invlenscale1 = 1/lenscale1+ invlenscale2 = 1/lenscale2+ invlenscale3 = 1/lenscale3+ invlenscale4 = 1/lenscale4+ in+ let pr = packetint sld (Ray neworig1 (vscale newdir1 invlenscale1)) + (Ray neworig2 (vscale newdir2 invlenscale2)) + (Ray neworig3 (vscale newdir3 invlenscale3)) + (Ray neworig4 (vscale newdir4 invlenscale4)) + (d*lenscale1) t+ PacketResult ri1 ri2 ri3 ri4 = pr + fix ri ils = + case ri of + RayMiss -> RayMiss+ RayHit depth pos n tex -> RayHit (depth*ils) + (xfm_point xfm pos) + (vnorm (invxfm_norm xfm n)) + tex+ in PacketResult (fix ri1 invlenscale1)+ (fix ri2 invlenscale2)+ (fix ri3 invlenscale3)+ (fix ri4 invlenscale4)++-- ugh, code duplication+rayint_debug_instance :: Instance -> Ray -> Flt -> Texture -> (Rayint,Int)+rayint_debug_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_debug sld (Ray neworig (vscale newdir invlenscale)) (d*lenscale) t) of+ (RayMiss, count) -> (RayMiss, count)+ (RayHit depth pos n tex, count) -> (RayHit (depth*invlenscale) + (xfm_point xfm pos) + (vnorm (invxfm_norm xfm n)) + tex, count)++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) ]++transform_leaf_instance :: Instance -> [Xfm] -> SolidItem+transform_leaf_instance (Instance s xfm2) xfm1 =+ transform_leaf 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) = + [SolidItem $ transform_leaf s [xfm]]+ -- group $ map (\x -> transform (flatten_transform x) [xfm]) (tolist s)++primcount_instance :: Instance -> Pcount+primcount_instance (Instance s xfm) = pcadd (primcount s) pcsinglexfm++instance Solid Instance where+ rayint = rayint_instance+ packetint = packetint_instance+ rayint_debug = rayint_debug_instance+ shadow = shadow_instance+ inside = inside_instance+ bound = bound_instance+ transform = transform_instance+ transform_leaf = transform_leaf_instance+ flatten_transform = flatten_transform_instance+ primcount = primcount_instance
+ Data/Glome/Spd.hs view
@@ -0,0 +1,250 @@+{-# LANGUAGE FlexibleInstances #-}++module Data.Glome.Spd where+import Data.Glome.Scene++-- NFF file format description:+-- http://tog.acm.org/resources/SPD/NFF.TXT++-- this would probably be much shorter if I used scanf instead of lex++lexignore s =+ let t = lex s + in+ case t of + [] -> [] -- newline+ [("",s1)] -> lexcr s1+ [(_,s1)] -> lexignore s1 ++lexcr s = + let t = lex s+ in + case t of+ [(s1,[])] -> t+ [("",s1)] -> lexcr s1+ [("#",s1)] -> lexignore s1+ _ -> t++data BgColor = BgColor(Color)++readsSpdVec :: ReadS Vec+readsSpdVec s = [((Vec x y z),s3) | (x,s1) <- reads s :: [(Flt,String)], + (y,s2) <- reads s1 :: [(Flt,String)],+ (z,s3) <- reads s2 :: [(Flt,String)] ]+instance Read Vec where+ readsPrec _ = readsSpdVec++readsSpdVecNorm :: ReadS (Vec,Vec)+readsSpdVecNorm s = [(((Vec x y z),(Vec nx ny nz)),s6) + | (x,s1) <- reads s :: [(Flt,String)], + (y,s2) <- reads s1 :: [(Flt,String)],+ (z,s3) <- reads s2 :: [(Flt,String)],+ (nx,s4) <- reads s3 :: [(Flt,String)],+ (ny,s5) <- reads s4 :: [(Flt,String)],+ (nz,s6) <- reads s5 :: [(Flt,String)] ]++instance Read (Vec,Vec) where+ readsPrec _ = readsSpdVecNorm+++-- if this seems intuitive, there's something wrong with you+readsSpdVecs :: ReadS [Vec]+readsSpdVecs s =+ let parses = reads s :: [(Vec,String)]+ in+ if null parses+ then [([],s)]+ else+ let (v,rest) = head parses+ (vs,returns) = head (readsSpdVecs rest)+ in [((v:vs),returns)]++instance Read [Vec] where+ readsPrec _ = readsSpdVecs++readsSpdVecsNorms :: ReadS [(Vec,Vec)]+readsSpdVecsNorms s = + let parses = readsSpdVecNorm s :: [((Vec,Vec),String)]+ in+ if null parses+ then [([],s)]+ else+ let (v,rest) = head parses+ (vs,returns) = head (readsSpdVecsNorms rest)+ in [((v:vs),returns)]++instance Read [(Vec,Vec)] where+ readsPrec _ = readsSpdVecsNorms++{- "v"+ "from" Fx Fy Fz+ "at" Ax Ay Az+ "up" Ux Uy Uz+ "angle" angle+ "hither" hither+ "resolution" xres yres -}+readsSpdCam :: ReadS Camera+readsSpdCam s = [ (camera from at up angle,s14) | ("v", s1) <- lexcr s,+ ("from", s2) <- lexcr s1,+ (from, s3) <- reads s2 :: [(Vec,String)],+ ("at", s4) <- lexcr s3,+ (at, s5) <- reads s4 :: [(Vec,String)],+ ("up", s6) <- lexcr s5,+ (up, s7) <- reads s6 :: [(Vec,String)],+ ("angle", s8) <- lexcr s7,+ (angle, s9) <- reads s8 :: [(Flt,String)],+ ("hither", s10)<- lexcr s9,+ (_,s11) <- lexcr s10,+ ("resolution", s12) <- lexcr s11,+ (_, s13) <- lexcr s12,+ (_, s14) <- lexcr s13 ]+instance Read Camera where+ readsPrec _ = readsSpdCam++++readsSpdClr :: ReadS Color+readsSpdClr s = [((Color r g b), s3) | (r, s1) <- reads s :: [(Flt,String)],+ (g, s2) <- reads s1 :: [(Flt,String)],+ (b, s3) <- reads s2 :: [(Flt,String)] ]+instance Read Color where+ readsPrec _ = readsSpdClr+++-- "b" R G B+readsSpdBackground :: ReadS BgColor+readsSpdBackground s = [((BgColor clr), s2) | ("b", s1) <- lexcr s,+ (clr, s2) <- reads s1 :: [(Color,String)] ]+instance Read BgColor where+ readsPrec _ = readsSpdBackground+++-- "l" X Y Z [R G B]+readsSpdLight :: ReadS Light+readsSpdLight s = [((Light pos clr),s3) | ("l", s1) <- lexcr s,+ (pos, s2) <- reads s1 :: [(Vec,String)],+ (clr, s3) <- reads s2 :: [(Color,String)] ]+ +++ [((Light pos (Color 1 1 1)),s2) | ("l", s1) <- lexcr s,+ (pos, s2) <- reads s1 :: [(Vec,String)] ]+instance Read Light where+ readsPrec _ = readsSpdLight++-- "f" red green blue Kd Ks Shine T index_of_refraction+-- data Material = Material {clr :: Color, reflect, refract, ior, kd, ks, shine :: Flt}+readsSpdFill :: ReadS Texture+readsSpdFill s = [(\ri->Material clr ks (1-trans) ior kd 0.5 shine, s7) | ("f", s1) <- lexcr s,+ (clr, s2) <- reads s1 :: [(Color,String)],+ (kd, s3) <- reads s2 :: [(Flt,String)],+ (ks, s4) <- reads s3 :: [(Flt,String)],+ (shine, s5) <- reads s4 :: [(Flt,String)],+ (trans, s6) <- reads s5 :: [(Flt,String)],+ (ior, s7) <- reads s6 :: [(Flt,String)] ]++instance Read (Rayint -> Material) where+ readsPrec _ = readsSpdFill+++-- "s" center.x center.y center.z radius++-- "c"+-- base.x base.y base.z base_radius+-- apex.x apex.y apex.z apex_radius++-- "p" total_vertices+-- vert1.x vert1.y vert1.z+-- [etc. for total_vertices vertices]++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)] ]+ +++ [((cone e1 r1 e2 r2),s5) | ("c",s1) <- lexcr s,+ (e1,s2) <- reads s1 :: [(Vec,String)],+ (r1,s3) <- reads s2 :: [(Flt,String)],+ (e2,s4) <- reads s3 :: [(Vec,String)],+ (r2,s5) <- reads s4 :: [(Flt,String)] ]+ +++ [(group (triangles verts),s3) | ("p",s1) <- lexcr s,+ (n,s2) <- reads s1 :: [(Int,String)],+ (verts,s3) <- readsSpdVecs s2 :: [([Vec],String)] ]+ +++ [(group (trianglesnorms (vns)),s3) | ("pp",s1) <- lexcr s,+ (n,s2) <- reads s1 :: [(Int,String)],+ (vns,s3) <- readsSpdVecsNorms s2 :: [([(Vec,Vec)],String)] ]+ {- +++ [(tex(Void,t),s1) | (t,s1) <- reads s :: [(Texture,String)]] -}+++-- instance Read Solid where+-- readsPrec _ = readsSpdSolid+++-- same as readSpdVecs, just different types+readsSpdPrims :: ReadS [SolidItem]+readsSpdPrims s =+ let parses = readsSpdSolid s :: [(SolidItem,String)]+ in+ if null parses+ then [([],s)]+ else+ let (v,rest) = head parses+ (vs,returns) = head (readsSpdPrims rest)+ in [((v:vs),returns)]++instance Read [SolidItem] where+ readsPrec _ = readsSpdPrims+++readsSpdTextureGroup :: ReadS SolidItem+readsSpdTextureGroup s =+ [((tex (bih prims) t),s2) | (t,s1) <- reads s :: [(Texture,String)],+ (prims,s2) <- readsSpdPrims s1 :: [([SolidItem],String)] ]+ +instance Read SolidItem where+ readsPrec _ = readsSpdTextureGroup++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 :: [(SolidItem,String)]+ lit = reads s :: [(Light,String)]+ bgc = reads s :: [(BgColor,String)]+ in+ if not $ null cam+ then+ let (c1,s1) = head cam+ in accum_rss (c1:cams) lights prims background s1+ else+ if not $ null sld+ then+ let (s2,s1) = head sld + in accum_rss cams lights (s2:prims) background s1+ else+ if not $ null lit+ then+ let (l1,s1) = head lit+ in accum_rss cams (l1:lights) prims background s1+ else+ if not $ null bgc+ then+ let (b1,s1) = head bgc+ in accum_rss cams lights prims (b1:background) s1+ else+ (cams,lights,prims,background,s)++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)]++-- | Read instance for scenes described in the Neutral File Format+-- (NFF) used by SPD, a collection of standard benchmark scenes put+-- together by Eric Haines. We support most NFF features, but not+-- all.+instance Read Scene where+ readsPrec _ = readsSpdScene
+ Data/Glome/Sphere.hs view
@@ -0,0 +1,88 @@+{-# OPTIONS_GHC -funbox-strict-fields #-}+{-# LANGUAGE BangPatterns #-}++module Data.Glome.Sphere (sphere) where+import Data.Glome.Vec+import Data.Glome.Solid++data Sphere = Sphere !Vec !Flt !Flt deriving Show+++-- | Construct a sphere given a center location and a radius.+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++packetint_sphere :: Sphere -> Ray -> Ray -> Ray -> Ray -> Flt -> Texture -> PacketResult+packetint_sphere s !r1 !r2 !r3 !r4 !d t =+ PacketResult (rayint_sphere s r1 d t)+ (rayint_sphere s r2 d t)+ (rayint_sphere s r3 d t)+ (rayint_sphere s r4 d t)++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+ packetint = packetint_sphere+ shadow = shadow_sphere+ inside = inside_sphere+ bound = bound_sphere
+ Data/Glome/Tex.hs view
@@ -0,0 +1,57 @@+module Data.Glome.Tex (tex) where+import Data.Glome.Vec+import Data.Glome.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++-- | Associate a texture with an object. For composite+-- objects, the shader uses the innermost texture.+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++rayint_debug_tex :: Tex -> Ray -> Flt -> Texture -> (Rayint,Int)+rayint_debug_tex (Tex s tex) r d t = rayint_debug 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++primcount_tex :: Tex -> Pcount+primcount_tex (Tex s _) = primcount s++instance Solid Tex where+ rayint = rayint_tex+ rayint_debug = rayint_debug_tex+ packetint = packetint_tex+ shadow = shadow_tex+ inside = inside_tex+ bound = bound_tex+ primcount = primcount_tex
+ Data/Glome/Trace.hs view
@@ -0,0 +1,174 @@+{-# OPTIONS_GHC -funbox-strict-fields #-}+{-# LANGUAGE BangPatterns #-}++module Data.Glome.Trace where+import Data.Glome.Scene+import Data.List++{-+We put lighting code in this file because it needs to be +mutually recursive with the trace function, for refraction+and reflection.+ -}++-- | Result of tracing a packet of 4 rays at once.+data PacketColor = PacketColor !Color !Color !Color !Color++{-+class (Show a) => Shader a where+ -- ray intersection, scene, recursion limit+ shade :: Rayint -> Ray -> Scene -> Int -> Color+ shadepacket :: PacketResult -> Ray -> Ray -> Ray -> Ray -> Scene -> Int -> PacketColor++ shadepacket (PacketResult ri1 ri2 ri3 ri4) r1 r2 r3 r4 scene recurs =+ PacketColor (shade ri1 r1 scene recurs)+ (shade ri2 r2 scene recurs)+ (shade ri3 r3 scene recurs)+ (shade ri4 r4 scene recurs)+-}++{-+simple_shade :: Rayint -> [Light] -> Solid -> Color -> Color+simple_shade ri lights s bg =+ case ri of+ (RayHit d p n t) ->+ let (Material clr refl refr ior kd shine) = t ri+ in cscale clr (vdot n (Vec 0.0 1.0 0.0))+ (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 (fabs $ nx/2) (fabs $ ny/2) (fabs $ nz/2))+ RayMiss -> bground scn++-- 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 ks shine) = t ri+ in clr++-- | This is the lighting routine that handles diffuse light, shadows, +-- specular highlights and reflection. Given a ray intersection, the ray,+-- a scene, and a recursion limit, return a color. "Debug" is a parameter+-- useful for debugging; sometimes we might want to tint the color by +-- the number of bounding boxes tested or something similar.+-- Todo: refraction+shade :: Rayint -- ^ ray intersection returned by rayint+ -> Ray -- ^ ray that resuted in the ray intersection+ -> Scene -- ^ scene we're rendering+ -> Int -- ^ recursion limit+ -> Int -- ^ debugging value (usualy not used)+ -> Color -- ^ computed color+shade ri (Ray o indir) scn recurs !debug = + case ri of+ (RayHit d p n t) ->+ let (Material clr refl_ refr ior kd ks shine) = t ri+ s = sld scn+ lights = lits scn+ direct = foldl' cadd c_black + (map (\ (Light lp lc) ->+ let eyedir = vinvert indir+ lvec = vsub lp p+ llen = vlen lvec+ ldir = vscale lvec (1.0/llen) + halfangle = bisect ldir eyedir+ ldotn = fmax 0 $ vdot ldir n+ -- blinn = fmax 0 ((vdot halfangle n)**(shine*3))+ blinn = fmax 0 $ ((vdot halfangle n) ** shine) * ldotn+ blinn_correct = if isNaN blinn then 0 else blinn+ -- indotn = fmax 0 $ vdot eyedir n+ intensity = 5.0 / (llen*llen)+ --intensity = 0.2+ in+ if vdot n lvec < 0 + then c_black+ else+ if not $ shadow s (Ray (vscaleadd p n delta) ldir) (llen-(2*delta))+ then+ cadd + -- diffuse+ --c_black+ (cmul clr $ cscale lc $ ldotn * intensity)+ -- blinn/torrance-sparrow highlight (pbrt p 440)+ (cscale lc $ blinn_correct * intensity * ks)+ -- c_black+ else + 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_+ else c_black+ refract = + if (refr > delta) && (recurs > 0)+ then c_black+ else c_black+ in+ cadd direct $ cadd reflect_ refract++ (RayMiss) -> bground scn++-- | Given a scene, a ray, a maximum distance, and a maximum+-- recursion depth, test the ray for intersection against +-- the object within the scene, then pass the ray intersection+-- to the shade routine (which may trace secondary rays of its +-- own), which returns a color. For most applications, this is+-- the entry point into the ray tracer.+trace :: Scene -> Ray -> Flt -> Int -> Color+trace scn ray depth recurs =+ let (Scene sld lights cam dtex bgcolor) = scn + in shade (rayint sld ray depth dtex) ray scn recurs 0+ +-- | Similar to trace, but return depth as well as color.+-- We might want the depth 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)++-- | Similar to trace, but 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 + 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)++-- | A trace function which returns some additional debugging+-- info, mainly for performance tuning.+trace_debug :: Scene -> Ray -> Flt -> Int -> Color+trace_debug scn ray depth recurs =+ let (Scene sld lights cam dtex bgcolor) = scn+ (ri,n) = rayint_debug sld ray depth dtex+ in + cadd (shade ri ray scn recurs 0) (Color 0 ((fromIntegral (Prelude.abs n)) * 0.01) 0)++-- | Trace a packet of four rays at a time. Sometimes, this+-- may be a performance advantage. However, ever since my +-- transition to typeclasses, this has not performed any better+-- than the mono-ray path.+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)
+ Data/Glome/Triangle.hs view
@@ -0,0 +1,141 @@+module Data.Glome.Triangle where+import Data.Glome.Vec+import Data.Glome.Solid++-- Simple triangles, and triangles with normal vectors+-- specified at each vertex.++data Triangle = Triangle Vec Vec Vec deriving Show+data TriangleNorm = TriangleNorm Vec Vec Vec Vec Vec Vec deriving Show++-- | Create a simple triangle from its 3 corners.+-- The normals are computed automatically.+triangle :: Vec -> Vec -> Vec -> SolidItem+triangle v1 v2 v3 =+ SolidItem (Triangle v1 v2 v3)++-- | Create a triangle fan from a list of verticies.+triangles :: [Vec] -> [SolidItem]+triangles (v1:vs) =+ zipWith (\v2 v3 -> triangle v1 v2 v3) vs (tail vs) ++-- | Create a triangle from a list of verticies, and +-- a list of normal vectors (one for each vertex).+trianglenorm v1 v2 v3 n1 n2 n3 =+ SolidItem (TriangleNorm v1 v2 v3 n1 n2 n3)++-- | Create a triangle fan from a list of verticies and normals.+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)++transform_triangle :: Triangle -> [Xfm] -> SolidItem+transform_triangle (Triangle p1 p2 p3) xfms =+ SolidItem $ Triangle (xfm_point (compose xfms) p1)+ (xfm_point (compose xfms) p2)+ (xfm_point (compose xfms) p3)++transform_trianglenorm :: TriangleNorm -> [Xfm] -> SolidItem+transform_trianglenorm (TriangleNorm p1 p2 p3 n1 n2 n3) xfms =+ SolidItem $ TriangleNorm (xfm_point (compose xfms) p1)+ (xfm_point (compose xfms) p2)+ (xfm_point (compose xfms) p3)+ (vnorm $ xfm_vec (compose xfms) n1)+ (vnorm $ xfm_vec (compose xfms) n2)+ (vnorm $ xfm_vec (compose xfms) n3)++instance Solid Triangle where+ rayint = rayint_triangle+ inside _ _ = False+ bound = bound_triangle+ transform = transform_triangle++instance Solid TriangleNorm where+ rayint = rayint_trianglenorm+ inside _ _ = False+ bound = bound_trianglenorm+ transform = transform_trianglenorm
+ GlomeTrace.cabal view
@@ -0,0 +1,34 @@+Name: GlomeTrace+Version: 0.1.1+Synopsis: Ray Tracing Library+Description: A ray tracing library with acceleration structure and many supported primitives.+License: GPL+License-file: LICENSE+Author: Jim Snow+Maintainer: Jim Snow <jsnow@cs.pdx.edu>+Copyright: Copyright 2008,2009 Jim Snow+Homepage: http://www.haskell.org/haskellwiki/Glome+Stability: experimental+Category: graphics+build-type: Simple+Cabal-Version: >= 1.2+extra-source-files:+ README.txt++library+ exposed-modules: Data.Glome.Trace+ Data.Glome.Scene+ Data.Glome.Clr+ Data.Glome.Solid+ Data.Glome.Spd+ Data.Glome.Bih+ Data.Glome.Bound+ Data.Glome.Box+ Data.Glome.Cone+ Data.Glome.Csg+ Data.Glome.Plane+ Data.Glome.Sphere+ Data.Glome.Tex+ Data.Glome.Triangle++ Build-Depends: base >= 3 && < 4, array, GlomeVec >= 0.1.1
+ LICENSE view
@@ -0,0 +1,14 @@+ This library, GlomeVec, is copyright 2008 Jim Snow++ This program is free software; you can redistribute it and/or modify+ it under the terms of version 2 of the GNU General Public License as + published by the Free Software Foundation;++ This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License for more details.++ You should have received a copy of the GNU General Public License+ along with this program; if not, write to the Free Software+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ README.txt view
@@ -0,0 +1,15 @@+This is GlomeTrace, a ray tracing library. Originally, it was part of Glome-hs, my haskell ray tracer. I decided to pull out the code to trace rays against objects and distribute it as a stand-alone library so that other projects can use it.++Glome-hs required HOpenGL, but only for display. Since this library is separated from any notion of display, it does not require HOpenGL, so it should run on more platforms than Glome-hs did.++A good source of documentation is the Haskell wiki. In particular, take a look at the tutorial I wrote for Glome-hs. A few things have changed, but most of the descriptions there are still valid.++I have begun to add haddock documentation to the actual code.++You will need to install GlomeVec first (or have cabal fetch it automatically).++http://www.haskell.org/haskellwiki/Glome++Direct all questions to:+Jim Snow+jsnow@cs.pdx.edu
+ Setup.hs view
@@ -0,0 +1,6 @@+module Main (main) where++import Distribution.Simple++main :: IO ()+main = defaultMain