GlomeVec-0.2: Data/Glome/Vec.hs
{-# 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
ray_ub :: Ray -> (# Flt, Flt, Flt, Flt, Flt, Flt #)
ray_ub (Ray (Vec !ox !oy !oz) (Vec !dx !dy !dz)) = (# ox, oy, oz, dx, dy, dz #)
-- | 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)
-}
bbclip_ub :: Ray -> Bbox -> (# Flt, Flt #)
bbclip_ub (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
(# !inx, !outx #) =
if dx > 0
then (# (p1x-ox)*dxrcp, (p2x-ox)*dxrcp #)
else (# (p2x-ox)*dxrcp, (p1x-ox)*dxrcp #)
(# !iny, !outy #) =
if dy > 0
then (# (p1y-oy)*dyrcp, (p2y-oy)*dyrcp #)
else (# (p2y-oy)*dyrcp, (p1y-oy)*dyrcp #)
(# !inz, !outz #) =
if dz > 0
then (# (p1z-oz)*dzrcp, (p2z-oz)*dzrcp #)
else (# (p2z-oz)*dzrcp, (p1z-oz)*dzrcp #)
in
(# fmax3 inx iny inz, fmin3 outx outy outz #)
bbclip :: Ray -> Bbox -> Interval
bbclip ray bbox =
let (# near, far #) = bbclip_ub ray bbox
in Interval near far