1 {-# LANGUAGE BangPatterns #-}
2 -- | The Grid module contains the Grid type, its tests, and the 'zoom'
3 -- function used to build the interpolation.
12 import qualified Data.Array.Repa as R
13 import qualified Data.Array.Repa.Operators.Traversal as R (unsafeTraverse)
14 import Test.HUnit (Assertion, assertEqual)
15 import Test.Framework (Test, testGroup)
16 import Test.Framework.Providers.HUnit (testCase)
17 import Test.Framework.Providers.QuickCheck2 (testProperty)
18 import Test.QuickCheck ((==>),
23 import Assertions (assertAlmostEqual, assertTrue)
24 import Comparisons ((~=))
25 import Cube (Cube(Cube),
26 find_containing_tetrahedron,
29 import Examples (trilinear, trilinear9x9x9, zeros)
30 import FunctionValues (make_values, value_at)
31 import Point (Point(..))
32 import ScaleFactor (ScaleFactor)
34 Tetrahedron(v0,v1,v2,v3),
38 import Values (Values3D, dims, empty3d, zoom_shape)
41 -- | Our problem is defined on a Grid. The grid size is given by the
42 -- positive number h, which we have defined to always be 1 for
43 -- performance reasons (and simplicity). The function values are the
44 -- values of the function at the grid points, which are distance h=1
45 -- from one another in each direction (x,y,z).
46 data Grid = Grid { function_values :: Values3D }
50 instance Arbitrary Grid where
52 fvs <- arbitrary :: Gen Values3D
57 -- | Takes a grid and a position as an argument and returns the cube
58 -- centered on that position. If there is no cube there, well, you
59 -- shouldn't have done that. The omitted "otherwise" case actually
60 -- does improve performance.
61 cube_at :: Grid -> Int -> Int -> Int -> Cube
63 Cube i j k fvs' tet_vol
65 fvs = function_values g
66 fvs' = make_values fvs i j k
70 -- The first cube along any axis covers (-1/2, 1/2). The second
71 -- covers (1/2, 3/2). The third, (3/2, 5/2), and so on.
73 -- We translate the (x,y,z) coordinates forward by 1/2 so that the
74 -- first covers (0, 1), the second covers (1, 2), etc. This makes
75 -- it easy to figure out which cube contains the given point.
76 calculate_containing_cube_coordinate :: Grid -> Double -> Int
77 calculate_containing_cube_coordinate g coord
78 -- Don't use a cube on the boundary if we can help it. This
79 -- returns cube #1 if we would have returned cube #0 and cube #1
82 | coord == offset && (xsize > 1 && ysize > 1 && zsize > 1) = 1
83 | otherwise = (ceiling (coord + offset)) - 1
85 (xsize, ysize, zsize) = dims (function_values g)
89 -- | Takes a 'Grid', and returns a 'Cube' containing the given 'Point'.
90 -- Since our grid is rectangular, we can figure this out without having
91 -- to check every cube.
92 find_containing_cube :: Grid -> Point -> Cube
93 find_containing_cube g (Point x y z) =
96 i = calculate_containing_cube_coordinate g x
97 j = calculate_containing_cube_coordinate g y
98 k = calculate_containing_cube_coordinate g z
101 zoom_lookup :: Values3D -> ScaleFactor -> a -> (R.DIM3 -> Double)
102 zoom_lookup v3d scale_factor _ =
103 zoom_result v3d scale_factor
106 zoom_result :: Values3D -> ScaleFactor -> R.DIM3 -> Double
107 zoom_result v3d (sfx, sfy, sfz) (R.Z R.:. m R.:. n R.:. o) =
112 m' = (fromIntegral m) / (fromIntegral sfx) - offset
113 n' = (fromIntegral n) / (fromIntegral sfy) - offset
114 o' = (fromIntegral o) / (fromIntegral sfz) - offset
116 cube = find_containing_cube g p
117 t = find_containing_tetrahedron cube p
121 zoom :: Values3D -> ScaleFactor -> Values3D
122 zoom v3d scale_factor
123 | xsize == 0 || ysize == 0 || zsize == 0 = empty3d
125 R.computeS $ R.unsafeTraverse v3d transExtent f
127 (xsize, ysize, zsize) = dims v3d
128 transExtent = zoom_shape scale_factor
129 f = zoom_lookup v3d scale_factor
132 -- | Check all coefficients of tetrahedron0 belonging to the cube
133 -- centered on (1,1,1) with a grid constructed from the trilinear
134 -- values. See example one in the paper.
136 -- We also verify that the four vertices on face0 of the cube are
137 -- in the correct location.
139 trilinear_c0_t0_tests :: Test.Framework.Test
140 trilinear_c0_t0_tests =
141 testGroup "trilinear c0 t0"
142 [testGroup "coefficients"
143 [testCase "c0030 is correct" test_trilinear_c0030,
144 testCase "c0003 is correct" test_trilinear_c0003,
145 testCase "c0021 is correct" test_trilinear_c0021,
146 testCase "c0012 is correct" test_trilinear_c0012,
147 testCase "c0120 is correct" test_trilinear_c0120,
148 testCase "c0102 is correct" test_trilinear_c0102,
149 testCase "c0111 is correct" test_trilinear_c0111,
150 testCase "c0210 is correct" test_trilinear_c0210,
151 testCase "c0201 is correct" test_trilinear_c0201,
152 testCase "c0300 is correct" test_trilinear_c0300,
153 testCase "c1020 is correct" test_trilinear_c1020,
154 testCase "c1002 is correct" test_trilinear_c1002,
155 testCase "c1011 is correct" test_trilinear_c1011,
156 testCase "c1110 is correct" test_trilinear_c1110,
157 testCase "c1101 is correct" test_trilinear_c1101,
158 testCase "c1200 is correct" test_trilinear_c1200,
159 testCase "c2010 is correct" test_trilinear_c2010,
160 testCase "c2001 is correct" test_trilinear_c2001,
161 testCase "c2100 is correct" test_trilinear_c2100,
162 testCase "c3000 is correct" test_trilinear_c3000],
164 testGroup "face0 vertices"
165 [testCase "v0 is correct" test_trilinear_f0_t0_v0,
166 testCase "v1 is correct" test_trilinear_f0_t0_v1,
167 testCase "v2 is correct" test_trilinear_f0_t0_v2,
168 testCase "v3 is correct" test_trilinear_f0_t0_v3]
172 cube = cube_at g 1 1 1
173 t = tetrahedron cube 0
175 test_trilinear_c0030 :: Assertion
176 test_trilinear_c0030 =
177 assertAlmostEqual "c0030 is correct" (c t 0 0 3 0) (17/8)
179 test_trilinear_c0003 :: Assertion
180 test_trilinear_c0003 =
181 assertAlmostEqual "c0003 is correct" (c t 0 0 0 3) (27/8)
183 test_trilinear_c0021 :: Assertion
184 test_trilinear_c0021 =
185 assertAlmostEqual "c0021 is correct" (c t 0 0 2 1) (61/24)
187 test_trilinear_c0012 :: Assertion
188 test_trilinear_c0012 =
189 assertAlmostEqual "c0012 is correct" (c t 0 0 1 2) (71/24)
191 test_trilinear_c0120 :: Assertion
192 test_trilinear_c0120 =
193 assertAlmostEqual "c0120 is correct" (c t 0 1 2 0) (55/24)
195 test_trilinear_c0102 :: Assertion
196 test_trilinear_c0102 =
197 assertAlmostEqual "c0102 is correct" (c t 0 1 0 2) (73/24)
199 test_trilinear_c0111 :: Assertion
200 test_trilinear_c0111 =
201 assertAlmostEqual "c0111 is correct" (c t 0 1 1 1) (8/3)
203 test_trilinear_c0210 :: Assertion
204 test_trilinear_c0210 =
205 assertAlmostEqual "c0210 is correct" (c t 0 2 1 0) (29/12)
207 test_trilinear_c0201 :: Assertion
208 test_trilinear_c0201 =
209 assertAlmostEqual "c0201 is correct" (c t 0 2 0 1) (11/4)
211 test_trilinear_c0300 :: Assertion
212 test_trilinear_c0300 =
213 assertAlmostEqual "c0300 is correct" (c t 0 3 0 0) (5/2)
215 test_trilinear_c1020 :: Assertion
216 test_trilinear_c1020 =
217 assertAlmostEqual "c1020 is correct" (c t 1 0 2 0) (8/3)
219 test_trilinear_c1002 :: Assertion
220 test_trilinear_c1002 =
221 assertAlmostEqual "c1002 is correct" (c t 1 0 0 2) (23/6)
223 test_trilinear_c1011 :: Assertion
224 test_trilinear_c1011 =
225 assertAlmostEqual "c1011 is correct" (c t 1 0 1 1) (13/4)
227 test_trilinear_c1110 :: Assertion
228 test_trilinear_c1110 =
229 assertAlmostEqual "c1110 is correct" (c t 1 1 1 0) (23/8)
231 test_trilinear_c1101 :: Assertion
232 test_trilinear_c1101 =
233 assertAlmostEqual "c1101 is correct" (c t 1 1 0 1) (27/8)
235 test_trilinear_c1200 :: Assertion
236 test_trilinear_c1200 =
237 assertAlmostEqual "c1200 is correct" (c t 1 2 0 0) 3
239 test_trilinear_c2010 :: Assertion
240 test_trilinear_c2010 =
241 assertAlmostEqual "c2010 is correct" (c t 2 0 1 0) (10/3)
243 test_trilinear_c2001 :: Assertion
244 test_trilinear_c2001 =
245 assertAlmostEqual "c2001 is correct" (c t 2 0 0 1) 4
247 test_trilinear_c2100 :: Assertion
248 test_trilinear_c2100 =
249 assertAlmostEqual "c2100 is correct" (c t 2 1 0 0) (7/2)
251 test_trilinear_c3000 :: Assertion
252 test_trilinear_c3000 =
253 assertAlmostEqual "c3000 is correct" (c t 3 0 0 0) 4
255 test_trilinear_f0_t0_v0 :: Assertion
256 test_trilinear_f0_t0_v0 =
257 assertEqual "v0 is correct" (v0 t) (Point 1 1 1)
259 test_trilinear_f0_t0_v1 :: Assertion
260 test_trilinear_f0_t0_v1 =
261 assertEqual "v1 is correct" (v1 t) (Point 0.5 1 1)
263 test_trilinear_f0_t0_v2 :: Assertion
264 test_trilinear_f0_t0_v2 =
265 assertEqual "v2 is correct" (v2 t) (Point 0.5 0.5 1.5)
267 test_trilinear_f0_t0_v3 :: Assertion
268 test_trilinear_f0_t0_v3 =
269 assertEqual "v3 is correct" (v3 t) (Point 0.5 1.5 1.5)
272 test_trilinear_reproduced :: Assertion
273 test_trilinear_reproduced =
274 assertTrue "trilinears are reproduced correctly" $
275 and [p (Point i' j' k') ~= value_at trilinear i j k
281 let p = polynomial t,
282 let i' = fromIntegral i,
283 let j' = fromIntegral j,
284 let k' = fromIntegral k]
287 cs = [ cube_at g ci cj ck | ci <- [0..2], cj <- [0..2], ck <- [0..2] ]
290 test_zeros_reproduced :: Assertion
291 test_zeros_reproduced =
292 assertTrue "the zero function is reproduced correctly" $
293 and [p (Point i' j' k') ~= value_at zeros i j k
297 let i' = fromIntegral i,
298 let j' = fromIntegral j,
299 let k' = fromIntegral k,
302 let p = polynomial t0 ]
305 cs = [ cube_at g ci cj ck | ci <- [0..2], cj <- [0..2], ck <- [0..2] ]
308 -- | Make sure we can reproduce a 9x9x9 trilinear from the 3x3x3 one.
309 test_trilinear9x9x9_reproduced :: Assertion
310 test_trilinear9x9x9_reproduced =
311 assertTrue "trilinear 9x9x9 is reproduced correctly" $
312 and [p (Point i' j' k') ~= value_at trilinear9x9x9 i j k
317 let p = polynomial t,
318 let i' = (fromIntegral i) * 0.5,
319 let j' = (fromIntegral j) * 0.5,
320 let k' = (fromIntegral k) * 0.5]
327 prop_cube_indices_never_go_out_of_bounds :: Grid -> Gen Bool
328 prop_cube_indices_never_go_out_of_bounds g =
330 let coordmin = negate (1/2)
332 let (xsize, ysize, zsize) = dims $ function_values g
333 let xmax = (fromIntegral xsize) - (1/2)
334 let ymax = (fromIntegral ysize) - (1/2)
335 let zmax = (fromIntegral zsize) - (1/2)
337 x <- choose (coordmin, xmax)
338 y <- choose (coordmin, ymax)
339 z <- choose (coordmin, zmax)
341 let idx_x = calculate_containing_cube_coordinate g x
342 let idx_y = calculate_containing_cube_coordinate g y
343 let idx_z = calculate_containing_cube_coordinate g z
347 idx_x <= xsize - 1 &&
349 idx_y <= ysize - 1 &&
354 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). Note that the
355 -- third and fourth indices of c-t10 have been switched. This is
356 -- because we store the triangles oriented such that their volume is
357 -- positive. If T and T-tilde share \<v1,v2,v3\> and v0,v0-tilde point
358 -- in opposite directions, one of them has to have negative volume!
359 prop_c0120_identity :: Grid -> Property
360 prop_c0120_identity g =
361 xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
362 c t0 0 1 2 0 ~= (c t0 1 0 2 0 + c t10 1 0 0 2) / 2
364 fvs = function_values g
365 (xsize, ysize, zsize) = dims fvs
366 cube0 = cube_at g 1 1 1
367 cube1 = cube_at g 0 1 1
368 t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
369 t10 = tetrahedron cube1 10
372 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
373 -- 'prop_c0120_identity'.
374 prop_c0111_identity :: Grid -> Property
375 prop_c0111_identity g =
376 xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
377 c t0 0 1 1 1 ~= (c t0 1 0 1 1 + c t10 1 0 1 1) / 2
379 fvs = function_values g
380 (xsize, ysize, zsize) = dims fvs
381 cube0 = cube_at g 1 1 1
382 cube1 = cube_at g 0 1 1
383 t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
384 t10 = tetrahedron cube1 10
387 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
388 -- 'prop_c0120_identity'.
389 prop_c0201_identity :: Grid -> Property
390 prop_c0201_identity g =
391 xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
392 c t0 0 2 0 1 ~= (c t0 1 1 0 1 + c t10 1 1 1 0) / 2
394 fvs = function_values g
395 (xsize, ysize, zsize) = dims fvs
396 cube0 = cube_at g 1 1 1
397 cube1 = cube_at g 0 1 1
398 t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
399 t10 = tetrahedron cube1 10
402 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
403 -- 'prop_c0120_identity'.
404 prop_c0102_identity :: Grid -> Property
405 prop_c0102_identity g =
406 xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
407 c t0 0 1 0 2 ~= (c t0 1 0 0 2 + c t10 1 0 2 0) / 2
409 fvs = function_values g
410 (xsize, ysize, zsize) = dims fvs
411 cube0 = cube_at g 1 1 1
412 cube1 = cube_at g 0 1 1
413 t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
414 t10 = tetrahedron cube1 10
417 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
418 -- 'prop_c0120_identity'.
419 prop_c0210_identity :: Grid -> Property
420 prop_c0210_identity g =
421 xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
422 c t0 0 2 1 0 ~= (c t0 1 1 1 0 + c t10 1 1 0 1) / 2
424 fvs = function_values g
425 (xsize, ysize, zsize) = dims fvs
426 cube0 = cube_at g 1 1 1
427 cube1 = cube_at g 0 1 1
428 t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
429 t10 = tetrahedron cube1 10
432 -- | Given in Sorokina and Zeilfelder, p. 80, (2.9). See
433 -- 'prop_c0120_identity'.
434 prop_c0300_identity :: Grid -> Property
435 prop_c0300_identity g =
436 xsize >= 3 && ysize >= 3 && zsize >= 3 ==>
437 c t0 0 3 0 0 ~= (c t0 1 2 0 0 + c t10 1 2 0 0) / 2
439 fvs = function_values g
440 (xsize, ysize, zsize) = dims fvs
441 cube0 = cube_at g 1 1 1
442 cube1 = cube_at g 0 1 1
443 t0 = tetrahedron cube0 0 -- These two tetrahedra share a face.
444 t10 = tetrahedron cube1 10
447 -- | All of the properties from Section (2.9), p. 80. These require a
448 -- grid since they refer to two adjacent cubes.
449 p80_29_properties :: Test.Framework.Test
451 testGroup "p. 80, Section (2.9) Properties" [
452 testProperty "c0120 identity" prop_c0120_identity,
453 testProperty "c0111 identity" prop_c0111_identity,
454 testProperty "c0201 identity" prop_c0201_identity,
455 testProperty "c0102 identity" prop_c0102_identity,
456 testProperty "c0210 identity" prop_c0210_identity,
457 testProperty "c0300 identity" prop_c0300_identity ]
460 grid_tests :: Test.Framework.Test
462 testGroup "Grid Tests" [
463 trilinear_c0_t0_tests,
465 testProperty "cube indices within bounds"
466 prop_cube_indices_never_go_out_of_bounds ]
469 -- Do the slow tests last so we can stop paying attention.
470 slow_tests :: Test.Framework.Test
472 testGroup "Slow Tests" [
473 testCase "trilinear reproduced" test_trilinear_reproduced,
474 testCase "trilinear9x9x9 reproduced" test_trilinear9x9x9_reproduced,
475 testCase "zeros reproduced" test_zeros_reproduced ]