snoc,
sum
)
-import Numeric.LinearAlgebra hiding (i, scale)
+
import Prelude hiding (LT)
import Test.Framework (Test, testGroup)
import Test.Framework.Providers.HUnit (testCase)
center (Tetrahedron _ v0' v1' v2' v3' _) =
(v0' + v1' + v2' + v3') `scale` (1/4)
- contains_point t p =
+ contains_point t p0 =
b0_unscaled `nearly_ge` 0 &&
b1_unscaled `nearly_ge` 0 &&
b2_unscaled `nearly_ge` 0 &&
-- would do if we used the regular b0,..b3 functions.
b0_unscaled :: Double
b0_unscaled = volume inner_tetrahedron
- where inner_tetrahedron = t { v0 = p }
+ where inner_tetrahedron = t { v0 = p0 }
b1_unscaled :: Double
b1_unscaled = volume inner_tetrahedron
- where inner_tetrahedron = t { v1 = p }
+ where inner_tetrahedron = t { v1 = p0 }
b2_unscaled :: Double
b2_unscaled = volume inner_tetrahedron
- where inner_tetrahedron = t { v2 = p }
+ where inner_tetrahedron = t { v2 = p0 }
b3_unscaled :: Double
b3_unscaled = volume inner_tetrahedron
- where inner_tetrahedron = t { v3 = p }
+ where inner_tetrahedron = t { v3 = p0 }
polynomial :: Tetrahedron -> (RealFunction Point)
--- | The matrix used in the tetrahedron volume calculation as given in
--- Lai & Schumaker, Definition 15.4, page 436.
-vol_matrix :: Tetrahedron -> Matrix Double
-vol_matrix t = (4><4)
- [1, 1, 1, 1,
- x1, x2, x3, x4,
- y1, y2, y3, y4,
- z1, z2, z3, z4 ]
- where
- (x1, y1, z1) = v0 t
- (x2, y2, z2) = v1 t
- (x3, y3, z3) = v2 t
- (x4, y4, z4) = v3 t
+-- | Compute the determinant of the 4x4 matrix,
+--
+-- [1]
+-- [x]
+-- [y]
+-- [z]
+--
+-- where [1] = [1, 1, 1, 1],
+-- [x] = [x1,x2,x3,x4],
+--
+-- et cetera.
+--
+-- The termX nonsense is an attempt to prevent Double overflow.
+-- which has been observed to happen with large coordinates.
+--
+det :: Point -> Point -> Point -> Point -> Double
+det p0 p1 p2 p3 =
+ term5 + term6
+ where
+ (x1, y1, z1) = p0
+ (x2, y2, z2) = p1
+ (x3, y3, z3) = p2
+ (x4, y4, z4) = p3
+ term1 = ((x2 - x4)*y1 - (x1 - x4)*y2 + (x1 - x2)*y4)*z3
+ term2 = ((x2 - x3)*y1 - (x1 - x3)*y2 + (x1 - x2)*y3)*z4
+ term3 = ((x3 - x4)*y2 - (x2 - x4)*y3 + (x2 - x3)*y4)*z1
+ term4 = ((x3 - x4)*y1 - (x1 - x4)*y3 + (x1 - x3)*y4)*z2
+ term5 = term1 - term2
+ term6 = term3 - term4
+
-- | Computed using the formula from Lai & Schumaker, Definition 15.4,
-- page 436.
volume :: Tetrahedron -> Double
volume t
- | (v0 t) == (v1 t) = 0
- | (v0 t) == (v2 t) = 0
- | (v0 t) == (v3 t) = 0
- | (v1 t) == (v2 t) = 0
- | (v1 t) == (v3 t) = 0
- | (v2 t) == (v3 t) = 0
- | otherwise = (1/6)*(det (vol_matrix t))
+ | v0' == v1' = 0
+ | v0' == v2' = 0
+ | v0' == v3' = 0
+ | v1' == v2' = 0
+ | v1' == v3' = 0
+ | v2' == v3' = 0
+ | otherwise = (1/6)*(det v0' v1' v2' v3')
+ where
+ v0' = v0 t
+ v1' = v1 t
+ v2' = v2 t
+ v3' = v3 t
-- | The barycentric coordinates of a point with respect to v0.