src/Cube.hs

   1 module Cube
   2 where
   3
   4 import Data.Maybe (fromJust)
   5 import qualified Data.Vector as V (
   6   Vector,
   7   findIndex,
   8   map,
   9   minimum,
  10   singleton,
  11   snoc,
  12   unsafeIndex
  13   )
  14 import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
  15
  16 import Cardinal
  17 import qualified Face (Face(Face, v0, v1, v2, v3))
  18 import FunctionValues
  19 import Point
  20 import Tetrahedron hiding (c)
  21 import ThreeDimensional
  22
  23 data Cube = Cube { h :: Double,
  24                    i :: Int,
  25                    j :: Int,
  26                    k :: Int,
  27                    fv :: FunctionValues,
  28                    tetrahedra_volume :: Double }
  29             deriving (Eq)
  30
  31
  32 instance Arbitrary Cube where
  33     arbitrary = do
  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)
  41         where
  42           coordmin = -268435456 -- -(2^29 / 2)
  43           coordmax = 268435456  -- +(2^29 / 2)
  44
  45
  46 instance Show Cube where
  47     show c =
  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"
  58         where
  59           subscript =
  60               (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c))
  61
  62
  63 -- | Returns an empty 'Cube'.
  64 empty_cube :: Cube
  65 empty_cube = Cube 0 0 0 0 empty_values 0
  66
  67
  68 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
  69 --   p. 76.
  70 xmin :: Cube -> Double
  71 xmin c = (2*i' - 1)*delta / 2
  72     where
  73       i' = fromIntegral (i c) :: Double
  74       delta = h c
  75
  76 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
  77 --   p. 76.
  78 xmax :: Cube -> Double
  79 xmax c = (2*i' + 1)*delta / 2
  80     where
  81       i' = fromIntegral (i c) :: Double
  82       delta = h c
  83
  84 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
  85 --   p. 76.
  86 ymin :: Cube -> Double
  87 ymin c = (2*j' - 1)*delta / 2
  88     where
  89       j' = fromIntegral (j c) :: Double
  90       delta = h c
  91
  92 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
  93 --   p. 76.
  94 ymax :: Cube -> Double
  95 ymax c = (2*j' + 1)*delta / 2
  96     where
  97       j' = fromIntegral (j c) :: Double
  98       delta = h c
  99
 100 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
 101 --   p. 76.
 102 zmin :: Cube -> Double
 103 zmin c = (2*k' - 1)*delta / 2
 104     where
 105       k' = fromIntegral (k c) :: Double
 106       delta = h c
 107
 108 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
 109 --   p. 76.
 110 zmax :: Cube -> Double
 111 zmax c = (2*k' + 1)*delta / 2
 112     where
 113       k' = fromIntegral (k c) :: Double
 114       delta = h c
 115
 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.
 119     center c = (x, y, z)
 120            where
 121              delta = h c
 122              i' = fromIntegral (i c) :: Double
 123              j' = fromIntegral (j c) :: Double
 124              k' = fromIntegral (k c) :: Double
 125              x = delta * i'
 126              y = delta * j'
 127              z = delta * k'
 128
 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
 138         | otherwise = True
 139
 140
 141
 142 -- Face stuff.
 143
 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'
 147     where
 148       delta = (1/2)*(h c)
 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)
 153
 154
 155
 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'
 159     where
 160       delta = (1/2)*(h c)
 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)
 165
 166
 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'
 170     where
 171       delta = (1/2)*(h c)
 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)
 176
 177
 178
 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'
 182     where
 183       delta = (1/2)*(h c)
 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)
 188
 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'
 192     where
 193       delta = (1/2)*(h c)
 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)
 198
 199
 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'
 203     where
 204       delta = (1/2)*(h c)
 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)
 209
 210
 211 tetrahedron :: Cube -> Int -> Tetrahedron
 212
 213 tetrahedron c 0 =
 214     Tetrahedron (Cube.fv c) v0' v1' v2' v3' vol 0
 215     where
 216       v0' = center c
 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
 221
 222 tetrahedron c 1 =
 223     Tetrahedron fv' v0' v1' v2' v3' vol 1
 224     where
 225       v0' = center c
 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
 231
 232 tetrahedron c 2 =
 233     Tetrahedron fv' v0' v1' v2' v3' vol 2
 234     where
 235       v0' = center c
 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
 241
 242 tetrahedron c 3 =
 243     Tetrahedron fv' v0' v1' v2' v3' vol 3
 244     where
 245       v0' = center c
 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
 251
 252 tetrahedron c 4 =
 253     Tetrahedron fv' v0' v1' v2' v3' vol 4
 254     where
 255       v0' = center c
 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
 261
 262 tetrahedron c 5 =
 263     Tetrahedron fv' v0' v1' v2' v3' vol 5
 264     where
 265       v0' = center c
 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
 271
 272 tetrahedron c 6 =
 273     Tetrahedron fv' v0' v1' v2' v3' vol 6
 274     where
 275       v0' = center c
 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
 280                        $ rotate cwz
 281                        $ Tetrahedron.fv (tetrahedron c 0)
 282       vol = tetrahedra_volume c
 283
 284 tetrahedron c 7 =
 285     Tetrahedron fv' v0' v1' v2' v3' vol 7
 286     where
 287       v0' = center c
 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
 293
 294 tetrahedron c 8 =
 295     Tetrahedron fv' v0' v1' v2' v3' vol 8
 296     where
 297       v0' = center c
 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
 303
 304 tetrahedron c 9 =
 305     Tetrahedron fv' v0' v1' v2' v3' vol 9
 306     where
 307       v0' = center c
 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
 312                        $ rotate cwx
 313                        $ Tetrahedron.fv (tetrahedron c 0)
 314       vol = tetrahedra_volume c
 315
 316 tetrahedron c 10 =
 317     Tetrahedron fv' v0' v1' v2' v3' vol 10
 318     where
 319       v0' = center c
 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
 324                        $ rotate cwx
 325                        $ rotate cwx
 326                        $ Tetrahedron.fv (tetrahedron c 0)
 327
 328       vol = tetrahedra_volume c
 329
 330 tetrahedron c 11 =
 331     Tetrahedron fv' v0' v1' v2' v3' vol 11
 332     where
 333       v0' = center c
 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
 338                        $ rotate ccwx
 339                        $ Tetrahedron.fv (tetrahedron c 0)
 340       vol = tetrahedra_volume c
 341
 342 tetrahedron c 12 =
 343     Tetrahedron fv' v0' v1' v2' v3' vol 12
 344     where
 345       v0' = center c
 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
 351
 352 tetrahedron c 13 =
 353     Tetrahedron fv' v0' v1' v2' v3' vol 13
 354     where
 355       v0' = center c
 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
 361
 362 tetrahedron c 14 =
 363     Tetrahedron fv' v0' v1' v2' v3' vol 14
 364     where
 365       v0' = center c
 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
 370                         $ rotate ccwz
 371                         $ Tetrahedron.fv (tetrahedron c 0)
 372       vol = tetrahedra_volume c
 373
 374 tetrahedron c 15 =
 375     Tetrahedron fv' v0' v1' v2' v3' vol 15
 376     where
 377       v0' = center c
 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
 383
 384 tetrahedron c 16 =
 385     Tetrahedron fv' v0' v1' v2' v3' vol 16
 386     where
 387       v0' = center c
 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
 393
 394 tetrahedron c 17 =
 395     Tetrahedron fv' v0' v1' v2' v3' vol 17
 396     where
 397       v0' = center c
 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
 403
 404 tetrahedron c 18 =
 405     Tetrahedron fv' v0' v1' v2' v3' vol 18
 406     where
 407       v0' = center c
 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
 412                         $ rotate cwy
 413                         $ Tetrahedron.fv (tetrahedron c 0)
 414       vol = tetrahedra_volume c
 415
 416 tetrahedron c 19 =
 417     Tetrahedron fv' v0' v1' v2' v3' vol 19
 418     where
 419       v0' = center c
 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
 426
 427 tetrahedron c 20 =
 428     Tetrahedron fv' v0' v1' v2' v3' vol 20
 429     where
 430       v0' = center c
 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
 436
 437 tetrahedron c 21 =
 438     Tetrahedron fv' v0' v1' v2' v3' vol 21
 439     where
 440       v0' = center c
 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
 446
 447 tetrahedron c 22 =
 448     Tetrahedron fv' v0' v1' v2' v3' vol 22
 449     where
 450       v0' = center c
 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
 455                        $ rotate ccwy
 456                        $ Tetrahedron.fv (tetrahedron c 0)
 457       vol = tetrahedra_volume c
 458
 459 tetrahedron c 23 =
 460     Tetrahedron fv' v0' v1' v2' v3' vol 23
 461     where
 462       v0' = center c
 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
 469
 470 -- Feels dirty, but whatever.
 471 tetrahedron _ _ = error "asked for a nonexistent tetrahedron"
 472
 473
 474 -- Only used in tests, so we don't need the added speed
 475 -- of Data.Vector.
 476 tetrahedra :: Cube -> [Tetrahedron]
 477 tetrahedra c = [ tetrahedron c n | n <- [0..23] ]
 478
 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`
 486     (tetrahedron c 21)
 487
 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`
 495     (tetrahedron c 21)
 496
 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`
 504     (tetrahedron c 19)
 505
 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`
 513     (tetrahedron c 19)
 514
 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`
 522     (tetrahedron c 21)
 523
 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`
 531     (tetrahedron c 23)
 532
 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`
 540     (tetrahedron c 17)
 541
 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`
 549     (tetrahedron c 18)
 550
 551 in_top_half :: Cube -> Point -> Bool
 552 in_top_half c (_,_,z) =
 553   distance_from_top <= distance_from_bottom
 554   where
 555     distance_from_top = abs $ (zmax c) - z
 556     distance_from_bottom = abs $ (zmin c) - z
 557
 558 in_front_half :: Cube -> Point -> Bool
 559 in_front_half c (x,_,_) =
 560     distance_from_front <= distance_from_back
 561   where
 562     distance_from_front = abs $ (xmin c) - x
 563     distance_from_back = abs $ (xmax c) - x
 564
 565
 566 in_left_half :: Cube -> Point -> Bool
 567 in_left_half c (_,y,_) =
 568     distance_from_left <= distance_from_right
 569   where
 570     distance_from_left = abs $ (ymin c) - y
 571     distance_from_right = abs $ (ymax c) - y
 572
 573
 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
 579 --   quickly.
 580 --
 581 --   This can throw an exception, but the use of 'head' might
 582 --   save us some unnecessary computations.
 583 --
 584 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
 585 find_containing_tetrahedron c p =
 586   candidates `V.unsafeIndex` (fromJust lucky_idx)
 587   where
 588     front_half = in_front_half c p
 589     top_half = in_top_half c p
 590     left_half = in_left_half c p
 591
 592     candidates =
 593       if front_half then
 594
 595         if left_half then
 596           if top_half then
 597             front_left_top_tetrahedra c
 598           else
 599             front_left_down_tetrahedra c
 600         else
 601           if top_half then
 602             front_right_top_tetrahedra c
 603           else
 604             front_right_down_tetrahedra c
 605
 606       else -- bottom half
 607
 608         if left_half then
 609           if top_half then
 610             back_left_top_tetrahedra c
 611           else
 612             back_left_down_tetrahedra c
 613         else
 614           if top_half then
 615             back_right_top_tetrahedra c
 616           else
 617             back_right_down_tetrahedra c
 618
 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)
 625                   candidates