Rosetta Code

Catmull–Clark subdivision surface: Difference between revisions

Page Discussion

Catmull–Clark subdivision surface (view source)

Revision as of 12:15, 12 November 2023

156,683 bytes added ,  6 months ago
m
→‎{{header|Wren}}: Minor tidy
(J: partial bugfix)
m (→‎{{header|Wren}}: Minor tidy)
(42 intermediate revisions by 18 users not shown)
Line 1:
{{task|3D}}[[Category:Graphics algorithms]]{{requires|Graphics}}
[[Category:Geometry]]
{{task|3D}}{{requires|Graphics}}
Implement the Catmull-Clark surface subdivision ([[wp:Catmull–Clark_subdivision_surface|description on Wikipedia]]), which is an algorithm that maps from a surface (described as a set of points and a set of polygons with vertices at those points) to another more refined surface. The resulting surface will always consist of a mesh of quadrilaterals.
 
Line 42 ⟶ 44:
 
For edges and vertices not next to a hole, the standard algorithm from above is used.
=={{header|C}}==
[[file:catmull-C.png|center]]
Only the subdivision part. The [[Catmull–Clark_subdivision_surface/C|full source]] is way too long to be shown here. Lots of macros, you'll have to see the full code to know what's what.
<syntaxhighlight lang="c">vertex face_point(face f)
{
int i;
vertex v;
 
if (!f->avg) {
f->avg = vertex_new();
foreach(i, v, f->v)
if (!i) f->avg->pos = v->pos;
else vadd(f->avg->pos, v->pos);
 
vdiv(f->avg->pos, len(f->v));
}
return f->avg;
}
 
#define hole_edge(e) (len(e->f)==1)
vertex edge_point(edge e)
{
int i;
face f;
 
if (!e->e_pt) {
e->e_pt = vertex_new();
e->avg = e->v[0]->pos;
vadd(e->avg, e->v[1]->pos);
e->e_pt->pos = e->avg;
 
if (!hole_edge(e)) {
foreach (i, f, e->f)
vadd(e->e_pt->pos, face_point(f)->pos);
vdiv(e->e_pt->pos, 4);
} else
vdiv(e->e_pt->pos, 2);
 
vdiv(e->avg, 2);
}
 
return e->e_pt;
}
 
#define hole_vertex(v) (len((v)->f) != len((v)->e))
vertex updated_point(vertex v)
{
int i, n = 0;
edge e;
face f;
coord_t sum = {0, 0, 0};
 
if (v->v_new) return v->v_new;
 
v->v_new = vertex_new();
if (hole_vertex(v)) {
v->v_new->pos = v->pos;
foreach(i, e, v->e) {
if (!hole_edge(e)) continue;
vadd(v->v_new->pos, edge_point(e)->pos);
n++;
}
vdiv(v->v_new->pos, n + 1);
} else {
n = len(v->f);
foreach(i, f, v->f)
vadd(sum, face_point(f)->pos);
foreach(i, e, v->e)
vmadd(sum, edge_point(e)->pos, 2, sum);
vdiv(sum, n);
vmadd(sum, v->pos, n - 3, sum);
vdiv(sum, n);
v->v_new->pos = sum;
}
 
return v->v_new;
}
 
model catmull(model m)
{
int i, j, a, b, c, d;
face f;
vertex v, x;
 
model nm = model_new();
foreach (i, f, m->f) {
foreach(j, v, f->v) {
_get_idx(a, updated_point(v));
_get_idx(b, edge_point(elem(f->e, (j + 1) % len(f->e))));
_get_idx(c, face_point(f));
_get_idx(d, edge_point(elem(f->e, j)));
model_add_face(nm, 4, a, b, c, d);
}
}
return nm;
}</syntaxhighlight>
=={{header|Go}}==
{{trans|Python}}
<br>
This just prints the new points and faces after 1 iteration of the Catmull-Clark algorithm, agreeing with what the Python results would have been had they been printed rather than plotted.
 
See the original version for full comments.
<syntaxhighlight lang="go">package main
 
import (
"fmt"
"sort"
)
 
type (
Point [3]float64
Face []int
 
Edge struct {
pn1 int // point number 1
pn2 int // point number 2
fn1 int // face number 1
fn2 int // face number 2
cp Point // center point
}
 
PointEx struct {
p Point
n int
}
)
 
func sumPoint(p1, p2 Point) Point {
sp := Point{}
for i := 0; i < 3; i++ {
sp[i] = p1[i] + p2[i]
}
return sp
}
 
func mulPoint(p Point, m float64) Point {
mp := Point{}
for i := 0; i < 3; i++ {
mp[i] = p[i] * m
}
return mp
}
 
func divPoint(p Point, d float64) Point {
return mulPoint(p, 1.0/d)
}
 
func centerPoint(p1, p2 Point) Point {
return divPoint(sumPoint(p1, p2), 2)
}
 
func getFacePoints(inputPoints []Point, inputFaces []Face) []Point {
facePoints := make([]Point, len(inputFaces))
for i, currFace := range inputFaces {
facePoint := Point{}
for _, cpi := range currFace {
currPoint := inputPoints[cpi]
facePoint = sumPoint(facePoint, currPoint)
}
facePoint = divPoint(facePoint, float64(len(currFace)))
facePoints[i] = facePoint
}
return facePoints
}
 
func getEdgesFaces(inputPoints []Point, inputFaces []Face) []Edge {
var edges [][3]int
for faceNum, face := range inputFaces {
numPoints := len(face)
for pointIndex := 0; pointIndex < numPoints; pointIndex++ {
pointNum1 := face[pointIndex]
var pointNum2 int
if pointIndex < numPoints-1 {
pointNum2 = face[pointIndex+1]
} else {
pointNum2 = face[0]
}
if pointNum1 > pointNum2 {
pointNum1, pointNum2 = pointNum2, pointNum1
}
edges = append(edges, [3]int{pointNum1, pointNum2, faceNum})
}
}
sort.Slice(edges, func(i, j int) bool {
if edges[i][0] == edges[j][0] {
if edges[i][1] == edges[j][1] {
return edges[i][2] < edges[j][2]
}
return edges[i][1] < edges[j][1]
}
return edges[i][0] < edges[j][0]
})
numEdges := len(edges)
eIndex := 0
var mergedEdges [][4]int
for eIndex < numEdges {
e1 := edges[eIndex]
if eIndex < numEdges-1 {
e2 := edges[eIndex+1]
if e1[0] == e2[0] && e1[1] == e2[1] {
mergedEdges = append(mergedEdges, [4]int{e1[0], e1[1], e1[2], e2[2]})
eIndex += 2
} else {
mergedEdges = append(mergedEdges, [4]int{e1[0], e1[1], e1[2], -1})
eIndex++
}
} else {
mergedEdges = append(mergedEdges, [4]int{e1[0], e1[1], e1[2], -1})
eIndex++
}
}
var edgesCenters []Edge
for _, me := range mergedEdges {
p1 := inputPoints[me[0]]
p2 := inputPoints[me[1]]
cp := centerPoint(p1, p2)
edgesCenters = append(edgesCenters, Edge{me[0], me[1], me[2], me[3], cp})
}
return edgesCenters
}
 
func getEdgePoints(inputPoints []Point, edgesFaces []Edge, facePoints []Point) []Point {
edgePoints := make([]Point, len(edgesFaces))
for i, edge := range edgesFaces {
cp := edge.cp
fp1 := facePoints[edge.fn1]
var fp2 Point
if edge.fn2 == -1 {
fp2 = fp1
} else {
fp2 = facePoints[edge.fn2]
}
cfp := centerPoint(fp1, fp2)
edgePoints[i] = centerPoint(cp, cfp)
}
return edgePoints
}
 
func getAvgFacePoints(inputPoints []Point, inputFaces []Face, facePoints []Point) []Point {
numPoints := len(inputPoints)
tempPoints := make([]PointEx, numPoints)
for faceNum := range inputFaces {
fp := facePoints[faceNum]
for _, pointNum := range inputFaces[faceNum] {
tp := tempPoints[pointNum].p
tempPoints[pointNum].p = sumPoint(tp, fp)
tempPoints[pointNum].n++
}
}
avgFacePoints := make([]Point, numPoints)
for i, tp := range tempPoints {
avgFacePoints[i] = divPoint(tp.p, float64(tp.n))
}
return avgFacePoints
}
 
func getAvgMidEdges(inputPoints []Point, edgesFaces []Edge) []Point {
numPoints := len(inputPoints)
tempPoints := make([]PointEx, numPoints)
for _, edge := range edgesFaces {
cp := edge.cp
for _, pointNum := range []int{edge.pn1, edge.pn2} {
tp := tempPoints[pointNum].p
tempPoints[pointNum].p = sumPoint(tp, cp)
tempPoints[pointNum].n++
}
}
avgMidEdges := make([]Point, len(tempPoints))
for i, tp := range tempPoints {
avgMidEdges[i] = divPoint(tp.p, float64(tp.n))
}
return avgMidEdges
}
 
func getPointsFaces(inputPoints []Point, inputFaces []Face) []int {
numPoints := len(inputPoints)
pointsFaces := make([]int, numPoints)
for faceNum := range inputFaces {
for _, pointNum := range inputFaces[faceNum] {
pointsFaces[pointNum]++
}
}
return pointsFaces
}
 
func getNewPoints(inputPoints []Point, pointsFaces []int, avgFacePoints, avgMidEdges []Point) []Point {
newPoints := make([]Point, len(inputPoints))
for pointNum := range inputPoints {
n := float64(pointsFaces[pointNum])
m1, m2, m3 := (n-3)/n, 1.0/n, 2.0/n
oldCoords := inputPoints[pointNum]
p1 := mulPoint(oldCoords, m1)
afp := avgFacePoints[pointNum]
p2 := mulPoint(afp, m2)
ame := avgMidEdges[pointNum]
p3 := mulPoint(ame, m3)
p4 := sumPoint(p1, p2)
newPoints[pointNum] = sumPoint(p4, p3)
}
return newPoints
}
 
func switchNums(pointNums [2]int) [2]int {
if pointNums[0] < pointNums[1] {
return pointNums
}
return [2]int{pointNums[1], pointNums[0]}
}
 
func cmcSubdiv(inputPoints []Point, inputFaces []Face) ([]Point, []Face) {
facePoints := getFacePoints(inputPoints, inputFaces)
edgesFaces := getEdgesFaces(inputPoints, inputFaces)
edgePoints := getEdgePoints(inputPoints, edgesFaces, facePoints)
avgFacePoints := getAvgFacePoints(inputPoints, inputFaces, facePoints)
avgMidEdges := getAvgMidEdges(inputPoints, edgesFaces)
pointsFaces := getPointsFaces(inputPoints, inputFaces)
newPoints := getNewPoints(inputPoints, pointsFaces, avgFacePoints, avgMidEdges)
var facePointNums []int
nextPointNum := len(newPoints)
for _, facePoint := range facePoints {
newPoints = append(newPoints, facePoint)
facePointNums = append(facePointNums, nextPointNum)
nextPointNum++
}
edgePointNums := make(map[[2]int]int)
for edgeNum := range edgesFaces {
pointNum1 := edgesFaces[edgeNum].pn1
pointNum2 := edgesFaces[edgeNum].pn2
edgePoint := edgePoints[edgeNum]
newPoints = append(newPoints, edgePoint)
edgePointNums[[2]int{pointNum1, pointNum2}] = nextPointNum
nextPointNum++
}
var newFaces []Face
for oldFaceNum, oldFace := range inputFaces {
if len(oldFace) == 4 {
a, b, c, d := oldFace[0], oldFace[1], oldFace[2], oldFace[3]
facePointAbcd := facePointNums[oldFaceNum]
edgePointAb := edgePointNums[switchNums([2]int{a, b})]
edgePointDa := edgePointNums[switchNums([2]int{d, a})]
edgePointBc := edgePointNums[switchNums([2]int{b, c})]
edgePointCd := edgePointNums[switchNums([2]int{c, d})]
newFaces = append(newFaces, Face{a, edgePointAb, facePointAbcd, edgePointDa})
newFaces = append(newFaces, Face{b, edgePointBc, facePointAbcd, edgePointAb})
newFaces = append(newFaces, Face{c, edgePointCd, facePointAbcd, edgePointBc})
newFaces = append(newFaces, Face{d, edgePointDa, facePointAbcd, edgePointCd})
}
}
return newPoints, newFaces
}
 
func main() {
inputPoints := []Point{
{-1.0, 1.0, 1.0},
{-1.0, -1.0, 1.0},
{1.0, -1.0, 1.0},
{1.0, 1.0, 1.0},
{1.0, -1.0, -1.0},
{1.0, 1.0, -1.0},
{-1.0, -1.0, -1.0},
{-1.0, 1.0, -1.0},
}
 
inputFaces := []Face{
{0, 1, 2, 3},
{3, 2, 4, 5},
{5, 4, 6, 7},
{7, 0, 3, 5},
{7, 6, 1, 0},
{6, 1, 2, 4},
}
 
outputPoints := make([]Point, len(inputPoints))
outputFaces := make([]Face, len(inputFaces))
copy(outputPoints, inputPoints)
copy(outputFaces, inputFaces)
iterations := 1
for i := 0; i < iterations; i++ {
outputPoints, outputFaces = cmcSubdiv(outputPoints, outputFaces)
}
for _, p := range outputPoints {
fmt.Printf("% .4f\n", p)
}
fmt.Println()
for _, f := range outputFaces {
fmt.Printf("%2d\n", f)
}
}</syntaxhighlight>
 
