Catmull Clark Surface Subdivider in C#

I wanted to play with stuff like this; http://www.flickr.com/photos/junearch/sets/72157602857326020/detail/ which meant I needed a surface subdivider. Since it’s what I do for fun, I implemented the method described here: http://rosettacode.org/wiki/Catmull%E2%80%93Clark_subdivision_surface in C#. You’ll find the code below.

Undivided

Notice how the box has a hole – you can of course divide solid manifolds as well.

undivided

Once Divided

once divided

Twice divided

twice divided

Thrice divided

image

Four times divided

Dividing the box 4 times results in 10460 triangles. The divisions took 245ms– but I haven’t spent a lot of time optimizing this.

image

The code

using System;
using System.Collections.Generic;
using System.Linq;

namespace Subdivision.Core
{
    public class CatmullClarkSubdivider : ISubdivider
    {
        public Shape Subdivide(Shape shape)
        {
            // http://rosettacode.org/wiki/Catmull%E2%80%93Clark_subdivision_surface 
            Shape subdivided = new Shape();

            //for each face, a face point is created which is the average of all the points of the face.
            CreateFacePoints(shape);

            //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.
            CreateEdgePoints(shape);

            //for each vertex point, its coordinates are updated from (new_coords):
            //    the old coordinates (old_coords),
            //    the average of the face points of the faces the point belongs to (avg_face_points),
            //    the average of the centers of edges the point belongs to (avg_mid_edges),
            //    how many faces a point belongs to (n), then use this formula: 
            //m1 = (n - 3) / n
            //m2 = 1 / n
            //m3 = 2 / n
            //new_coords = (m1 * old_coords)
            //           + (m2 * avg_face_points)
            //           + (m3 * avg_mid_edges)
            CreateVertexPoints(shape);

            //for a triangle face (a,b,c): 
            //   (a, edge_pointab, face_pointabc, edge_pointca)
            //   (b, edge_pointbc, face_pointabc, edge_pointab)
            //   (c, edge_pointca, face_pointabc, edge_pointbc)

            //for a quad face (a,b,c,d): 
            //   (a, edge_pointab, face_pointabcd, edge_pointda)
            //   (b, edge_pointbc, face_pointabcd, edge_pointab)
            //   (c, edge_pointcd, face_pointabcd, edge_pointbc)
            //   (d, edge_pointda, face_pointabcd, edge_pointcd)
            CreateFaces(shape, subdivided);

            return subdivided;
        }

        private void CreateFacePoints(Shape shape)
        {
            //for each face, a face point is created which is the average of all the points of the face.
            foreach (Face face in shape.Faces)
            {
                List points = face.AllPoints;
                face.FacePoint = new Point(Average(points));
            }
        }

        private void CreateEdgePoints(Shape shape)
        {
            //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.
            List edges = shape.AllEdges;
            foreach (Edge edge in edges)
            {
                if (edge.IsOnBorderOfHole)
                {
                    Vector3 position =
                        Average(
                            edge.Points[0],
                            edge.Points[1]);

                    edge.EdgePoint = new Point(position);
                }
                else
                {
                    Vector3 position =
                        Average(
                            edge.Points[0],
                            edge.Points[1],
                            edge.Faces[0].FacePoint,
                            edge.Faces[1].FacePoint);

                    edge.EdgePoint = new Point(position);
                }
            }
        }

        private void CreateVertexPoints(Shape shape)
        {
            //for each vertex point, its coordinates are updated from (new_coords):
            //    the old coordinates (old_coords),
            //    the average of the centers of edges the point belongs to (avg_mid_edges),
            //    the average of the face points of the faces the point belongs to (avg_face_points),
            //    how many faces a point belongs to (n), then use this formula: 
            //m1 = (n - 3) / n
            //m2 = 1 / n
            //m3 = 2 / n
            //new_coords = (m1 * old_coords)
            //           + (m2 * avg_face_points)
            //           + (m3 * avg_mid_edges) 
            List allPoints = shape.AllPoints;
            List allEdges = shape.AllEdges;

            foreach (Point oldPoint in allPoints)
            {
                if (oldPoint.IsOnBorderOfHole)
                {
                    oldPoint.Successor = CreateVertexPointForBorderPoint(oldPoint);
                }
                else
                {
                    oldPoint.Successor = CreateVertexPoint(allEdges, oldPoint);
                }
            }
        }

