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
-- | The Face module just contains the definition of the 'Face' data
-- type and its two typeclass instances.
-module Face
- ( Face(..) )
+module Face (
+ Face(..),
+ center
+ )
where
import Point
-import ThreeDimensional
data Face = Face { v0 :: !Point,
v1 :: !Point,
deriving (Eq)
instance Show Face where
- show f = "Face:\n" ++
- " v0: " ++ (show (v0 f)) ++ "\n" ++
- " v1: " ++ (show (v1 f)) ++ "\n" ++
- " v2: " ++ (show (v2 f)) ++ "\n" ++
- " v3: " ++ (show (v3 f)) ++ "\n"
+ show (Face v0' v1' v2' v3') =
+ "Face:\n" ++
+ " v0: " ++ (show v0') ++ "\n" ++
+ " v1: " ++ (show v1') ++ "\n" ++
+ " v2: " ++ (show v2') ++ "\n" ++
+ " v3: " ++ (show v3') ++ "\n"
--- | The 'Face' type is an instance of 'ThreeDimensional' so that we
--- can call the 'center' function on it. This is useful because the
--- center of a face is always a vertex of a tetrahedron.
-instance ThreeDimensional Face where
- -- | Since a face is square, we can just average the four vertices
- -- to find the center.
- center f = ((v0 f) + (v1 f) + (v2 f) + (v3 f)) `scale` (1/4)
-
- -- | It's possible to implement this, but it hasn't been done
- -- yet. A face will contain a point if the point lies in the same
- -- plane as the vertices of the face, and if it falls on the
- -- correct side of the four sides of the face.
- contains_point _ _ = False
+-- | Returns the center of the given face. Since a face is just
+-- square, we can average the four vertices to find its center. This
+-- is useful because the center of a face is always a vertex of a
+-- tetrahedron.
+center :: Face -> Point
+center (Face v0' v1' v2' v3') =
+ (v0' + v1' + v2' + v3') `scale` (1/4)
import FunctionValues (make_values, value_at)
import Point (Point(..))
import ScaleFactor (ScaleFactor)
-import Tetrahedron (Tetrahedron, c, polynomial, v0, v1, v2, v3)
-import ThreeDimensional (ThreeDimensional(..))
+import Tetrahedron (
+ Tetrahedron(v0,v1,v2,v3),
+ c,
+ contains_point,
+ polynomial,
+ )
import Values (Values3D, dims, empty3d, zoom_shape)
b1, -- Cube test
b2, -- Cube test
b3, -- Cube test
+ barycenter,
c,
+ contains_point,
polynomial,
tetrahedron_properties,
tetrahedron_tests,
import Misc (factorial)
import Point (Point(..), scale)
import RealFunction (RealFunction, cmult, fexp)
-import ThreeDimensional (ThreeDimensional(..))
data Tetrahedron =
Tetrahedron { function_values :: FunctionValues,
" v3: " ++ (show (v3 t)) ++ "\n"
-instance ThreeDimensional Tetrahedron where
- center (Tetrahedron _ v0' v1' v2' v3' _) =
- (v0' + v1' + v2' + v3') `scale` (1/4)
-
- -- contains_point is only used in tests.
- contains_point t p0 =
- b0_unscaled `nearly_ge` 0 &&
- b1_unscaled `nearly_ge` 0 &&
- b2_unscaled `nearly_ge` 0 &&
- b3_unscaled `nearly_ge` 0
+-- | Find the barycenter of the given tetrahedron.
+-- We just average the four vertices.
+barycenter :: Tetrahedron -> Point
+barycenter (Tetrahedron _ v0' v1' v2' v3' _) =
+ (v0' + v1' + v2' + v3') `scale` (1/4)
+
+-- | A point is internal to a tetrahedron if all of its barycentric
+-- coordinates with respect to that tetrahedron are non-negative.
+contains_point :: Tetrahedron -> Point -> Bool
+contains_point t p0 =
+ b0_unscaled `nearly_ge` 0 &&
+ b1_unscaled `nearly_ge` 0 &&
+ b2_unscaled `nearly_ge` 0 &&
+ b3_unscaled `nearly_ge` 0
+ where
+ -- Drop the useless division and volume calculation that we
+ -- would do if we used the regular b0,..b3 functions.
+ b0_unscaled :: Double
+ b0_unscaled = volume inner_tetrahedron
where
- -- Drop the useless division and volume calculation that we
- -- would do if we used the regular b0,..b3 functions.
- b0_unscaled :: Double
- b0_unscaled = volume inner_tetrahedron
- where inner_tetrahedron = t { v0 = p0 }
-
- b1_unscaled :: Double
- b1_unscaled = volume inner_tetrahedron
- where inner_tetrahedron = t { v1 = p0 }
-
- b2_unscaled :: Double
- b2_unscaled = volume inner_tetrahedron
- where inner_tetrahedron = t { v2 = p0 }
-
- b3_unscaled :: Double
- b3_unscaled = volume inner_tetrahedron
- where inner_tetrahedron = t { v3 = p0 }
+ inner_tetrahedron = t { v0 = p0 }
+
+ b1_unscaled :: Double
+ b1_unscaled = volume inner_tetrahedron
+ where inner_tetrahedron = t { v1 = p0 }
+
+ b2_unscaled :: Double
+ b2_unscaled = volume inner_tetrahedron
+ where inner_tetrahedron = t { v2 = p0 }
+
+ b3_unscaled :: Double
+ b3_unscaled = volume inner_tetrahedron
+ where inner_tetrahedron = t { v3 = p0 }
+
{-# INLINE polynomial #-}
polynomial :: Tetrahedron -> (RealFunction Point)
+++ /dev/null
-module ThreeDimensional
- ( ThreeDimensional(..) )
-where
-
-import Point(Point)
-
-class ThreeDimensional a where
- center :: a -> Point
- contains_point :: a -> Point -> Bool
projects / spline3.git / commitdiff