{{out}}
<pre>
[-0.5556 0.5556 0.5556]
[-0.5556 -0.5556 0.5556]
[ 0.5556 -0.5556 0.5556]
[ 0.5556 0.5556 0.5556]
[ 0.5556 -0.5556 -0.5556]
[ 0.5556 0.5556 -0.5556]
[-0.5556 -0.5556 -0.5556]
[-0.5556 0.5556 -0.5556]
[ 0.0000 0.0000 1.0000]
[ 1.0000 0.0000 0.0000]
[ 0.0000 0.0000 -1.0000]
[ 0.0000 1.0000 0.0000]
[-1.0000 0.0000 0.0000]
[ 0.0000 -1.0000 0.0000]
[-0.7500 0.0000 0.7500]
[ 0.0000 0.7500 0.7500]
[-0.7500 0.7500 0.0000]
[ 0.0000 -0.7500 0.7500]
[-0.7500 -0.7500 0.0000]
[ 0.7500 0.0000 0.7500]
[ 0.7500 -0.7500 0.0000]
[ 0.7500 0.7500 0.0000]
[ 0.7500 0.0000 -0.7500]
[ 0.0000 -0.7500 -0.7500]
[ 0.0000 0.7500 -0.7500]
[-0.7500 0.0000 -0.7500]
 
[ 0 14 8 15]
[ 1 17 8 14]
[ 2 19 8 17]
[ 3 15 8 19]
[ 3 19 9 21]
[ 2 20 9 19]
[ 4 22 9 20]
[ 5 21 9 22]
[ 5 22 10 24]
[ 4 23 10 22]
[ 6 25 10 23]
[ 7 24 10 25]
[ 7 16 11 24]
[ 0 15 11 16]
[ 3 21 11 15]
[ 5 24 11 21]
[ 7 25 12 16]
[ 6 18 12 25]
[ 1 14 12 18]
[ 0 16 12 14]
[ 6 18 13 23]
[ 1 17 13 18]
[ 2 20 13 17]
[ 4 23 13 20]
</pre>
=={{header|Haskell}}==
<syntaxhighlight 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</syntaxhighlight>
{{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}}==
 
<syntaxhighlight lang="j">avg=: +/ % #
FIXME: result points do not correspond to sample data on talk page (I have factors of 3/4 where sample has 2/3).
 
havePoints=: e."1/~ i.@#
 
