4 import Data.List ( (\\) )
5 import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
8 import qualified Face (Face(Face, v0, v1, v2, v3))
11 import Tetrahedron hiding (c)
12 import ThreeDimensional
14 data Cube = Cube { h :: Double,
19 tetrahedra_volume :: Double }
23 instance Arbitrary Cube where
25 (Positive h') <- arbitrary :: Gen (Positive Double)
26 i' <- choose (coordmin, coordmax)
27 j' <- choose (coordmin, coordmax)
28 k' <- choose (coordmin, coordmax)
29 fv' <- arbitrary :: Gen FunctionValues
30 (Positive tet_vol) <- arbitrary :: Gen (Positive Double)
31 return (Cube h' i' j' k' fv' tet_vol)
33 coordmin = -268435456 -- -(2^29 / 2)
34 coordmax = 268435456 -- +(2^29 / 2)
37 instance Show Cube where
39 "Cube_" ++ subscript ++ "\n" ++
40 " h: " ++ (show (h c)) ++ "\n" ++
41 " Center: " ++ (show (center c)) ++ "\n" ++
42 " xmin: " ++ (show (xmin c)) ++ "\n" ++
43 " xmax: " ++ (show (xmax c)) ++ "\n" ++
44 " ymin: " ++ (show (ymin c)) ++ "\n" ++
45 " ymax: " ++ (show (ymax c)) ++ "\n" ++
46 " zmin: " ++ (show (zmin c)) ++ "\n" ++
47 " zmax: " ++ (show (zmax c)) ++ "\n" ++
48 " fv: " ++ (show (Cube.fv c)) ++ "\n"
51 (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c))
54 -- | Returns an empty 'Cube'.
56 empty_cube = Cube 0 0 0 0 empty_values 0
59 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
61 xmin :: Cube -> Double
62 xmin c = (2*i' - 1)*delta / 2
64 i' = fromIntegral (i c) :: Double
67 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
69 xmax :: Cube -> Double
70 xmax c = (2*i' + 1)*delta / 2
72 i' = fromIntegral (i c) :: Double
75 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
77 ymin :: Cube -> Double
78 ymin c = (2*j' - 1)*delta / 2
80 j' = fromIntegral (j c) :: Double
83 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
85 ymax :: Cube -> Double
86 ymax c = (2*j' + 1)*delta / 2
88 j' = fromIntegral (j c) :: Double
91 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
93 zmin :: Cube -> Double
94 zmin c = (2*k' - 1)*delta / 2
96 k' = fromIntegral (k c) :: Double
99 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
101 zmax :: Cube -> Double
102 zmax c = (2*k' + 1)*delta / 2
104 k' = fromIntegral (k c) :: Double
107 instance ThreeDimensional Cube where
108 -- | The center of Cube_ijk coincides with v_ijk at
109 -- (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
113 i' = fromIntegral (i c) :: Double
114 j' = fromIntegral (j c) :: Double
115 k' = fromIntegral (k c) :: Double
120 -- | It's easy to tell if a point is within a cube; just make sure
121 -- that it falls on the proper side of each of the cube's faces.
122 contains_point c (x, y, z)
123 | x < (xmin c) = False
124 | x > (xmax c) = False
125 | y < (ymin c) = False
126 | y > (ymax c) = False
127 | z < (zmin c) = False
128 | z > (zmax c) = False
135 -- | The top (in the direction of z) face of the cube.
136 top_face :: Cube -> Face.Face
137 top_face c = Face.Face v0' v1' v2' v3'
140 v0' = (center c) + (delta, -delta, delta)
141 v1' = (center c) + (delta, delta, delta)
142 v2' = (center c) + (-delta, delta, delta)
143 v3' = (center c) + (-delta, -delta, delta)
147 -- | The back (in the direction of x) face of the cube.
148 back_face :: Cube -> Face.Face
149 back_face c = Face.Face v0' v1' v2' v3'
152 v0' = (center c) + (delta, -delta, -delta)
153 v1' = (center c) + (delta, delta, -delta)
154 v2' = (center c) + (delta, delta, delta)
155 v3' = (center c) + (delta, -delta, delta)
158 -- The bottom face (in the direction of -z) of the cube.
159 down_face :: Cube -> Face.Face
160 down_face c = Face.Face v0' v1' v2' v3'
163 v0' = (center c) + (-delta, -delta, -delta)
164 v1' = (center c) + (-delta, delta, -delta)
165 v2' = (center c) + (delta, delta, -delta)
166 v3' = (center c) + (delta, -delta, -delta)
170 -- | The front (in the direction of -x) face of the cube.
171 front_face :: Cube -> Face.Face
172 front_face c = Face.Face v0' v1' v2' v3'
175 v0' = (center c) + (-delta, -delta, delta)
176 v1' = (center c) + (-delta, delta, delta)
177 v2' = (center c) + (-delta, delta, -delta)
178 v3' = (center c) + (-delta, -delta, -delta)
180 -- | The left (in the direction of -y) face of the cube.
181 left_face :: Cube -> Face.Face
182 left_face c = Face.Face v0' v1' v2' v3'
185 v0' = (center c) + (delta, -delta, delta)
186 v1' = (center c) + (-delta, -delta, delta)
187 v2' = (center c) + (-delta, -delta, -delta)
188 v3' = (center c) + (delta, -delta, -delta)
191 -- | The right (in the direction of y) face of the cube.
192 right_face :: Cube -> Face.Face
193 right_face c = Face.Face v0' v1' v2' v3'
196 v0' = (center c) + (-delta, delta, delta)
197 v1' = (center c) + (delta, delta, delta)
198 v2' = (center c) + (delta, delta, -delta)
199 v3' = (center c) + (-delta, delta, -delta)
202 tetrahedron0 :: Cube -> Tetrahedron
204 Tetrahedron (Cube.fv c) v0' v1' v2' v3' vol
207 v1' = center (front_face c)
208 v2' = Face.v0 (front_face c)
209 v3' = Face.v1 (front_face c)
210 vol = tetrahedra_volume c
212 tetrahedron1 :: Cube -> Tetrahedron
214 Tetrahedron fv' v0' v1' v2' v3' vol
217 v1' = center (front_face c)
218 v2' = Face.v1 (front_face c)
219 v3' = Face.v2 (front_face c)
220 fv' = rotate ccwx (Cube.fv c)
221 vol = tetrahedra_volume c
223 tetrahedron2 :: Cube -> Tetrahedron
225 Tetrahedron fv' v0' v1' v2' v3' vol
228 v1' = center (front_face c)
229 v2' = Face.v2 (front_face c)
230 v3' = Face.v3 (front_face c)
231 fv' = rotate ccwx $ rotate ccwx $ Cube.fv c
232 vol = tetrahedra_volume c
234 tetrahedron3 :: Cube -> Tetrahedron
236 Tetrahedron fv' v0' v1' v2' v3' vol
239 v1' = center (front_face c)
240 v2' = Face.v3 (front_face c)
241 v3' = Face.v0 (front_face c)
242 fv' = rotate cwx (Cube.fv c)
243 vol = tetrahedra_volume c
245 tetrahedron4 :: Cube -> Tetrahedron
247 Tetrahedron fv' v0' v1' v2' v3' vol
250 v1' = center (top_face c)
251 v2' = Face.v0 (top_face c)
252 v3' = Face.v1 (top_face c)
253 fv' = rotate cwy (Cube.fv c)
254 vol = tetrahedra_volume c
256 tetrahedron5 :: Cube -> Tetrahedron
258 Tetrahedron fv' v0' v1' v2' v3' vol
261 v1' = center (top_face c)
262 v2' = Face.v1 (top_face c)
263 v3' = Face.v2 (top_face c)
264 fv' = rotate cwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
265 vol = tetrahedra_volume c
267 tetrahedron6 :: Cube -> Tetrahedron
269 Tetrahedron fv' v0' v1' v2' v3' vol
272 v1' = center (top_face c)
273 v2' = Face.v2 (top_face c)
274 v3' = Face.v3 (top_face c)
275 fv' = rotate cwy $ rotate cwz
277 $ Tetrahedron.fv (tetrahedron0 c)
278 vol = tetrahedra_volume c
280 tetrahedron7 :: Cube -> Tetrahedron
282 Tetrahedron fv' v0' v1' v2' v3' vol
285 v1' = center (top_face c)
286 v2' = Face.v3 (top_face c)
287 v3' = Face.v0 (top_face c)
288 fv' = rotate cwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
289 vol = tetrahedra_volume c
291 tetrahedron8 :: Cube -> Tetrahedron
293 Tetrahedron fv' v0' v1' v2' v3' vol
296 v1' = center (back_face c)
297 v2' = Face.v0 (back_face c)
298 v3' = Face.v1 (back_face c)
299 fv' = rotate cwy $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
300 vol = tetrahedra_volume c
302 tetrahedron9 :: Cube -> Tetrahedron
304 Tetrahedron fv' v0' v1' v2' v3' vol
307 v1' = center (back_face c)
308 v2' = Face.v1 (back_face c)
309 v3' = Face.v2 (back_face c)
310 fv' = rotate cwy $ rotate cwy
312 $ Tetrahedron.fv (tetrahedron0 c)
313 vol = tetrahedra_volume c
315 tetrahedron10 :: Cube -> Tetrahedron
317 Tetrahedron fv' v0' v1' v2' v3' vol
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 (tetrahedron0 c)
328 vol = tetrahedra_volume c
330 tetrahedron11 :: Cube -> Tetrahedron
332 Tetrahedron fv' v0' v1' v2' v3' vol
335 v1' = center (back_face c)
336 v2' = Face.v3 (back_face c)
337 v3' = Face.v0 (back_face c)
338 fv' = rotate cwy $ rotate cwy
340 $ Tetrahedron.fv (tetrahedron0 c)
341 vol = tetrahedra_volume c
344 tetrahedron12 :: Cube -> Tetrahedron
346 Tetrahedron fv' v0' v1' v2' v3' vol
349 v1' = center (down_face c)
350 v2' = Face.v0 (down_face c)
351 v3' = Face.v1 (down_face c)
352 fv' = rotate ccwy (Tetrahedron.fv (tetrahedron0 c))
353 vol = tetrahedra_volume c
356 tetrahedron13 :: Cube -> Tetrahedron
358 Tetrahedron fv' v0' v1' v2' v3' vol
361 v1' = center (down_face c)
362 v2' = Face.v1 (down_face c)
363 v3' = Face.v2 (down_face c)
364 fv' = rotate ccwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
365 vol = tetrahedra_volume c
368 tetrahedron14 :: Cube -> Tetrahedron
370 Tetrahedron fv' v0' v1' v2' v3' vol
373 v1' = center (down_face c)
374 v2' = Face.v2 (down_face c)
375 v3' = Face.v3 (down_face c)
376 fv' = rotate ccwy $ rotate ccwz
378 $ Tetrahedron.fv (tetrahedron0 c)
379 vol = tetrahedra_volume c
382 tetrahedron15 :: Cube -> Tetrahedron
384 Tetrahedron fv' v0' v1' v2' v3' vol
387 v1' = center (down_face c)
388 v2' = Face.v3 (down_face c)
389 v3' = Face.v0 (down_face c)
390 fv' = rotate ccwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
391 vol = tetrahedra_volume c
394 tetrahedron16 :: Cube -> Tetrahedron
396 Tetrahedron fv' v0' v1' v2' v3' vol
399 v1' = center (right_face c)
400 v2' = Face.v0 (right_face c)
401 v3' = Face.v1 (right_face c)
402 fv' = rotate ccwz (Tetrahedron.fv (tetrahedron0 c))
403 vol = tetrahedra_volume c
406 tetrahedron17 :: Cube -> Tetrahedron
408 Tetrahedron fv' v0' v1' v2' v3' vol
411 v1' = center (right_face c)
412 v2' = Face.v1 (right_face c)
413 v3' = Face.v2 (right_face c)
414 fv' = rotate ccwz $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
415 vol = tetrahedra_volume c
418 tetrahedron18 :: Cube -> Tetrahedron
420 Tetrahedron fv' v0' v1' v2' v3' vol
423 v1' = center (right_face c)
424 v2' = Face.v2 (right_face c)
425 v3' = Face.v3 (right_face c)
426 fv' = rotate ccwz $ rotate cwy
428 $ Tetrahedron.fv (tetrahedron0 c)
429 vol = tetrahedra_volume c
432 tetrahedron19 :: Cube -> Tetrahedron
434 Tetrahedron fv' v0' v1' v2' v3' vol
437 v1' = center (right_face c)
438 v2' = Face.v3 (right_face c)
439 v3' = Face.v0 (right_face c)
440 fv' = rotate ccwz $ rotate ccwy
441 $ Tetrahedron.fv (tetrahedron0 c)
442 vol = tetrahedra_volume c
445 tetrahedron20 :: Cube -> Tetrahedron
447 Tetrahedron fv' v0' v1' v2' v3' vol
450 v1' = center (left_face c)
451 v2' = Face.v0 (left_face c)
452 v3' = Face.v1 (left_face c)
453 fv' = rotate cwz (Tetrahedron.fv (tetrahedron0 c))
454 vol = tetrahedra_volume c
457 tetrahedron21 :: Cube -> Tetrahedron
459 Tetrahedron fv' v0' v1' v2' v3' vol
462 v1' = center (left_face c)
463 v2' = Face.v1 (left_face c)
464 v3' = Face.v2 (left_face c)
465 fv' = rotate cwz $ rotate ccwy $ Tetrahedron.fv (tetrahedron0 c)
466 vol = tetrahedra_volume c
469 tetrahedron22 :: Cube -> Tetrahedron
471 Tetrahedron fv' v0' v1' v2' v3' vol
474 v1' = center (left_face c)
475 v2' = Face.v2 (left_face c)
476 v3' = Face.v3 (left_face c)
477 fv' = rotate cwz $ rotate ccwy
479 $ Tetrahedron.fv (tetrahedron0 c)
480 vol = tetrahedra_volume c
483 tetrahedron23 :: Cube -> Tetrahedron
485 Tetrahedron fv' v0' v1' v2' v3' vol
488 v1' = center (left_face c)
489 v2' = Face.v3 (left_face c)
490 v3' = Face.v0 (left_face c)
491 fv' = rotate cwz $ rotate cwy
492 $ Tetrahedron.fv (tetrahedron0 c)
493 vol = tetrahedra_volume c
496 tetrahedra :: Cube -> [Tetrahedron]
523 -- | All completely contained in the front half of the cube.
524 front_half_tetrahedra :: Cube -> [Tetrahedron]
525 front_half_tetrahedra c =
536 -- | All tetrahedra completely contained in the top half of the cube.
537 top_half_tetrahedra :: Cube -> [Tetrahedron]
538 top_half_tetrahedra c =
549 -- | All tetrahedra completely contained in the back half of the cube.
550 back_half_tetrahedra :: Cube -> [Tetrahedron]
551 back_half_tetrahedra c =
562 -- | All tetrahedra completely contained in the down half of the cube.
563 down_half_tetrahedra :: Cube -> [Tetrahedron]
564 down_half_tetrahedra c =
575 -- | All tetrahedra completely contained in the right half of the cube.
576 right_half_tetrahedra :: Cube -> [Tetrahedron]
577 right_half_tetrahedra c =
588 -- | All tetrahedra completely contained in the left half of the cube.
589 left_half_tetrahedra :: Cube -> [Tetrahedron]
590 left_half_tetrahedra c =
601 in_top_half :: Cube -> Point -> Bool
602 in_top_half c (_,_,z) =
603 distance_from_top <= distance_from_bottom
605 distance_from_top = abs $ (zmax c) - z
606 distance_from_bottom = abs $ (zmin c) - z
608 in_front_half :: Cube -> Point -> Bool
609 in_front_half c (x,_,_) =
610 distance_from_front <= distance_from_back
612 distance_from_front = abs $ (xmin c) - x
613 distance_from_back = abs $ (xmax c) - x
616 in_left_half :: Cube -> Point -> Bool
617 in_left_half c (_,y,_) =
618 distance_from_left <= distance_from_right
620 distance_from_left = abs $ (ymin c) - y
621 distance_from_right = abs $ (ymax c) - y
624 -- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
625 -- contain the given 'Point'. This should be faster than checking
626 -- every tetrahedron individually, since we determine which half
627 -- (hemisphere?) of the cube the point lies in three times: once in
628 -- each dimension. This allows us to eliminate non-candidates
631 -- This can throw an exception, but the use of 'head' might
632 -- save us some unnecessary computations.
634 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
635 find_containing_tetrahedron c p =
636 head containing_tetrahedra
638 candidates = tetrahedra c
640 if (in_front_half c p) then
641 back_half_tetrahedra c
643 front_half_tetrahedra c
645 candidates_x = candidates \\ non_candidates_x
648 if (in_left_half c p) then
649 right_half_tetrahedra c
651 left_half_tetrahedra c
653 candidates_xy = candidates_x \\ non_candidates_y
656 if (in_top_half c p) then
657 down_half_tetrahedra c
659 top_half_tetrahedra c
661 candidates_xyz = candidates_xy \\ non_candidates_z
663 contains_our_point = flip contains_point p
664 containing_tetrahedra = filter contains_our_point candidates_xyz