4 import Data.Maybe (fromJust)
5 import qualified Data.Vector as V (
14 import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
17 import qualified Face (Face(Face, v0, v1, v2, v3))
20 import Tetrahedron hiding (c)
21 import ThreeDimensional
23 data Cube = Cube { h :: Double,
28 tetrahedra_volume :: Double }
32 instance Arbitrary Cube where
34 (Positive h') <- arbitrary :: Gen (Positive Double)
35 i' <- choose (coordmin, coordmax)
36 j' <- choose (coordmin, coordmax)
37 k' <- choose (coordmin, coordmax)
38 fv' <- arbitrary :: Gen FunctionValues
39 (Positive tet_vol) <- arbitrary :: Gen (Positive Double)
40 return (Cube h' i' j' k' fv' tet_vol)
42 coordmin = -268435456 -- -(2^29 / 2)
43 coordmax = 268435456 -- +(2^29 / 2)
46 instance Show Cube where
48 "Cube_" ++ subscript ++ "\n" ++
49 " h: " ++ (show (h c)) ++ "\n" ++
50 " Center: " ++ (show (center c)) ++ "\n" ++
51 " xmin: " ++ (show (xmin c)) ++ "\n" ++
52 " xmax: " ++ (show (xmax c)) ++ "\n" ++
53 " ymin: " ++ (show (ymin c)) ++ "\n" ++
54 " ymax: " ++ (show (ymax c)) ++ "\n" ++
55 " zmin: " ++ (show (zmin c)) ++ "\n" ++
56 " zmax: " ++ (show (zmax c)) ++ "\n" ++
57 " fv: " ++ (show (Cube.fv c)) ++ "\n"
60 (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c))
63 -- | Returns an empty 'Cube'.
65 empty_cube = Cube 0 0 0 0 empty_values 0
68 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
70 xmin :: Cube -> Double
71 xmin c = (2*i' - 1)*delta / 2
73 i' = fromIntegral (i c) :: Double
76 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
78 xmax :: Cube -> Double
79 xmax c = (2*i' + 1)*delta / 2
81 i' = fromIntegral (i c) :: Double
84 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
86 ymin :: Cube -> Double
87 ymin c = (2*j' - 1)*delta / 2
89 j' = fromIntegral (j c) :: Double
92 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
94 ymax :: Cube -> Double
95 ymax c = (2*j' + 1)*delta / 2
97 j' = fromIntegral (j c) :: Double
100 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
102 zmin :: Cube -> Double
103 zmin c = (2*k' - 1)*delta / 2
105 k' = fromIntegral (k c) :: Double
108 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
110 zmax :: Cube -> Double
111 zmax c = (2*k' + 1)*delta / 2
113 k' = fromIntegral (k c) :: Double
116 instance ThreeDimensional Cube where
117 -- | The center of Cube_ijk coincides with v_ijk at
118 -- (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
122 i' = fromIntegral (i c) :: Double
123 j' = fromIntegral (j c) :: Double
124 k' = fromIntegral (k c) :: Double
129 -- | It's easy to tell if a point is within a cube; just make sure
130 -- that it falls on the proper side of each of the cube's faces.
131 contains_point c (x, y, z)
132 | x < (xmin c) = False
133 | x > (xmax c) = False
134 | y < (ymin c) = False
135 | y > (ymax c) = False
136 | z < (zmin c) = False
137 | z > (zmax c) = False
144 -- | The top (in the direction of z) face of the cube.
145 top_face :: Cube -> Face.Face
146 top_face c = Face.Face v0' v1' v2' v3'
149 v0' = (center c) + (delta, -delta, delta)
150 v1' = (center c) + (delta, delta, delta)
151 v2' = (center c) + (-delta, delta, delta)
152 v3' = (center c) + (-delta, -delta, delta)
156 -- | The back (in the direction of x) face of the cube.
157 back_face :: Cube -> Face.Face
158 back_face c = Face.Face v0' v1' v2' v3'
161 v0' = (center c) + (delta, -delta, -delta)
162 v1' = (center c) + (delta, delta, -delta)
163 v2' = (center c) + (delta, delta, delta)
164 v3' = (center c) + (delta, -delta, delta)
167 -- The bottom face (in the direction of -z) of the cube.
168 down_face :: Cube -> Face.Face
169 down_face c = Face.Face v0' v1' v2' v3'
172 v0' = (center c) + (-delta, -delta, -delta)
173 v1' = (center c) + (-delta, delta, -delta)
174 v2' = (center c) + (delta, delta, -delta)
175 v3' = (center c) + (delta, -delta, -delta)
179 -- | The front (in the direction of -x) face of the cube.
180 front_face :: Cube -> Face.Face
181 front_face c = Face.Face v0' v1' v2' v3'
184 v0' = (center c) + (-delta, -delta, delta)
185 v1' = (center c) + (-delta, delta, delta)
186 v2' = (center c) + (-delta, delta, -delta)
187 v3' = (center c) + (-delta, -delta, -delta)
189 -- | The left (in the direction of -y) face of the cube.
190 left_face :: Cube -> Face.Face
191 left_face c = Face.Face v0' v1' v2' v3'
194 v0' = (center c) + (delta, -delta, delta)
195 v1' = (center c) + (-delta, -delta, delta)
196 v2' = (center c) + (-delta, -delta, -delta)
197 v3' = (center c) + (delta, -delta, -delta)
200 -- | The right (in the direction of y) face of the cube.
201 right_face :: Cube -> Face.Face
202 right_face c = Face.Face v0' v1' v2' v3'
205 v0' = (center c) + (-delta, delta, delta)
206 v1' = (center c) + (delta, delta, delta)
207 v2' = (center c) + (delta, delta, -delta)
208 v3' = (center c) + (-delta, delta, -delta)
211 tetrahedron :: Cube -> Int -> Tetrahedron
214 Tetrahedron (Cube.fv c) v0' v1' v2' v3' vol 0
217 v1' = center (front_face c)
218 v2' = Face.v0 (front_face c)
219 v3' = Face.v1 (front_face c)
220 vol = tetrahedra_volume c
223 Tetrahedron fv' v0' v1' v2' v3' vol 1
226 v1' = center (front_face c)
227 v2' = Face.v1 (front_face c)
228 v3' = Face.v2 (front_face c)
229 fv' = rotate ccwx (Cube.fv c)
230 vol = tetrahedra_volume c
233 Tetrahedron fv' v0' v1' v2' v3' vol 2
236 v1' = center (front_face c)
237 v2' = Face.v2 (front_face c)
238 v3' = Face.v3 (front_face c)
239 fv' = rotate ccwx $ rotate ccwx $ Cube.fv c
240 vol = tetrahedra_volume c
243 Tetrahedron fv' v0' v1' v2' v3' vol 3
246 v1' = center (front_face c)
247 v2' = Face.v3 (front_face c)
248 v3' = Face.v0 (front_face c)
249 fv' = rotate cwx (Cube.fv c)
250 vol = tetrahedra_volume c
253 Tetrahedron fv' v0' v1' v2' v3' vol 4
256 v1' = center (top_face c)
257 v2' = Face.v0 (top_face c)
258 v3' = Face.v1 (top_face c)
259 fv' = rotate cwy (Cube.fv c)
260 vol = tetrahedra_volume c
263 Tetrahedron fv' v0' v1' v2' v3' vol 5
266 v1' = center (top_face c)
267 v2' = Face.v1 (top_face c)
268 v3' = Face.v2 (top_face c)
269 fv' = rotate cwy $ rotate cwz $ Tetrahedron.fv (tetrahedron c 0)
270 vol = tetrahedra_volume c
273 Tetrahedron fv' v0' v1' v2' v3' vol 6
276 v1' = center (top_face c)
277 v2' = Face.v2 (top_face c)
278 v3' = Face.v3 (top_face c)
279 fv' = rotate cwy $ rotate cwz
281 $ Tetrahedron.fv (tetrahedron c 0)
282 vol = tetrahedra_volume c
285 Tetrahedron fv' v0' v1' v2' v3' vol 7
288 v1' = center (top_face c)
289 v2' = Face.v3 (top_face c)
290 v3' = Face.v0 (top_face c)
291 fv' = rotate cwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron c 0)
292 vol = tetrahedra_volume c
295 Tetrahedron fv' v0' v1' v2' v3' vol 8
298 v1' = center (back_face c)
299 v2' = Face.v0 (back_face c)
300 v3' = Face.v1 (back_face c)
301 fv' = rotate cwy $ rotate cwy $ Tetrahedron.fv (tetrahedron c 0)
302 vol = tetrahedra_volume c
305 Tetrahedron fv' v0' v1' v2' v3' vol 9
308 v1' = center (back_face c)
309 v2' = Face.v1 (back_face c)
310 v3' = Face.v2 (back_face c)
311 fv' = rotate cwy $ rotate cwy
313 $ Tetrahedron.fv (tetrahedron c 0)
314 vol = tetrahedra_volume c
317 Tetrahedron fv' v0' v1' v2' v3' vol 10
320 v1' = center (back_face c)
321 v2' = Face.v2 (back_face c)
322 v3' = Face.v3 (back_face c)
323 fv' = rotate cwy $ rotate cwy
326 $ Tetrahedron.fv (tetrahedron c 0)
328 vol = tetrahedra_volume c
331 Tetrahedron fv' v0' v1' v2' v3' vol 11
334 v1' = center (back_face c)
335 v2' = Face.v3 (back_face c)
336 v3' = Face.v0 (back_face c)
337 fv' = rotate cwy $ rotate cwy
339 $ Tetrahedron.fv (tetrahedron c 0)
340 vol = tetrahedra_volume c
343 Tetrahedron fv' v0' v1' v2' v3' vol 12
346 v1' = center (down_face c)
347 v2' = Face.v0 (down_face c)
348 v3' = Face.v1 (down_face c)
349 fv' = rotate ccwy (Tetrahedron.fv (tetrahedron c 0))
350 vol = tetrahedra_volume c
353 Tetrahedron fv' v0' v1' v2' v3' vol 13
356 v1' = center (down_face c)
357 v2' = Face.v1 (down_face c)
358 v3' = Face.v2 (down_face c)
359 fv' = rotate ccwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron c 0)
360 vol = tetrahedra_volume c
363 Tetrahedron fv' v0' v1' v2' v3' vol 14
366 v1' = center (down_face c)
367 v2' = Face.v2 (down_face c)
368 v3' = Face.v3 (down_face c)
369 fv' = rotate ccwy $ rotate ccwz
371 $ Tetrahedron.fv (tetrahedron c 0)
372 vol = tetrahedra_volume c
375 Tetrahedron fv' v0' v1' v2' v3' vol 15
378 v1' = center (down_face c)
379 v2' = Face.v3 (down_face c)
380 v3' = Face.v0 (down_face c)
381 fv' = rotate ccwy $ rotate cwz $ Tetrahedron.fv (tetrahedron c 0)
382 vol = tetrahedra_volume c
385 Tetrahedron fv' v0' v1' v2' v3' vol 16
388 v1' = center (right_face c)
389 v2' = Face.v0 (right_face c)
390 v3' = Face.v1 (right_face c)
391 fv' = rotate ccwz (Tetrahedron.fv (tetrahedron c 0))
392 vol = tetrahedra_volume c
395 Tetrahedron fv' v0' v1' v2' v3' vol 17
398 v1' = center (right_face c)
399 v2' = Face.v1 (right_face c)
400 v3' = Face.v2 (right_face c)
401 fv' = rotate ccwz $ rotate cwy $ Tetrahedron.fv (tetrahedron c 0)
402 vol = tetrahedra_volume c
405 Tetrahedron fv' v0' v1' v2' v3' vol 18
408 v1' = center (right_face c)
409 v2' = Face.v2 (right_face c)
410 v3' = Face.v3 (right_face c)
411 fv' = rotate ccwz $ rotate cwy
413 $ Tetrahedron.fv (tetrahedron c 0)
414 vol = tetrahedra_volume c
417 Tetrahedron fv' v0' v1' v2' v3' vol 19
420 v1' = center (right_face c)
421 v2' = Face.v3 (right_face c)
422 v3' = Face.v0 (right_face c)
423 fv' = rotate ccwz $ rotate ccwy
424 $ Tetrahedron.fv (tetrahedron c 0)
425 vol = tetrahedra_volume c
428 Tetrahedron fv' v0' v1' v2' v3' vol 20
431 v1' = center (left_face c)
432 v2' = Face.v0 (left_face c)
433 v3' = Face.v1 (left_face c)
434 fv' = rotate cwz (Tetrahedron.fv (tetrahedron c 0))
435 vol = tetrahedra_volume c
438 Tetrahedron fv' v0' v1' v2' v3' vol 21
441 v1' = center (left_face c)
442 v2' = Face.v1 (left_face c)
443 v3' = Face.v2 (left_face c)
444 fv' = rotate cwz $ rotate ccwy $ Tetrahedron.fv (tetrahedron c 0)
445 vol = tetrahedra_volume c
448 Tetrahedron fv' v0' v1' v2' v3' vol 22
451 v1' = center (left_face c)
452 v2' = Face.v2 (left_face c)
453 v3' = Face.v3 (left_face c)
454 fv' = rotate cwz $ rotate ccwy
456 $ Tetrahedron.fv (tetrahedron c 0)
457 vol = tetrahedra_volume c
460 Tetrahedron fv' v0' v1' v2' v3' vol 23
463 v1' = center (left_face c)
464 v2' = Face.v3 (left_face c)
465 v3' = Face.v0 (left_face c)
466 fv' = rotate cwz $ rotate cwy
467 $ Tetrahedron.fv (tetrahedron c 0)
468 vol = tetrahedra_volume c
470 -- Feels dirty, but whatever.
471 tetrahedron _ _ = error "asked for a nonexistent tetrahedron"
474 -- Only used in tests, so we don't need the added speed
476 tetrahedra :: Cube -> [Tetrahedron]
477 tetrahedra c = [ tetrahedron c n | n <- [0..23] ]
479 front_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
480 front_left_top_tetrahedra c =
481 V.singleton (tetrahedron c 0) `V.snoc`
482 (tetrahedron c 3) `V.snoc`
483 (tetrahedron c 6) `V.snoc`
484 (tetrahedron c 7) `V.snoc`
485 (tetrahedron c 20) `V.snoc`
488 front_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
489 front_left_down_tetrahedra c =
490 V.singleton (tetrahedron c 0) `V.snoc`
491 (tetrahedron c 2) `V.snoc`
492 (tetrahedron c 3) `V.snoc`
493 (tetrahedron c 12) `V.snoc`
494 (tetrahedron c 15) `V.snoc`
497 front_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
498 front_right_top_tetrahedra c =
499 V.singleton (tetrahedron c 0) `V.snoc`
500 (tetrahedron c 1) `V.snoc`
501 (tetrahedron c 5) `V.snoc`
502 (tetrahedron c 6) `V.snoc`
503 (tetrahedron c 16) `V.snoc`
506 front_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
507 front_right_down_tetrahedra c =
508 V.singleton (tetrahedron c 1) `V.snoc`
509 (tetrahedron c 2) `V.snoc`
510 (tetrahedron c 12) `V.snoc`
511 (tetrahedron c 13) `V.snoc`
512 (tetrahedron c 18) `V.snoc`
515 back_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
516 back_left_top_tetrahedra c =
517 V.singleton (tetrahedron c 0) `V.snoc`
518 (tetrahedron c 3) `V.snoc`
519 (tetrahedron c 6) `V.snoc`
520 (tetrahedron c 7) `V.snoc`
521 (tetrahedron c 20) `V.snoc`
524 back_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
525 back_left_down_tetrahedra c =
526 V.singleton (tetrahedron c 8) `V.snoc`
527 (tetrahedron c 11) `V.snoc`
528 (tetrahedron c 14) `V.snoc`
529 (tetrahedron c 15) `V.snoc`
530 (tetrahedron c 22) `V.snoc`
533 back_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
534 back_right_top_tetrahedra c =
535 V.singleton (tetrahedron c 4) `V.snoc`
536 (tetrahedron c 5) `V.snoc`
537 (tetrahedron c 9) `V.snoc`
538 (tetrahedron c 10) `V.snoc`
539 (tetrahedron c 16) `V.snoc`
542 back_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
543 back_right_down_tetrahedra c =
544 V.singleton (tetrahedron c 8) `V.snoc`
545 (tetrahedron c 9) `V.snoc`
546 (tetrahedron c 13) `V.snoc`
547 (tetrahedron c 14) `V.snoc`
548 (tetrahedron c 17) `V.snoc`
551 in_top_half :: Cube -> Point -> Bool
552 in_top_half c (_,_,z) =
553 distance_from_top <= distance_from_bottom
555 distance_from_top = abs $ (zmax c) - z
556 distance_from_bottom = abs $ (zmin c) - z
558 in_front_half :: Cube -> Point -> Bool
559 in_front_half c (x,_,_) =
560 distance_from_front <= distance_from_back
562 distance_from_front = abs $ (xmin c) - x
563 distance_from_back = abs $ (xmax c) - x
566 in_left_half :: Cube -> Point -> Bool
567 in_left_half c (_,y,_) =
568 distance_from_left <= distance_from_right
570 distance_from_left = abs $ (ymin c) - y
571 distance_from_right = abs $ (ymax c) - y
574 -- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
575 -- contain the given 'Point'. This should be faster than checking
576 -- every tetrahedron individually, since we determine which half
577 -- (hemisphere?) of the cube the point lies in three times: once in
578 -- each dimension. This allows us to eliminate non-candidates
581 -- This can throw an exception, but the use of 'head' might
582 -- save us some unnecessary computations.
584 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
585 find_containing_tetrahedron c p =
586 candidates `V.unsafeIndex` (fromJust lucky_idx)
588 front_half = in_front_half c p
589 top_half = in_top_half c p
590 left_half = in_left_half c p
597 front_left_top_tetrahedra c
599 front_left_down_tetrahedra c
602 front_right_top_tetrahedra c
604 front_right_down_tetrahedra c
610 back_left_top_tetrahedra c
612 back_left_down_tetrahedra c
615 back_right_top_tetrahedra c
617 back_right_down_tetrahedra c
619 -- Use the dot product instead of 'distance' here to save a
620 -- sqrt(). So, "distances" below really means "distances squared."
621 distances = V.map ((dot p) . center) candidates
622 shortest_distance = V.minimum distances
623 lucky_idx = V.findIndex
624 (\t -> (center t) `dot` p == shortest_distance)