diff --git a/Data/Glome/Bih.hs b/Data/Glome/Bih.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Bih.hs @@ -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 diff --git a/Data/Glome/Bound.hs b/Data/Glome/Bound.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Bound.hs @@ -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 diff --git a/Data/Glome/Box.hs b/Data/Glome/Box.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Box.hs @@ -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 diff --git a/Data/Glome/Clr.hs b/Data/Glome/Clr.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Clr.hs @@ -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) diff --git a/Data/Glome/Cone.hs b/Data/Glome/Cone.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Cone.hs @@ -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 diff --git a/Data/Glome/Csg.hs b/Data/Glome/Csg.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Csg.hs @@ -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 diff --git a/Data/Glome/Plane.hs b/Data/Glome/Plane.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Plane.hs @@ -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 diff --git a/Data/Glome/Scene.hs b/Data/Glome/Scene.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Scene.hs @@ -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) +-} diff --git a/Data/Glome/Solid.hs b/Data/Glome/Solid.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Solid.hs @@ -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 diff --git a/Data/Glome/Spd.hs b/Data/Glome/Spd.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Spd.hs @@ -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 diff --git a/Data/Glome/Sphere.hs b/Data/Glome/Sphere.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Sphere.hs @@ -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 diff --git a/Data/Glome/Tex.hs b/Data/Glome/Tex.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Tex.hs @@ -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 diff --git a/Data/Glome/Trace.hs b/Data/Glome/Trace.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Trace.hs @@ -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) diff --git a/Data/Glome/Triangle.hs b/Data/Glome/Triangle.hs new file mode 100644 --- /dev/null +++ b/Data/Glome/Triangle.hs @@ -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 diff --git a/GlomeTrace.cabal b/GlomeTrace.cabal new file mode 100644 --- /dev/null +++ b/GlomeTrace.cabal @@ -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 +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 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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 diff --git a/README.txt b/README.txt new file mode 100644 --- /dev/null +++ b/README.txt @@ -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 diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,6 @@ +module Main (main) where + +import Distribution.Simple + +main :: IO () +main = defaultMain