GlomeVec 0.1.1 → 0.1.2

raw patch · 5 files changed

+867/−867 lines, 5 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Data/Glome/Texture.hs +126/−0
• Data/Glome/Vec.hs +738/−0
• Data/GlomeTexture.hs +0/−126
• Data/GlomeVec.hs +0/−738
• GlomeVec.cabal +3/−3

@@ -0,0 +1,126 @@+{-+module SolidTexture (square_wave, triangle_wave, sine_wave, +                     stripe, noise, turbulence+                    ) where -}++module Data.Glome.Texture where+import Data.Glome.Vec+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

@@ -0,0 +1,738 @@+{-# OPTIONS_GHC -fexcess-precision #-}+{-# OPTIONS_GHC -funbox-strict-fields #-}+{-# LANGUAGE BangPatterns #-}++module Data.Glome.Vec where++-- | Performance is pretty similar with Floats or Doubles.+-- Todo: make separate Float and Double instances of this library.+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 to native format (noop).+rad :: Flt -> Flt+rad !x = x++-- | Convert from rotations to native format.  (rot 1 == deg 360)+rot :: Flt -> Flt+rot !x = x*3.1415926535897*2++-- | Trig with degrees instead of radians.+dcos :: Flt -> Flt+dcos d = cos $ deg d++-- | Force a value to be within a range.  Usage: clamp min x max+clamp :: Flt -> Flt -> Flt -> Flt+clamp !min !x !max+ | x < min = min+ | x > max = max+ | otherwise = x++-- | Tuning parameter.+delta = 0.0001 :: Flt++-- | Non-polymorphic fmin; this speeds+-- things up in ocaml, not sure about haskell.+fmin :: Flt -> Flt -> Flt+fmin !a !b = if a > b then b else a++-- | Non-polymorphic fmax.+fmax :: Flt -> Flt -> Flt+fmax !a !b = if a > b then a else b++-- | Non-polymorphic min of 3 values.+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++-- | Non-polymorphic max of 3 values.+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++-- | Min of 4 values.+fmin4 :: Flt -> Flt -> Flt -> Flt -> Flt+fmin4 !a !b !c !d = fmin (fmin a b) (fmin c d)++-- | Max of 4 values.+fmax4 :: Flt -> Flt -> Flt -> Flt -> Flt+fmax4 !a !b !c !d = fmax (fmax a b) (fmax c d)++-- | Non-polymorphic absolute value.+fabs :: Flt -> Flt+fabs !a = + if a < 0 then (-a) else a++-- | Non-polymorphic integer absolute value.+iabs :: Int -> Int+iabs !a =+ if a < 0 then (-a) else a++-- | Force user to use fabs or iabs, for performance reasons.  Not sure if+-- this really helps, though.+abs a = error "use non-polymorphic version, fabs"++-- | Approximate equality for Flt.  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)++-- | 3d type represented as a record of unboxed floats.+data Vec = Vec !Flt !Flt !Flt deriving Show++-- | A Ray is made up of an origin and direction Vec.+data Ray = Ray {origin, dir :: !Vec} deriving Show+--data Plane = Plane {norm :: !Vec, offset :: !Flt} deriving Show++-- | Vec constructor.+vec :: Flt -> Flt -> Flt -> Vec+vec !x !y !z = (Vec x y z)++-- | Zero Vec.+vzero :: Vec+vzero = Vec 0.0 0.0 0.0++-- | For when we need a unit vector, but we +-- don't care where it points.+vunit :: Vec+vunit = vx++-- | Unit X vector.+vx :: Vec+vx  = Vec 1 0 0++-- | Unit y vector.+vy :: Vec+vy  = Vec 0 1 0++-- | Unit z vector.+vz :: Vec+vz  = Vec 0 0 1++-- | Negative x vector.+nvx :: Vec+nvx = Vec (-1) 0 0++-- | Negative y vector.+nvy :: Vec+nvy = Vec 0 (-1) 0++-- | Negative z vector.+nvz :: Vec+nvz = Vec 0 0 (-1)++-- Extract x coordinate.+x :: Vec -> Flt+x (Vec x_ _ _) = x_++-- Extract y coordinate.+y :: Vec -> Flt+y (Vec _ y_ _) = y_++-- Extract z coordinate.+z :: Vec -> Flt+z (Vec _ _ z_) = z_++-- | Access the Vec as if it were an array indexed from 0..2.+-- Note: this actually accounts for a noticeable amount of cpu +-- time in the Glome ray tracer.+va :: Vec -> Int -> Flt+va !(Vec x y z) !n = + case n of+  0 -> x+  1 -> y+  2 -> z++-- | Create a new Vec with the Nth field overwritten by new value.+-- I could have used record update syntax.+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++-- | Dot product of 2 vectors.  We use this all the time.  Dot product of 2+-- normal vectors is the cosine of the angle between them.+vdot :: Vec -> Vec -> Flt+vdot !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ (x1*x2)+(y1*y2)+(z1*z2)++-- | Cross product of 2 vectors.  Produces a vector perpendicular +-- to the given vectors.  We use this for things like making the forward,+-- up, and right camera vectors orthogonal.  If the input vectors are+-- normalized, the output vector will be as well.+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))++-- | Apply a unary Flt operator to each field of the Vec.+vmap :: (Flt -> Flt) -> Vec -> Vec+vmap f !v1 = + Vec (f (x v1)) (f (y v1)) (f (z v1))++-- | Apply a binary Flt operator to pairs of fields from 2 Vecs.+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))++-- | Reverse the direction of a Vec.+vinvert :: Vec -> Vec+vinvert !(Vec x1 y1 z1) =+ Vec (-x1) (-y1) (-z1)++-- | Get the length of a Vec squared.  We use this to avoid a slow sqrt. +vlensqr :: Vec -> Flt+vlensqr !v1 = vdot v1 v1++-- | Get the length of a Vec.  This is expensive because sqrt is slow.+vlen :: Vec -> Flt+vlen !v1 = sqrt (vdot v1 v1)++-- | Add 2 vectors.+vadd :: Vec -> Vec -> Vec+vadd !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ Vec (x1 + x2)+     (y1 + y2)+     (z1 + z2)++-- | Add 3 vectors.+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)++-- | Subtract vectors.  "vsub b a" is the vector from a to b.+vsub :: Vec -> Vec -> Vec+vsub !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ Vec (x1 - x2)+     (y1 - y2)+     (z1 - z2)++-- | Multiply corresponding fields.  Rarely useful.+vmul :: Vec -> Vec -> Vec+vmul !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ Vec (x1 * x2)+     (y1 * y2)+     (z1 * z2)++-- | Add a value to all the fields of a Vec.  Useful, for instance, to get+-- one corner of the bounding box around a sphere.+vinc :: Vec -> Flt -> Vec+vinc !(Vec x y z) !n =+ Vec (x + n)+     (y + n)+     (z + n)++-- | Subtract a value from all fields of a Vec.+vdec :: Vec -> Flt -> Vec+vdec !(Vec x y z) !n =+ Vec (x - n)+     (y - n)+     (z - n)++-- | Get the maximum of all corresponding fields between 2 Vecs.+vmax :: Vec -> Vec -> Vec+vmax !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ Vec (fmax x1 x2)+     (fmax y1 y2)+     (fmax z1 z2)++-- | Get the minimum of all corresponding fields between 2 Vecs.+vmin :: Vec -> Vec -> Vec+vmin !(Vec x1 y1 z1) !(Vec x2 y2 z2) =+ Vec (fmin x1 x2)+     (fmin y1 y2)+     (fmin z1 z2)++-- | Return the largest axis.  Often used with "va".+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++-- | Scale a Vec by some value.+vscale :: Vec -> Flt -> Vec+vscale !(Vec x y z) !fac =+ Vec (x * fac)+     (y * fac)+     (z * fac)++-- | Take the first Vec, and add to it the second Vec scaled by some amount.+-- This is used quite a lot in Glome.+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 of a Vec just a little shorter.+vnudge :: Vec -> Vec+vnudge x = vscale x (1-delta)++-- | Normalize a vector.  Division is expensive, so we compute the reciprocol +-- of the length and multiply by that.  The sqrt is also expensive.+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)++-- | Throw an exception if a vector hasn't been normalized.+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++-- | Get the victor bisecting two other vectors (which ought to be the same+-- length).+bisect :: Vec -> Vec -> Vec+bisect !v1 !v2 = vnorm (vadd v1 v2)++-- | Distance between 2 vectors.+vdist :: Vec -> Vec -> Flt+vdist v1 v2 = + let d = vsub v2 v1 in vlen d++-- | Reflect a vector "v" off of a surface with normal "norm".+reflect :: Vec -> Vec -> Vec+reflect !v !norm =+  -- vadd v $ vscale norm $ (-2) * (vdot v norm)+  vscaleadd v norm $ (-2) * (vdot v norm)++-- | Reciprocol of all fields of a Vec.+vrcp :: Vec -> Vec+vrcp !(Vec x y z) =+ Vec (1/x) (1/y) (1/z)++-- | Test Vecs for approximate 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)++-- | Test Vecs for matching sign on all fields.  Returns false if any value is+-- zero.  Used by packet tracing.+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 a pair of orthogonal vectors to the one given.+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 "p" and a normal "norm".+-- (Ray does not need to 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++-- | Find the distance along a ray until it intersects with a plane defined+-- by a point "p" and normal "norm".+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 --++-- | 3x4 Transformation matrix.  These are described in most graphics texts.+data Matrix = Matrix !Flt !Flt !Flt !Flt  +                     !Flt !Flt !Flt !Flt  +                     !Flt !Flt !Flt !Flt deriving Show++-- | A transformation.  Inverting a matrix is expensive, so we keep a forward+-- transformation matrix and a reverse transformation matrix.+-- Note: This can be made a little faster if the matricies are non-strict.+data Xfm = Xfm Matrix Matrix deriving Show++-- | Identity matrix.  Transforming a vector by this matrix does nothing.+ident_matrix :: Matrix+ident_matrix = (Matrix 1 0 0 0  0 1 0 0  0 0 1 0)++-- | Identity transformation.+ident_xfm :: Xfm+ident_xfm = Xfm ident_matrix ident_matrix++-- | Multiply two matricies.  This is unrolled for efficiency, and it's also+-- a little bit easier (in my opinion) to see what's going on.+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)++-- | Multiply two tranformations.  This just multiplies the forward and +-- reverse transformations.+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 --++-- | There is a seemingly-magical property of transformation matricies, that+-- we can combine the effects of any number of transformations into a single+-- transformation just by multiplying them together in reverse order.  For +-- instance, we could move a point, then rotate it about the origin by some +-- angle around some vector, then move it again, and this can all be done by +-- a single transformation.+-- This function combines transformations in this way, though it reverses the+-- list first so the transformations take effect in their expected order.+compose :: [Xfm] -> Xfm+compose xfms = check_xfm $ foldr xfm_mult ident_xfm (reverse xfms)++-- | Make sure a transformation is valid.  Multipy the forward and reverse+-- matrix and verify that the result is the identity matrix.+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)) ++-- | Complex transformations: Rotate point (or vector) "pt" about ray by +-- angle c.  The angle is in radians,+-- but using the angle conversion routines "deg", "rad" and "rot" is +-- recommended.+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++-- | Transform a point.  The 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)++-- | Inverse transform a point.+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)++-- | Transform a vector.  The 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)++-- | Inverse transform a vector.+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)++-- | Inverse transform a normal.  This one is tricky: we need to 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)++-- | Transform a Ray.+xfm_ray :: Xfm -> Ray -> Ray+xfm_ray !xfm !(Ray orig dir) =+ Ray (xfm_point xfm orig) (vnorm (xfm_vec xfm dir))++-- | Inverse transform a Ray.+invxfm_ray :: Xfm -> Ray -> Ray+invxfm_ray !xfm !(Ray orig dir) =+ Ray (invxfm_point xfm orig) (vnorm (invxfm_vec xfm dir))++-- BASIC TRANSFORMS --+-- | Basic transforms: move by some displacement vector.+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))++-- | Basic transforms: stretch along the three axes, by the amount+-- in the given vector.  (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)++-- | Basic transforms: rotate about a given axis by some 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)++-- | Basic transforms: Convert coordinate system from canonical xyz +-- 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)++-- | Basic transforms: Convert from uvw coordinates back to normal xyz +-- coordinates.+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 --+-- | Axis-aligned Bounding Box (AABB), defined by opposite corners.  P1 is the+-- min values, p2 has the max values.+data Bbox = Bbox {p1 :: !Vec, p2 :: !Vec} deriving Show++-- | A near-far pair of distances.  Basically just a tuple.+data Interval = Interval !Flt !Flt deriving Show -- used instead of a tuple++-- | Bounding box that encloses two bounding boxes.+bbjoin :: Bbox -> Bbox -> Bbox+bbjoin (Bbox p1a p2a) (Bbox p1b p2b) =+ (Bbox (vmin p1a p1b) (vmax p2a p2b))++-- | Find the overlap of two bounding boxes.+bboverlap :: Bbox -> Bbox -> Bbox+bboverlap (Bbox p1a p2a) (Bbox p1b p2b) =+ (Bbox (vmax p1a p1b) (vmin p2a p2b))++-- | Test if a Vec is inside the bounding box.+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, given an axis and offset.  Throw exception+-- if the offset isn't inside the bounding box.+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 minimum bounding box that encloses 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 of a bounding box.  Useful for cost heuristics when attempting+-- to build optimal bounding box heirarchies.  Undefined for degenerate bounding+-- boxes.+bbsa :: Bbox -> Flt+bbsa (Bbox p1 p2) =+ let Vec dx dy dz = vsub p2 p1 + in dx*dy + dx*dz + dy*dz++-- | Volume of a bounding box.  Undefined for degenerate bounding boxes.+bbvol :: Bbox -> Flt+bbvol (Bbox p1 p2) =+ let (Vec dx dy dz) = vsub p2 p1+ in dx*dy*dz++-- | Degenerate bounding box that contains an empty volume.+empty_bbox :: Bbox+empty_bbox = + Bbox (Vec infinity infinity infinity) +      (Vec (-infinity) (-infinity) (-infinity))++-- | "Infinite" bounding box.+everything_bbox :: Bbox+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)++

@@ -1,126 +0,0 @@-{--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

@@ -1,738 +0,0 @@-{-# OPTIONS_GHC -fexcess-precision #-}-{-# OPTIONS_GHC -funbox-strict-fields #-}-{-# LANGUAGE BangPatterns #-}--module Data.GlomeVec where---- | Performance is pretty similar with Floats or Doubles.--- Todo: make separate Float and Double instances of this library.-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 to native format (noop).-rad :: Flt -> Flt-rad !x = x---- | Convert from rotations to native format.  (rot 1 == deg 360)-rot :: Flt -> Flt-rot !x = x*3.1415926535897*2---- | Trig with degrees instead of radians.-dcos :: Flt -> Flt-dcos d = cos $ deg d---- | Force a value to be within a range.  Usage: clamp min x max-clamp :: Flt -> Flt -> Flt -> Flt-clamp !min !x !max- | x < min = min- | x > max = max- | otherwise = x---- | Tuning parameter.-delta = 0.0001 :: Flt---- | Non-polymorphic fmin; this speeds--- things up in ocaml, not sure about haskell.-fmin :: Flt -> Flt -> Flt-fmin !a !b = if a > b then b else a---- | Non-polymorphic fmax.-fmax :: Flt -> Flt -> Flt-fmax !a !b = if a > b then a else b---- | Non-polymorphic min of 3 values.-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---- | Non-polymorphic max of 3 values.-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---- | Min of 4 values.-fmin4 :: Flt -> Flt -> Flt -> Flt -> Flt-fmin4 !a !b !c !d = fmin (fmin a b) (fmin c d)---- | Max of 4 values.-fmax4 :: Flt -> Flt -> Flt -> Flt -> Flt-fmax4 !a !b !c !d = fmax (fmax a b) (fmax c d)---- | Non-polymorphic absolute value.-fabs :: Flt -> Flt-fabs !a = - if a < 0 then (-a) else a---- | Non-polymorphic integer absolute value.-iabs :: Int -> Int-iabs !a =- if a < 0 then (-a) else a---- | Force user to use fabs or iabs, for performance reasons.  Not sure if--- this really helps, though.-abs a = error "use non-polymorphic version, fabs"---- | Approximate equality for Flt.  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)---- | 3d type represented as a record of unboxed floats.-data Vec = Vec !Flt !Flt !Flt deriving Show---- | A Ray is made up of an origin and direction Vec.-data Ray = Ray {origin, dir :: !Vec} deriving Show---data Plane = Plane {norm :: !Vec, offset :: !Flt} deriving Show---- | Vec constructor.-vec :: Flt -> Flt -> Flt -> Vec-vec !x !y !z = (Vec x y z)---- | Zero Vec.-vzero :: Vec-vzero = Vec 0.0 0.0 0.0---- | For when we need a unit vector, but we --- don't care where it points.-vunit :: Vec-vunit = vx---- | Unit X vector.-vx :: Vec-vx  = Vec 1 0 0---- | Unit y vector.-vy :: Vec-vy  = Vec 0 1 0---- | Unit z vector.-vz :: Vec-vz  = Vec 0 0 1---- | Negative x vector.-nvx :: Vec-nvx = Vec (-1) 0 0---- | Negative y vector.-nvy :: Vec-nvy = Vec 0 (-1) 0---- | Negative z vector.-nvz :: Vec-nvz = Vec 0 0 (-1)---- Extract x coordinate.-x :: Vec -> Flt-x (Vec x_ _ _) = x_---- Extract y coordinate.-y :: Vec -> Flt-y (Vec _ y_ _) = y_---- Extract z coordinate.-z :: Vec -> Flt-z (Vec _ _ z_) = z_---- | Access the Vec as if it were an array indexed from 0..2.--- Note: this actually accounts for a noticeable amount of cpu --- time in the Glome ray tracer.-va :: Vec -> Int -> Flt-va !(Vec x y z) !n = - case n of-  0 -> x-  1 -> y-  2 -> z---- | Create a new Vec with the Nth field overwritten by new value.--- I could have used record update syntax.-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---- | Dot product of 2 vectors.  We use this all the time.  Dot product of 2--- normal vectors is the cosine of the angle between them.-vdot :: Vec -> Vec -> Flt-vdot !(Vec x1 y1 z1) !(Vec x2 y2 z2) =- (x1*x2)+(y1*y2)+(z1*z2)---- | Cross product of 2 vectors.  Produces a vector perpendicular --- to the given vectors.  We use this for things like making the forward,--- up, and right camera vectors orthogonal.  If the input vectors are--- normalized, the output vector will be as well.-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))---- | Apply a unary Flt operator to each field of the Vec.-vmap :: (Flt -> Flt) -> Vec -> Vec-vmap f !v1 = - Vec (f (x v1)) (f (y v1)) (f (z v1))---- | Apply a binary Flt operator to pairs of fields from 2 Vecs.-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))---- | Reverse the direction of a Vec.-vinvert :: Vec -> Vec-vinvert !(Vec x1 y1 z1) =- Vec (-x1) (-y1) (-z1)---- | Get the length of a Vec squared.  We use this to avoid a slow sqrt. -vlensqr :: Vec -> Flt-vlensqr !v1 = vdot v1 v1---- | Get the length of a Vec.  This is expensive because sqrt is slow.-vlen :: Vec -> Flt-vlen !v1 = sqrt (vdot v1 v1)---- | Add 2 vectors.-vadd :: Vec -> Vec -> Vec-vadd !(Vec x1 y1 z1) !(Vec x2 y2 z2) =- Vec (x1 + x2)-     (y1 + y2)-     (z1 + z2)---- | Add 3 vectors.-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)---- | Subtract vectors.  "vsub b a" is the vector from a to b.-vsub :: Vec -> Vec -> Vec-vsub !(Vec x1 y1 z1) !(Vec x2 y2 z2) =- Vec (x1 - x2)-     (y1 - y2)-     (z1 - z2)---- | Multiply corresponding fields.  Rarely useful.-vmul :: Vec -> Vec -> Vec-vmul !(Vec x1 y1 z1) !(Vec x2 y2 z2) =- Vec (x1 * x2)-     (y1 * y2)-     (z1 * z2)---- | Add a value to all the fields of a Vec.  Useful, for instance, to get--- one corner of the bounding box around a sphere.-vinc :: Vec -> Flt -> Vec-vinc !(Vec x y z) !n =- Vec (x + n)-     (y + n)-     (z + n)---- | Subtract a value from all fields of a Vec.-vdec :: Vec -> Flt -> Vec-vdec !(Vec x y z) !n =- Vec (x - n)-     (y - n)-     (z - n)---- | Get the maximum of all corresponding fields between 2 Vecs.-vmax :: Vec -> Vec -> Vec-vmax !(Vec x1 y1 z1) !(Vec x2 y2 z2) =- Vec (fmax x1 x2)-     (fmax y1 y2)-     (fmax z1 z2)---- | Get the minimum of all corresponding fields between 2 Vecs.-vmin :: Vec -> Vec -> Vec-vmin !(Vec x1 y1 z1) !(Vec x2 y2 z2) =- Vec (fmin x1 x2)-     (fmin y1 y2)-     (fmin z1 z2)---- | Return the largest axis.  Often used with "va".-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---- | Scale a Vec by some value.-vscale :: Vec -> Flt -> Vec-vscale !(Vec x y z) !fac =- Vec (x * fac)-     (y * fac)-     (z * fac)---- | Take the first Vec, and add to it the second Vec scaled by some amount.--- This is used quite a lot in Glome.-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 of a Vec just a little shorter.-vnudge :: Vec -> Vec-vnudge x = vscale x (1-delta)---- | Normalize a vector.  Division is expensive, so we compute the reciprocol --- of the length and multiply by that.  The sqrt is also expensive.-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)---- | Throw an exception if a vector hasn't been normalized.-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---- | Get the victor bisecting two other vectors (which ought to be the same--- length).-bisect :: Vec -> Vec -> Vec-bisect !v1 !v2 = vnorm (vadd v1 v2)---- | Distance between 2 vectors.-vdist :: Vec -> Vec -> Flt-vdist v1 v2 = - let d = vsub v2 v1 in vlen d---- | Reflect a vector "v" off of a surface with normal "norm".-reflect :: Vec -> Vec -> Vec-reflect !v !norm =-  -- vadd v $ vscale norm $ (-2) * (vdot v norm)-  vscaleadd v norm $ (-2) * (vdot v norm)---- | Reciprocol of all fields of a Vec.-vrcp :: Vec -> Vec-vrcp !(Vec x y z) =- Vec (1/x) (1/y) (1/z)---- | Test Vecs for approximate 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)---- | Test Vecs for matching sign on all fields.  Returns false if any value is--- zero.  Used by packet tracing.-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 a pair of orthogonal vectors to the one given.-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 "p" and a normal "norm".--- (Ray does not need to 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---- | Find the distance along a ray until it intersects with a plane defined--- by a point "p" and normal "norm".-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 ------ | 3x4 Transformation matrix.  These are described in most graphics texts.-data Matrix = Matrix !Flt !Flt !Flt !Flt  -                     !Flt !Flt !Flt !Flt  -                     !Flt !Flt !Flt !Flt deriving Show---- | A transformation.  Inverting a matrix is expensive, so we keep a forward--- transformation matrix and a reverse transformation matrix.--- Note: This can be made a little faster if the matricies are non-strict.-data Xfm = Xfm Matrix Matrix deriving Show---- | Identity matrix.  Transforming a vector by this matrix does nothing.-ident_matrix :: Matrix-ident_matrix = (Matrix 1 0 0 0  0 1 0 0  0 0 1 0)---- | Identity transformation.-ident_xfm :: Xfm-ident_xfm = Xfm ident_matrix ident_matrix---- | Multiply two matricies.  This is unrolled for efficiency, and it's also--- a little bit easier (in my opinion) to see what's going on.-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)---- | Multiply two tranformations.  This just multiplies the forward and --- reverse transformations.-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 ------ | There is a seemingly-magical property of transformation matricies, that--- we can combine the effects of any number of transformations into a single--- transformation just by multiplying them together in reverse order.  For --- instance, we could move a point, then rotate it about the origin by some --- angle around some vector, then move it again, and this can all be done by --- a single transformation.--- This function combines transformations in this way, though it reverses the--- list first so the transformations take effect in their expected order.-compose :: [Xfm] -> Xfm-compose xfms = check_xfm $ foldr xfm_mult ident_xfm (reverse xfms)---- | Make sure a transformation is valid.  Multipy the forward and reverse--- matrix and verify that the result is the identity matrix.-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)) ---- | Complex transformations: Rotate point (or vector) "pt" about ray by --- angle c.  The angle is in radians,--- but using the angle conversion routines "deg", "rad" and "rot" is --- recommended.-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---- | Transform a point.  The 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)---- | Inverse transform a point.-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)---- | Transform a vector.  The 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)---- | Inverse transform a vector.-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)---- | Inverse transform a normal.  This one is tricky: we need to 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)---- | Transform a Ray.-xfm_ray :: Xfm -> Ray -> Ray-xfm_ray !xfm !(Ray orig dir) =- Ray (xfm_point xfm orig) (vnorm (xfm_vec xfm dir))---- | Inverse transform a Ray.-invxfm_ray :: Xfm -> Ray -> Ray-invxfm_ray !xfm !(Ray orig dir) =- Ray (invxfm_point xfm orig) (vnorm (invxfm_vec xfm dir))---- BASIC TRANSFORMS ----- | Basic transforms: move by some displacement vector.-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))---- | Basic transforms: stretch along the three axes, by the amount--- in the given vector.  (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)---- | Basic transforms: rotate about a given axis by some 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)---- | Basic transforms: Convert coordinate system from canonical xyz --- 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)---- | Basic transforms: Convert from uvw coordinates back to normal xyz --- coordinates.-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 ----- | Axis-aligned Bounding Box (AABB), defined by opposite corners.  P1 is the--- min values, p2 has the max values.-data Bbox = Bbox {p1 :: !Vec, p2 :: !Vec} deriving Show---- | A near-far pair of distances.  Basically just a tuple.-data Interval = Interval !Flt !Flt deriving Show -- used instead of a tuple---- | Bounding box that encloses two bounding boxes.-bbjoin :: Bbox -> Bbox -> Bbox-bbjoin (Bbox p1a p2a) (Bbox p1b p2b) =- (Bbox (vmin p1a p1b) (vmax p2a p2b))---- | Find the overlap of two bounding boxes.-bboverlap :: Bbox -> Bbox -> Bbox-bboverlap (Bbox p1a p2a) (Bbox p1b p2b) =- (Bbox (vmax p1a p1b) (vmin p2a p2b))---- | Test if a Vec is inside the bounding box.-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, given an axis and offset.  Throw exception--- if the offset isn't inside the bounding box.-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 minimum bounding box that encloses 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 of a bounding box.  Useful for cost heuristics when attempting--- to build optimal bounding box heirarchies.  Undefined for degenerate bounding--- boxes.-bbsa :: Bbox -> Flt-bbsa (Bbox p1 p2) =- let Vec dx dy dz = vsub p2 p1 - in dx*dy + dx*dz + dy*dz---- | Volume of a bounding box.  Undefined for degenerate bounding boxes.-bbvol :: Bbox -> Flt-bbvol (Bbox p1 p2) =- let (Vec dx dy dz) = vsub p2 p1- in dx*dy*dz---- | Degenerate bounding box that contains an empty volume.-empty_bbox :: Bbox-empty_bbox = - Bbox (Vec infinity infinity infinity) -      (Vec (-infinity) (-infinity) (-infinity))---- | "Infinite" bounding box.-everything_bbox :: Bbox-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)--

@@ -1,5 +1,5 @@ Name:                GlomeVec-Version:             0.1.1+Version:             0.1.2 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@@ -14,7 +14,7 @@ Cabal-Version: >= 1.2 extra-source-files:   README.txt- library-  exposed-modules:   Data.GlomeVec, Data.GlomeTexture+  exposed-modules:   Data.Glome.Vec+                     Data.Glome.Texture   Build-Depends:     base >= 3 && < 4, array