Rosetta Code

Catmull–Clark subdivision surface: Difference between revisions

Page Discussion

Catmull–Clark subdivision surface (view source)

Revision as of 05:23, 5 October 2017

13,061 bytes added ,  6 years ago
Add Haskell
mNo edit summary
(Add Haskell)
Line 139:
return nm;
}</lang>
 
=={{header|Haskell}}==
<lang haskell>{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}
 
import Data.Array
import Data.Foldable (length, concat, sum)
import Data.List (genericLength)
import Data.Maybe (mapMaybe)
import Prelude hiding (length, concat, sum)
import qualified Data.Map.Strict as Map
 
{-
A SimpleMesh consists of only vertices and faces that refer to them.
A Mesh extends the SimpleMesh to contain edges as well as references to
adjoining mesh components for each other component, such as a vertex
also contains what faces it belongs to.
An isolated edge can be represented as a degenerate face with 2 vertices.
Faces with 0 or 1 vertices can be thrown out, as they do not contribute to
the result (they can also propagate NaNs).
-}
 
newtype VertexId = VertexId { getVertexId :: Int } deriving (Ix, Ord, Eq, Show)
newtype EdgeId = EdgeId { getEdgeId :: Int } deriving (Ix, Ord, Eq, Show)
newtype FaceId = FaceId { getFaceId :: Int } deriving (Ix, Ord, Eq, Show)
 
data Vertex a = Vertex
{ vertexPoint :: a
, vertexEdges :: [EdgeId]
, vertexFaces :: [FaceId]
} deriving Show
 
data Edge = Edge
{ edgeVertexA :: VertexId
, edgeVertexB :: VertexId
, edgeFaces :: [FaceId]
} deriving Show
 
data Face = Face
{ faceVertices :: [VertexId]
, faceEdges :: [EdgeId]
} deriving Show
 
type VertexArray a = Array VertexId (Vertex a)
type EdgeArray = Array EdgeId Edge
type FaceArray = Array FaceId Face
 
data Mesh a = Mesh
{ meshVertices :: VertexArray a
, meshEdges :: EdgeArray
, meshFaces :: FaceArray
} deriving Show
 
data SimpleVertex a = SimpleVertex { sVertexPoint :: a } deriving Show
data SimpleFace = SimpleFace { sFaceVertices :: [VertexId] } deriving Show
 
type SimpleVertexArray a = Array VertexId (SimpleVertex a)
type SimpleFaceArray = Array FaceId SimpleFace
 
data SimpleMesh a = SimpleMesh
{ sMeshVertices :: SimpleVertexArray a
, sMeshFaces :: SimpleFaceArray
} deriving Show
 
-- Generic helpers.
fmap1 :: Functor f => (t -> a -> b) -> (t -> f a) -> t -> f b
fmap1 g h x = fmap (g x) (h x)
 
aZipWith :: Ix i1 => (a -> b -> e) -> Array i1 a -> Array i b -> Array i1 e
aZipWith f a b = listArray (bounds a) $ zipWith f (elems a) (elems b)
 
average :: (Foldable f, Fractional a) => f a -> a
average xs = (sum xs) / (fromIntegral $ length xs)
 
-- Intermediary point types for ultimately converting into a point `a`.
newtype FacePoint a = FacePoint { getFacePoint :: a } deriving Show
newtype EdgeCenterPoint a = EdgeCenterPoint { getEdgeCenterPoint :: a } deriving Show
newtype EdgePoint a = EdgePoint { getEdgePoint :: a } deriving Show
newtype VertexPoint a = VertexPoint { getVertexPoint :: a } deriving Show
 
type FacePointArray a = Array FaceId (FacePoint a)
type EdgePointArray a = Array EdgeId (EdgePoint a)
type EdgeCenterPointArray a = Array EdgeId (EdgeCenterPoint a)
type IsEdgeHoleArray = Array EdgeId Bool
type VertexPointArray a = Array VertexId (VertexPoint a)
 
-- Subdivision helpers.
facePoint :: Fractional a => Mesh a -> Face -> FacePoint a
facePoint mesh = FacePoint . average . (fmap $ vertexPointById mesh) . faceVertices
 
allFacePoints :: Fractional a => Mesh a -> FacePointArray a
allFacePoints = fmap1 facePoint meshFaces
 
vertexPointById :: Mesh a -> VertexId -> a
vertexPointById mesh = vertexPoint . (meshVertices mesh !)
 
edgeCenterPoint :: Fractional a => Mesh a -> Edge -> EdgeCenterPoint a
edgeCenterPoint mesh (Edge ea eb _)
= EdgeCenterPoint . average $ fmap (vertexPointById mesh) [ea, eb]
 
allEdgeCenterPoints :: Fractional a => Mesh a -> EdgeCenterPointArray a
allEdgeCenterPoints = fmap1 edgeCenterPoint meshEdges
 
allIsEdgeHoles :: Mesh a -> IsEdgeHoleArray
allIsEdgeHoles = fmap ((< 2) . length . edgeFaces) . meshEdges
 
edgePoint :: Fractional a => Edge -> FacePointArray a -> EdgeCenterPoint a -> EdgePoint a
edgePoint (Edge _ _ [_]) _ (EdgeCenterPoint ecp) = EdgePoint ecp
edgePoint (Edge _ _ faceIds) facePoints (EdgeCenterPoint ecp)
= EdgePoint $ average [ecp, average $ fmap (getFacePoint . (facePoints !)) faceIds]
 
allEdgePoints :: Fractional a => Mesh a -> FacePointArray a -> EdgeCenterPointArray a -> EdgePointArray a
allEdgePoints mesh fps ecps = aZipWith (\e ecp -> edgePoint e fps ecp) (meshEdges mesh) ecps
 
vertexPoint' :: Fractional a => Vertex a -> FacePointArray a -> EdgeCenterPointArray a -> IsEdgeHoleArray -> VertexPoint a
vertexPoint' vertex facePoints ecps iehs
| length faceIds == length edgeIds = VertexPoint newCoords
| otherwise = VertexPoint avgHoleEcps
where
newCoords = (oldCoords * m1) + (avgFacePoints * m2) + (avgMidEdges * m3)
oldCoords = vertexPoint vertex
avgFacePoints = average $ fmap (getFacePoint . (facePoints !)) faceIds
avgMidEdges = average $ fmap (getEdgeCenterPoint . (ecps !)) edgeIds
m1 = (n - 3) / n
m2 = 1 / n
m3 = 2 / n
n = genericLength faceIds
faceIds = vertexFaces vertex
edgeIds = vertexEdges vertex
avgHoleEcps = average . (oldCoords:) . fmap (getEdgeCenterPoint . (ecps !)) $ filter (iehs !) edgeIds
 
