X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FTetrahedron.hs;h=2342ba139d2b79cf64c2b5dbac60c3fa0de55b3a;hb=279936ce0f59cbc2b5c1c3c748ae9fae03ee1146;hp=b8d4fe67943ef30db4b14f00f87c7bc5bee96b2c;hpb=3a954903101eca7594a65824868517b9758e188d;p=spline3.git

diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs
index b8d4fe6..2342ba1 100644
--- a/src/Tetrahedron.hs
+++ b/src/Tetrahedron.hs
@@ -17,7 +17,7 @@ import qualified Data.Vector as V (
   snoc,
   sum
   )
-import Numeric.LinearAlgebra hiding (i, scale)
+
 import Prelude hiding (LT)
 import Test.Framework (Test, testGroup)
 import Test.Framework.Providers.HUnit (testCase)
@@ -72,7 +72,7 @@ instance ThreeDimensional Tetrahedron where
     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 &&
@@ -82,19 +82,19 @@ instance ThreeDimensional Tetrahedron where
         -- 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)
@@ -281,31 +281,47 @@ c _ _ _ _ _ = error "coefficient index out of bounds"
 
 
 
--- | 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.
+--
+det :: Point -> Point -> Point -> Point -> Double
+det p0 p1 p2 p3 =
+  x1*y2*z4 - x1*y2*z3 + x1*y3*z2 - x1*y3*z4 - x1*y4*z2 + x1*y4*z3 +
+  x2*y1*z3 - x2*y1*z4 - x2*y3*z1 + x2*y3*z4 + x2*y4*z1 + x3*y1*z4 +
+  x3*y2*z1 - x3*y2*z4 - x3*y4*z1 - x2*y4*z3 - x3*y1*z2 + x3*y4*z2 +
+  x4*y1*z2 - x4*y1*z3 - x4*y2*z1 + x4*y2*z3 + x4*y3*z1 - x4*y3*z2
+  where
+    (x1, y1, z1) = p0
+    (x2, y2, z2) = p1
+    (x3, y3, z3) = p2
+    (x4, y4, z4) = p3
+
 
 -- | 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.