<lang j>havePoints=: e."1/~ i.@#
catmullclark=:3 :0
'mesh points'=. y
face_point=:. avg"2 mesh{points
point_face=:. |: mesh havePoints points
avg_face_points=:. point_face avg@#"1 2 face_point
edges=:. ~.,/ meshEdges=:. mesh /:~@,"+1|."1 mesh
edge_face=:. *./"2 edges e."0 1/ mesh
edge_center=:. avg"2 edges{points
edge_point=:. (0.5*edge_center) + 0.25 * edge_face +/ .* face_point
point_edge=:. |: edges havePoints points
avg_mid_edges=:. point_edge avg@#"1 2 edge_center
n=:. +/"1 point_edge
'm3 m2 m1'=:. (2,1,:n-3)%"1 n
new_coords=:. (m1 * points) + (m2 * avg_face_points) + (m3 * avg_mid_edges)
pts=:. face_point,edge_point,new_coords
c0=:. (#edge_point)+ e0=:. #face_point
msh=:. (,c0+mesh),.(,e0+edgeedges i. meshEdges),.((#i.)~/$mesh),.,e0+_1|."1 edgeedges i. meshEdges
msh;pts
)</langsyntaxhighlight>
 
Example use:
 
<langsyntaxhighlight lang="j">NB.cube
points=: _1+2*#:i.8
mesh=: 1 A."1 I.(,1-|.)8&$@#&0 1">4 2 1
 
catmullclark mesh;points
Line 102 ⟶ 912:
│ │ 0.555556 0.555556 _0.555556│
│ │ 0.555556 0.555556 0.555556│
└──────────┴─────────────────────────────┘</langsyntaxhighlight>
=={{header|Julia}}==
<syntaxhighlight lang="julia">using Makie, Statistics
 
# Point3f0 is a 3-tuple of 32-bit floats for 3-dimensional space, and all Points are 3D.
=={{header|Mathematica}}==
Point = Point3f0
 
# a Face is defined by the points that are its vertices, in order.
{{improve|Mathematica|Reformat the code to be more readable; to make its structure more visible.}}
Face = Vector{Point}
 
# an Edge is a line segment where the points are sorted
struct Edge
p1::Point
p2::Point
Edge(a, b) = new(min(a, b), max(a, b))
end
 
edgemidpoint(edge) = (edge.p1 + edge.p2) / 2.0
facesforpoint(p, faces) = [f for f in faces if p in f]
facesforedge(e, faces) = [f for f in faces if (e.p1 in f) && (e.p2 in f)]
nexttohole(edge, faces) = length(facesforedge(edge, faces)) < 2
 
function newedgepoint(edge, faces)
f = facesforedge(edge, faces)
p1, p2, len = edge.p1, edge.p2, length(f)
if len == 2
return (p1 + p2 + mean(f[1]) + mean(f[2])) / 4.0
elseif len == 1
return (p1 + p2 + mean(f[1])) / 3.0
end
return (p1 + p2) / 2.0
end
 
function edgesforface(face)
ret, indices = Vector{Edge}(), collect(1:length(face))
for i in 1:length(face)-1
push!(ret, Edge(face[indices[1]], face[indices[2]]))
indices .= circshift(indices, 1)
end
ret
end
 
function edgesforpoint(p, faces)
f = filter(x -> p in x, faces)
return filter(e -> p == e.p1 || p == e.p2, mapreduce(edgesforface, vcat, f))
end
 
function adjacentpoints(point, face)
a = indexin([point], face)
if a[1] != nothing
adjacent = (a[1] == 1) ? [face[end], face[2]] :
a[1] == length(face) ? [face[end-1], face[1]] :
[face[a[1] - 1], face[a[1] + 1]]
return sort(adjacent)
else
throw("point $point not in face $face")
end
end
 
adjacentedges(point, face) = [Edge(point, x) for x in adjacentpoints(point, face)]
facewrapped(face) = begin f = deepcopy(face); push!(f, f[1]); f end
drawface(face, colr) = lines(facewrapped(face); color=colr)
drawface!(face, colr) = lines!(facewrapped(face); color=colr)
drawfaces!(faces, colr) = for f in faces drawface!(f, colr) end
const colors = [:red, :green, :blue, :gold]
 
function drawfaces(faces, colr)
scene = drawface(faces[1], colr)
if length(faces) > 1
for f in faces[2:end]
drawface!(f, colr)
end
end
scene
end
 
function catmullclarkstep(faces)
d, E, dprime = Set(reduce(vcat, faces)), Dict{Vector, Point}(), Dict{Point, Point}()
for face in faces, (i, p) in enumerate(face)
edge = (p == face[end]) ? Edge(p, face[1]) : Edge(p, face[i + 1])
E[[edge, face]] = newedgepoint(edge, faces)
end
for p in d
F = mean([mean(face) for face in facesforpoint(p, faces)])
pe = edgesforpoint(p, faces)
R = mean(map(edgemidpoint, pe))
n = length(pe)
dprime[p] = (F + 2 * R + p * (n - 3)) / n
end
newfaces = Vector{Face}()
for face in faces
v = mean(face)
for point in face
fp1, fp2 = map(x -> E[[x, face]], adjacentedges(point, face))
push!(newfaces, [fp1, dprime[point], fp2, v])
end
end
return newfaces
end
 
"""
catmullclark(faces, iters, scene)
 
Perform a multistep Catmull-Clark subdivision of a surface. See Wikipedia or page 53
of http://graphics.stanford.edu/courses/cs468-10-fall/LectureSlides/10_Subdivision.pdf
Plots the iterations, with colors for each iteration as set in the colors array.
Uses a Makie Scene of scene to plot the iters iterations.
"""
function catmullclark(faces, iters, scene)
nextfaces = deepcopy(faces)
for i in 1:iters
nextfaces = catmullclarkstep(nextfaces)
drawfaces!(nextfaces, colors[i])
display(scene)
sleep(1)
end
end
 
const inputpoints = [
[-1.0, -1.0, -1.0],
[-1.0, -1.0, 1.0],
[-1.0, 1.0, -1.0],
[-1.0, 1.0, 1.0],
[1.0, -1.0, -1.0],
[1.0, -1.0, 1.0],
[1.0, 1.0, -1.0],
[1.0, 1.0, 1.0]
]
 
const inputfaces = [
[0, 4, 5, 1],
[4, 6, 7, 5],
[6, 2, 3, 7],
[2, 0, 1, 3],
[1, 5, 7, 3],
[0, 2, 6, 4]
]
 
const faces = [map(x -> Point3f0(inputpoints[x]), p .+ 1) for p in inputfaces]
 
scene = drawfaces(faces, :black)
display(scene)
sleep(1)
 
catmullclark(faces, 4, scene)
 
println("Press Enter to continue", readline())
</syntaxhighlight>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
This implementation supports tris, quads, and higher polys, as well as surfaces with holes.
The function relies on three externally defined general functionality functions:
 
<syntaxhighlight lang="mathematica">subSetQ[large_,small_] := MemberQ[large,small]
subSetQ[large_,small_List] := And@@(MemberQ[large,#]&/@small)
 
containing[groupList_,item_]:= Flatten[Position[groupList,group_/;subSetQ[group,item]]]
 
ReplaceFace[face_]:=Transpose[Prepend[Transpose[{#[[1]],face,#[[2]]}&/@Transpose[Partition[face,2,1,1]//{#,RotateRight[#]}&]],face]]</syntaxhighlight>
subSetQ[small,large] is a boolean test for whether small is a subset of large. Note this is not a general purpose implimentation and only serves this purpose under the constrictions of the following program.
 
containing[{obj1,obj2,...},item] Will return a list of indices of the objects containing item, where objects are faces or edges and item is edges or vertexes.
faces containing a given vertex, faces containing a given edge, edges containing a given point.
It is used for each such purpose in the code called via infix notation, the specific usage is easily distinguised by variable names. For example faces~containing~edge would be a list of the indices for the faces containing the given edge.
 
ReplaceFace[face] replaces the face with a list of descriptions for the new faces. It will return a list containing mixed objects, vertexes, edges and faces where edges and faces referes to the new vertexes to be generated by the code. When the new vertexes have been appended to the updated old vertexes, these mixed objects will be recalcluated into correct indices into the new vertex list by the later defined function newIndex[].
 
 
<syntaxhighlight lang="mathematica">CatMullClark[{Points_,faces_}]:=Block[{avgFacePoints,avgEdgePoints,updatedPoints,newEdgePoints,newPoints,edges,newFaces,weights,pointUpdate,edgeUpdate,newIndex},
edges = DeleteDuplicates[Flatten[Partition[#,2,1,-1]&/@faces,1],Sort[#1]==Sort[#2]&];
avgFacePoints=Mean[Points[[#]]] &/@ faces;
avgEdgePoints=Mean[Points[[#]]] &/@ edges;
 
weights[vertex_]:= Count[faces,vertex,2]//{(#-3),1,2}/#&;
pointUpdate[vertex_]:=
If[Length[faces~containing~vertex]!=Length[edges~containing~vertex],
Mean[avgEdgePoints[[Select[edges~containing~vertex,holeQ[edges[[#]],faces]&]]]],
Total[weights[vertex]{ Points[[vertex]], Mean[avgFacePoints[[faces~containing~vertex]]], Mean[avgEdgePoints[[edges~containing~vertex]]]}]
];
 
edgeUpdate[edge_]:=
If[Length[faces~containing~edge]==1,
Mean[Points[[edge]]],
Mean[Points[[Flatten[{edge, faces[[faces~containing~edge]]}]]]]
];
 
updatedPoints = pointUpdate/@Range[1,Length[Points]];
newEdgePoints = edgeUpdate/@edges;
newPoints = Join[updatedPoints,avgFacePoints,newEdgePoints];
 
newIndex[edge_/;Length[edge]==2] := Length[Points]+Length[faces]+Position[Sort/@edges,Sort@edge][[1,1]]
newIndex[face_] := Length[Points]+Position[faces,face][[1,1]]
 
newFaces = Flatten[Map[newIndex[#,{Points,edges,faces}]&,ReplaceFace/@faces,{-2}],1];
{newPoints,newFaces}
]</syntaxhighlight>
 
The implimentation can be tested with polygons with and without holes by using the polydata
<syntaxhighlight lang="mathematica">{points,faces}=PolyhedronData["Cube",{"VertexCoordinates","FaceIndices"}];
 
Function[iteration,
Graphics3D[(Polygon[iteration[[1]][[#]]]&/@iteration[[2]])]
]/@NestList[CatMullClark,{points,faces},3]//GraphicsRow</syntaxhighlight>
[[File:CAM noholes 1.png]]
 
For a surface with holes the resulting iterative subdivision will be:
<syntaxhighlight lang="mathematica">faces = Delete[faces, 6];
Function[iteration, Graphics3D[
(Polygon[iteration[[1]][[#]]] & /@ iteration[[2]])
]] /@ NestList[CatMullClark, {points, faces}, 3] // GraphicsRow</syntaxhighlight>
[[File:CAM holes 1.png]]
 
This code was written in Mathematica 8.
=={{header|Nim}}==
{{trans|Go}}
{{trans|Python}}
<syntaxhighlight lang="nim">import algorithm
import tables
 
const None = -1 # Index number used to indicate no data.
 
type
 
Point = array[3, float]
Face = seq[int]
 
Edge = object
pn1: int # Point number 1.
pn2: int # Point number 2.
fn1: int # Face number 1.
fn2: int # Face number 2.
cp: Point # Center point.
 
PointEx = object
p: Point
n: int
 
PointNums = tuple[pn1, pn2: int]
 
 
####################################################################################################
# Point operations.
 
func `+`(p1, p2: Point): Point =
## Adds points p1 and p2.
for i in 0..Point.high:
result[i] = p1[i] + p2[i]
 
#---------------------------------------------------------------------------------------------------
 
func `*`(p: Point; m: float): Point =
## Multiply point p by m.
for i in 0..Point.high:
result[i] = p[i] * m
 
#---------------------------------------------------------------------------------------------------
 
func `/`(p: Point; d: float): Point =
## Divide point p by d.
for i in 0..Point.high:
result[i] = p[i] / d
 
#---------------------------------------------------------------------------------------------------
 
func centerPoint(p1, p2: Point): Point =
## Return a point in the center of the segment ended by points p1 and p2.
for i in 0..Point.high:
result[i] = (p1[i] + p2[i]) / 2
 
 
####################################################################################################
 
func getFacePoints(inputPoints: seq[Point]; inputFaces: seq[Face]): seq[Point] =
## For each face, a face point is created which is the average of all the points of the face.
 
result.setLen(inputFaces.len)
for i, currFace in inputFaces:
var facePoint: Point
for currPointIndex in currFace:
# Add currPoint to facePoint. Will divide later.
facePoint = facePoint + inputPoints[currPointIndex]
result[i] = facePoint / currFace.len.toFloat
 
#---------------------------------------------------------------------------------------------------
 
func getEdgesFaces(inputPoints: seq[Point]; inputFaces: seq[Face]): seq[Edge] =
## Get list of edges and the one or two adjacent faces in a list.
## Also get center point of edge.
 
# Get edges for each face.
var edges: seq[array[3, int]]
for faceNum, face in inputFaces:
# Loop over index into face.
for pointIndex in 0..face.high:
# If not last point then edge is current point and next point
# and for last point, edge is current point and first point.
var pointNum1 = face[pointIndex]
var pointNum2 = if pointIndex < face.high: face[pointIndex + 1] else: face[0]
# Order points in edge by lowest point number.
if pointNum1 > pointNum2:
swap pointNum1, pointNum2
edges &= [pointNum1, pointNum2, faceNum]
 
# Sort edges by pointNum1, pointNum2, faceNum.
edges.sort(proc(a1, a2: array[3, int]): int =
result = cmp(a1[0], a2[0])
if result == 0:
result = cmp(a1[1], a2[1])
if result == 0:
result = cmp(a1[2], a2[2]))
 
# Merge edges with 2 adjacent faces:
# [pointNum1, pointNum2, faceNum1, faceNum2] or
# [pointNum1, pointNum2, faceNum1, None]
var eIndex = 0
var mergedEdges: seq[array[4, int]]
while eIndex < edges.len:
let e1 = edges[eIndex]
# Check if not last edge.
if eIndex < edges.high:
let e2 = edges[eIndex + 1]
if e1[0] == e2[0] and e1[1] == e2[1]:
mergedEdges &= [e1[0], e1[1], e1[2], e2[2]]
inc eIndex, 2
else:
mergedEdges &= [e1[0], e1[1], e1[2], None]
inc eIndex
else:
mergedEdges &= [e1[0], e1[1], e1[2], None]
inc eIndex
 
# Add edge centers.
for me in mergedEdges:
let p1 = inputPoints[me[0]]
let p2 = inputPoints[me[1]]
let cp = centerPoint(p1, p2)
result.add(Edge(pn1: me[0], pn2: me[1], fn1: me[2], fn2: me[3], cp: cp))
 
#---------------------------------------------------------------------------------------------------
 
func getEdgePoints(inputPoints: seq[Point];
edgesFaces: seq[Edge];
facePoints: seq[Point]): seq[Point] =
## For each edge, an edge point is created which is the average between the center of the
## edge and the center of the segment made with the face points of the two adjacent faces.
 
result.setLen(edgesFaces.len)
for i, edge in edgesFaces:
# Get center of two facepoints.
let fp1 = facePoints[edge.fn1]
# If not two faces, just use one facepoint (should not happen for solid like a cube).
let fp2 = if edge.fn2 == None: fp1 else: facePoints[edge.fn2]
let cfp = centerPoint(fp1, fp2)
# Get average between center of edge and center of facePoints.
result[i] = centerPoint(edge.cp, cfp)
 
#---------------------------------------------------------------------------------------------------
 
func getAvgFacePoints(inputPoints: seq[Point];
inputFaces: seq[Face];
facePoints: seq[Point]): seq[Point] =
## For each point calculate the average of the face points of the faces the
## point belongs to (avgFacePoints), create a list of lists of two numbers
## [facePointSum, numPoints] by going through the points in all the faces then
## create the avgFacePoints list of point by dividing pointSum(x, y, z) by numPoints
 
var tempPoints = newSeq[PointEx](inputPoints.len)
 
# Loop through faces, updating tempPoints.
for faceNum, pointNums in inputFaces:
let fp = facePoints[faceNum]
for pointNum in pointNums:
let tp = tempPoints[pointNum].p
tempPoints[pointNum].p = tp + fp
inc tempPoints[pointNum].n
 
# Divide to build the result.
result.setLen(inputPoints.len)
for i, tp in tempPoints:
result[i] = tp.p / tp.n.toFloat
 
#---------------------------------------------------------------------------------------------------
 
func getAvgMidEdges(inputPoints: seq[Point]; edgesFaces: seq[Edge]): seq[Point] =
## Return the average of the centers of edges the point belongs to (avgMidEdges).
## Create list with entry for each point. Each entry has two elements. One is a point
## that is the sum of the centers of the edges and the other is the number of edges.
## After going through all edges divide by number of edges.
 
var tempPoints = newSeq[PointEx](inputPoints.len)
 
# Go through edgesFaces using center updating each point.
for edge in edgesFaces:
for pointNum in [edge.pn1, edge.pn2]:
let tp = tempPoints[pointNum].p
tempPoints[pointNum].p = tp + edge.cp
inc tempPoints[pointNum].n
 
# Divide out number of points to get average.
result.setLen(inputPoints.len)
for i, tp in tempPoints:
result[i] = tp.p / tp.n.toFloat
 
#---------------------------------------------------------------------------------------------------
 
func getPointsFaces(inputPoints: seq[Point]; inputFaces: seq[Face]): seq[int] =
# Return the number of faces for each point.
 
result.setLen(inputPoints.len)
 
# Loop through faces updating the result.
for pointNums in inputFaces:
for pointNum in pointNums:
inc result[pointNum]
 
#---------------------------------------------------------------------------------------------------
 
func getNewPoints(inputPoints: seq[Point]; pointsFaces: seq[int];
avgFacePoints, avgMidEdges: seq[Point]): seq[Point] =
## m1 = (n - 3.0) / n
## m2 = 1.0 / n
## m3 = 2.0 / n
## newCoords = (m1 * oldCoords) + (m2 * avgFacePoints) + (m3 * avgMidEdges)
 
result.setLen(inputPoints.len)
for pointNum, oldCoords in inputPoints:
let n = pointsFaces[pointNum].toFloat
let (m1, m2, m3) = ((n - 3) / n, 1 / n, 2 / n)
let p1 = oldCoords * m1
let afp = avgFacePoints[pointNum]
let p2 = afp * m2
let ame = avgMidEdges[pointNum]
let p3 = ame * m3
let p4 = p1 + p2
result[pointNum] = p4 + p3
 
#---------------------------------------------------------------------------------------------------
 
func switchNums(pointNums: PointNums): PointNums =
## Return tuple of point numbers sorted least to most.
 
if pointNums.pn1 < pointNums.pn2: pointNums
else: (pointNums.pn2, pointNums.pn1)
 
#---------------------------------------------------------------------------------------------------
 
func cmcSubdiv(inputPoints: seq[Point];
inputFaces: seq[Face]): tuple[p: seq[Point], f: seq[Face]] =
## For each face, a face point is created which is the average of all the points of the face.
## Each entry in the returned list is a point (x, y, z).
 
let facePoints = getFacePoints(inputPoints, inputFaces)
let edgesFaces = getEdgesFaces(inputPoints, inputFaces)
let edgePoints = getEdgePoints(inputPoints, edgesFaces, facePoints)
let avgFacePoints = getAvgFacePoints(inputPoints, inputFaces, facePoints)
let avgMidEdges = getAvgMidEdges(inputPoints, edgesFaces)
let pointsFaces = getPointsFaces(inputPoints, inputFaces)
var newPoints = getNewPoints(inputPoints, pointsFaces, avgFacePoints, avgMidEdges)
 
#[Then each face is replaced by new faces made with the new points,
 
for a triangle face (a,b,c):
(a, edgePoint ab, facePoint abc, edgePoint ca)
(b, edgePoint bc, facePoint abc, edgePoint ab)
(c, edgePoint ca, facePoint abc, edgePoint bc)
 
for a quad face (a,b,c,d):
(a, edgePoint ab, facePoint abcd, edgePoint da)
(b, edgePoint bc, facePoint abcd, edgePoint ab)
(c, edgePoint cd, facePoint abcd, edgePoint bc)
(d, edgePoint da, facePoint abcd, edgePoint cd)
 
facePoints is a list indexed by face number so that is easy to get.
 
edgePoints is a list indexed by the edge number which is an index into edgesFaces.
 
Need to add facePoints and edgePoints to newPoints and get index into each.
 
Then create two new structures:
 
facePointNums: list indexes by faceNum whose value is the index into newPoints.
 
edgePointNums: dictionary with key (pointNum1, pointNum2) and value is index into newPoints.
]#
 
# Add face points to newPoints.
var facePointNums: seq[int]
var nextPointNum = newPoints.len # PointNum after next append to newPoints.
for facePoint in facePoints:
newPoints.add(facePoint)
facePointNums.add(nextPointNum)
inc nextPointNum
 
# Add edge points to newPoints.
var edgePointNums: Table[tuple[pn1, pn2: int], int]
for edgeNum, edgesFace in edgesFaces:
let pointNum1 = edgesFace.pn1
let pointNum2 = edgesFace.pn2
newPoints.add(edgePoints[edgeNum])
edgePointNums[(pointNum1, pointNum2)] = nextPointNum
inc nextPointNum
 
# newPoints now has the points to output. Need new faces.
 
#[Just doing this case for now:
 
for a quad face (a,b,c,d):
(a, edgePoint ab, facePoint abcd, edgePoint da)
(b, edgePoint bc, facePoint abcd, edgePoint ab)
(c, edgePoint cd, facePoint abcd, edgePoint bc)
(d, edgePoint da, facePoint abcd, edgePoint cd)
 
newFaces will be a list of lists where the elements are like this:
[pointNum1, pointNum2, pointNum3, pointNum4]
]#
 
var newFaces: seq[Face]
for oldFaceNum, oldFace in inputFaces:
# 4 points face.
if oldFace.len == 4:
let (a, b, c, d) = (oldFace[0], oldface[1], oldface[2], oldface[3])
let facePointAbcd = facePointNums[oldFaceNum]
let edgePointAb = edgePointNums[switchNums((a, b))]
let edgePointDa = edgePointNums[switchNums((d, a))]
let edgePointBc = edgePointNums[switchNums((b, c))]
let edgePointCd = edgePointNums[switchNums((c, d))]
newFaces &= @[a, edgePointAb, facePointAbcd, edgePointDa]
newFaces &= @[b, edgePointBc, facePointAbcd, edgePointAb]
newFaces &= @[c, edgePointCd, facePointAbcd, edgePointBc]
newFaces &= @[d, edgePointDa, facePointAbcd, edgePointCd]
 
result = (newPoints, newFaces)
 
#———————————————————————————————————————————————————————————————————————————————————————————————————
 
when isMainModule:
 
import strformat, strutils
 
let inputPoints = @[[-1.0, 1.0, 1.0],
[-1.0, -1.0, 1.0],
[1.0, -1.0, 1.0],
[1.0, 1.0, 1.0],
[1.0, -1.0, -1.0],
[1.0, 1.0, -1.0],
[-1.0, -1.0, -1.0],
[-1.0, 1.0, -1.0]]
 
let inputFaces = @[@[0, 1, 2, 3],
@[3, 2, 4, 5],
@[5, 4, 6, 7],
@[7, 0, 3, 5],
@[7, 6, 1, 0],
@[6, 1, 2, 4]]
 
var outputPoints = inputPoints
var outputFaces = inputFaces
const Iterations = 1
 
for i in 1..Iterations:
(outputPoints, outputFaces) = cmcSubdiv(outputPoints, outputFaces)
 
for p in outputPoints:
<lang Mathematica>CatmullClark[{v_, i_}] := Block[{e, vc, fp, ep, vp},
echo fmt"[{p[0]: .4f}, {p[1]: .4f}, {p[2]: .4f}]"
e = Function[a, {a, Select[Transpose[{i, Range@Length@i}],
echo ""
Length@Intersection[#[[1]], a] == 2 &][[All, 2]]}] /@
for nums in outputFaces:
Union[Sort /@ Flatten[Partition[#, 2, 1, 1] & /@ i, 1]];
var s = "["
vc = Table[{n, Select[Transpose[{i, Range@Length@i}], MemberQ[#[[1]], n] &][[All, 2]],
for n in nums:
Select[Transpose[{e[[All, 1]], Range@Length@e[[All, 1]]}], MemberQ[#[[1]], n] &][[All, 2]]}, {n, Length@v}];
s.addSep(", ", 1)
fp = Mean[v[[#]]] & /@ i;
s.add(fmt"{n:2d}")
ep = If[Length[#[[2]]] == 1, Mean[v[[#[[1]]]]], Mean@Join[v[[#[[1]]]], fp[[#[[2]]]]]] & /@ e;
s.add(']')
vp = If[Length[#[[2]]] != Length[#[[3]]], Mean@Join[{v[[#[[1]]]]}, ep[[Select[#[[3]],
echo s</syntaxhighlight>
Length[e[[#, 2]]] != 2 &]]]], ((Length@#[[2]] - 3) v[[#[[1]]]] + Mean@fp[[#[[2]]]] +
2 Mean@ep[[#[[3]]]])/Length@#[[2]]] & /@ vc;
{Join[vp, ep, fp], Flatten[Function[a, Function[
b, {a[[1]], #[[1]] + Length[vc], b + Length[vc] + Length[e], #[[2]] + Length[vc]} &@
Sort[Select[Transpose[{e, Range@Length@e}], MemberQ[#[[1, 1]], a[[1]]] && MemberQ[#[[1, 2]], b] &],
With[{f = i[[Intersection[#[[1, 2]], #2[[1, 2]]][[1]]]],
n = Intersection[#[[1, 1]], #2[[1, 1]]][[1]]},
Xor[Abs[#] == 1, # < 0] &@(Position[f, Complement[#[[1, 1]], {n}][[1]]] -
Position[f, n])[[1, 1]]] &][[All, 2]]] /@ a[[2]]] /@ vc, 1]}]
v = PolyhedronData["Cube", "VertexCoordinates"] // N
i = PolyhedronData["Cube", "FaceIndices"]
NestList[CatmullClark, {v, i}, 4];
Graphics3D[{FaceForm[{Opacity[0.3]}, {Opacity[0.1]}], GraphicsComplex[#[[1]], Polygon[#[[2]]]]}] & /@ %
Graphics3D[{EdgeForm[], FaceForm[White, Black],
GraphicsComplex[#[[1]], Polygon[#[[2]]], VertexNormals -> #[[1]]]}, Boxed -> False] & /@ %%
</lang>
 
{{out}}
The last few lines, after the function definition, do a test by using the built-in polyhedron data to generate the vertices and face indices. Then it repeatedly applies the method and graphs the results. Note that this was written in Mathematica 7, although it should be easy enough to port to maybe v5.2.
The output is the same as that of Go.
<pre>[-0.5556, 0.5556, 0.5556]
[-0.5556, -0.5556, 0.5556]
[ 0.5556, -0.5556, 0.5556]
[ 0.5556, 0.5556, 0.5556]
[ 0.5556, -0.5556, -0.5556]
[ 0.5556, 0.5556, -0.5556]
[-0.5556, -0.5556, -0.5556]
[-0.5556, 0.5556, -0.5556]
[ 0.0000, 0.0000, 1.0000]
[ 1.0000, 0.0000, 0.0000]
[ 0.0000, 0.0000, -1.0000]
[ 0.0000, 1.0000, 0.0000]
[-1.0000, 0.0000, 0.0000]
[ 0.0000, -1.0000, 0.0000]
[-0.7500, 0.0000, 0.7500]
[ 0.0000, 0.7500, 0.7500]
[-0.7500, 0.7500, 0.0000]
[ 0.0000, -0.7500, 0.7500]
[-0.7500, -0.7500, 0.0000]
[ 0.7500, 0.0000, 0.7500]
[ 0.7500, -0.7500, 0.0000]
[ 0.7500, 0.7500, 0.0000]
[ 0.7500, 0.0000, -0.7500]
[ 0.0000, -0.7500, -0.7500]
[ 0.0000, 0.7500, -0.7500]
[-0.7500, 0.0000, -0.7500]
 
[ 0, 14, 8, 15]
[ 1, 17, 8, 14]
[ 2, 19, 8, 17]
[ 3, 15, 8, 19]
[ 3, 19, 9, 21]
[ 2, 20, 9, 19]
[ 4, 22, 9, 20]
[ 5, 21, 9, 22]
[ 5, 22, 10, 24]
[ 4, 23, 10, 22]
[ 6, 25, 10, 23]
[ 7, 24, 10, 25]
[ 7, 16, 11, 24]
[ 0, 15, 11, 16]
[ 3, 21, 11, 15]
[ 5, 24, 11, 21]
[ 7, 25, 12, 16]
[ 6, 18, 12, 25]
[ 1, 14, 12, 18]
[ 0, 16, 12, 14]
[ 6, 18, 13, 23]
[ 1, 17, 13, 18]
[ 2, 20, 13, 17]
[ 4, 23, 13, 20]</pre>
=={{header|OCaml}}==
{{incorrect|OCaml|wrong output data}}
 
The implementation below only supports quad faces, but it does handle surfaces with holes.
Line 146 ⟶ 1,546:
In the [[Catmull–Clark subdivision surface/OCaml|sub-page]] there is also a program in OCaml+OpenGL which displays a cube subdivided 2 times with this algorithm.
 
<langsyntaxhighlight lang="ocaml">open Dynar
 
let add3 (x1, y1, z1) (x2, y2, z2) (x3, y3, z3) =
Line 281 ⟶ 1,681:
(Dynar.to_array da_points,
Dynar.to_array new_faces)
;;</langsyntaxhighlight>
=== Another implementation ===
Another implementation which should work with holes, but has only been tested on a cube
{{works with|OCaml|4.02+}}
<syntaxhighlight lang="ocaml">type point = { x: float; y : float; z : float }
let zero = { x = 0.0; y = 0.0; z = 0.0 }
let add a b = { x = a.x+.b.x; y = a.y+.b.y; z = a.z+.b.z }
let mul a k = { x = a.x*.k; y = a.y*.k; z= a.z*.k }
let div p k = mul p (1.0/.k)
 
type face = Face of point list
type edge = Edge of point*point
 
let make_edge a b = Edge (min a b, max a b)
let make_face a b c d = Face [a;b;c;d]
 
let centroid plist = div (List.fold_left add zero plist) (float (List.length plist))
let mid_edge (Edge (p1,p2)) = div (add p1 p2) 2.0
let face_point (Face pl) = centroid pl
let point_in_face p (Face pl) = List.mem p pl
let point_in_edge p (Edge (p1,p2)) = p = p1 || p = p2
let edge_in_face (Edge (p1,p2)) f = point_in_face p1 f && point_in_face p2 f
 
let border_edge faces e =
List.length (List.filter (edge_in_face e) faces) < 2
 
let edge_point faces e =
if border_edge faces e then mid_edge e else
let adjacent = List.filter (edge_in_face e) faces in
let fps = List.map face_point adjacent in
centroid [mid_edge e; centroid fps]
 
let mod_vertex faces edges p =
let v_edges = List.filter (point_in_edge p) edges in
let v_faces = List.filter (point_in_face p) faces in
let n = List.length v_faces in
let is_border = n <> (List.length v_edges) in
if is_border then
let border_mids = List.map mid_edge (List.filter (border_edge faces) v_edges) in
(* description ambiguity: average (border+p) or average(average(border),p) ?? *)
centroid (p :: border_mids)
else
let avg_face = centroid (List.map face_point v_faces) in
let avg_mid = centroid (List.map mid_edge v_edges) in
div (add (add (mul p (float(n-3))) avg_face) (mul avg_mid 2.0)) (float n)
 
let edges_of_face (Face pl) =
let rec next acc = function
| [] -> invalid_arg "empty face"
| a :: [] -> List.rev (make_edge a (List.hd pl) :: acc)
| a :: (b :: _ as xs) -> next (make_edge a b :: acc) xs in
next [] pl
 
let catmull_clark faces =
let module EdgeSet = Set.Make(struct type t = edge let compare = compare end) in
let edges = EdgeSet.elements (EdgeSet.of_list (List.concat (List.map edges_of_face faces))) in
let mod_face ((Face pl) as face) =
let fp = face_point face in
let ep = List.map (edge_point faces) (edges_of_face face) in
let e_tl = List.hd (List.rev ep) in
let vl = List.map (mod_vertex faces edges) pl in
let add_facet (e', acc) v e = e, (make_face e' v e fp :: acc) in
let _, new_faces = List.fold_left2 add_facet (e_tl, []) vl ep in
List.rev new_faces in
List.concat (List.map mod_face faces)
 
let show_faces fl =
let pr_point p = Printf.printf " (%.4f, %.4f, %.4f)" p.x p.y p.z in
let pr_face (Face pl) = print_string "Face:"; List.iter pr_point pl; print_string "\n" in
(print_string "surface {\n"; List.iter pr_face fl; print_string "}\n")
 
let c p q r = let s i = if i = 0 then -1.0 else 1.0 in { x = s p; y = s q; z = s r } ;;
let cube = [
Face [c 0 0 0; c 0 0 1; c 0 1 1; c 0 1 0]; Face [c 1 0 0; c 1 0 1; c 1 1 1; c 1 1 0];
Face [c 0 0 0; c 1 0 0; c 1 0 1; c 0 0 1]; Face [c 0 1 0; c 1 1 0; c 1 1 1; c 0 1 1];
Face [c 0 0 0; c 0 1 0; c 1 1 0; c 1 0 0]; Face [c 0 0 1; c 0 1 1; c 1 1 1; c 1 0 1] ] in
show_faces cube;
show_faces (catmull_clark cube)</syntaxhighlight>
with output:
<pre>surface {
Face: (-1.0000, -1.0000, -1.0000) (-1.0000, -1.0000, 1.0000) (-1.0000, 1.0000, 1.0000) (-1.0000, 1.0000, -1.0000)
Face: (1.0000, -1.0000, -1.0000) (1.0000, -1.0000, 1.0000) (1.0000, 1.0000, 1.0000) (1.0000, 1.0000, -1.0000)
Face: (-1.0000, -1.0000, -1.0000) (1.0000, -1.0000, -1.0000) (1.0000, -1.0000, 1.0000) (-1.0000, -1.0000, 1.0000)
Face: (-1.0000, 1.0000, -1.0000) (1.0000, 1.0000, -1.0000) (1.0000, 1.0000, 1.0000) (-1.0000, 1.0000, 1.0000)
Face: (-1.0000, -1.0000, -1.0000) (-1.0000, 1.0000, -1.0000) (1.0000, 1.0000, -1.0000) (1.0000, -1.0000, -1.0000)
Face: (-1.0000, -1.0000, 1.0000) (-1.0000, 1.0000, 1.0000) (1.0000, 1.0000, 1.0000) (1.0000, -1.0000, 1.0000)
}
surface {
Face: (-0.7500, 0.0000, -0.7500) (-0.5556, -0.5556, -0.5556) (-0.7500, -0.7500, 0.0000) (-1.0000, 0.0000, 0.0000)
Face: (-0.7500, -0.7500, 0.0000) (-0.5556, -0.5556, 0.5556) (-0.7500, 0.0000, 0.7500) (-1.0000, 0.0000, 0.0000)
Face: (-0.7500, 0.0000, 0.7500) (-0.5556, 0.5556, 0.5556) (-0.7500, 0.7500, 0.0000) (-1.0000, 0.0000, 0.0000)
Face: (-0.7500, 0.7500, 0.0000) (-0.5556, 0.5556, -0.5556) (-0.7500, 0.0000, -0.7500) (-1.0000, 0.0000, 0.0000)
Face: (0.7500, 0.0000, -0.7500) (0.5556, -0.5556, -0.5556) (0.7500, -0.7500, 0.0000) (1.0000, 0.0000, 0.0000)
Face: (0.7500, -0.7500, 0.0000) (0.5556, -0.5556, 0.5556) (0.7500, 0.0000, 0.7500) (1.0000, 0.0000, 0.0000)
Face: (0.7500, 0.0000, 0.7500) (0.5556, 0.5556, 0.5556) (0.7500, 0.7500, 0.0000) (1.0000, 0.0000, 0.0000)
Face: (0.7500, 0.7500, 0.0000) (0.5556, 0.5556, -0.5556) (0.7500, 0.0000, -0.7500) (1.0000, 0.0000, 0.0000)
Face: (-0.7500, -0.7500, 0.0000) (-0.5556, -0.5556, -0.5556) (0.0000, -0.7500, -0.7500) (0.0000, -1.0000, 0.0000)
Face: (0.0000, -0.7500, -0.7500) (0.5556, -0.5556, -0.5556) (0.7500, -0.7500, 0.0000) (0.0000, -1.0000, 0.0000)
Face: (0.7500, -0.7500, 0.0000) (0.5556, -0.5556, 0.5556) (0.0000, -0.7500, 0.7500) (0.0000, -1.0000, 0.0000)
Face: (0.0000, -0.7500, 0.7500) (-0.5556, -0.5556, 0.5556) (-0.7500, -0.7500, 0.0000) (0.0000, -1.0000, 0.0000)
Face: (-0.7500, 0.7500, 0.0000) (-0.5556, 0.5556, -0.5556) (0.0000, 0.7500, -0.7500) (0.0000, 1.0000, 0.0000)
Face: (0.0000, 0.7500, -0.7500) (0.5556, 0.5556, -0.5556) (0.7500, 0.7500, 0.0000) (0.0000, 1.0000, 0.0000)
Face: (0.7500, 0.7500, 0.0000) (0.5556, 0.5556, 0.5556) (0.0000, 0.7500, 0.7500) (0.0000, 1.0000, 0.0000)
Face: (0.0000, 0.7500, 0.7500) (-0.5556, 0.5556, 0.5556) (-0.7500, 0.7500, 0.0000) (0.0000, 1.0000, 0.0000)
Face: (0.0000, -0.7500, -0.7500) (-0.5556, -0.5556, -0.5556) (-0.7500, 0.0000, -0.7500) (0.0000, 0.0000, -1.0000)
Face: (-0.7500, 0.0000, -0.7500) (-0.5556, 0.5556, -0.5556) (0.0000, 0.7500, -0.7500) (0.0000, 0.0000, -1.0000)
Face: (0.0000, 0.7500, -0.7500) (0.5556, 0.5556, -0.5556) (0.7500, 0.0000, -0.7500) (0.0000, 0.0000, -1.0000)
Face: (0.7500, 0.0000, -0.7500) (0.5556, -0.5556, -0.5556) (0.0000, -0.7500, -0.7500) (0.0000, 0.0000, -1.0000)
Face: (0.0000, -0.7500, 0.7500) (-0.5556, -0.5556, 0.5556) (-0.7500, 0.0000, 0.7500) (0.0000, 0.0000, 1.0000)
Face: (-0.7500, 0.0000, 0.7500) (-0.5556, 0.5556, 0.5556) (0.0000, 0.7500, 0.7500) (0.0000, 0.0000, 1.0000)
Face: (0.0000, 0.7500, 0.7500) (0.5556, 0.5556, 0.5556) (0.7500, 0.0000, 0.7500) (0.0000, 0.0000, 1.0000)
Face: (0.7500, 0.0000, 0.7500) (0.5556, -0.5556, 0.5556) (0.0000, -0.7500, 0.7500) (0.0000, 0.0000, 1.0000)
}</pre>
=={{header|Phix}}==
{{trans|Go}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #000080;font-style:italic;">-- demo\rosetta\Catmull_Clark_subdivision_surface.exw</span>
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">newPoint</span><span style="color: #0000FF;">()</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">newPointEx</span><span style="color: #0000FF;">()</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">newPoint</span><span style="color: #0000FF;">(),</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">centerPoint</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p2</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p2</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">getFacePoints</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">facePoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">currFace</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">f</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">facePoint</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">currFace</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]]</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">currFace</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">facePoint</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">facePoint</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">currFace</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">facePoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">f</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">facePoint</span><span style="color: #0000FF;">,</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">currFace</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">facePoints</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">getEdgesFaces</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">edges</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">faceNum</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">face</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">faceNum</span><span style="color: #0000FF;">]</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">numPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">face</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">pointIndex</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">numPoints</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">pointNum1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">face</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointIndex</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">pointNum2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">face</span><span style="color: #0000FF;">[</span><span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pointIndex</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">numPoints</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">?</span><span style="color: #000000;">pointIndex</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">:</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)]</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">pointNum1</span> <span style="color: #0000FF;">></span> <span style="color: #000000;">pointNum2</span> <span style="color: #008080;">then</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">pointNum1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pointNum2</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">pointNum2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pointNum1</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">edges</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">edges</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">pointNum1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pointNum2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">faceNum</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">edges</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sort</span><span style="color: #0000FF;">(</span><span style="color: #000000;">edges</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">numEdges</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">edges</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">eIndex</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">mergedEdges</span> <span style="color: #0000FF;">={}</span>
<span style="color: #008080;">while</span> <span style="color: #000000;">eIndex</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">numEdges</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">e1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">edges</span><span style="color: #0000FF;">[</span><span style="color: #000000;">eIndex</span><span style="color: #0000FF;">])</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">e4</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">eIndex</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">numEdges</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span>
<span style="color: #008080;">and</span> <span style="color: #000000;">e1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]==</span><span style="color: #000000;">edges</span><span style="color: #0000FF;">[</span><span style="color: #000000;">eIndex</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">][</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">e4</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">edges</span><span style="color: #0000FF;">[</span><span style="color: #000000;">eIndex</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">][</span><span style="color: #000000;">3</span><span style="color: #0000FF;">])</span>
<span style="color: #000000;">eIndex</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">e1</span> <span style="color: #0000FF;">&=</span> <span style="color: #000000;">e4</span>
<span style="color: #000000;">mergedEdges</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mergedEdges</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">e1</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">eIndex</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">edgesCenters</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mergedEdges</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">me</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mergedEdges</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]),</span>
<span style="color: #000000;">cp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">centerPoint</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">me</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]],</span>
<span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">me</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]])</span>
<span style="color: #000000;">me</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">me</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cp</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">edgesCenters</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">edgesCenters</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">me</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">edgesCenters</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">getEdgePoints</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgesFaces</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">facePoints</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">edgePoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">edgesFaces</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">edgesFaces</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">edge</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">edgesFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">cp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">edge</span><span style="color: #0000FF;">[</span><span style="color: #000000;">5</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">fp1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">facePoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">edge</span><span style="color: #0000FF;">[</span><span style="color: #000000;">3</span><span style="color: #0000FF;">]],</span>
<span style="color: #000000;">fp2</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">edge</span><span style="color: #0000FF;">[</span><span style="color: #000000;">4</span><span style="color: #0000FF;">]=-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">?</span><span style="color: #000000;">fp1</span><span style="color: #0000FF;">:</span><span style="color: #000000;">facePoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">edge</span><span style="color: #0000FF;">[</span><span style="color: #000000;">4</span><span style="color: #0000FF;">]]),</span>
<span style="color: #000000;">cfp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">centerPoint</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fp1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fp2</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">edgePoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">centerPoint</span><span style="color: #0000FF;">(</span><span style="color: #000000;">cp</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">cfp</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">edgePoints</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">getAvgFacePoints</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">facePoints</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">numPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">tempPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">newPointEx</span><span style="color: #0000FF;">(),</span><span style="color: #000000;">numPoints</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">faceNum</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">fp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">facePoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">faceNum</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">faceNum</span><span style="color: #0000FF;">])</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">pointNum</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">faceNum</span><span style="color: #0000FF;">][</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">tp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tempPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">][</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">tempPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">][</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tp</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fp</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">tempPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">][</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">avgFacePoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">numPoints</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tempPoints</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">tp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tempPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">avgFacePoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tp</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">tp</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">avgFacePoints</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">getAvgMidEdges</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgesFaces</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">numPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">tempPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">newPointEx</span><span style="color: #0000FF;">(),</span> <span style="color: #000000;">numPoints</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">edgesFaces</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">edge</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">edgesFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">cp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">edge</span><span style="color: #0000FF;">[</span><span style="color: #000000;">5</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">edx</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">2</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">pointNum</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">edge</span><span style="color: #0000FF;">[</span><span style="color: #000000;">edx</span><span style="color: #0000FF;">]</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">tp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tempPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">][</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">tempPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">][</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tp</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">cp</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">tempPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">][</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">avgMidEdges</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tempPoints</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tempPoints</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">tp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tempPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">avgMidEdges</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tp</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">tp</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">avgMidEdges</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">getPointsFaces</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">numPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pointsFaces</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">numPoints</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">faceNum</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">faceNum</span><span style="color: #0000FF;">])</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">pointNum</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">faceNum</span><span style="color: #0000FF;">][</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">pointsFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">pointsFaces</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">getNewPoints</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pointsFaces</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">avgFacePoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">avgMidEdges</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">newPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">pointNum</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pointsFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">]</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">p1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">],</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">p2</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">avgFacePoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">p3</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">avgMidEdges</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">/</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">newPoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">pointNum</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p2</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">p3</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">newPoints</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">switchNums</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">pointNums</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">pointNums</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;"><</span> <span style="color: #000000;">pointNums</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">then</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">pointNums</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">pointNums</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">pointNums</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">cmcSubdiv</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">facePoints</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">getFacePoints</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">edgesFaces</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">getEdgesFaces</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">edgePoints</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">getEdgePoints</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgesFaces</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">facePoints</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">avgFacePoints</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">getAvgFacePoints</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">facePoints</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">avgMidEdges</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">getAvgMidEdges</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgesFaces</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">pointsFaces</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">getPointsFaces</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">newPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">getNewPoints</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pointsFaces</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">avgFacePoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">avgMidEdges</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">facePointNums</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">nextPointNum</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">newPoints</span><span style="color: #0000FF;">)+</span><span style="color: #000000;">1</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">facePoints</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">facePoint</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">facePoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">newPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">newPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">facePoint</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">facePointNums</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">facePointNums</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nextPointNum</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">nextPointNum</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">edgePointNums</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">new_dict</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">edgeNum</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">edgesFaces</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">pointNum1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pointNum2</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">edgesFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">edgeNum</span><span style="color: #0000FF;">]</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">edgePoint</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">edgePoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">edgeNum</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">newPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">newPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgePoint</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">setd</span><span style="color: #0000FF;">({</span><span style="color: #000000;">pointNum1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pointNum2</span><span style="color: #0000FF;">},</span><span style="color: #000000;">nextPointNum</span><span style="color: #0000FF;">,</span><span style="color: #000000;">edgePointNums</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">nextPointNum</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">newFaces</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">oldFaceNum</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">oldFace</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">[</span><span style="color: #000000;">oldFaceNum</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">oldFace</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">==</span> <span style="color: #000000;">4</span> <span style="color: #008080;">then</span>
<span style="color: #004080;">integer</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">oldFace</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">facePointAbcd</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">facePointNums</span><span style="color: #0000FF;">[</span><span style="color: #000000;">oldFaceNum</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">edgePointAb</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">getd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">switchNums</span><span style="color: #0000FF;">({</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">}),</span><span style="color: #000000;">edgePointNums</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">edgePointDa</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">getd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">switchNums</span><span style="color: #0000FF;">({</span><span style="color: #000000;">d</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">}),</span><span style="color: #000000;">edgePointNums</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">edgePointBc</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">getd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">switchNums</span><span style="color: #0000FF;">({</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">}),</span><span style="color: #000000;">edgePointNums</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">edgePointCd</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">getd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">switchNums</span><span style="color: #0000FF;">({</span><span style="color: #000000;">c</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">}),</span><span style="color: #000000;">edgePointNums</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">newFaces</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">newFaces</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgePointAb</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">facePointAbcd</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgePointDa</span><span style="color: #0000FF;">})</span>
<span style="color: #000000;">newFaces</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">newFaces</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgePointBc</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">facePointAbcd</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgePointAb</span><span style="color: #0000FF;">})</span>
<span style="color: #000000;">newFaces</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">newFaces</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">c</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgePointCd</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">facePointAbcd</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgePointBc</span><span style="color: #0000FF;">})</span>
<span style="color: #000000;">newFaces</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">newFaces</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">d</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgePointDa</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">facePointAbcd</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">edgePointCd</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">newPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">newFaces</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">inputPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}},</span>
<span style="color: #000000;">inputFaces</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">}}</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">outputPoints</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">inputPoints</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">outputFaces</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inputFaces</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">iterations</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">iterations</span> <span style="color: #008080;">do</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">outputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">outputFaces</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">cmcSubdiv</span><span style="color: #0000FF;">(</span><span style="color: #000000;">outputPoints</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">outputFaces</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #7060A8;">ppOpt</span><span style="color: #0000FF;">({</span><span style="color: #004600;">pp_Nest</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_FltFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%7.4f"</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_IntFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%7.4f"</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_IntCh</span><span style="color: #0000FF;">,</span><span style="color: #004600;">false</span><span style="color: #0000FF;">})</span>
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">outputPoints</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">outputFaces</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),{</span><span style="color: #004600;">pp_IntFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%2d"</span><span style="color: #0000FF;">})</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
{{-0.5556, 0.5556, 0.5556},
{-0.5556,-0.5556, 0.5556},
{ 0.5556,-0.5556, 0.5556},
{ 0.5556, 0.5556, 0.5556},
{ 0.5556,-0.5556,-0.5556},
{ 0.5556, 0.5556,-0.5556},
{-0.5556,-0.5556,-0.5556},
{-0.5556, 0.5556,-0.5556},
{ 0.0000, 0.0000, 1.0000},
{ 1.0000, 0.0000, 0.0000},
{ 0.0000, 0.0000,-1.0000},
{ 0.0000, 1.0000, 0.0000},
{-1.0000, 0.0000, 0.0000},
{ 0.0000,-1.0000, 0.0000},
{-0.7500, 0.0000, 0.7500},
{ 0.0000, 0.7500, 0.7500},
{-0.7500, 0.7500, 0.0000},
{ 0.0000,-0.7500, 0.7500},
{-0.7500,-0.7500, 0.0000},
{ 0.7500, 0.0000, 0.7500},
{ 0.7500,-0.7500, 0.0000},
{ 0.7500, 0.7500, 0.0000},
{ 0.7500, 0.0000,-0.7500},
{ 0.0000,-0.7500,-0.7500},
{ 0.0000, 0.7500,-0.7500},
{-0.7500, 0.0000,-0.7500}}
{{ 0,14, 8,15},
{ 1,17, 8,14},
{ 2,19, 8,17},
{ 3,15, 8,19},
{ 3,19, 9,21},
{ 2,20, 9,19},
{ 4,22, 9,20},
{ 5,21, 9,22},
{ 5,22,10,24},
{ 4,23,10,22},
{ 6,25,10,23},
{ 7,24,10,25},
{ 7,16,11,24},
{ 0,15,11,16},
{ 3,21,11,15},
{ 5,24,11,21},
{ 7,25,12,16},
{ 6,18,12,25},
{ 1,14,12,18},
{ 0,16,12,14},
{ 6,18,13,23},
{ 1,17,13,18},
{ 2,20,13,17},
{ 4,23,13,20}}
</pre>
=={{header|Python}}==
<syntaxhighlight lang="python">
"""
 
Input and output are assumed to be in this form based on the talk
page for the task:
 
input_points = [
[-1.0, 1.0, 1.0],
[-1.0, -1.0, 1.0],
[ 1.0, -1.0, 1.0],
[ 1.0, 1.0, 1.0],
[ 1.0, -1.0, -1.0],
[ 1.0, 1.0, -1.0],
[-1.0, -1.0, -1.0],
[-1.0, 1.0, -1.0]
]
 
input_faces = [
[0, 1, 2, 3],
[3, 2, 4, 5],
[5, 4, 6, 7],
[7, 0, 3, 5],
[7, 6, 1, 0],
[6, 1, 2, 4],
]
 
So, the graph is a list of points and a list of faces.
Each face is a list of indexes into the points list.
"""
 
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
import sys
 
def center_point(p1, p2):
"""
returns a point in the center of the
segment ended by points p1 and p2
"""
cp = []
for i in range(3):
cp.append((p1[i]+p2[i])/2)
return cp
def sum_point(p1, p2):
"""
adds points p1 and p2
"""
sp = []
for i in range(3):
sp.append(p1[i]+p2[i])
return sp
 
def div_point(p, d):
"""
divide point p by d
"""
sp = []
for i in range(3):
sp.append(p[i]/d)
return sp
def mul_point(p, m):
"""
multiply point p by m
"""
sp = []
for i in range(3):
sp.append(p[i]*m)
return sp
 
def get_face_points(input_points, input_faces):
"""
From http://rosettacode.org/wiki/Catmull%E2%80%93Clark_subdivision_surface
1. for each face, a face point is created which is the average of all the points of the face.
"""
 
# 3 dimensional space
NUM_DIMENSIONS = 3
# face_points will have one point for each face
face_points = []
for curr_face in input_faces:
face_point = [0.0, 0.0, 0.0]
for curr_point_index in curr_face:
curr_point = input_points[curr_point_index]
# add curr_point to face_point
# will divide later
for i in range(NUM_DIMENSIONS):
face_point[i] += curr_point[i]
# divide by number of points for average
num_points = len(curr_face)
for i in range(NUM_DIMENSIONS):
face_point[i] /= num_points
face_points.append(face_point)
return face_points
def get_edges_faces(input_points, input_faces):
"""
Get list of edges and the one or two adjacent faces in a list.
also get center point of edge
Each edge would be [pointnum_1, pointnum_2, facenum_1, facenum_2, center]
"""
# will have [pointnum_1, pointnum_2, facenum]
edges = []
# get edges from each face
for facenum in range(len(input_faces)):
face = input_faces[facenum]
num_points = len(face)
# loop over index into face
for pointindex in range(num_points):
# if not last point then edge is curr point and next point
if pointindex < num_points - 1:
pointnum_1 = face[pointindex]
pointnum_2 = face[pointindex+1]
else:
# for last point edge is curr point and first point
pointnum_1 = face[pointindex]
pointnum_2 = face[0]
# order points in edge by lowest point number
if pointnum_1 > pointnum_2:
temp = pointnum_1
pointnum_1 = pointnum_2
pointnum_2 = temp
edges.append([pointnum_1, pointnum_2, facenum])
# sort edges by pointnum_1, pointnum_2, facenum
edges = sorted(edges)
# merge edges with 2 adjacent faces
# [pointnum_1, pointnum_2, facenum_1, facenum_2] or
# [pointnum_1, pointnum_2, facenum_1, None]
num_edges = len(edges)
eindex = 0
merged_edges = []
while eindex < num_edges:
e1 = edges[eindex]
# check if not last edge
if eindex < num_edges - 1:
e2 = edges[eindex+1]
if e1[0] == e2[0] and e1[1] == e2[1]:
merged_edges.append([e1[0],e1[1],e1[2],e2[2]])
eindex += 2
else:
merged_edges.append([e1[0],e1[1],e1[2],None])
eindex += 1
else:
merged_edges.append([e1[0],e1[1],e1[2],None])
eindex += 1
# add edge centers
edges_centers = []
for me in merged_edges:
p1 = input_points[me[0]]
p2 = input_points[me[1]]
cp = center_point(p1, p2)
edges_centers.append(me+[cp])
return edges_centers
def get_edge_points(input_points, edges_faces, face_points):
"""
for each edge, an edge point is created which is the average
between the center of the edge and the center of the segment made
with the face points of the two adjacent faces.
"""
edge_points = []
for edge in edges_faces:
# get center of edge
cp = edge[4]
# get center of two facepoints
fp1 = face_points[edge[2]]
# if not two faces just use one facepoint
# should not happen for solid like a cube
if edge[3] == None:
fp2 = fp1
else:
fp2 = face_points[edge[3]]
cfp = center_point(fp1, fp2)
# get average between center of edge and
# center of facepoints
edge_point = center_point(cp, cfp)
edge_points.append(edge_point)
return edge_points
def get_avg_face_points(input_points, input_faces, face_points):
"""
for each point calculate
the average of the face points of the faces the point belongs to (avg_face_points)
create a list of lists of two numbers [facepoint_sum, num_points] by going through the
points in all the faces.
then create the avg_face_points list of point by dividing point_sum (x, y, z) by num_points
"""
# initialize list with [[0.0, 0.0, 0.0], 0]
num_points = len(input_points)
temp_points = []
for pointnum in range(num_points):
temp_points.append([[0.0, 0.0, 0.0], 0])
# loop through faces updating temp_points
for facenum in range(len(input_faces)):
fp = face_points[facenum]
for pointnum in input_faces[facenum]:
tp = temp_points[pointnum][0]
temp_points[pointnum][0] = sum_point(tp,fp)
temp_points[pointnum][1] += 1
# divide to create avg_face_points
avg_face_points = []
for tp in temp_points:
afp = div_point(tp[0], tp[1])
avg_face_points.append(afp)
return avg_face_points
def get_avg_mid_edges(input_points, edges_faces):
"""
the average of the centers of edges the point belongs to (avg_mid_edges)
create list with entry for each point
each entry has two elements. one is a point that is the sum of the centers of the edges
and the other is the number of edges. after going through all edges divide by
number of edges.
"""
# initialize list with [[0.0, 0.0, 0.0], 0]
num_points = len(input_points)
temp_points = []
for pointnum in range(num_points):
temp_points.append([[0.0, 0.0, 0.0], 0])
# go through edges_faces using center updating each point
for edge in edges_faces:
cp = edge[4]
for pointnum in [edge[0], edge[1]]:
tp = temp_points[pointnum][0]
temp_points[pointnum][0] = sum_point(tp,cp)
temp_points[pointnum][1] += 1
# divide out number of points to get average
avg_mid_edges = []
for tp in temp_points:
ame = div_point(tp[0], tp[1])
avg_mid_edges.append(ame)
return avg_mid_edges
 
def get_points_faces(input_points, input_faces):
# initialize list with 0
num_points = len(input_points)
points_faces = []
for pointnum in range(num_points):
points_faces.append(0)
# loop through faces updating points_faces
for facenum in range(len(input_faces)):
for pointnum in input_faces[facenum]:
points_faces[pointnum] += 1
return points_faces
 
def get_new_points(input_points, points_faces, avg_face_points, avg_mid_edges):
"""
m1 = (n - 3.0) / n
m2 = 1.0 / n
m3 = 2.0 / n
new_coords = (m1 * old_coords)
+ (m2 * avg_face_points)
+ (m3 * avg_mid_edges)
"""
new_points =[]
for pointnum in range(len(input_points)):
n = points_faces[pointnum]
m1 = (n - 3.0) / n
m2 = 1.0 / n
m3 = 2.0 / n
old_coords = input_points[pointnum]
p1 = mul_point(old_coords, m1)
afp = avg_face_points[pointnum]
p2 = mul_point(afp, m2)
ame = avg_mid_edges[pointnum]
p3 = mul_point(ame, m3)
p4 = sum_point(p1, p2)
new_coords = sum_point(p4, p3)
new_points.append(new_coords)
return new_points
def switch_nums(point_nums):
"""
Returns tuple of point numbers
sorted least to most
"""
if point_nums[0] < point_nums[1]:
return point_nums
else:
return (point_nums[1], point_nums[0])
 
def cmc_subdiv(input_points, input_faces):
# 1. for each face, a face point is created which is the average of all the points of the face.
# each entry in the returned list is a point (x, y, z).
face_points = get_face_points(input_points, input_faces)
# get list of edges with 1 or 2 adjacent faces
# [pointnum_1, pointnum_2, facenum_1, facenum_2, center] or
# [pointnum_1, pointnum_2, facenum_1, None, center]
edges_faces = get_edges_faces(input_points, input_faces)
# get edge points, a list of points
edge_points = get_edge_points(input_points, edges_faces, face_points)
# the average of the face points of the faces the point belongs to (avg_face_points)
avg_face_points = get_avg_face_points(input_points, input_faces, face_points)
# the average of the centers of edges the point belongs to (avg_mid_edges)
avg_mid_edges = get_avg_mid_edges(input_points, edges_faces)
# how many faces a point belongs to
points_faces = get_points_faces(input_points, input_faces)
"""
m1 = (n - 3) / n
m2 = 1 / n
m3 = 2 / n
new_coords = (m1 * old_coords)
+ (m2 * avg_face_points)
+ (m3 * avg_mid_edges)
"""
new_points = get_new_points(input_points, points_faces, avg_face_points, avg_mid_edges)
"""
Then each face is replaced by new faces made with the new points,
for a triangle face (a,b,c):
(a, edge_point ab, face_point abc, edge_point ca)
(b, edge_point bc, face_point abc, edge_point ab)
(c, edge_point ca, face_point abc, edge_point bc)
for a quad face (a,b,c,d):
(a, edge_point ab, face_point abcd, edge_point da)
(b, edge_point bc, face_point abcd, edge_point ab)
(c, edge_point cd, face_point abcd, edge_point bc)
(d, edge_point da, face_point abcd, edge_point cd)
face_points is a list indexed by face number so that is
easy to get.
edge_points is a list indexed by the edge number
which is an index into edges_faces.
need to add face_points and edge points to
new_points and get index into each.
then create two new structures
face_point_nums - list indexes by facenum
whose value is the index into new_points
edge_point num - dictionary with key (pointnum_1, pointnum_2)
and value is index into new_points
"""
# add face points to new_points
face_point_nums = []
# point num after next append to new_points
next_pointnum = len(new_points)
for face_point in face_points:
new_points.append(face_point)
face_point_nums.append(next_pointnum)
next_pointnum += 1
# add edge points to new_points
edge_point_nums = dict()
for edgenum in range(len(edges_faces)):
pointnum_1 = edges_faces[edgenum][0]
pointnum_2 = edges_faces[edgenum][1]
edge_point = edge_points[edgenum]
new_points.append(edge_point)
edge_point_nums[(pointnum_1, pointnum_2)] = next_pointnum
next_pointnum += 1
# new_points now has the points to output. Need new
# faces
"""
just doing this case for now:
for a quad face (a,b,c,d):
(a, edge_point ab, face_point abcd, edge_point da)
(b, edge_point bc, face_point abcd, edge_point ab)
(c, edge_point cd, face_point abcd, edge_point bc)
(d, edge_point da, face_point abcd, edge_point cd)
new_faces will be a list of lists where the elements are like this:
[pointnum_1, pointnum_2, pointnum_3, pointnum_4]
"""
new_faces =[]
for oldfacenum in range(len(input_faces)):
oldface = input_faces[oldfacenum]
# 4 point face
if len(oldface) == 4:
a = oldface[0]
b = oldface[1]
c = oldface[2]
d = oldface[3]
face_point_abcd = face_point_nums[oldfacenum]
edge_point_ab = edge_point_nums[switch_nums((a, b))]
edge_point_da = edge_point_nums[switch_nums((d, a))]
edge_point_bc = edge_point_nums[switch_nums((b, c))]
edge_point_cd = edge_point_nums[switch_nums((c, d))]
new_faces.append((a, edge_point_ab, face_point_abcd, edge_point_da))
new_faces.append((b, edge_point_bc, face_point_abcd, edge_point_ab))
new_faces.append((c, edge_point_cd, face_point_abcd, edge_point_bc))
new_faces.append((d, edge_point_da, face_point_abcd, edge_point_cd))
return new_points, new_faces
 
def graph_output(output_points, output_faces):
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
"""
Plot each face
"""
for facenum in range(len(output_faces)):
curr_face = output_faces[facenum]
xcurr = []
ycurr = []
zcurr = []
for pointnum in range(len(curr_face)):
xcurr.append(output_points[curr_face[pointnum]][0])
ycurr.append(output_points[curr_face[pointnum]][1])
zcurr.append(output_points[curr_face[pointnum]][2])
xcurr.append(output_points[curr_face[0]][0])
ycurr.append(output_points[curr_face[0]][1])
zcurr.append(output_points[curr_face[0]][2])
ax.plot(xcurr,ycurr,zcurr,color='b')
plt.show()
 
 
# cube
 
input_points = [
[-1.0, 1.0, 1.0],
[-1.0, -1.0, 1.0],
[ 1.0, -1.0, 1.0],
[ 1.0, 1.0, 1.0],
[ 1.0, -1.0, -1.0],
[ 1.0, 1.0, -1.0],
[-1.0, -1.0, -1.0],
[-1.0, 1.0, -1.0]
]
 
input_faces = [
[0, 1, 2, 3],
[3, 2, 4, 5],
[5, 4, 6, 7],
[7, 0, 3, 5],
[7, 6, 1, 0],
[6, 1, 2, 4],
]
 
if len(sys.argv) != 2:
print("Should have one argument integer number of iterations")
sys.exit()
else:
iterations = int(sys.argv[1])
output_points, output_faces = input_points, input_faces
for i in range(iterations):
output_points, output_faces = cmc_subdiv(output_points, output_faces)
graph_output(output_points, output_faces)
</syntaxhighlight>
=={{header|Rust}}==
{{trans|Python}}
<syntaxhighlight lang="rust">
pub struct Vector3 {pub x: f64, pub y: f64, pub z: f64, pub w: f64}
 
pub struct Triangle {pub r: [usize; 3], pub(crate) col: [f32; 4], pub(crate) p: [Vector3; 3], n: Vector3, pub t: [Vector2; 3]}
pub struct Mesh{pub v: Vec<Vector3>, pub f: Vec<Triangle>}
 
impl Triangle{
pub fn new() -> Triangle {return Triangle {r: [0, 0, 0], col: [0.0; 4], p: [Vector3::new(0.0, 0.0, 0.0), Vector3::new(0.0, 0.0, 0.0), Vector3::new(0.0, 0.0, 0.0)], n: Vector3::new(0.0, 0.0, 0.0), t: [Vector2::new(0.0, 0.0), Vector2::new(0.0, 0.0), Vector2::new(0.0, 0.0)]}}
pub fn copy(&self) -> Triangle {return Triangle {r: self.r.clone(), col: self.col, p: [self.p[0].copy(), self.p[1].copy(), self.p[2].copy()], n: self.n.copy(), t: [self.t[0].copy(), self.t[1].copy(), self.t[2].copy()]}}
}
 
 
impl Vector3 {
pub fn new(x: f64, y: f64, z: f64) -> Vector3 {return Vector3 {x, y, z, w: 1.0}}
pub fn normalize(&mut self) {
let l = (self.x * self.x + self.y * self.y + self.z * self.z).sqrt();
self.x /= l;
self.y /= l;
self.z /= l;
}
pub fn dot_product(v1: &Vector3, v2: &Vector3) -> f64 {return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z}
pub fn cross_product(v1: &Vector3, v2: &Vector3) -> Vector3 {return Vector3::new(v1.y * v2.z - v1.z * v2.y, v1.z * v2.x - v1.x * v2.z, v1.x * v2.y - v1.y * v2.x)}
pub fn intersect_plane(plane_n: &Vector3, plane_p: &Vector3, line_start: &Vector3, line_end: &Vector3, mut t: f64) -> Vector3 {
let mut p_n = plane_n.copy();
p_n.normalize();
let plane_d = -Vector3::dot_product(&p_n, plane_p);
let ad = Vector3::dot_product(line_start, &p_n);
let bd = Vector3::dot_product(line_end, &p_n);
t = (-plane_d - ad) / (bd - ad);
let line = line_end.copy() - line_start;
let line_i = line * t;
return line_start.copy() + &line_i;
}
pub fn copy(&self) -> Vector3 {return Vector3 {x: self.x, y: self.y, z: self.z, w: self.w}}
}
 
impl Mesh {
pub fn get_face_points(&self) -> Vec<Vector3> {
let mut face_points: Vec<Vector3> = Vec::new();
 
for curr_face in &self.f {
let mut face_point = Vector3::new(0.0, 0.0, 0.0);
for curr_point_index in curr_face.r {
let curr_point = &self.v[curr_point_index];
face_point += &curr_point
}
 
face_point /= curr_face.r.len() as f64;
face_points.push(face_point.copy());
}
return face_points;
}
pub fn get_edges_faces(&self) -> Vec<[f64; 7]> {
let mut edges: Vec<[usize; 3]> = Vec::new();
 
for face_num in 0..self.f.len() {
let face: Triangle = self.f[face_num].copy();
let num_points = face.p.len();
for point_index in 0..num_points {
let mut point_num_1 = 0;
let mut point_num_2 = 0;
if point_index < num_points - 1 {
point_num_1 = face.r[point_index];
point_num_2 = face.r[point_index + 1];
} else {
point_num_1 = face.r[point_index];
point_num_2 = face.r[0];
}
if point_num_1 > point_num_2 {
let temp = point_num_1;
point_num_1 = point_num_2;
point_num_2 = temp;
}
edges.push([point_num_1, point_num_2, face_num]);
}
}
edges.sort();
 
let num_edges = edges.len();
let mut index = 0;
let mut merged_edges: Vec<[f64; 4]> = Vec::new();
 
while index < num_edges {
let e1 = edges[index];
if index < num_edges - 1 {
let e2 = edges[index + 1];
if e1[0] == e2[0] && e1[1] == e2[1] {
merged_edges.push([e1[0] as f64, e1[1] as f64, e1[2] as f64, e2[2] as f64]);
index += 2;
} else {
merged_edges.push([e1[0] as f64, e1[1] as f64, e1[2] as f64, -1.0]);
index += 1;
}
} else {
merged_edges.push([e1[0] as f64, e1[1] as f64, e1[2] as f64, -1.0]);
index += 1
}
}
 
let mut edges_centers = Vec::new();
 
for me in merged_edges {
let p1 = self.v[me[0] as usize].copy();
let p2 = self.v[me[1] as usize].copy();
let cp: Vector3 = Mesh::center_point(&p1, &p2);
edges_centers.push([me[0] as f64, me[1] as f64, me[2] as f64, me[3] as f64, cp.x, cp.y, cp.z]);
}
return edges_centers;
}
pub fn get_edge_points(&self, edges_faces: &Vec<[f64; 7]>, face_points: &Vec<Vector3>) -> Vec<Vector3> {
let mut edge_points = Vec::new();
 
for edge in edges_faces {
let cp = Vector3::new(edge[4], edge[5], edge[6]);
let fp1: Vector3 = face_points[edge[2] as usize].copy();
let mut fp2: Vector3 = fp1.copy();
if edge[3] != -1.0 { fp2 = face_points[edge[3] as usize].copy() };
let cfp = Mesh::center_point(&fp1, &fp2);
let edge_point = Mesh::center_point(&cp, &cfp);
edge_points.push(edge_point);
}
 
return edge_points
}
pub fn get_avg_face_points(&self, face_points: &Vec<Vector3>) -> Vec<Vector3> {
let num_points = self.v.len();
let mut temp_points = Vec::new();
let mut div: Vec<i32> = Vec::new();
 
for _ in 0..num_points {
temp_points.push(Vector3::new(0.0, 0.0, 0.0));
div.push(0)
};
 
for face_num in 0..self.f.len() {
let mut fp = face_points[face_num].copy();
for point_num in self.f[face_num].r {
let tp = temp_points[point_num].copy();
temp_points[point_num] = tp + &fp;
div[point_num] += 1;
}
}
 
let mut avg_face_points: Vec<Vector3> = Vec::new();
for i in 0..temp_points.len() {
let tp: Vector3 = temp_points[i].copy();
let t = tp / (div[i] as f64);
avg_face_points.push(t.copy());
}
 
return avg_face_points;
}
pub fn get_avg_mid_edges(&self, edges_faces: &Vec<[f64; 7]>) -> Vec<Vector3> {
let num_points = self.v.len();
let mut temp_points = Vec::new();
let mut div: Vec<i32> = Vec::new();
 
for point_num in 0..num_points{ temp_points.push(Vector3::new(0.0, 0.0, 0.0)); div.push(0)}
for edge in edges_faces {
let cp = Vector3::new(edge[4], edge[5], edge[6]);
for point_num in [edge[0] as usize, edge[1] as usize] {
let tp = temp_points[point_num].copy();
temp_points[point_num] = tp + &cp;
div[point_num] += 1
}
}
 
let mut avg_mid_edges: Vec<Vector3> = Vec::new();
 
for i in 0..temp_points.len(){
let ame: Vector3 = temp_points[i].copy() / (div[i] as f64);
avg_mid_edges.push(ame)}
 
return avg_mid_edges
}
pub fn get_points_faces(&self) -> Vec<i32> {
let num_points = self.v.len();
let mut points_faces: Vec<i32> = Vec::new();
 
for point_num in 0..num_points{points_faces.push(0)}
 
for face_num in 0..self.f.len() {
for point_num in self.f[face_num].r {
points_faces[point_num] += 1;
}
}
return points_faces
}
pub fn get_new_points(&self, points_faces: &Vec<i32>, avg_face_points: &Vec<Vector3>, avg_mid_edges: &Vec<Vector3>) -> Vec<Vector3> {
let mut new_points: Vec<Vector3> = Vec::new();
 
for point_num in 0..self.v.len() {
let n = points_faces[point_num] as f64;
let m1 = (n - 3.0) / n;
let m2 = 1.0 / n;
let m3 = 2.0 / n;
let old_coords = self.v[point_num].copy();
let p1 = old_coords * m1;
let afp = avg_face_points[point_num].copy();
let p2 = afp * m2;
let ame = avg_mid_edges[point_num].copy();
let p3 = ame * m3;
let p4 = p1 + &p2;
let new_coords = p4 + &p3;
 
new_points.push(new_coords);
}
 
return new_points;
}
 
pub fn switch_nums(point_nums: [f64; 2]) -> [f64; 2] {
return if point_nums[0] < point_nums[1] { point_nums } else {[point_nums[1], point_nums[0]]}
}
 
pub fn get_key(points: [f64; 2]) -> String {
return points[0].to_string() + ";" + &*points[1].to_string();
}
 
pub fn subdivide(&mut self) {
let face_points = self.get_face_points();
let edges_faces = self.get_edges_faces();
let edge_points = self.get_edge_points(&edges_faces, &face_points);
let avg_face_points = self.get_avg_face_points(&face_points);
let avg_mid_edges = self.get_avg_mid_edges(&edges_faces);
let points_faces = self.get_points_faces();
let mut new_points = self.get_new_points(&points_faces, &avg_face_points, &avg_mid_edges);
 
let mut face_point_nums = Vec::new();
let mut next_point_num = new_points.len();
 
for face_point in face_points {
new_points.push(face_point);
face_point_nums.push(next_point_num);
next_point_num += 1;
}
 
let mut edge_point_nums: HashMap<String, usize> = HashMap::new();
 
for edge_num in 0..edges_faces.len() {
let point_num_1 = edges_faces[edge_num][0];
let point_num_2 = edges_faces[edge_num][1];
let edge_point = edge_points[edge_num].copy();
new_points.push(edge_point);
edge_point_nums.insert(Mesh::get_key([point_num_1, point_num_2]), next_point_num);
next_point_num += 1;
}
 
let mut new_faces = Vec::new();
 
for old_face_num in 0..self.f.len() {
let old_face = self.f[old_face_num].copy();
let a = old_face.r[0] as f64;
let b = old_face.r[1] as f64;
let c = old_face.r[2] as f64;
let face_point_abc = face_point_nums[old_face_num];
let edge_point_ab = *edge_point_nums.get(&*Mesh::get_key(Mesh::switch_nums([a, b]))).unwrap();
let edge_point_bc = *edge_point_nums.get(&*Mesh::get_key(Mesh::switch_nums([b, c]))).unwrap();
let edge_point_ca = *edge_point_nums.get(&*Mesh::get_key(Mesh::switch_nums([c, a]))).unwrap();
new_faces.push([a, edge_point_ab as f64, face_point_abc as f64, edge_point_ca as f64]);
new_faces.push([b, edge_point_bc as f64, face_point_abc as f64, edge_point_ab as f64]);
new_faces.push([c, edge_point_ca as f64, face_point_abc as f64, edge_point_bc as f64]);
}
self.f = Vec::new();
 
for face_num in 0..new_faces.len() {
let curr_face = new_faces[face_num];
let mut t1: Triangle = Triangle::new();
let mut t2: Triangle = Triangle::new();
t1.p[0] = Vector3::new(new_points[curr_face[0] as usize].x, new_points[curr_face[0] as usize].y, new_points[curr_face[0] as usize].z);
t1.p[1] = Vector3::new(new_points[curr_face[1] as usize].x, new_points[curr_face[1] as usize].y, new_points[curr_face[1] as usize].z);
t1.p[2] = Vector3::new(new_points[curr_face[2] as usize].x, new_points[curr_face[2] as usize].y, new_points[curr_face[2] as usize].z);
t2.p[0] = Vector3::new(new_points[curr_face[2] as usize].x, new_points[curr_face[2] as usize].y, new_points[curr_face[2] as usize].z);
t2.p[1] = Vector3::new(new_points[curr_face[3] as usize].x, new_points[curr_face[3] as usize].y, new_points[curr_face[3] as usize].z);
t2.p[2] = Vector3::new(new_points[curr_face[0] as usize].x, new_points[curr_face[0] as usize].y, new_points[curr_face[0] as usize].z);
t1.r = [curr_face[0] as usize, curr_face[1] as usize, curr_face[2] as usize];
t2.r = [curr_face[2] as usize, curr_face[3] as usize, curr_face[0] as usize];
self.f.push(t1);
self.f.push(t2);
}
self.v = new_points;
}
}
</syntaxhighlight>
=={{header|Tcl}}==
This code handles both holes and arbitrary polygons in the input data.
<langsyntaxhighlight lang="tcl">package require Tcl 8.5
 
# Use math functions and operators as commands (Lisp-like).
Line 365 ⟶ 2,995:
foreach e $edges {
set ep [selectFrom $points $e]
set mid [centroid $ep]
if {[llength $fp4e($e)] > 1} {
lappendset edgepointsmid [centroid $ep $fp4e($e)]
} else {
lappendset edgepointsmid [centroid $midep]
foreach p $ep {
lappend ep_heavy($p) $mid
}
}
lappend edgepoints $mid
set en4e($e) $i
foreach p $ep {
Line 380 ⟶ 3,013:
# Generate the new vertex points with our lookup tables.
foreach p $points {
set nif {[llength $fp4p($p)] >= 4} {
if {$n == set n [llength $ep4pfp4p($p)]} {
lappend newPoints [add3 [mulC [/ [- $n 3.0] $n] $p] \
[mulC [/ 1.0 $n] [centroid $fp4p($p)]] \
[mulC [/ 2.0 $n] [centroid $ep4p($p)]]]
} else {
# Update a point on the edge of a hole. This formula is not
lappend newPoints [centroid [list $p] $ep4p($p)]
# described on the WP page, but produces a nice result.
}
lappend newPoints [centroid $ep_heavy($p) [list $p $p]]
}
}
 
Line 401 ⟶ 3,036:
 
list [concat $newPoints $facepoints $edgepoints] $newFaces
}</langsyntaxhighlight>
 
The [[/Tcl Test Code|test code]] for this solution is available as well. The example there produces the following partial toroid output image:
 
<center>[[File:Tcl-Catmull.png]]</center>
=={{header|Wren}}==
{{trans|Go}}
{{libheader|Wren-dynamic}}
{{libheader|Wren-sort}}
{{libheader|Wren-math}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./dynamic" for Tuple, Struct
import "./sort" for Sort
import "./math" for Int
import "./fmt" for Fmt
 
var Point = Tuple.create("Point", ["x", "y", "z"])
var fields = [
"pn1", // point number 1
"pn2", // point number 2
"fn1", // face number 1
"fn2", // face number 2
"cp" // center point
]
var Edge = Tuple.create("Edge", fields)
var PointEx = Struct.create("PointEx", ["p", "n"])
 
var sumPoint = Fn.new { |p1, p2| Point.new(p1.x + p2.x, p1.y + p2.y, p1.z + p2.z) }
 
var mulPoint = Fn.new { |p, m| Point.new(p.x * m, p.y * m, p.z * m) }
 
var divPoint = Fn.new { |p, d| mulPoint.call(p, 1/d) }
 
var centerPoint = Fn.new { |p1, p2| divPoint.call(sumPoint.call(p1, p2), 2) }
 
var getFacePoints = Fn.new { |inputPoints, inputFaces|
var facePoints = List.filled(inputFaces.count, null)
var i = 0
for (currFace in inputFaces) {
var facePoint = Point.new(0, 0, 0)
for (cpi in currFace) {
var currPoint = inputPoints[cpi]
facePoint = sumPoint.call(facePoint, currPoint)
}
facePoints[i] = divPoint.call(facePoint, currFace.count)
i = i + 1
}
return facePoints
}
 
var getEdgesFaces = Fn.new { |inputPoints, inputFaces|
var edges = []
var faceNum = 0
for (face in inputFaces) {
var numPoints = face.count
for (pointIndex in 0...numPoints) {
var pointNum1 = face[pointIndex]
var pointNum2 = (pointIndex < numPoints-1) ? face[pointIndex+1] : face[0]
if (pointNum1 > pointNum2) {
var t = pointNum1
pointNum1 = pointNum2
pointNum2 = t
}
edges.add([pointNum1, pointNum2, faceNum])
}
faceNum = faceNum + 1
}
var cmp = Fn.new { |e1, e2|
if (e1[0] == e2[0]) {
if (e1[1] == e2[1]) return (e1[2] - e2[2]).sign
return (e1[1] - e2[1]).sign
}
return (e1[0] - e2[0]).sign
}
var numEdges = edges.count
Sort.quick(edges, 0, numEdges-1, cmp)
var eIndex = 0
var mergedEdges = []
while (eIndex < numEdges) {
var e1 = edges[eIndex]
if (eIndex < numEdges-1) {
var e2 = edges[eIndex+1]
if (e1[0] == e2[0] && e1[1] == e2[1]) {
mergedEdges.add([e1[0], e1[1], e1[2], e2[2]])
eIndex = eIndex + 2
} else {
mergedEdges.add([e1[0], e1[1], e1[2], -1])
eIndex = eIndex + 1
}
} else {
mergedEdges.add([e1[0], e1[1], e1[2], -1])
eIndex = eIndex + 1
}
}
var edgesCenters = []
for (me in mergedEdges) {
var p1 = inputPoints[me[0]]
var p2 = inputPoints[me[1]]
var cp = centerPoint.call(p1, p2)
edgesCenters.add(Edge.new(me[0], me[1], me[2], me[3], cp))
}
return edgesCenters
}
 
var getEdgePoints = Fn.new { |inputPoints, edgesFaces, facePoints|
var edgePoints = List.filled(edgesFaces.count, null)
var i = 0
for (edge in edgesFaces) {
var cp = edge.cp
var fp1 = facePoints[edge.fn1]
var fp2 = (edge.fn2 == -1) ? fp1 : facePoints[edge.fn2]
var cfp = centerPoint.call(fp1, fp2)
edgePoints[i] = centerPoint.call(cp, cfp)
i = i + 1
}
return edgePoints
}
 
var getAvgFacePoints = Fn.new { |inputPoints, inputFaces, facePoints|
var numPoints = inputPoints.count
var tempPoints = List.filled(numPoints, null)
for (i in 0...numPoints) tempPoints[i] = PointEx.new(Point.new(0, 0, 0), 0)
for (faceNum in 0...inputFaces.count) {
var fp = facePoints[faceNum]
for (pointNum in inputFaces[faceNum]) {
var tp = tempPoints[pointNum].p
tempPoints[pointNum].p = sumPoint.call(tp, fp)
tempPoints[pointNum].n = tempPoints[pointNum].n + 1
}
}
var avgFacePoints = List.filled(numPoints, null)
var i = 0
for (tp in tempPoints) {
avgFacePoints[i] = divPoint.call(tp.p, tp.n)
i = i + 1
}
return avgFacePoints
}
 
var getAvgMidEdges = Fn.new { |inputPoints, edgesFaces|
var numPoints = inputPoints.count
var tempPoints = List.filled(numPoints, null)
for (i in 0...numPoints) tempPoints[i] = PointEx.new(Point.new(0, 0, 0), 0)
for (edge in edgesFaces) {
var cp = edge.cp
for (pointNum in [edge.pn1, edge.pn2]) {
var tp = tempPoints[pointNum].p
tempPoints[pointNum].p = sumPoint.call(tp, cp)
tempPoints[pointNum].n = tempPoints[pointNum].n + 1
}
}
var avgMidEdges = List.filled(tempPoints.count, null)
var i = 0
for (tp in tempPoints) {
avgMidEdges[i] = divPoint.call(tp.p, tp.n)
i = i + 1
}
return avgMidEdges
}
 
var getPointsFaces = Fn.new { |inputPoints, inputFaces|
var numPoints = inputPoints.count
var pointsFaces = List.filled(numPoints, 0)
for (faceNum in 0...inputFaces.count) {
for (pointNum in inputFaces[faceNum]) {
pointsFaces[pointNum] = pointsFaces[pointNum] + 1
}
}
return pointsFaces
}
 
var getNewPoints = Fn.new { |inputPoints, pointsFaces, avgFacePoints, avgMidEdges|
var newPoints = List.filled(inputPoints.count, null)
for (pointNum in 0...inputPoints.count) {
var n = pointsFaces[pointNum]
var m1 = (n-3) / n
var m2 = 1 / n
var m3 = 2 / n
var oldCoords = inputPoints[pointNum]
var p1 = mulPoint.call(oldCoords, m1)
var afp = avgFacePoints[pointNum]
var p2 = mulPoint.call(afp, m2)
var ame = avgMidEdges[pointNum]
var p3 = mulPoint.call(ame, m3)
var p4 = sumPoint.call(p1, p2)
newPoints[pointNum] = sumPoint.call(p4, p3)
}
return newPoints
}
 
var switchNums = Fn.new { |pointNums|
if (pointNums[0] < pointNums[1]) return pointNums
return [pointNums[1], pointNums[0]]
}
 
var cmcSubdiv = Fn.new { |inputPoints, inputFaces|
var facePoints = getFacePoints.call(inputPoints, inputFaces)
var edgesFaces = getEdgesFaces.call(inputPoints, inputFaces)
var edgePoints = getEdgePoints.call(inputPoints, edgesFaces, facePoints)
var avgFacePoints = getAvgFacePoints.call(inputPoints, inputFaces, facePoints)
var avgMidEdges = getAvgMidEdges.call(inputPoints, edgesFaces)
var pointsFaces = getPointsFaces.call(inputPoints, inputFaces)
var newPoints = getNewPoints.call(inputPoints, pointsFaces, avgFacePoints, avgMidEdges)
var facePointNums = []
var nextPointNum = newPoints.count
for (facePoint in facePoints) {
newPoints.add(facePoint)
facePointNums.add(nextPointNum)
nextPointNum = nextPointNum + 1
}
var edgePointNums = {}
for (edgeNum in 0...edgesFaces.count) {
var pointNum1 = edgesFaces[edgeNum].pn1
var pointNum2 = edgesFaces[edgeNum].pn2
var edgePoint = edgePoints[edgeNum]
newPoints.add(edgePoint)
edgePointNums[Int.cantorPair(pointNum1, pointNum2)] = nextPointNum
nextPointNum = nextPointNum + 1
}
var newFaces = []
var oldFaceNum = 0
for (oldFace in inputFaces) {
if (oldFace.count == 4) {
var a = oldFace[0]
var b = oldFace[1]
var c = oldFace[2]
var d = oldFace[3]
var facePointAbcd = facePointNums[oldFaceNum]
var p = switchNums.call([a, b])
var edgePointAb = edgePointNums[Int.cantorPair(p[0], p[1])]
p = switchNums.call([d, a])
var edgePointDa = edgePointNums[Int.cantorPair(p[0], p[1])]
p = switchNums.call([b, c])
var edgePointBc = edgePointNums[Int.cantorPair(p[0], p[1])]
p = switchNums.call([c, d])
var edgePointCd = edgePointNums[Int.cantorPair(p[0], p[1])]
newFaces.add([a, edgePointAb, facePointAbcd, edgePointDa])
newFaces.add([b, edgePointBc, facePointAbcd, edgePointAb])
newFaces.add([c, edgePointCd, facePointAbcd, edgePointBc])
newFaces.add([d, edgePointDa, facePointAbcd, edgePointCd])
}
oldFaceNum = oldFaceNum + 1
}
return [newPoints, newFaces]
}
 
var inputPoints = [
Point.new(-1, 1, 1),
Point.new(-1, -1, 1),
Point.new( 1, -1, 1),
Point.new( 1, 1, 1),
Point.new( 1, -1, -1),
Point.new( 1, 1, -1),
Point.new(-1, -1, -1),
Point.new(-1, 1, -1)
]
 
var inputFaces = [
[0, 1, 2, 3],
[3, 2, 4, 5],
[5, 4, 6, 7],
[7, 0, 3, 5],
[7, 6, 1, 0],
[6, 1, 2, 4]
]
 
var outputPoints = inputPoints.toList
var outputFaces = inputFaces.toList
var iterations = 1
for (i in 0...iterations) {
var res = cmcSubdiv.call(outputPoints, outputFaces)
outputPoints = res[0]
outputFaces = res[1]
}
for (p in outputPoints) {
Fmt.aprint([p.x, p.y, p.z], 7, 4, "[]")
}
System.print()
for (f in outputFaces) {
Fmt.aprint(f, 2, 0, "[]")
}</syntaxhighlight>
 
{{out}}
<pre>
[-0.5556 0.5556 0.5556]
[-0.5556 -0.5556 0.5556]
[ 0.5556 -0.5556 0.5556]
[ 0.5556 0.5556 0.5556]
[ 0.5556 -0.5556 -0.5556]
[ 0.5556 0.5556 -0.5556]
[-0.5556 -0.5556 -0.5556]
[-0.5556 0.5556 -0.5556]
[ 0.0000 0.0000 1.0000]
[ 1.0000 0.0000 0.0000]
[ 0.0000 0.0000 -1.0000]
[ 0.0000 1.0000 0.0000]
[-1.0000 0.0000 0.0000]
[ 0.0000 -1.0000 0.0000]
[-0.7500 0.0000 0.7500]
[ 0.0000 0.7500 0.7500]
[-0.7500 0.7500 0.0000]
[ 0.0000 -0.7500 0.7500]
[-0.7500 -0.7500 0.0000]
[ 0.7500 0.0000 0.7500]
[ 0.7500 -0.7500 0.0000]
[ 0.7500 0.7500 0.0000]
[ 0.7500 0.0000 -0.7500]
[ 0.0000 -0.7500 -0.7500]
[ 0.0000 0.7500 -0.7500]
[-0.7500 0.0000 -0.7500]
 
[ 0 14 8 15]
[ 1 17 8 14]
[ 2 19 8 17]
[ 3 15 8 19]
[ 3 19 9 21]
[ 2 20 9 19]
[ 4 22 9 20]
[ 5 21 9 22]
[ 5 22 10 24]
[ 4 23 10 22]
[ 6 25 10 23]
[ 7 24 10 25]
[ 7 16 11 24]
[ 0 15 11 16]
[ 3 21 11 15]
[ 5 24 11 21]
[ 7 25 12 16]
[ 6 18 12 25]
[ 1 14 12 18]
[ 0 16 12 14]
[ 6 18 13 23]
[ 1 17 13 18]
[ 2 20 13 17]
[ 4 23 13 20]
</pre>
PureFox
9,476

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/136015...349994"