        private Point CreateVertexPoint(List allEdges, Point oldPoint)
        {
            //    the average of the face points of the faces the point belongs to (avg_face_points),
            //    how many faces a point belongs to (n), then use this formula: 

            //    the average of the centers of edges the point belongs to (avg_mid_edges),
            Vector3 avgMidEdges = Vector3.Average(oldPoint.Edges.Select(e => e.Middle));

            //    the average of the face points of the faces the point belongs to (avg_face_points),
            List pointFaces = oldPoint.Edges.SelectMany(e => e.Faces).Distinct().ToList();
            Vector3 avgFacePoints = Average(pointFaces.Select(pf => pf.FacePoint));

            int faceCount = pointFaces.Count;

            double m1 = (faceCount - 3f) / faceCount;
            double m2 = 1f / faceCount;
            double m3 = 2f / faceCount;

            Point newPoint = new Point(m1 * oldPoint.Position + m2 * avgFacePoints + m3 * avgMidEdges);
            return newPoint;
        }

        private Point CreateVertexPointForBorderPoint(Point oldPoint)
        {
            //for the vertex points that are on the border of a hole, the new coordinates are calculated as follows:
            // in all the edges the point belongs to, only take in account the middles of the edges that are on the border of the hole
            // calculate the average between these points (on the hole boundary) and the old coordinates (also on the hole boundary). 
            List positions = oldPoint.Edges.Where(e => e.IsOnBorderOfHole).Select(e => e.Middle).ToList();
            positions.Add(oldPoint.Position);

            return new Point(Vector3.Average(positions));
        }

        private void CreateFaces(Shape shape, Shape subdivided)
        {
            List faces = shape.Faces;
            List existingEdges = new List();
            foreach (Face face in faces)
            {
                if (face.AllPoints.Count() == 3)
                {
                    CreateTriangleFace(existingEdges, subdivided, face);
                }
                else if (face.AllPoints.Count() == 4)
                {
                    CreateQuadFace(existingEdges, subdivided, face);
                }
                else
                {
                    throw new InvalidOperationException(string.Format("Unhandled facetype (point count={0})!", face.AllPoints.Count()));
                }
            }
            SubdivisionUtilities.VerifyThatThereAreNoEdgeDuplicates(existingEdges);
        }

        private void CreateTriangleFace(List existingEdges, Shape subdivided, Face face)
        {
            List points = face.AllPoints;
            Point a = points[0];
            Point b = points[1];
            Point c = points[2];

            //for a triangle face (a,b,c): 
            //   (a, edge_pointab, face_pointabc, edge_pointca)
            //   (b, edge_pointbc, face_pointabc, edge_pointab)
            //   (c, edge_pointca, face_pointabc, edge_pointbc)
            Point facePoint = face.FacePoint;

            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, a, face.Edges[0].EdgePoint, facePoint, face.Edges[2].EdgePoint));
            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, b, face.Edges[1].EdgePoint, facePoint, face.Edges[0].EdgePoint));
            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, c, face.Edges[2].EdgePoint, facePoint, face.Edges[1].EdgePoint));

            SubdivisionUtilities.VerifyThatThereAreNoEdgeDuplicates(existingEdges);
        }

        private void CreateQuadFace(List existingEdges, Shape subdivided, Face face)
        {
            //                  0 1 2 -> 3 
            //for a quad face (a,b,c,d): 
            //   (a, edge_pointab, face_pointabcd, edge_pointda)
            //   (b, edge_pointbc, face_pointabcd, edge_pointab)
            //   (c, edge_pointcd, face_pointabcd, edge_pointbc)
            //   (d, edge_pointda, face_pointabcd, edge_pointcd)
            List points = face.AllPoints;
            Point a = points[0].Successor;
            Point b = points[1].Successor;
            Point c = points[2].Successor;
            Point d = points[3].Successor;

            Point facePoint = face.FacePoint;

            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, a, face.Edges[0].EdgePoint, facePoint, face.Edges[3].EdgePoint));
            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, b, face.Edges[1].EdgePoint, facePoint, face.Edges[0].EdgePoint));
            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, c, face.Edges[2].EdgePoint, facePoint, face.Edges[1].EdgePoint));
            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, d, face.Edges[3].EdgePoint, facePoint, face.Edges[2].EdgePoint));

            SubdivisionUtilities.VerifyThatThereAreNoEdgeDuplicates(existingEdges);
        }

        private Vector3 Average(IEnumerable points)
        {
            return Vector3.Average(points.Select(p => p.Position));
        }

        private Vector3 Average(params Point[] points)
        {
            return Vector3.Average(points.Select(p => p.Position));
        }
    }
}

You can download the complete source code here (be warned, it’s 17MB).

About mfagerlund
Writes code in my sleep - and sometimes it even compiles!

2 Responses to Catmull Clark Surface Subdivider in C#

  1. Alex Torbin says:

    Hello! I’m working on a project where I plan to use the same algorithm for smoothing procedurally generated mesh. And I initially operate quads, since triangles give poor results when smoothing (distortion of surface shape). I want to have quads at the input and the output. So that further smoothing was unproblematic, and increased the number of triangles exactly 4 times. I was going to implement it myself, but would love to have a ready-made reference.Can you share the full code?

  2. mfagerlund says:

    Email me at mattias.fagerlund@carretera.se and I’ll email you the code.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: