Begin fixing value_at outside of the grid.

[spline3.git] / src / FunctionValues.hs

diff --git a/src/FunctionValues.hs b/src/FunctionValues.hs
index 00bb0a828d3072f3d4d4b2cdc51c8c4e2ccd2549..9d5232334bc5ad805f449faa3d51249644689b4e 100644 (file)
--- a/src/FunctionValues.hs
+++ b/src/FunctionValues.hs
@@ -7,6 +7,7 @@ module FunctionValues (
   make_values,
   rotate,
   function_values_tests,
+  function_values_properties,
   value_at
   )
 where
@@ -15,10 +16,11 @@ import Prelude hiding (LT)
 import Test.HUnit
 import Test.Framework (Test, testGroup)
 import Test.Framework.Providers.HUnit (testCase)
+import Test.Framework.Providers.QuickCheck2 (testProperty)
 import Test.QuickCheck (Arbitrary(..), choose)
 
 import Assertions (assertTrue)
-import Cardinal ( Cardinal(..) )
+import Cardinal ( Cardinal(..), cwx, cwy, cwz )
 import Examples (trilinear)
 import Values (Values3D, dims, idx)
 
@@ -175,8 +177,11 @@ eval f (Quotient x y) = (eval f x) / (eval f y)
 
 -- | Takes a three-dimensional list of 'Double' and a set of 3D
 --   coordinates (i,j,k), and returns the value at (i,j,k) in the
---   supplied list. If there is no such value, we choose a nearby
---   point and use its value.
+--   supplied list. If there is no such value, we calculate one
+--   according to Sorokina and Zeilfelder, remark 7.3, p. 99.
+--
+--   We specifically do not consider values more than one unit away
+--   from our grid.
 --
 --   Examples:
 --
@@ -184,7 +189,7 @@ eval f (Quotient x y) = (eval f x) / (eval f y)
 --   1.0
 --
 --   >>> value_at Examples.trilinear (-1) 0 0
---   1.0
+--   0.0
 --
 --   >>> value_at Examples.trilinear 0 0 4
 --   1.0
@@ -194,13 +199,55 @@ eval f (Quotient x y) = (eval f x) / (eval f y)
 --
 value_at :: Values3D -> Int -> Int -> Int -> Double
 value_at v3d i j k
-         | i < 0 = value_at v3d 0 j k
-         | j < 0 = value_at v3d i 0 k
-         | k < 0 = value_at v3d i j 0
-         | xsize <= i = value_at v3d (xsize - 1) j k
-         | ysize <= j = value_at v3d i (ysize - 1) k
-         | zsize <= k = value_at v3d i j (zsize - 1)
-         | otherwise = idx v3d i j k
+  -- Put the most common case first!
+  | (i >= 0) && (j >= 0) && (k >= 0) =
+      idx v3d i j k
+
+  -- The next three are from the first line in (7.3).
+  | (i == -1) && (j >= 0) && (k >= 0) =
+      2*(value_at v3d 0 j k) - (value_at v3d 1 j k)
+
+  | (i >= 0) && (j == -1) && (k >= 0) =
+      2*(value_at v3d i 0 k) - (value_at v3d i 1 k)
+
+  | (i >= 0) && (j >= 0) && (k == -1) =
+      2*(value_at v3d i j 0) - (value_at v3d i j 1)
+
+  -- The next two are from the second line in (7.3).
+  | (i == -1) && (j == -1) && (k >= 0) =
+      2*(value_at v3d i 0 k) - (value_at v3d i 1 k)
+
+  | (i == -1) && (j == ysize) && (k >= 0) =
+      2*(value_at v3d i (ysize - 1) k) - (value_at v3d i (ysize - 2) k)
+
+  -- The next two are from the third line in (7.3).
+  | (i == -1) && (j >= 0) && (k == -1) =
+      2*(value_at v3d i j 0) - (value_at v3d i j 1)
+
+  | (i == -1) && (j >= 0) && (k == zsize) =
+      2*(value_at v3d i j (zsize - 1)) - (value_at v3d i j (zsize - 2))
+
+  -- Repeat the above (j and k) cases for i >= 0.
+  | (i >= 0) && (j == -1) && (k == -1) =
+      2*(value_at v3d i j 0) - (value_at v3d i j 1)
+
+  | (i == xsize) && (j == -1) && (k >= 0) =
+      2*(value_at v3d (xsize - 1) j k) - (value_at v3d (xsize - 2) j k)
+
+  -- These two cases I made up.
+  | (i == -1) && (j == -1) && (k == -1) =
+      2*(value_at v3d i j 0) - (value_at v3d i j 1)
+
+  | (i == xsize) && (j == ysize) && (k == zsize) =
+      2*(value_at v3d i j (zsize - 1)) - (value_at v3d i j (zsize - 2))
+
+  | otherwise =
+      let istr = show i
+          jstr = show j
+          kstr = show k
+          coordstr = "(" ++ istr ++ "," ++ jstr ++ "," ++ kstr ++ ")"
+      in
+        error $ "value_at called outside of domain: " ++ coordstr
   where
     (xsize, ysize, zsize) = dims v3d
 
@@ -314,3 +361,68 @@ function_values_tests :: Test.Framework.Test
 function_values_tests =
     testGroup "FunctionValues Tests"
                   [ testCase "test directions" test_directions ]
+
+
+prop_x_rotation_doesnt_affect_front :: FunctionValues -> Bool
+prop_x_rotation_doesnt_affect_front fv0 =
+    expr1 == expr2
+    where
+      fv1 = rotate cwx fv0
+      expr1 = front fv0
+      expr2 = front fv1
+
+prop_x_rotation_doesnt_affect_back :: FunctionValues -> Bool
+prop_x_rotation_doesnt_affect_back fv0 =
+    expr1 == expr2
+    where
+      fv1 = rotate cwx fv0
+      expr1 = back fv0
+      expr2 = back fv1
+
+
+prop_y_rotation_doesnt_affect_left :: FunctionValues -> Bool
+prop_y_rotation_doesnt_affect_left fv0 =
+    expr1 == expr2
+    where
+      fv1 = rotate cwy fv0
+      expr1 = left fv0
+      expr2 = left fv1
+
+prop_y_rotation_doesnt_affect_right :: FunctionValues -> Bool
+prop_y_rotation_doesnt_affect_right fv0 =
+    expr1 == expr2
+    where
+      fv1 = rotate cwy fv0
+      expr1 = right fv0
+      expr2 = right fv1
+
+
+prop_z_rotation_doesnt_affect_down :: FunctionValues -> Bool
+prop_z_rotation_doesnt_affect_down fv0 =
+    expr1 == expr2
+    where
+      fv1 = rotate cwz fv0
+      expr1 = down fv0
+      expr2 = down fv1
+
+
+prop_z_rotation_doesnt_affect_top :: FunctionValues -> Bool
+prop_z_rotation_doesnt_affect_top fv0 =
+    expr1 == expr2
+    where
+      fv1 = rotate cwz fv0
+      expr1 = top fv0
+      expr2 = top fv1
+
+
+function_values_properties :: Test.Framework.Test
+function_values_properties =
+  let tp = testProperty
+  in
+    testGroup "FunctionValues Properties" [
+      tp "x rotation doesn't affect front" prop_x_rotation_doesnt_affect_front,
+      tp "x rotation doesn't affect back" prop_x_rotation_doesnt_affect_back,
+      tp "y rotation doesn't affect left" prop_y_rotation_doesnt_affect_left,
+      tp "y rotation doesn't affect right" prop_y_rotation_doesnt_affect_right,
+      tp "z rotation doesn't affect top" prop_z_rotation_doesnt_affect_top,
+      tp "z rotation doesn't affect down" prop_z_rotation_doesnt_affect_down ]