allVertexPoints :: Fractional a => Mesh a -> FacePointArray a -> EdgeCenterPointArray a -> IsEdgeHoleArray -> VertexPointArray a
allVertexPoints mesh fps ecps iehs = fmap (\v -> vertexPoint' v fps ecps iehs) (meshVertices mesh)
 
-- For each vertex in a face, generate a set of new faces from it with its vertex point,
-- neighbor edge points, and face point. The new faces will refer to vertices in the
-- combined vertex array.
newFaces :: Face -> FaceId -> Int -> Int -> [SimpleFace]
newFaces (Face vertexIds edgeIds) faceId epOffset vpOffset
= take (genericLength vertexIds)
$ zipWith3 newFace (cycle vertexIds) (cycle edgeIds) (drop 1 (cycle edgeIds))
where
f = VertexId . (+ epOffset) . getEdgeId
newFace vid epA epB = SimpleFace
[ VertexId . (+ vpOffset) $ getVertexId vid
, f epA
, VertexId $ getFaceId faceId
, f epB]
 
subdivide :: Fractional a => SimpleMesh a -> SimpleMesh a
subdivide simpleMesh
= SimpleMesh combinedVertices (listArray (FaceId 0, FaceId (genericLength faces - 1)) faces)
where
mesh = makeComplexMesh simpleMesh
fps = allFacePoints mesh
ecps = allEdgeCenterPoints mesh
eps = allEdgePoints mesh fps ecps
iehs = allIsEdgeHoles mesh
vps = allVertexPoints mesh fps ecps iehs
edgePointOffset = length fps
vertexPointOffset = edgePointOffset + length eps
combinedVertices
= listArray (VertexId 0, VertexId (vertexPointOffset + length vps - 1))
. fmap SimpleVertex
$ concat [ fmap getFacePoint $ elems fps
, fmap getEdgePoint $ elems eps
, fmap getVertexPoint $ elems vps]
faces
= concat $ zipWith (\face fid -> newFaces face fid edgePointOffset vertexPointOffset)
(elems $ meshFaces mesh) (fmap FaceId [0..])
 
-- Transform to a Mesh by filling in the missing references and generating edges.
-- Faces can be updated with their edges, but must be ordered.
-- Edge and face order does not matter for vertices.
-- TODO: Discard degenerate faces (ones with 0 to 2 vertices/edges),
-- or we could transform these into single edges or vertices.
makeComplexMesh :: forall a. SimpleMesh a -> Mesh a
makeComplexMesh (SimpleMesh sVertices sFaces) = Mesh vertices edges faces
where
makeEdgesFromFace :: SimpleFace -> FaceId -> [Edge]
makeEdgesFromFace (SimpleFace vertexIds) fid
= take (genericLength vertexIds)
$ zipWith (\a b -> Edge a b [fid]) verts (drop 1 verts)
where
verts = cycle vertexIds
 
edgeKey :: VertexId -> VertexId -> (VertexId, VertexId)
edgeKey a b = (min a b, max a b)
 
sFacesList :: [SimpleFace]
sFacesList = elems sFaces
 
fids :: [FaceId]
fids = fmap FaceId [0..]
 
eids :: [EdgeId]
eids = fmap EdgeId [0..]
 
faceEdges :: [[Edge]]
faceEdges = zipWith makeEdgesFromFace sFacesList fids
 
edgeMap :: Map.Map (VertexId, VertexId) Edge
edgeMap
= Map.fromListWith (\(Edge a b fidsA) (Edge _ _ fidsB) -> Edge a b (fidsA ++ fidsB))
. fmap (\edge@(Edge a b _) -> (edgeKey a b, edge))
$ concat faceEdges
 
edges :: EdgeArray
edges = listArray (EdgeId 0, EdgeId $ (Map.size edgeMap) - 1) $ Map.elems edgeMap
 
edgeIdMap :: Map.Map (VertexId, VertexId) EdgeId
edgeIdMap = Map.fromList $ zipWith (\(Edge a b _) eid -> ((edgeKey a b), eid)) (elems edges) eids
 
faceEdgeIds :: [[EdgeId]]
faceEdgeIds = fmap (mapMaybe (\(Edge a b _) -> Map.lookup (edgeKey a b) edgeIdMap)) faceEdges
 
faces :: FaceArray
faces
= listArray (FaceId 0, FaceId $ (length sFaces) - 1)
$ zipWith (\(SimpleFace verts) edgeIds -> Face verts edgeIds) sFacesList faceEdgeIds
 
vidsToFids :: Map.Map VertexId [FaceId]
vidsToFids
= Map.fromListWith (++)
. concat
$ zipWith (\(SimpleFace vertexIds) fid -> fmap (\vid -> (vid, [fid])) vertexIds) sFacesList fids
 
vidsToEids :: Map.Map VertexId [EdgeId]
vidsToEids
= Map.fromListWith (++)
. concat
$ zipWith (\(Edge a b _) eid -> [(a, [eid]), (b, [eid])]) (elems edges) eids
 
simpleToComplexVert :: SimpleVertex a -> VertexId -> Vertex a
simpleToComplexVert (SimpleVertex point) vid
= Vertex point
(Map.findWithDefault [] vid vidsToEids)
(Map.findWithDefault [] vid vidsToFids)
 
vertices :: VertexArray a
vertices
= listArray (bounds sVertices)
$ zipWith simpleToComplexVert (elems sVertices) (fmap VertexId [0..])
 
pShowSimpleMesh :: Show a => SimpleMesh a -> String
pShowSimpleMesh (SimpleMesh vertices faces)
= "Vertices:\n" ++ (arrShow vertices sVertexPoint)
++ "Faces:\n" ++ (arrShow faces (fmap getVertexId . sFaceVertices))
where
arrShow a f = concatMap ((++ "\n") . show . (\(i, e) -> (i, f e))) . zip [0 :: Int ..] $ elems a
 
-- Testing types.
data Point a = Point a a a deriving (Show)
 
instance Functor Point where
fmap f (Point x y z) = Point (f x) (f y) (f z)
 
zipPoint :: (a -> b -> c) -> Point a -> Point b -> Point c
zipPoint f (Point x y z) (Point x' y' z') = Point (f x x') (f y y') (f z z')
 
instance Num a => Num (Point a) where
(+) = zipPoint (+)
(-) = zipPoint (-)
(*) = zipPoint (*)
negate = fmap negate
abs = fmap abs
signum = fmap signum
fromInteger i = let i' = fromInteger i in Point i' i' i'
 
instance Fractional a => Fractional (Point a) where
recip = fmap recip
fromRational r = let r' = fromRational r in Point r' r' r'
 
testCube :: SimpleMesh (Point Double)
testCube = SimpleMesh vertices faces
where
vertices = listArray (VertexId 0, VertexId 7)
$ fmap SimpleVertex
[ Point (-1) (-1) (-1)
, Point (-1) (-1) 1
, Point (-1) 1 (-1)
, Point (-1) 1 1
, Point 1 (-1) (-1)
, Point 1 (-1) 1
, Point 1 1 (-1)
, Point 1 1 1]
faces = listArray (FaceId 0, FaceId 5)
$ fmap (SimpleFace . (fmap VertexId))
[ [0, 4, 5, 1]
, [4, 6, 7, 5]
, [6, 2, 3, 7]
, [2, 0, 1, 3]
, [1, 5, 7, 3]
, [0, 2, 6, 4]]
 
testCubeWithHole :: SimpleMesh (Point Double)
testCubeWithHole
= SimpleMesh (sMeshVertices testCube) (ixmap (FaceId 0, FaceId 4) id (sMeshFaces testCube))
 
testTriangle :: SimpleMesh (Point Double)
testTriangle = SimpleMesh vertices faces
where
vertices = listArray (VertexId 0, VertexId 2)
$ fmap SimpleVertex
[ Point 0 0 0
, Point 0 0 1
, Point 0 1 0]
faces = listArray (FaceId 0, FaceId 0)
$ fmap (SimpleFace . (fmap VertexId))
[ [0, 1, 2]]
 
main :: IO ()
main = putStr . pShowSimpleMesh $ subdivide testCube</lang>
{{out}}
<pre>Vertices:
(0,Point 0.0 (-1.0) 0.0)
(1,Point 1.0 0.0 0.0)
(2,Point 0.0 1.0 0.0)
(3,Point (-1.0) 0.0 0.0)
(4,Point 0.0 0.0 1.0)
(5,Point 0.0 0.0 (-1.0))
(6,Point (-0.75) (-0.75) 0.0)
(7,Point (-0.75) 0.0 (-0.75))
(8,Point 0.0 (-0.75) (-0.75))
(9,Point (-0.75) 0.0 0.75)
(10,Point 0.0 (-0.75) 0.75)
(11,Point (-0.75) 0.75 0.0)
(12,Point 0.0 0.75 (-0.75))
(13,Point 0.0 0.75 0.75)
(14,Point 0.75 (-0.75) 0.0)
(15,Point 0.75 0.0 (-0.75))
(16,Point 0.75 0.0 0.75)
(17,Point 0.75 0.75 0.0)
(18,Point (-0.5555555555555556) (-0.5555555555555556) (-0.5555555555555556))
(19,Point (-0.5555555555555556) (-0.5555555555555556) 0.5555555555555556)
(20,Point (-0.5555555555555556) 0.5555555555555556 (-0.5555555555555556))
(21,Point (-0.5555555555555556) 0.5555555555555556 0.5555555555555556)
(22,Point 0.5555555555555556 (-0.5555555555555556) (-0.5555555555555556))
(23,Point 0.5555555555555556 (-0.5555555555555556) 0.5555555555555556)
(24,Point 0.5555555555555556 0.5555555555555556 (-0.5555555555555556))
(25,Point 0.5555555555555556 0.5555555555555556 0.5555555555555556)
Faces:
(0,[18,8,0,14])
(1,[22,14,0,10])
(2,[23,10,0,6])
(3,[19,6,0,8])
(4,[22,15,1,17])
(5,[24,17,1,16])
(6,[25,16,1,14])
(7,[23,14,1,15])
(8,[24,12,2,11])
(9,[20,11,2,13])
(10,[21,13,2,17])
(11,[25,17,2,12])
(12,[20,7,3,6])
(13,[18,6,3,9])
(14,[19,9,3,11])
(15,[21,11,3,7])
(16,[19,10,4,16])
(17,[23,16,4,13])
(18,[25,13,4,9])
(19,[21,9,4,10])
(20,[18,7,5,12])
(21,[20,12,5,15])
(22,[24,15,5,8])
(23,[22,8,5,7])</pre>
 
=={{header|J}}==
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/136040"