X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FCube.hs;h=b0b153e782e3fd6ab98b4856d7f0cf589ee60c40;hb=15ccf02b51c72a7bb6eb098c157005a78cd3f84d;hp=6e31f1911cdef2ed70fe7fcdf836d0bed541c3d8;hpb=e5151f5050e80027f69640813db618aaea54946e;p=spline3.git

diff --git a/src/Cube.hs b/src/Cube.hs
index 6e31f19..b0b153e 100644
--- a/src/Cube.hs
+++ b/src/Cube.hs
@@ -1,14 +1,23 @@
 module Cube
 where
 
-import Data.List ( (\\) )
+import Data.Maybe (fromJust)
+import qualified Data.Vector as V (
+  Vector,
+  findIndex,
+  map,
+  minimum,
+  singleton,
+  snoc,
+  unsafeIndex
+  )
 import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
 
 import Cardinal
 import qualified Face (Face(Face, v0, v1, v2, v3))
 import FunctionValues
 import Point
-import Tetrahedron hiding (c)
+import Tetrahedron hiding (c, fv)
 import ThreeDimensional
 
 data Cube = Cube { h :: Double,
@@ -202,7 +211,7 @@ right_face c = Face.Face v0' v1' v2' v3'
 tetrahedron :: Cube -> Int -> Tetrahedron
 
 tetrahedron c 0 =
-    Tetrahedron (Cube.fv c) v0' v1' v2' v3' vol 0
+    Tetrahedron (Cube.fv c) v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (front_face c)
@@ -211,7 +220,7 @@ tetrahedron c 0 =
       vol = tetrahedra_volume c
 
 tetrahedron c 1 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 1
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (front_face c)
@@ -221,7 +230,7 @@ tetrahedron c 1 =
       vol = tetrahedra_volume c
 
 tetrahedron c 2 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 2
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (front_face c)
@@ -231,7 +240,7 @@ tetrahedron c 2 =
       vol = tetrahedra_volume c
 
 tetrahedron c 3 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 3
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (front_face c)
@@ -241,7 +250,7 @@ tetrahedron c 3 =
       vol = tetrahedra_volume c
 
 tetrahedron c 4 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 4
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (top_face c)
@@ -251,17 +260,17 @@ tetrahedron c 4 =
       vol = tetrahedra_volume c
 
 tetrahedron c 5 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 5
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (top_face c)
       v2' = Face.v1 (top_face c)
       v3' = Face.v2 (top_face c)
-      fv' = rotate cwy $ rotate cwz $ Tetrahedron.fv (tetrahedron c 0)
+      fv' = rotate cwy $ rotate cwz $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 6 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 6
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (top_face c)
@@ -269,31 +278,31 @@ tetrahedron c 6 =
       v3' = Face.v3 (top_face c)
       fv' = rotate cwy $ rotate cwz
                        $ rotate cwz
-                       $ Tetrahedron.fv (tetrahedron c 0)
+                       $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 7 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 7
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (top_face c)
       v2' = Face.v3 (top_face c)
       v3' = Face.v0 (top_face c)
-      fv' = rotate cwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron c 0)
+      fv' = rotate cwy $ rotate ccwz $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 8 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 8
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (back_face c)
       v2' = Face.v0 (back_face c)
       v3' = Face.v1 (back_face c)
-      fv' = rotate cwy $ rotate cwy $ Tetrahedron.fv (tetrahedron c 0)
+      fv' = rotate cwy $ rotate cwy $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 9 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 9
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (back_face c)
@@ -301,11 +310,11 @@ tetrahedron c 9 =
       v3' = Face.v2 (back_face c)
       fv' = rotate cwy $ rotate cwy
                        $ rotate cwx
-                       $ Tetrahedron.fv (tetrahedron c 0)
+                       $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 10 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 10
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (back_face c)
@@ -314,12 +323,12 @@ tetrahedron c 10 =
       fv' = rotate cwy $ rotate cwy
                        $ rotate cwx
                        $ rotate cwx
-                       $ Tetrahedron.fv (tetrahedron c 0)
+                       $ fv c
 
       vol = tetrahedra_volume c
 
 tetrahedron c 11 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 11
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (back_face c)
@@ -327,31 +336,31 @@ tetrahedron c 11 =
       v3' = Face.v0 (back_face c)
       fv' = rotate cwy $ rotate cwy
                        $ rotate ccwx
-                       $ Tetrahedron.fv (tetrahedron c 0)
+                       $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 12 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 12
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (down_face c)
       v2' = Face.v0 (down_face c)
       v3' = Face.v1 (down_face c)
-      fv' = rotate ccwy (Tetrahedron.fv (tetrahedron c 0))
+      fv' = rotate ccwy $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 13 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 13
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (down_face c)
       v2' = Face.v1 (down_face c)
       v3' = Face.v2 (down_face c)
-      fv' = rotate ccwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron c 0)
+      fv' = rotate ccwy $ rotate ccwz $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 14 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 14
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (down_face c)
@@ -359,41 +368,41 @@ tetrahedron c 14 =
       v3' = Face.v3 (down_face c)
       fv' = rotate ccwy $ rotate ccwz
                         $ rotate ccwz
-                        $ Tetrahedron.fv (tetrahedron c 0)
+                        $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 15 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 15
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (down_face c)
       v2' = Face.v3 (down_face c)
       v3' = Face.v0 (down_face c)
-      fv' = rotate ccwy $ rotate cwz $ Tetrahedron.fv (tetrahedron c 0)
+      fv' = rotate ccwy $ rotate cwz $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 16 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 16
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (right_face c)
       v2' = Face.v0 (right_face c)
       v3' = Face.v1 (right_face c)
-      fv' = rotate ccwz (Tetrahedron.fv (tetrahedron c 0))
+      fv' = rotate ccwz $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 17 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 17
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (right_face c)
       v2' = Face.v1 (right_face c)
       v3' = Face.v2 (right_face c)
-      fv' = rotate ccwz $ rotate cwy $ Tetrahedron.fv (tetrahedron c 0)
+      fv' = rotate ccwz $ rotate cwy $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 18 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 18
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (right_face c)
@@ -401,42 +410,42 @@ tetrahedron c 18 =
       v3' = Face.v3 (right_face c)
       fv' = rotate ccwz $ rotate cwy
                         $ rotate cwy
-                        $ Tetrahedron.fv (tetrahedron c 0)
+                        $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 19 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 19
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (right_face c)
       v2' = Face.v3 (right_face c)
       v3' = Face.v0 (right_face c)
       fv' = rotate ccwz $ rotate ccwy
-                        $ Tetrahedron.fv (tetrahedron c 0)
+                        $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 20 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 20
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (left_face c)
       v2' = Face.v0 (left_face c)
       v3' = Face.v1 (left_face c)
-      fv' = rotate cwz (Tetrahedron.fv (tetrahedron c 0))
+      fv' = rotate cwz $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 21 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 21
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (left_face c)
       v2' = Face.v1 (left_face c)
       v3' = Face.v2 (left_face c)
-      fv' = rotate cwz $ rotate ccwy $ Tetrahedron.fv (tetrahedron c 0)
+      fv' = rotate cwz $ rotate ccwy $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 22 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 22
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (left_face c)
@@ -444,57 +453,100 @@ tetrahedron c 22 =
       v3' = Face.v3 (left_face c)
       fv' = rotate cwz $ rotate ccwy
                        $ rotate ccwy
-                       $ Tetrahedron.fv (tetrahedron c 0)
+                       $ fv c
       vol = tetrahedra_volume c
 
 tetrahedron c 23 =
-    Tetrahedron fv' v0' v1' v2' v3' vol 23
+    Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
       v1' = center (left_face c)
       v2' = Face.v3 (left_face c)
       v3' = Face.v0 (left_face c)
       fv' = rotate cwz $ rotate cwy
-                       $ Tetrahedron.fv (tetrahedron c 0)
+                       $ fv c
       vol = tetrahedra_volume c
 
 -- Feels dirty, but whatever.
 tetrahedron _ _ = error "asked for a nonexistent tetrahedron"
 
 
+-- Only used in tests, so we don't need the added speed
+-- of Data.Vector.
 tetrahedra :: Cube -> [Tetrahedron]
-tetrahedra c =
-    [ tetrahedron c n | n <- [0..23] ]
-
--- | All completely contained in the front half of the cube.
-front_half_tetrahedra :: Cube -> [Tetrahedron]
-front_half_tetrahedra c =
-  [ tetrahedron c n | n <- [0,1,2,3,6,12,19,21] ]
-
--- | All tetrahedra completely contained in the top half of the cube.
-top_half_tetrahedra :: Cube -> [Tetrahedron]
-top_half_tetrahedra c =
-  [ tetrahedron c n | n <- [4,5,6,7,0,10,16,20] ]
-
--- | All tetrahedra completely contained in the back half of the cube.
-back_half_tetrahedra :: Cube -> [Tetrahedron]
-back_half_tetrahedra c =
-  [ tetrahedron c n | n <- [8,9,10,11,4,14,17,23] ]
-
--- | All tetrahedra completely contained in the down half of the cube.
-down_half_tetrahedra :: Cube -> [Tetrahedron]
-down_half_tetrahedra c =
-  [ tetrahedron c n | n <- [12,13,14,15,2,8,18,22] ]
-
--- | All tetrahedra completely contained in the right half of the cube.
-right_half_tetrahedra :: Cube -> [Tetrahedron]
-right_half_tetrahedra c =
-  [ tetrahedron c n | n <- [16,17,18,19,1,5,9,13] ]
-
--- | All tetrahedra completely contained in the left half of the cube.
-left_half_tetrahedra :: Cube -> [Tetrahedron]
-left_half_tetrahedra c =
-  [ tetrahedron c n | n <- [20,21,22,23,3,7,11,15] ]
+tetrahedra c = [ tetrahedron c n | n <- [0..23] ]
+
+front_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
+front_left_top_tetrahedra  c =
+  V.singleton (tetrahedron c 0) `V.snoc`
+    (tetrahedron c 3) `V.snoc`
+    (tetrahedron c 6) `V.snoc`
+    (tetrahedron c 7) `V.snoc`
+    (tetrahedron c 20) `V.snoc`
+    (tetrahedron c 21)
+
+front_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
+front_left_down_tetrahedra  c =
+  V.singleton (tetrahedron c 0) `V.snoc`
+    (tetrahedron c 2) `V.snoc`
+    (tetrahedron c 3) `V.snoc`
+    (tetrahedron c 12) `V.snoc`
+    (tetrahedron c 15) `V.snoc`
+    (tetrahedron c 21)
+
+front_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
+front_right_top_tetrahedra  c =
+  V.singleton (tetrahedron c 0) `V.snoc`
+    (tetrahedron c 1) `V.snoc`
+    (tetrahedron c 5) `V.snoc`
+    (tetrahedron c 6) `V.snoc`
+    (tetrahedron c 16) `V.snoc`
+    (tetrahedron c 19)
+
+front_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
+front_right_down_tetrahedra  c =
+  V.singleton (tetrahedron c 1) `V.snoc`
+    (tetrahedron c 2) `V.snoc`
+    (tetrahedron c 12) `V.snoc`
+    (tetrahedron c 13) `V.snoc`
+    (tetrahedron c 18) `V.snoc`
+    (tetrahedron c 19)
+
+back_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
+back_left_top_tetrahedra  c =
+  V.singleton (tetrahedron c 0) `V.snoc`
+    (tetrahedron c 3) `V.snoc`
+    (tetrahedron c 6) `V.snoc`
+    (tetrahedron c 7) `V.snoc`
+    (tetrahedron c 20) `V.snoc`
+    (tetrahedron c 21)
+
+back_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
+back_left_down_tetrahedra  c =
+  V.singleton (tetrahedron c 8) `V.snoc`
+    (tetrahedron c 11) `V.snoc`
+    (tetrahedron c 14) `V.snoc`
+    (tetrahedron c 15) `V.snoc`
+    (tetrahedron c 22) `V.snoc`
+    (tetrahedron c 23)
+
+back_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
+back_right_top_tetrahedra  c =
+  V.singleton (tetrahedron c 4) `V.snoc`
+    (tetrahedron c 5) `V.snoc`
+    (tetrahedron c 9) `V.snoc`
+    (tetrahedron c 10) `V.snoc`
+    (tetrahedron c 16) `V.snoc`
+    (tetrahedron c 17)
+
+back_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
+back_right_down_tetrahedra  c =
+  V.singleton (tetrahedron c 8) `V.snoc`
+    (tetrahedron c 9) `V.snoc`
+    (tetrahedron c 13) `V.snoc`
+    (tetrahedron c 14) `V.snoc`
+    (tetrahedron c 17) `V.snoc`
+    (tetrahedron c 18)
 
 in_top_half :: Cube -> Point -> Bool
 in_top_half c (_,_,z) =
@@ -531,33 +583,43 @@ in_left_half c (_,y,_) =
 --
 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
 find_containing_tetrahedron c p =
-  head containing_tetrahedra
+  candidates `V.unsafeIndex` (fromJust lucky_idx)
   where
-    candidates = tetrahedra c
-    non_candidates_x =
-        if (in_front_half c p) then
-          back_half_tetrahedra c
+    front_half = in_front_half c p
+    top_half = in_top_half c p
+    left_half = in_left_half c p
+
+    candidates =
+      if front_half then
+
+        if left_half then
+          if top_half then
+            front_left_top_tetrahedra c
+          else
+            front_left_down_tetrahedra c
         else
-          front_half_tetrahedra c
-
-    candidates_x = candidates \\ non_candidates_x
-
-    non_candidates_y =
-      if (in_left_half c p) then
-        right_half_tetrahedra c
-      else
-        left_half_tetrahedra c
-
-    candidates_xy = candidates_x \\ non_candidates_y
-
-    non_candidates_z =
-      if (in_top_half c p) then
-        down_half_tetrahedra c
-      else
-        top_half_tetrahedra c
-
-    candidates_xyz = candidates_xy \\ non_candidates_z
-
-    contains_our_point = flip contains_point p
-    containing_tetrahedra = filter contains_our_point candidates_xyz
-
+          if top_half then
+            front_right_top_tetrahedra c
+          else
+            front_right_down_tetrahedra c
+
+      else -- bottom half
+
+        if left_half then
+          if top_half then
+            back_left_top_tetrahedra c
+          else
+            back_left_down_tetrahedra c
+        else
+          if top_half then
+            back_right_top_tetrahedra c
+          else
+            back_right_down_tetrahedra c
+
+    -- Use the dot product instead of 'distance' here to save a
+    -- sqrt(). So, "distances" below really means "distances squared."
+    distances = V.map ((dot p) . center) candidates
+    shortest_distance = V.minimum distances
+    lucky_idx = V.findIndex
+                  (\t -> (center t) `dot` p == shortest_distance)
+                  candidates