GlomeVec (empty) → 0.1
raw patch · 6 files changed
+803/−0 lines, 6 filesdep +arraydep +basesetup-changed
Dependencies added: array, base
Files
- Data/GlomeTexture.hs +126/−0
- Data/GlomeVec.hs +617/−0
- GlomeVec.cabal +20/−0
- LICENSE +14/−0
- README.txt +20/−0
- Setup.hs +6/−0
+ Data/GlomeTexture.hs view
@@ -0,0 +1,126 @@+{-+module SolidTexture (square_wave, triangle_wave, sine_wave, + stripe, noise, turbulence+ ) where -}++module Data.GlomeTexture where+import Data.GlomeVec+import Data.Array.IArray++-- INTERPOLATION FUNCTIONS --+square_wave :: Flt -> Flt+square_wave x =+ let offset = x - (fromIntegral (floor x))+ in if offset < 0.5 then 0 else 1++triangle_wave :: Flt -> Flt+triangle_wave x =+ let offset = x - (fromIntegral (floor x))+ in if offset < 0.5 + then (offset*2)+ else (2-(offset*2))++sine_wave :: Flt -> Flt+sine_wave x = (sin (x*2*pi))*0.5 + 0.5+++lump_wave :: Flt -> Flt+lump_wave x = 1 - x*x*x++-- SCALAR TEXTURE FUNCTIONS --++-- These are simple solid texture functions that take a +-- point as argument and return a number 0 < n < 1++stripe :: Vec -> (Flt -> Flt) -> (Vec -> Flt)+stripe axis interp =+ let len = vlen axis + in+ (\pos -> let offset = vdot pos axis + in interp offset)+++-- PERLIN NOISE --++-- (-6 t^5 + 15 t^4 - 10t^3 +1)+-- "realistic ray tracing 2nd edition" inconsistent +-- on whether it should be t^5 or t^6,+-- but t^5 works and t^6 doesn't.+omega :: Flt -> Flt+omega t_ = + let t = fabs t_+ tsqr = t*t+ tcube = tsqr*t+ in (-6)*tcube*tsqr + 15*tcube*t - 10*tcube + 1++-- questionably random+phi :: Array Int Int+phi = listArray (0,11) [3,0,2,7,4,1,5,11,8,10,9,6]++grad :: Array Int Vec+grad = listArray (0,11) + $ filter (\x -> let l = vlen x in l < 1.5 && l > 1.1) + [Vec x y z | x <- [(-1),0,1],+ y <- [(-1),0,1],+ z <- [(-1),0,1]] ++gamma :: Int -> Int -> Int -> Vec+gamma i j k =+ let a = phi!(mod (iabs k) 12)+ b = phi!(mod (iabs (j+a)) 12)+ c = phi!(mod (iabs (i+b)) 12)+ in grad!c++knot :: Int -> Int -> Int -> Vec -> Flt+knot i j k v =+ let Vec x y z = v+ in (omega x) * (omega y) * (omega z) * (vdot (gamma i j k) v)++intGamma :: Int -> Int -> Int+intGamma i j =+ let a = phi!(mod (iabs j) 16)+ b = phi!(mod (iabs (i+a)) 16)+ in b++turbulence :: Vec -> Int -> Flt+turbulence p 1 = fabs(noise(p))+turbulence p n =+ let newp = vscale p 0.5+ t = fabs (noise p)+ in t + (0.5 * (turbulence newp (n-1)))++noise :: Vec -> Flt +noise (Vec x y z) =+ let i = floor x+ j = floor y+ k = floor z+ u = x-(fromIntegral i)+ v = y-(fromIntegral j)+ w = z-(fromIntegral k)+ in knot i j k (Vec u v w) ++ knot (i+1) j k (Vec (u-1) v w) ++ knot i (j+1) k (Vec u (v-1) w) ++ knot i j (k+1) (Vec u v (w-1)) ++ knot (i+1) (j+1) k (Vec (u-1) (v-1) w) ++ knot (i+1) j (k+1) (Vec (u-1) v (w-1)) ++ knot i (j+1) (k+1) (Vec u (v-1) (w-1)) ++ knot (i+1) (j+1) (k+1) (Vec (u-1) (v-1) (w-1))++perlin :: Vec -> Flt+perlin v =+ let p = ((noise v)+1)*0.5+ in if p > 1 + then error $ "perlin noise error, 1 < " ++ (show p)+ else if p < 0 + then error $ "perlin noise error, 0 > " ++ (show p)+ else p++--untested+perlin_turb :: Vec -> Int -> Flt+perlin_turb v l =+ let p = turbulence v l+ in if p > 1 + then error $ "perlin turbulence error, 1 < " ++ (show p)+ else if p < 0 + then error $ "perlin turbulence error, 0 > " ++ (show p)+ else p
+ Data/GlomeVec.hs view
@@ -0,0 +1,617 @@+{-# OPTIONS_GHC -fexcess-precision #-}+{-# OPTIONS_GHC -funbox-strict-fields #-}+{-# LANGUAGE BangPatterns #-}++module Data.GlomeVec where++-- Performance is pretty similar with Floats or Doubles+-- best performance seems to be doubles with -fvia-C+type Flt = Double++-- maybe this is defined somewhere?+infinity :: Flt+--infinity = 1.0 / 0.0+infinity = 1000000.0++-- convert from degrees to native angle format (radians)+deg :: Flt -> Flt+deg !x = (x*3.1415926535897)/180++-- convert from radians (noop)+rad :: Flt -> Flt+rad !x = x++-- convert from rotations+rot :: Flt -> Flt+rot !x = x*3.1415926535897*2++-- trig with degrees +dcos :: Flt -> Flt+dcos d = cos $ deg d++-- force a value to be within a range+clamp :: Flt -> Flt -> Flt -> Flt+clamp !min !x !max+ | x < min = min+ | x > max = max+ | otherwise = x++-- delta = 0.00001 :: Flt+delta = 0.0001 :: Flt++-- non-polymorphic versions; this speeds+-- things up in ocaml, not sure about haskell+fmin :: Flt -> Flt -> Flt+fmin !a !b = if a > b then b else a++fmax :: Flt -> Flt -> Flt+fmax !a !b = if a > b then a else b++fmin3 :: Flt -> Flt -> Flt -> Flt+fmin3 !a !b !c = if a > b + then if b > c + then c+ else b+ else if a > c+ then c+ else a++fmax3 :: Flt -> Flt -> Flt -> Flt+fmax3 !a !b !c = if a > b+ then if a > c+ then a+ else c+ else if b > c+ then b+ else c++fmin4 :: Flt -> Flt -> Flt -> Flt -> Flt+fmin4 !a !b !c !d = fmin (fmin a b) (fmin c d)++fmax4 :: Flt -> Flt -> Flt -> Flt -> Flt+fmax4 !a !b !c !d = fmax (fmax a b) (fmax c d)++fabs :: Flt -> Flt+fabs !a = + if a < 0 then (-a) else a++iabs :: Int -> Int+iabs !a =+ if a < 0 then (-a) else a++abs a = error "use non-polymorphic version, fabs"++-- true if a and b are "almost" equal+-- the (abs $ a-b) test doesn't work if+-- a and b are large+about_equal :: Flt -> Flt -> Bool+about_equal !a !b =+ if a > 1 + then+ fabs (1 - (a/b)) < (delta*10) + else+ (fabs $ a - b) < (delta*10)+++data Vec = Vec !Flt !Flt !Flt deriving Show++data Ray = Ray {origin, dir :: !Vec} deriving Show+--data Plane = Plane {norm :: !Vec, offset :: !Flt} deriving Show++vec !x !y !z = (Vec x y z)+vzero = Vec 0.0 0.0 0.0++-- for when we need a unit vector, but we +-- don't care where it points+vunit = vx++-- unit axis vectors+vx = Vec 1 0 0+vy = Vec 0 1 0+vz = Vec 0 0 1+nvx = Vec (-1) 0 0+nvy = Vec 0 (-1) 0+nvz = Vec 0 0 (-1)++x (Vec x_ _ _) = x_+y (Vec _ y_ _) = y_+z (Vec _ _ z_) = z_++-- this actually accounts for a+-- noticeable amount of cpu time+va :: Vec -> Int -> Flt+va !(Vec x y z) !n = + case n of+ 0 -> x+ 1 -> y+ 2 -> z++vset :: Vec -> Int -> Flt -> Vec+vset !(Vec x y z) !i !f =+ case i of+ 0 -> Vec f y z+ 1 -> Vec x f z+ 2 -> Vec x y f++vdot :: Vec -> Vec -> Flt+vdot !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ (x1*x2)+(y1*y2)+(z1*z2)++vcross :: Vec -> Vec -> Vec+vcross !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ Vec + ((y1 * z2) - (z1 * y2))+ ((z1 * x2) - (x1 * z2))+ ((x1 * y2) - (y1 * x2))++vmap :: (Flt -> Flt) -> Vec -> Vec+vmap f !v1 = + Vec (f (x v1)) (f (y v1)) (f (z v1))++vmap2 :: (Flt -> Flt -> Flt) -> Vec -> Vec -> Vec+vmap2 f !v1 !v2 =+ Vec (f (x v1) (x v2)) + (f (y v1) (y v2)) + (f (z v1) (z v2))++vinvert :: Vec -> Vec+vinvert !(Vec x1 y1 z1) =+ Vec (-x1) (-y1) (-z1)++vlensqr :: Vec -> Flt+vlensqr !v1 = vdot v1 v1++vlen :: Vec -> Flt+vlen !v1 = sqrt (vdot v1 v1)++vadd :: Vec -> Vec -> Vec+vadd !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ Vec (x1 + x2)+ (y1 + y2)+ (z1 + z2)++vadd3 :: Vec -> Vec -> Vec -> Vec+vadd3 !(Vec x1 y1 z1) !(Vec x2 y2 z2) !(Vec x3 y3 z3) =+ Vec (x1 + x2 + x3)+ (y1 + y2 + y3)+ (z1 + z2 + z3)++vsub :: Vec -> Vec -> Vec+vsub !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ Vec (x1 - x2)+ (y1 - y2)+ (z1 - z2)++vmul :: Vec -> Vec -> Vec+vmul !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ Vec (x1 * x2)+ (y1 * y2)+ (z1 * z2)++vinc :: Vec -> Flt -> Vec+vinc !(Vec x y z) !n =+ Vec (x + n)+ (y + n)+ (z + n)++vdec :: Vec -> Flt -> Vec+vdec !(Vec x y z) !n =+ Vec (x - n)+ (y - n)+ (z - n)++vmax :: Vec -> Vec -> Vec+vmax !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ Vec (fmax x1 x2)+ (fmax y1 y2)+ (fmax z1 z2)++vmin :: Vec -> Vec -> Vec+vmin !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ Vec (fmin x1 x2)+ (fmin y1 y2)+ (fmin z1 z2)++vmaxaxis :: Vec -> Int+vmaxaxis !(Vec x y z) =+ if (x > y) + then if (x > z) + then 0+ else 2+ else if (y > z) + then 1+ else 2++vscale :: Vec -> Flt -> Vec+vscale !(Vec x y z) !fac =+ Vec (x * fac)+ (y * fac)+ (z * fac)++vscaleadd :: Vec -> Vec -> Flt -> Vec+vscaleadd !(Vec x1 y1 z1) !(Vec x2 y2 z2) fac =+ Vec (x1 + (x2 * fac))+ (y1 + (y2 * fac))+ (z1 + (z2 * fac))+ +-- make the length just a little shorter+vnudge :: Vec -> Vec+vnudge x = vscale x (1-delta)++vnorm :: Vec -> Vec+vnorm !(Vec x1 y1 z1) = + let !invlen = 1.0 / (sqrt ((x1*x1)+(y1*y1)+(z1*z1))) in+ Vec (x1*invlen) (y1*invlen) (z1*invlen)++assert_norm :: Vec -> Vec+assert_norm v =+ let l = vdot v v+ in if l > (1+delta) + then error $ "vector too long" ++ (show v)+ else if l < (1-delta)+ then error $ "vector too short: " ++ (show v)+ else v++bisect :: Vec -> Vec -> Vec+bisect !v1 !v2 = vnorm (vadd v1 v2)++vdist :: Vec -> Vec -> Flt+vdist v1 v2 = + let d = vsub v2 v1 in vlen d++reflect :: Vec -> Vec -> Vec+reflect !v !norm =+ -- vadd v $ vscale norm $ (-2) * (vdot v norm)+ vscaleadd v norm $ (-2) * (vdot v norm)++vrcp :: Vec -> Vec+vrcp !(Vec x y z) =+ Vec (1/x) (1/y) (1/z)++-- test for equality+veq :: Vec -> Vec -> Bool+veq !(Vec ax ay az) !(Vec bx by bz) =+ (about_equal ax bx) && (about_equal ay by) && (about_equal az bz)++--returns false on zero value+veqsign :: Vec -> Vec -> Bool+veqsign !(Vec ax ay az) !(Vec bx by bz) =+ ax*bx > 0 && ay*by > 0 && az*bz > 0++-- translate a ray's origin in ray's direction by d amount+ray_move :: Ray -> Flt -> Ray+ray_move !(Ray orig dir) !d =+ (Ray (vscaleadd orig dir d) dir)++-- find orthogonal vectors+orth :: Vec -> (Vec,Vec)+orth v1 =+ if about_equal (vdot v1 v1) 1+ then+ let x = (Vec 1 0 0)+ y = (Vec 0 1 0)+ dvx = vdot v1 x+ v2 = if dvx < 0.8 && dvx > (-0.8) -- don't want to cross with a+ then vnorm $ vcross v1 x -- vector that's too similar+ else vnorm $ vcross v1 y+ v3 = vcross v1 v2+ in (v2,v3)+ else error $ "orth: unnormalized vector" ++ (show v1)++-- intersect a ray with a plane +-- defined by a point and a normal+-- (ray need not be normalized)+plane_int :: Ray -> Vec -> Vec -> Vec+plane_int !(Ray orig dir) !p !norm =+ let newo = vsub orig p+ dist = -(vdot norm newo) / (vdot norm dir)+ in vscaleadd orig dir dist++plane_int_dist :: Ray -> Vec -> Vec -> Flt+plane_int_dist !(Ray orig dir) !p !norm =+ let newo = vsub orig p+ in -(vdot norm newo) / (vdot norm dir)++-- find intersection with plane+-- from graphics gems -- an efficient ray-polygon intersection+-- it seems that the ray need not be normalized+-- let plane_intersect ray (n,d) =+-- let t = -.((d +. (vdot n ray.origin)) /. (vdot n ray.dir))+-- in vadd ray.origin (vscale ray.dir t)+++-- TRANSFORMATIONS --++data Matrix = Matrix !Flt !Flt !Flt !Flt + !Flt !Flt !Flt !Flt + !Flt !Flt !Flt !Flt deriving Show++-- this is a little faster if the matricies are non-strict+data Xfm = Xfm Matrix Matrix deriving Show++ident_matrix = (Matrix 1 0 0 0 0 1 0 0 0 0 1 0)+ident_xfm = Xfm ident_matrix ident_matrix++mat_mult :: Matrix -> Matrix -> Matrix+mat_mult (Matrix a00 a01 a02 a03 a10 a11 a12 a13 a20 a21 a22 a23)+ (Matrix b00 b01 b02 b03 b10 b11 b12 b13 b20 b21 b22 b23) =+ Matrix+ (a00*b00 + a01*b10 + a02*b20)+ (a00*b01 + a01*b11 + a02*b21)+ (a00*b02 + a01*b12 + a02*b22)+ (a00*b03 + a01*b13 + a02*b23 + a03)++ (a10*b00 + a11*b10 + a12*b20)+ (a10*b01 + a11*b11 + a12*b21)+ (a10*b02 + a11*b12 + a12*b22)+ (a10*b03 + a11*b13 + a12*b23 + a13)++ (a20*b00 + a21*b10 + a22*b20)+ (a20*b01 + a21*b11 + a22*b21)+ (a20*b02 + a21*b12 + a22*b22)+ (a20*b03 + a21*b13 + a22*b23 + a23)++xfm_mult :: Xfm -> Xfm -> Xfm+xfm_mult (Xfm a inva) (Xfm b invb) =+ Xfm (mat_mult a b) (mat_mult invb inva)++-- TRANSFORM UTILITY FUNCTIONS --++-- If we multiply two transformation matricies, we get+-- a transformation matrix equivalent to applying the +-- second then the first.++-- By reversing the list, the transforms are applied in the expected order.+compose :: [Xfm] -> Xfm+compose xfms = check_xfm $ foldr xfm_mult ident_xfm (reverse xfms)++check_xfm :: Xfm -> Xfm+check_xfm (Xfm m i) = + let (Matrix m00 m01 m02 m03 + m10 m11 m12 m13 + m20 m21 m22 m23) = mat_mult m i+ ae = about_equal+ in+ if ae m00 1 && ae m01 0 && ae m02 0 && ae m03 0 &&+ ae m10 0 && ae m11 1 && ae m12 0 && ae m13 0 &&+ ae m20 0 && ae m21 0 && ae m22 1 && ae m23 0+ then (Xfm m i)+ else error $ "corrupt matrix " ++ (show (Xfm m i)) ++ "\n" ++ (show (mat_mult m i)) ++-- rotate point (or vector) a about ray b by angle c+vrotate :: Vec -> Ray -> Flt -> Vec+vrotate pt (Ray orig axis_) angle =+ let axis = assert_norm axis_+ transform = compose [ translate (vinvert orig)+ , rotate axis angle+ , translate orig+ ]+ new_pt = xfm_point transform pt+ in if about_equal (vlen (vsub orig pt)) (vlen (vsub orig new_pt))+ then new_pt+ else error $ "something is wrong with vrotate" ++ + (show $ vlen (vsub orig pt)) ++ " " ++ + (show $ vlen (vsub orig new_pt))+++-- TRANSFORM APPLICATION --+-- these need to be fast++-- point is treated as (x y z 1)+xfm_point :: Xfm -> Vec -> Vec+xfm_point !(Xfm (Matrix m00 m01 m02 m03 + m10 m11 m12 m13 + m20 m21 m22 m23) inv) + !(Vec x y z) =+ Vec (m00*x + m01*y + m02*z + m03)+ (m10*x + m11*y + m12*z + m13)+ (m20*x + m21*y + m22*z + m23)++invxfm_point :: Xfm -> Vec -> Vec+invxfm_point !(Xfm fwd (Matrix i00 i01 i02 i03 + i10 i11 i12 i13 + i20 i21 i22 i23)) + !(Vec x y z) =+ Vec (i00*x + i01*y + i02*z + i03)+ (i10*x + i11*y + i12*z + i13)+ (i20*x + i21*y + i22*z + i23)++-- vector is treated as (x y z 0)+xfm_vec :: Xfm -> Vec -> Vec+xfm_vec !(Xfm (Matrix m00 m01 m02 m03 + m10 m11 m12 m13 + m20 m21 m22 m23) inv) + !(Vec x y z) =+ Vec (m00*x + m01*y + m02*z)+ (m10*x + m11*y + m12*z)+ (m20*x + m21*y + m22*z)++invxfm_vec :: Xfm -> Vec -> Vec+invxfm_vec !(Xfm fwd (Matrix i00 i01 i02 i03 + i10 i11 i12 i13 + i20 i21 i22 i23)) + !(Vec x y z) =+ Vec (i00*x + i01*y + i02*z)+ (i10*x + i11*y + i12*z)+ (i20*x + i21*y + i22*z)++-- this one is tricky+-- we transform by the inverse transpose+invxfm_norm :: Xfm -> Vec -> Vec+invxfm_norm !(Xfm fwd (Matrix i00 i01 i02 i03 + i10 i11 i12 i13 + i20 i21 i22 i23)) + !(Vec x y z) =+ Vec (i00*x + i10*y + i20*z)+ (i01*x + i11*y + i21*z)+ (i02*x + i12*y + i22*z)++xfm_ray :: Xfm -> Ray -> Ray+xfm_ray !xfm !(Ray orig dir) =+ Ray (xfm_point xfm orig) (vnorm (xfm_vec xfm dir))++invxfm_ray !xfm !(Ray orig dir) =+ Ray (invxfm_point xfm orig) (vnorm (invxfm_vec xfm dir))++-- BASIC TRANSFORMS --+-- move+translate :: Vec -> Xfm+translate (Vec x y z) =+ check_xfm $ Xfm (Matrix 1 0 0 x 0 1 0 y 0 0 1 z) + (Matrix 1 0 0 (-x) 0 1 0 (-y) 0 0 1 (-z))++-- strectch along three axes (if x==y==z, then it's uniform scaling)+scale :: Vec -> Xfm+scale (Vec x y z) =+ check_xfm $ Xfm (Matrix x 0 0 0 0 y 0 0 0 0 z 0)+ (Matrix (1/x) 0 0 0 0 (1/y) 0 0 0 0 (1/z) 0)++-- rotate about an arbitrary axis and angle+rotate :: Vec -> Flt -> Xfm+rotate v@(Vec x y z) angle =+ if not $ (vlen v) `about_equal` 1+ then error $ "please use a normalized vector for rotation: " ++ (show (vlen v))+ else + let s = sin angle+ c = cos angle ++ m00 = ((x*x)+((1-(x*x))*c)) + m01 = (((x*y)*(1-c))-(z*s)) + m02 = ((x*z*(1-c))+(y*s))++ m10 = (((x*y)*(1-c))+(z*s))+ m11 = ((y*y)+((1-(y*y))*c))+ m12 = ((y*z*(1-c))-(x*s))++ m20 = ((x*z*(1-c))-(y*s))+ m21 = ((y*z*(1-c))+(x*s))+ m22 = ((z*z)+((1-(z*z))*c))+ in+ check_xfm $ Xfm (Matrix m00 m01 m02 0 m10 m11 m12 0 m20 m21 m22 0)+ (Matrix m00 m10 m20 0 m01 m11 m21 0 m02 m12 m22 0)++-- convert canonical coordinates to uvw coordinates+xyz_to_uvw :: Vec -> Vec -> Vec -> Xfm+xyz_to_uvw u v w =+ let Vec ux uy uz = u+ Vec vx vy vz = v+ Vec wx wy wz = w+ in if (vdot u u) `about_equal` 1+ then+ if (vdot v v) `about_equal` 1+ then+ if (vdot w w) `about_equal` 1+ then + if ((vdot u v) `about_equal` 0) && + ((vdot u w) `about_equal` 0) && + ((vdot v w) `about_equal` 0)+ then+ check_xfm $ Xfm (Matrix ux vx wx 0 uy vy wy 0 uz vz wz 0)+ (Matrix ux uy uz 0 vx vy vz 0 wx wy wz 0)+ else error "vectors aren't orthogonal"+ else error $ "unnormalized w " ++ (show w)+ else error $ "unnormalized v " ++ (show v)+ else error $ "unnormalized u " ++ (show u)++uvw_to_xyz :: Vec -> Vec -> Vec -> Xfm+uvw_to_xyz (Vec ux uy uz) (Vec vx vy vz) (Vec wx wy wz) =+ check_xfm $ Xfm (Matrix ux uy uz 0 vx vy vz 0 wx wy wz 0)+ (Matrix ux vx wx 0 uy vy wy 0 uz vz wz 0)++++-- TRIANGLE UTILITY FUNCTIONS --++-- given a side, angle, and side of a triangle, produce the length of the opposite side+sas2s :: Flt -> Flt -> Flt -> Flt+sas2s s1 a s2 =+ sqrt (((s1 * s1) + (s2 * s2)) - ((2 * s1 * s2 * (dcos a))))++++-- BOUNDING BOXES --+data Bbox = Bbox {p1 :: !Vec, p2 :: !Vec} deriving Show+data Interval = Interval !Flt !Flt deriving Show -- used instead of a tuple++--union of two bounding boxes+bbjoin :: Bbox -> Bbox -> Bbox+bbjoin (Bbox p1a p2a) (Bbox p1b p2b) =+ (Bbox (vmin p1a p1b) (vmax p2a p2b))++--overlap of two bounding boxes+bboverlap :: Bbox -> Bbox -> Bbox+bboverlap (Bbox p1a p2a) (Bbox p1b p2b) =+ (Bbox (vmax p1a p1b) (vmin p2a p2b))++bbinside :: Bbox -> Vec -> Bool+bbinside (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) (Vec x y z) =+ p1x <= x && x <= p2x && p1y <= y && y <= p2y && p1z <= z && z <= p2z++--split a bounding box into two+bbsplit :: Bbox -> Int -> Flt -> (Bbox,Bbox)+bbsplit (Bbox p1 p2) axis offset =+ if (offset < (va p1 axis)) || (offset > (va p2 axis))+ then error "degenerate bounding box split"+ else ((Bbox p1 (vset p2 axis offset)),+ (Bbox (vset p1 axis offset) p2))++-- generate a bounding box from a list of points+bbpts :: [Vec] -> Bbox+bbpts [] = empty_bbox+bbpts ((Vec x y z):[]) =+ Bbox (Vec (x-delta) (y-delta) (z-delta)) + (Vec (x+delta) (y+delta) (z+delta))++bbpts ((Vec x y z):pts) =+ let (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bbpts pts+ minx = fmin (x-delta) p1x+ miny = fmin (y-delta) p1y+ minz = fmin (z-delta) p1z+ maxx = fmax (x+delta) p2x+ maxy = fmax (y+delta) p2y+ maxz = fmax (z+delta) p2z in+ Bbox (Vec minx miny minz) (Vec maxx maxy maxz)++-- surface area, volume of bounding boxes+bbsa :: Bbox -> Flt+bbsa (Bbox p1 p2) =+ let Vec dx dy dz = vsub p2 p1 + in dx*dy + dx*dz + dy*dz++bbvol :: Bbox -> Flt+bbvol (Bbox p1 p2) =+ let (Vec dx dy dz) = vsub p2 p1+ in dx*dy*dz++empty_bbox = + Bbox (Vec infinity infinity infinity) + (Vec (-infinity) (-infinity) (-infinity))++everything_bbox =+ Bbox (Vec (-infinity) (-infinity) (-infinity))+ (Vec infinity infinity infinity)++-- Find a ray's entrance and exit from a bounding +-- box. If last entrance is before the first exit,+-- we hit. Otherwise, we miss. (It's up to the +-- caller to figure that out.)++bbclip :: Ray -> Bbox -> Interval+bbclip (Ray (Vec ox oy oz) (Vec dx dy dz)) + (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) =+ let dxrcp = 1/dx+ dyrcp = 1/dy+ dzrcp = 1/dz+ Interval inx outx = if dx > 0 + then Interval ((p1x-ox)*dxrcp) ((p2x-ox)*dxrcp)+ else Interval ((p2x-ox)*dxrcp) ((p1x-ox)*dxrcp)+ Interval iny outy = if dy > 0+ then Interval ((p1y-oy)*dyrcp) ((p2y-oy)*dyrcp)+ else Interval ((p2y-oy)*dyrcp) ((p1y-oy)*dyrcp)+ Interval inz outz = if dz > 0+ then Interval ((p1z-oz)*dzrcp) ((p2z-oz)*dzrcp)+ else Interval ((p2z-oz)*dzrcp) ((p1z-oz)*dzrcp)+ in+ Interval (fmax3 inx iny inz) (fmin3 outx outy outz)++
+ GlomeVec.cabal view
@@ -0,0 +1,20 @@+Name: GlomeVec+Version: 0.1+Synopsis: Simple 3D vector library+Description: A simple library for dealing with 3D vectors, suitable for graphics projects. A small texture library with Perlin noise is included as well.+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.GlomeVec, Data.GlomeTexture+ Build-Depends: base >= 3 && < 4, array
+ 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,20 @@+This is the vector library used by the Glome raytracer. It has been separated +out of the main Glome distribution as the beginning of an effort to separate +all of the separately-useful modules in Glome. (I have not yet removed the+Vec library from the main Glome distribution.)++This library may prove useful for graphics and computational geometry+algorithms. It includes basic operations such as dot product, adding vectors,+etc, but it also includes transformations matricies and some useful operations+on axis-aligned bounding boxes such as clipping a ray to an AABB. ++See the Glome tutorial on the haskell wiki for details.++This was one of the first things I wrote in Haskell. As such, it has a few+rough edges.++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