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,
18 fv :: FunctionValues }
22 instance Arbitrary Cube where
24 (Positive h') <- arbitrary :: Gen (Positive Double)
25 i' <- choose (coordmin, coordmax)
26 j' <- choose (coordmin, coordmax)
27 k' <- choose (coordmin, coordmax)
28 fv' <- arbitrary :: Gen FunctionValues
29 return (Cube h' i' j' k' fv')
31 coordmin = -268435456 -- -(2^29 / 2)
32 coordmax = 268435456 -- +(2^29 / 2)
35 instance Show Cube where
37 "Cube_" ++ subscript ++ "\n" ++
38 " h: " ++ (show (h c)) ++ "\n" ++
39 " Center: " ++ (show (center c)) ++ "\n" ++
40 " xmin: " ++ (show (xmin c)) ++ "\n" ++
41 " xmax: " ++ (show (xmax c)) ++ "\n" ++
42 " ymin: " ++ (show (ymin c)) ++ "\n" ++
43 " ymax: " ++ (show (ymax c)) ++ "\n" ++
44 " zmin: " ++ (show (zmin c)) ++ "\n" ++
45 " zmax: " ++ (show (zmax c)) ++ "\n" ++
46 " fv: " ++ (show (Cube.fv c)) ++ "\n"
49 (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c))
52 -- | Returns an empty 'Cube'.
54 empty_cube = Cube 0 0 0 0 empty_values
57 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
59 xmin :: Cube -> Double
60 xmin c = (2*i' - 1)*delta / 2
62 i' = fromIntegral (i c) :: Double
65 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
67 xmax :: Cube -> Double
68 xmax c = (2*i' + 1)*delta / 2
70 i' = fromIntegral (i c) :: Double
73 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
75 ymin :: Cube -> Double
76 ymin c = (2*j' - 1)*delta / 2
78 j' = fromIntegral (j c) :: Double
81 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
83 ymax :: Cube -> Double
84 ymax c = (2*j' + 1)*delta / 2
86 j' = fromIntegral (j c) :: Double
89 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
91 zmin :: Cube -> Double
92 zmin c = (2*k' - 1)*delta / 2
94 k' = fromIntegral (k c) :: Double
97 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
99 zmax :: Cube -> Double
100 zmax c = (2*k' + 1)*delta / 2
102 k' = fromIntegral (k c) :: Double
105 instance ThreeDimensional Cube where
106 -- | The center of Cube_ijk coincides with v_ijk at
107 -- (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
111 i' = fromIntegral (i c) :: Double
112 j' = fromIntegral (j c) :: Double
113 k' = fromIntegral (k c) :: Double
118 -- | It's easy to tell if a point is within a cube; just make sure
119 -- that it falls on the proper side of each of the cube's faces.
120 contains_point c (x, y, z)
121 | x < (xmin c) = False
122 | x > (xmax c) = False
123 | y < (ymin c) = False
124 | y > (ymax c) = False
125 | z < (zmin c) = False
126 | z > (zmax c) = False
133 -- | The top (in the direction of z) face of the cube.
134 top_face :: Cube -> Face.Face
135 top_face c = Face.Face v0' v1' v2' v3'
138 v0' = (center c) + (delta, -delta, delta)
139 v1' = (center c) + (delta, delta, delta)
140 v2' = (center c) + (-delta, delta, delta)
141 v3' = (center c) + (-delta, -delta, delta)
145 -- | The back (in the direction of x) face of the cube.
146 back_face :: Cube -> Face.Face
147 back_face c = Face.Face v0' v1' v2' v3'
150 v0' = (center c) + (delta, -delta, -delta)
151 v1' = (center c) + (delta, delta, -delta)
152 v2' = (center c) + (delta, delta, delta)
153 v3' = (center c) + (delta, -delta, delta)
156 -- The bottom face (in the direction of -z) of the cube.
157 down_face :: Cube -> Face.Face
158 down_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)
168 -- | The front (in the direction of -x) face of the cube.
169 front_face :: Cube -> Face.Face
170 front_face c = Face.Face v0' v1' v2' v3'
173 v0' = (center c) + (-delta, -delta, delta)
174 v1' = (center c) + (-delta, delta, delta)
175 v2' = (center c) + (-delta, delta, -delta)
176 v3' = (center c) + (-delta, -delta, -delta)
178 -- | The left (in the direction of -y) face of the cube.
179 left_face :: Cube -> Face.Face
180 left_face c = Face.Face v0' v1' v2' v3'
183 v0' = (center c) + (delta, -delta, delta)
184 v1' = (center c) + (-delta, -delta, delta)
185 v2' = (center c) + (-delta, -delta, -delta)
186 v3' = (center c) + (delta, -delta, -delta)
189 -- | The right (in the direction of y) face of the cube.
190 right_face :: Cube -> Face.Face
191 right_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 tetrahedron0 :: Cube -> Tetrahedron
202 Tetrahedron (Cube.fv c) v0' v1' v2' v3'
205 v1' = center (front_face c)
206 v2' = Face.v0 (front_face c)
207 v3' = Face.v1 (front_face c)
209 tetrahedron1 :: Cube -> Tetrahedron
211 Tetrahedron fv' v0' v1' v2' v3'
214 v1' = center (front_face c)
215 v2' = Face.v1 (front_face c)
216 v3' = Face.v2 (front_face c)
217 fv' = rotate ccwx (Cube.fv c)
219 tetrahedron2 :: Cube -> Tetrahedron
221 Tetrahedron fv' v0' v1' v2' v3'
224 v1' = center (front_face c)
225 v2' = Face.v2 (front_face c)
226 v3' = Face.v3 (front_face c)
227 fv' = rotate ccwx $ rotate ccwx $ Cube.fv c
229 tetrahedron3 :: Cube -> Tetrahedron
231 Tetrahedron fv' v0' v1' v2' v3'
234 v1' = center (front_face c)
235 v2' = Face.v3 (front_face c)
236 v3' = Face.v0 (front_face c)
237 fv' = rotate cwx (Cube.fv c)
239 tetrahedron4 :: Cube -> Tetrahedron
241 Tetrahedron fv' v0' v1' v2' v3'
244 v1' = center (top_face c)
245 v2' = Face.v0 (top_face c)
246 v3' = Face.v1 (top_face c)
247 fv' = rotate cwy (Cube.fv c)
249 tetrahedron5 :: Cube -> Tetrahedron
251 Tetrahedron fv' v0' v1' v2' v3'
254 v1' = center (top_face c)
255 v2' = Face.v1 (top_face c)
256 v3' = Face.v2 (top_face c)
257 fv' = rotate cwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
259 tetrahedron6 :: Cube -> Tetrahedron
261 Tetrahedron fv' v0' v1' v2' v3'
264 v1' = center (top_face c)
265 v2' = Face.v2 (top_face c)
266 v3' = Face.v3 (top_face c)
267 fv' = rotate cwy $ rotate cwz
269 $ Tetrahedron.fv (tetrahedron0 c)
271 tetrahedron7 :: Cube -> Tetrahedron
273 Tetrahedron fv' v0' v1' v2' v3'
276 v1' = center (top_face c)
277 v2' = Face.v3 (top_face c)
278 v3' = Face.v0 (top_face c)
279 fv' = rotate cwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
281 tetrahedron8 :: Cube -> Tetrahedron
283 Tetrahedron fv' v0' v1' v2' v3'
286 v1' = center (back_face c)
287 v2' = Face.v0 (back_face c)
288 v3' = Face.v1 (back_face c)
289 fv' = rotate cwy $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
291 tetrahedron9 :: Cube -> Tetrahedron
293 Tetrahedron fv' v0' v1' v2' v3'
296 v1' = center (back_face c)
297 v2' = Face.v1 (back_face c)
298 v3' = Face.v2 (back_face c)
299 fv' = rotate cwy $ rotate cwy
301 $ Tetrahedron.fv (tetrahedron0 c)
303 tetrahedron10 :: Cube -> Tetrahedron
305 Tetrahedron fv' v0' v1' v2' v3'
308 v1' = center (back_face c)
309 v2' = Face.v2 (back_face c)
310 v3' = Face.v3 (back_face c)
311 fv' = rotate cwy $ rotate cwy
314 $ Tetrahedron.fv (tetrahedron0 c)
317 tetrahedron11 :: Cube -> Tetrahedron
319 Tetrahedron fv' v0' v1' v2' v3'
322 v1' = center (back_face c)
323 v2' = Face.v3 (back_face c)
324 v3' = Face.v0 (back_face c)
325 fv' = rotate cwy $ rotate cwy
327 $ Tetrahedron.fv (tetrahedron0 c)
330 tetrahedron12 :: Cube -> Tetrahedron
332 Tetrahedron fv' v0' v1' v2' v3'
335 v1' = center (down_face c)
336 v2' = Face.v0 (down_face c)
337 v3' = Face.v1 (down_face c)
338 fv' = rotate ccwy (Tetrahedron.fv (tetrahedron0 c))
341 tetrahedron13 :: Cube -> Tetrahedron
343 Tetrahedron fv' v0' v1' v2' v3'
346 v1' = center (down_face c)
347 v2' = Face.v1 (down_face c)
348 v3' = Face.v2 (down_face c)
349 fv' = rotate ccwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
352 tetrahedron14 :: Cube -> Tetrahedron
354 Tetrahedron fv' v0' v1' v2' v3'
357 v1' = center (down_face c)
358 v2' = Face.v2 (down_face c)
359 v3' = Face.v3 (down_face c)
360 fv' = rotate ccwy $ rotate ccwz
362 $ Tetrahedron.fv (tetrahedron0 c)
365 tetrahedron15 :: Cube -> Tetrahedron
367 Tetrahedron fv' v0' v1' v2' v3'
370 v1' = center (down_face c)
371 v2' = Face.v3 (down_face c)
372 v3' = Face.v0 (down_face c)
373 fv' = rotate ccwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
376 tetrahedron16 :: Cube -> Tetrahedron
378 Tetrahedron fv' v0' v1' v2' v3'
381 v1' = center (right_face c)
382 v2' = Face.v0 (right_face c)
383 v3' = Face.v1 (right_face c)
384 fv' = rotate ccwz (Tetrahedron.fv (tetrahedron0 c))
387 tetrahedron17 :: Cube -> Tetrahedron
389 Tetrahedron fv' v0' v1' v2' v3'
392 v1' = center (right_face c)
393 v2' = Face.v1 (right_face c)
394 v3' = Face.v2 (right_face c)
395 fv' = rotate ccwz $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
398 tetrahedron18 :: Cube -> Tetrahedron
400 Tetrahedron fv' v0' v1' v2' v3'
403 v1' = center (right_face c)
404 v2' = Face.v2 (right_face c)
405 v3' = Face.v3 (right_face c)
406 fv' = rotate ccwz $ rotate cwy
408 $ Tetrahedron.fv (tetrahedron0 c)
411 tetrahedron19 :: Cube -> Tetrahedron
413 Tetrahedron fv' v0' v1' v2' v3'
416 v1' = center (right_face c)
417 v2' = Face.v3 (right_face c)
418 v3' = Face.v0 (right_face c)
419 fv' = rotate ccwz $ rotate ccwy
420 $ Tetrahedron.fv (tetrahedron0 c)
423 tetrahedron20 :: Cube -> Tetrahedron
425 Tetrahedron fv' v0' v1' v2' v3'
428 v1' = center (left_face c)
429 v2' = Face.v0 (left_face c)
430 v3' = Face.v1 (left_face c)
431 fv' = rotate cwz (Tetrahedron.fv (tetrahedron0 c))
434 tetrahedron21 :: Cube -> Tetrahedron
436 Tetrahedron fv' v0' v1' v2' v3'
439 v1' = center (left_face c)
440 v2' = Face.v1 (left_face c)
441 v3' = Face.v2 (left_face c)
442 fv' = rotate cwz $ rotate ccwy $ Tetrahedron.fv (tetrahedron0 c)
445 tetrahedron22 :: Cube -> Tetrahedron
447 Tetrahedron fv' v0' v1' v2' v3'
450 v1' = center (left_face c)
451 v2' = Face.v2 (left_face c)
452 v3' = Face.v3 (left_face c)
453 fv' = rotate cwz $ rotate ccwy
455 $ Tetrahedron.fv (tetrahedron0 c)
458 tetrahedron23 :: Cube -> Tetrahedron
460 Tetrahedron fv' v0' v1' v2' v3'
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 (tetrahedron0 c)
470 tetrahedra :: Cube -> [Tetrahedron]
497 -- | All completely contained in the front half of the cube.
498 front_half_tetrahedra :: Cube -> [Tetrahedron]
499 front_half_tetrahedra c =
510 -- | All tetrahedra completely contained in the top half of the cube.
511 top_half_tetrahedra :: Cube -> [Tetrahedron]
512 top_half_tetrahedra c =
523 -- | All tetrahedra completely contained in the back half of the cube.
524 back_half_tetrahedra :: Cube -> [Tetrahedron]
525 back_half_tetrahedra c =
536 -- | All tetrahedra completely contained in the down half of the cube.
537 down_half_tetrahedra :: Cube -> [Tetrahedron]
538 down_half_tetrahedra c =
549 -- | All tetrahedra completely contained in the right half of the cube.
550 right_half_tetrahedra :: Cube -> [Tetrahedron]
551 right_half_tetrahedra c =
562 -- | All tetrahedra completely contained in the left half of the cube.
563 left_half_tetrahedra :: Cube -> [Tetrahedron]
564 left_half_tetrahedra c =
575 in_top_half :: Cube -> Point -> Bool
576 in_top_half c (_,_,z) =
577 distance_from_top <= distance_from_bottom
579 distance_from_top = abs $ (zmax c) - z
580 distance_from_bottom = abs $ (zmin c) - z
582 in_front_half :: Cube -> Point -> Bool
583 in_front_half c (x,_,_) =
584 distance_from_front <= distance_from_back
586 distance_from_front = abs $ (xmin c) - x
587 distance_from_back = abs $ (xmax c) - x
590 in_left_half :: Cube -> Point -> Bool
591 in_left_half c (_,y,_) =
592 distance_from_left <= distance_from_right
594 distance_from_left = abs $ (ymin c) - y
595 distance_from_right = abs $ (ymax c) - y
598 -- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
599 -- contain the given 'Point'. This should be faster than checking
600 -- every tetrahedron individually, since we determine which half
601 -- (hemisphere?) of the cube the point lies in three times: once in
602 -- each dimension. This allows us to eliminate non-candidates
605 -- This can throw an exception, but the use of 'head' might
606 -- save us some unnecessary computations.
608 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
609 find_containing_tetrahedron c p =
610 head containing_tetrahedra
612 candidates = tetrahedra c
614 if (in_front_half c p) then
615 back_half_tetrahedra c
617 front_half_tetrahedra c
619 candidates_x = candidates \\ non_candidates_x
622 if (in_left_half c p) then
623 right_half_tetrahedra c
625 left_half_tetrahedra c
627 candidates_xy = candidates_x \\ non_candidates_y
630 if (in_top_half c p) then
631 down_half_tetrahedra c
633 top_half_tetrahedra c
635 candidates_xyz = candidates_xy \\ non_candidates_z
637 contains_our_point = flip contains_point p
638 containing_tetrahedra = filter contains_our_point candidates_xyz