import Cardinal
import Comparisons ((~=), (~~=))
-import qualified Face (Face(Face, v0, v1, v2, v3))
+import qualified Face (Face(..), center)
import FunctionValues (FunctionValues, eval, rotate)
import Misc (all_equal, disjoint)
import Point (Point(..), dot)
-import Tetrahedron (Tetrahedron(..), c, volume)
-import ThreeDimensional
+import Tetrahedron (Tetrahedron(..), barycenter, c, volume)
data Cube = Cube { h :: !Double,
i :: !Int,
k' = fromIntegral (k cube) :: Double
delta = h cube
-instance ThreeDimensional Cube where
- -- | The center of Cube_ijk coincides with v_ijk at
- -- (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
- center cube = Point x y z
- where
- delta = h cube
- i' = fromIntegral (i cube) :: Double
- j' = fromIntegral (j cube) :: Double
- k' = fromIntegral (k cube) :: Double
- x = delta * i'
- y = delta * j'
- z = delta * k'
-
- -- | It's easy to tell if a point is within a cube; just make sure
- -- that it falls on the proper side of each of the cube's faces.
- contains_point cube (Point x y z)
- | x < (xmin cube) = False
- | x > (xmax cube) = False
- | y < (ymin cube) = False
- | y > (ymax cube) = False
- | z < (zmin cube) = False
- | z > (zmax cube) = False
- | otherwise = True
+-- | The center of Cube_ijk coincides with v_ijk at
+-- (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
+center :: Cube -> Point
+center cube =
+ Point x y z
+ where
+ delta = h cube
+ i' = fromIntegral (i cube) :: Double
+ j' = fromIntegral (j cube) :: Double
+ k' = fromIntegral (k cube) :: Double
+ x = delta * i'
+ y = delta * j'
+ z = delta * k'
-- Face stuff.
where
v0' = center cube
ff = front_face cube
- v1' = center ff
+ v1' = Face.center ff
v2' = Face.v0 ff
v3' = Face.v1 ff
vol = tetrahedra_volume cube
where
v0' = center cube
ff = front_face cube
- v1' = center ff
+ v1' = Face.center ff
v2' = Face.v1 ff
v3' = Face.v2 ff
fv' = rotate ccwx (fv cube)
where
v0' = center cube
ff = front_face cube
- v1' = center ff
+ v1' = Face.center ff
v2' = Face.v2 ff
v3' = Face.v3 ff
fv' = rotate ccwx $ rotate ccwx $ fv cube
where
v0' = center cube
ff = front_face cube
- v1' = center ff
+ v1' = Face.center ff
v2' = Face.v3 ff
v3' = Face.v0 ff
fv' = rotate cwx (fv cube)
where
v0' = center cube
tf = top_face cube
- v1' = center tf
+ v1' = Face.center tf
v2' = Face.v0 tf
v3' = Face.v1 tf
fv' = rotate cwy (fv cube)
where
v0' = center cube
tf = top_face cube
- v1' = center tf
+ v1' = Face.center tf
v2' = Face.v1 tf
v3' = Face.v2 tf
fv' = rotate cwy $ rotate cwz $ fv cube
where
v0' = center cube
tf = top_face cube
- v1' = center tf
+ v1' = Face.center tf
v2' = Face.v2 tf
v3' = Face.v3 tf
fv' = rotate cwy $ rotate cwz
where
v0' = center cube
tf = top_face cube
- v1' = center tf
+ v1' = Face.center tf
v2' = Face.v3 tf
v3' = Face.v0 tf
fv' = rotate cwy $ rotate ccwz $ fv cube
where
v0' = center cube
bf = back_face cube
- v1' = center bf
+ v1' = Face.center bf
v2' = Face.v0 bf
v3' = Face.v1 bf
fv' = rotate cwy $ rotate cwy $ fv cube
where
v0' = center cube
bf = back_face cube
- v1' = center bf
+ v1' = Face.center bf
v2' = Face.v1 bf
v3' = Face.v2 bf
fv' = rotate cwy $ rotate cwy
where
v0' = center cube
bf = back_face cube
- v1' = center bf
+ v1' = Face.center bf
v2' = Face.v2 bf
v3' = Face.v3 bf
fv' = rotate cwy $ rotate cwy
where
v0' = center cube
bf = back_face cube
- v1' = center bf
+ v1' = Face.center bf
v2' = Face.v3 bf
v3' = Face.v0 bf
fv' = rotate cwy $ rotate cwy
where
v0' = center cube
df = down_face cube
- v1' = center df
+ v1' = Face.center df
v2' = Face.v0 df
v3' = Face.v1 df
fv' = rotate ccwy $ fv cube
where
v0' = center cube
df = down_face cube
- v1' = center df
+ v1' = Face.center df
v2' = Face.v1 df
v3' = Face.v2 df
fv' = rotate ccwy $ rotate ccwz $ fv cube
where
v0' = center cube
df = down_face cube
- v1' = center df
+ v1' = Face.center df
v2' = Face.v2 df
v3' = Face.v3 df
fv' = rotate ccwy $ rotate ccwz
where
v0' = center cube
df = down_face cube
- v1' = center df
+ v1' = Face.center df
v2' = Face.v3 df
v3' = Face.v0 df
fv' = rotate ccwy $ rotate cwz $ fv cube
where
v0' = center cube
rf = right_face cube
- v1' = center rf
+ v1' = Face.center rf
v2' = Face.v0 rf
v3' = Face.v1 rf
fv' = rotate ccwz $ fv cube
where
v0' = center cube
rf = right_face cube
- v1' = center rf
+ v1' = Face.center rf
v2' = Face.v1 rf
v3' = Face.v2 rf
fv' = rotate ccwz $ rotate cwy $ fv cube
where
v0' = center cube
rf = right_face cube
- v1' = center rf
+ v1' = Face.center rf
v2' = Face.v2 rf
v3' = Face.v3 rf
fv' = rotate ccwz $ rotate cwy
where
v0' = center cube
rf = right_face cube
- v1' = center rf
+ v1' = Face.center rf
v2' = Face.v3 rf
v3' = Face.v0 rf
fv' = rotate ccwz $ rotate ccwy
where
v0' = center cube
lf = left_face cube
- v1' = center lf
+ v1' = Face.center lf
v2' = Face.v0 lf
v3' = Face.v1 lf
fv' = rotate cwz $ fv cube
where
v0' = center cube
lf = left_face cube
- v1' = center lf
+ v1' = Face.center lf
v2' = Face.v1 lf
v3' = Face.v2 lf
fv' = rotate cwz $ rotate ccwy $ fv cube
where
v0' = center cube
lf = left_face cube
- v1' = center lf
+ v1' = Face.center lf
v2' = Face.v2 lf
v3' = Face.v3 lf
fv' = rotate cwz $ rotate ccwy
where
v0' = center cube
lf = left_face cube
- v1' = center lf
+ v1' = Face.center lf
v2' = Face.v3 lf
v3' = Face.v0 lf
fv' = rotate cwz $ rotate cwy
top_half = in_top_half cube p
left_half = in_left_half cube p
+ candidates :: V.Vector Tetrahedron
candidates =
if front_half then
-- Use the dot product instead of Euclidean distance here to save
-- a sqrt(). So, "distances" below really means "distances
-- squared."
- distances = V.map ((dot p) . center) candidates
+ distances :: V.Vector Double
+ distances = V.map ((dot p) . barycenter) candidates
+
+ shortest_distance :: Double
shortest_distance = V.minimum distances
+
+ -- Compute the index of the tetrahedron with the center closest to
+ -- p. This is a bad algorithm, but don't change it! If you make it
+ -- smarter by finding the index of shortest_distance in distances
+ -- (this should give the same answer and avoids recomputing the
+ -- dot product), the program gets slower. Seriously!
+ lucky_idx :: Maybe Int
lucky_idx = V.findIndex
- (\t -> (center t) `dot` p == shortest_distance)
+ (\t -> (barycenter t) `dot` p == shortest_distance)
candidates
where
delta = h cube
--- | All tetrahedron should have their v0 located at the center of the cube.
+-- | All tetrahedron should have their v0 located at the center of the
+-- cube.
prop_v0_all_equal :: Cube -> Bool
prop_v0_all_equal cube = (v0 t0) == (v0 t1)
where
projects / spline3.git / blobdiff