X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FCube.hs;h=d863c290f2084ef8c79c0dc944dbb7c0c5c6c191;hb=5f01596d42cca3ec2b8236d697adb468cfcdb055;hp=9f931431b93bafb4ca6e0a13c0940bc036a3e01a;hpb=1bf996325008f79215a607d765adb042026f7444;p=spline3.git

diff --git a/src/Cube.hs b/src/Cube.hs
index 9f93143..d863c29 100644
--- a/src/Cube.hs
+++ b/src/Cube.hs
@@ -24,12 +24,11 @@ import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
 
 import Cardinal
 import Comparisons ((~=), (~~=))
-import qualified Face (Face(Face, v0, v1, v2, v3))
+import qualified Face (Face(..), center)
 import FunctionValues (FunctionValues, eval, rotate)
 import Misc (all_equal, disjoint)
 import Point (Point(..), dot)
-import Tetrahedron (Tetrahedron(..), c, volume)
-import ThreeDimensional
+import Tetrahedron (Tetrahedron(..), barycenter, c, volume)
 
 data Cube = Cube { h  :: !Double,
                    i  :: !Int,
@@ -122,30 +121,20 @@ zmax cube = (k' + 1/2)*delta
       k' = fromIntegral (k cube) :: Double
       delta = h cube
 
-instance ThreeDimensional Cube where
-    -- | The center of Cube_ijk coincides with v_ijk at
-    --   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
-    center cube = Point x y z
-           where
-             delta = h cube
-             i' = fromIntegral (i cube) :: Double
-             j' = fromIntegral (j cube) :: Double
-             k' = fromIntegral (k cube) :: Double
-             x = delta * i'
-             y = delta * j'
-             z = delta * k'
-
-    -- | It's easy to tell if a point is within a cube; just make sure
-    --   that it falls on the proper side of each of the cube's faces.
-    contains_point cube (Point x y z)
-        | x < (xmin cube) = False
-        | x > (xmax cube) = False
-        | y < (ymin cube) = False
-        | y > (ymax cube) = False
-        | z < (zmin cube) = False
-        | z > (zmax cube) = False
-        | otherwise = True
 
+-- | The center of Cube_ijk coincides with v_ijk at
+--   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
+center :: Cube -> Point
+center cube =
+  Point x y z
+  where
+    delta = h cube
+    i' = fromIntegral (i cube) :: Double
+    j' = fromIntegral (j cube) :: Double
+    k' = fromIntegral (k cube) :: Double
+    x = delta * i'
+    y = delta * j'
+    z = delta * k'
 
 
 -- Face stuff.
@@ -230,7 +219,7 @@ tetrahedron cube 0 =
     where
       v0' = center cube
       ff  = front_face cube
-      v1' = center ff
+      v1' = Face.center ff
       v2' = Face.v0 ff
       v3' = Face.v1 ff
       vol = tetrahedra_volume cube
@@ -240,7 +229,7 @@ tetrahedron cube 1 =
     where
       v0' = center cube
       ff  = front_face cube
-      v1' = center ff
+      v1' = Face.center ff
       v2' = Face.v1 ff
       v3' = Face.v2 ff
       fv' = rotate ccwx (fv cube)
@@ -251,7 +240,7 @@ tetrahedron cube 2 =
     where
       v0' = center cube
       ff  = front_face cube
-      v1' = center ff
+      v1' = Face.center ff
       v2' = Face.v2 ff
       v3' = Face.v3 ff
       fv' = rotate ccwx $ rotate ccwx $ fv cube
@@ -262,7 +251,7 @@ tetrahedron cube 3 =
     where
       v0' = center cube
       ff  = front_face cube
-      v1' = center ff
+      v1' = Face.center ff
       v2' = Face.v3 ff
       v3' = Face.v0 ff
       fv' = rotate cwx (fv cube)
@@ -273,7 +262,7 @@ tetrahedron cube 4 =
     where
       v0' = center cube
       tf  = top_face cube
-      v1' = center tf
+      v1' = Face.center tf
       v2' = Face.v0 tf
       v3' = Face.v1 tf
       fv' = rotate cwy (fv cube)
@@ -284,7 +273,7 @@ tetrahedron cube 5 =
     where
       v0' = center cube
       tf  = top_face cube
-      v1' = center tf
+      v1' = Face.center tf
       v2' = Face.v1 tf
       v3' = Face.v2 tf
       fv' = rotate cwy $ rotate cwz $ fv cube
@@ -295,7 +284,7 @@ tetrahedron cube 6 =
     where
       v0' = center cube
       tf  = top_face cube
-      v1' = center tf
+      v1' = Face.center tf
       v2' = Face.v2 tf
       v3' = Face.v3 tf
       fv' = rotate cwy $ rotate cwz
@@ -308,7 +297,7 @@ tetrahedron cube 7 =
     where
       v0' = center cube
       tf  = top_face cube
-      v1' = center tf
+      v1' = Face.center tf
       v2' = Face.v3 tf
       v3' = Face.v0 tf
       fv' = rotate cwy $ rotate ccwz $ fv cube
@@ -319,7 +308,7 @@ tetrahedron cube 8 =
     where
       v0' = center cube
       bf  = back_face cube
-      v1' = center bf
+      v1' = Face.center bf
       v2' = Face.v0 bf
       v3' = Face.v1 bf
       fv' = rotate cwy $ rotate cwy $ fv cube
@@ -330,7 +319,7 @@ tetrahedron cube 9 =
     where
       v0' = center cube
       bf  = back_face cube
-      v1' = center bf
+      v1' = Face.center bf
       v2' = Face.v1 bf
       v3' = Face.v2 bf
       fv' = rotate cwy $ rotate cwy
@@ -343,7 +332,7 @@ tetrahedron cube 10 =
     where
       v0' = center cube
       bf  = back_face cube
-      v1' = center bf
+      v1' = Face.center bf
       v2' = Face.v2 bf
       v3' = Face.v3 bf
       fv' = rotate cwy $ rotate cwy
@@ -358,7 +347,7 @@ tetrahedron cube 11 =
     where
       v0' = center cube
       bf  = back_face cube
-      v1' = center bf
+      v1' = Face.center bf
       v2' = Face.v3 bf
       v3' = Face.v0 bf
       fv' = rotate cwy $ rotate cwy
@@ -371,7 +360,7 @@ tetrahedron cube 12 =
     where
       v0' = center cube
       df  = down_face cube
-      v1' = center df
+      v1' = Face.center df
       v2' = Face.v0 df
       v3' = Face.v1 df
       fv' = rotate ccwy $ fv cube
@@ -382,7 +371,7 @@ tetrahedron cube 13 =
     where
       v0' = center cube
       df  = down_face cube
-      v1' = center df
+      v1' = Face.center df
       v2' = Face.v1 df
       v3' = Face.v2 df
       fv' = rotate ccwy $ rotate ccwz $ fv cube
@@ -393,7 +382,7 @@ tetrahedron cube 14 =
     where
       v0' = center cube
       df  = down_face cube
-      v1' = center df
+      v1' = Face.center df
       v2' = Face.v2 df
       v3' = Face.v3 df
       fv' = rotate ccwy $ rotate ccwz
@@ -406,7 +395,7 @@ tetrahedron cube 15 =
     where
       v0' = center cube
       df  = down_face cube
-      v1' = center df
+      v1' = Face.center df
       v2' = Face.v3 df
       v3' = Face.v0 df
       fv' = rotate ccwy $ rotate cwz $ fv cube
@@ -417,7 +406,7 @@ tetrahedron cube 16 =
     where
       v0' = center cube
       rf  = right_face cube
-      v1' = center rf
+      v1' = Face.center rf
       v2' = Face.v0 rf
       v3' = Face.v1 rf
       fv' = rotate ccwz $ fv cube
@@ -428,7 +417,7 @@ tetrahedron cube 17 =
     where
       v0' = center cube
       rf  = right_face cube
-      v1' = center rf
+      v1' = Face.center rf
       v2' = Face.v1 rf
       v3' = Face.v2 rf
       fv' = rotate ccwz $ rotate cwy $ fv cube
@@ -439,7 +428,7 @@ tetrahedron cube 18 =
     where
       v0' = center cube
       rf  = right_face cube
-      v1' = center rf
+      v1' = Face.center rf
       v2' = Face.v2 rf
       v3' = Face.v3 rf
       fv' = rotate ccwz $ rotate cwy
@@ -452,7 +441,7 @@ tetrahedron cube 19 =
     where
       v0' = center cube
       rf  = right_face cube
-      v1' = center rf
+      v1' = Face.center rf
       v2' = Face.v3 rf
       v3' = Face.v0 rf
       fv' = rotate ccwz $ rotate ccwy
@@ -464,7 +453,7 @@ tetrahedron cube 20 =
     where
       v0' = center cube
       lf  = left_face cube
-      v1' = center lf
+      v1' = Face.center lf
       v2' = Face.v0 lf
       v3' = Face.v1 lf
       fv' = rotate cwz $ fv cube
@@ -475,7 +464,7 @@ tetrahedron cube 21 =
     where
       v0' = center cube
       lf  = left_face cube
-      v1' = center lf
+      v1' = Face.center lf
       v2' = Face.v1 lf
       v3' = Face.v2 lf
       fv' = rotate cwz $ rotate ccwy $ fv cube
@@ -486,7 +475,7 @@ tetrahedron cube 22 =
     where
       v0' = center cube
       lf  = left_face cube
-      v1' = center lf
+      v1' = Face.center lf
       v2' = Face.v2 lf
       v3' = Face.v3 lf
       fv' = rotate cwz $ rotate ccwy
@@ -499,7 +488,7 @@ tetrahedron cube 23 =
     where
       v0' = center cube
       lf  = left_face cube
-      v1' = center lf
+      v1' = Face.center lf
       v2' = Face.v3 lf
       v3' = Face.v0 lf
       fv' = rotate cwz $ rotate cwy
@@ -626,6 +615,7 @@ find_containing_tetrahedron cube p =
     top_half = in_top_half cube p
     left_half = in_left_half cube p
 
+    candidates :: V.Vector Tetrahedron
     candidates =
       if front_half then
 
@@ -656,10 +646,20 @@ find_containing_tetrahedron cube p =
     -- Use the dot product instead of Euclidean distance here to save
     -- a sqrt(). So, "distances" below really means "distances
     -- squared."
-    distances = V.map ((dot p) . center) candidates
+    distances :: V.Vector Double
+    distances = V.map ((dot p) . barycenter) candidates
+
+    shortest_distance :: Double
     shortest_distance = V.minimum distances
+
+    -- Compute the index of the tetrahedron with the center closest to
+    -- p. This is a bad algorithm, but don't change it! If you make it
+    -- smarter by finding the index of shortest_distance in distances
+    -- (this should give the same answer and avoids recomputing the
+    -- dot product), the program gets slower. Seriously!
+    lucky_idx :: Maybe Int
     lucky_idx = V.findIndex
-                  (\t -> (center t) `dot` p == shortest_distance)
+                  (\t -> (barycenter t) `dot` p == shortest_distance)
                   candidates
 
 
@@ -716,7 +716,8 @@ prop_all_volumes_exact cube =
     where
       delta = h cube
 
--- | All tetrahedron should have their v0 located at the center of the cube.
+-- | All tetrahedron should have their v0 located at the center of the
+--   cube.
 prop_v0_all_equal :: Cube -> Bool
 prop_v0_all_equal cube = (v0 t0) == (v0 t1)
     where