Matt Rajca

Modeling Parametric Surfaces with SceneKit

May 14, 2014

In a previous blog post, I wrote about algorithmically generating images with Core Image. Computer-generated art is a fascinating subject, and things get even more interesting when we enter the third dimension. In fact, mathematical equations can be used to parameterize 3D surfaces such as the one pictured below:

The seashell surface, shown above, can be parameterized by the following equations:

\[x(u,v) = \frac{5}{4}(1 - \frac{v}{2 \pi})\cos(2v)(1 + \cos(u)) + \cos(2v)\] \[y(u,v) = \frac{5}{4}(1 - \frac{v}{2 \pi})\sin(2v)(1 + \cos(u)) + \sin(2v)\] \[z(u,v) = \frac{10 v}{2 \pi} + \frac{5}{4}(1 - \frac{v}{2 \pi}) \sin(u) + 15\]

Note that $ x, y, z $ are functions of $ u $ and $ v $, where $ 0 \leq u \leq 2 \pi $ and $ -2 \pi \leq v \leq 2 \pi $. In other words, any point in 2D $ uv $-space maps to a 3D point in our model. If we were to generate these points at infinitesimally small intervals, we would get a “point cloud” that resembles the seashell rendered above. That’s all for the math; now let’s translate this to code.

In OS X “Mountain Lion”, Apple introduced SceneKit, a new Objective-C framework for rendering 3D graphics. While SceneKit is technically a wrapper over OpenGL, it provides an object-oriented (rather than procedural) model of programming. More specifically, a scene in SceneKit consists of a hierarchy of nodes to which we can attach geometry, lights, cameras, and so on. SceneKit’s node hierarchy will feel familiar to anyone who has worked with Core Animation layers. In fact, a lot of the concepts and terminology used in Core Animation carry over to SceneKit, providing a low learning curve once you get used to thinking in 3D space. I won’t focus on the basics of SceneKit here; Apple’s 2012 and 2013 WWDC sessions are good starting points.

Like any good framework, SceneKit makes common tasks easy and complex tasks possible. We’ll focus on the latter here. Going beyond SceneKit’s primitive geometries – cubes, cylinders, spheres – we’ll explore how we can programmatically construct the geometry for the seashell pictured above.

The three building blocks of geometry in SceneKit are vertices, vertex normals, and triangles. As you probably guessed, the $ x, y, z $ equations listed above give us our model’s vertices. Unfortunately, we can’t pass these equations along to SceneKit. Rather, we need to come up with a set of discrete vertices corresponding to points in our $ uv $-space. We can do this by subdividing our $ uv $-space into an evenly-spaced grid of points. As the number of subdivisions grows larger, the distances between the vertices gets smaller, ultimately resulting in smoother surfaces.

The code below does just this.

typedef struct {
	float x, y, z;
	float nx, ny, nz;
} Vertex;

#define SUBDIVISIONS 250

- (SCNGeometry *)newSeashell {
	// Allocate enough space for our vertices
	NSInteger vertexCount = (SUBDIVISIONS + 1) * (SUBDIVISIONS + 1);
	Vertex *vertices = malloc(sizeof(Vertex) * vertexCount);
	
	// Calculate the uv step interval given the number of subdivisions
	float uStep = 2.0f * M_PI / SUBDIVISIONS; // (2pi - 0) / subdivisions
	float vStep = 4.0f * M_PI / SUBDIVISIONS; // (2pi - -2pi) / subdivisions
	
	Vertex *curr = vertices;
	float u = 0.0f;
	
	// Loop through our uv-space, generating 3D vertices
	for (NSInteger i = 0; i <= SUBDIVISIONS; ++i, u += uStep) {
		float v = -2 * M_PI;
		
		for (NSInteger j = 0; j <= SUBDIVISIONS; ++j, v += vStep, ++curr) {
			curr->x = 5/4.0f * (1-v/(2*M_PI)) * cos(2*v) * (1 + cos(u)) + cos(2*v);
			curr->y = 5/4.0f * (1-v/(2*M_PI)) * sin(2*v) * (1 + cos(u)) + sin(2*v);
			curr->z = 5*v / M_PI + 5/4.0f * (1 - v/(2*M_PI)) * sin(u) + 15;
			
			// STEP 2: add normal equations here
		}
	}
	
	// STEP 3: generate indices
	// STEP 4: return geometry
}

First, we define a custom Vertex type that stores $ x, y, z $ coordinates and $ nx, ny, nz $ normal coordinates. Then we generate a 2D matrix of points in $ uv $-space. You can tweak the SUBDIVISIONS constant to reach a good balance between rendering performance and detail.

Next, we need to compute vertex normals for each of the vertices. Among other things, vertex normals ensure light is properly reflected off the model during rendering. We can compute a normal vector for a vertex by taking the cross product of the partial derivatives $ \frac{dr}{du} $ and $ \frac{dr}{dv} $ of the parameterization above. The resulting equations, obtained with Mathematica, are quite complex:

$$ nx = (-5 * (2 \pi - v0 * (2 * (20 + 18 \pi - 5v 0 * \cos(u - 2v) + 5 * (2 \pi - v0 * \cos(2 * (u - v0 ) + 20 * \pi * \cos(2*v) - 10v * \cos(2v) + 10 * \pi * \cos(2 * (u + v0 ) - 5v * \cos(2 * (u + v0 ) - 40 \cos(u + 2v) + 36 \pi * \cos(u + 2v) - 10v * \cos(u + 2v) + 5 \sin(u - 2v) - 10 \sin(2v) - 5 \sin(u + 2v) 0 0 / (128 * \pi^{2} 0 0 $$ $$ ny = (-5 * (2 \pi - v0 * (5 \cos(v)^{2} 0 * (1 + \cos(u) - 8\sin(u) 0 - 5 \sin(v) * (4v \cos(u)^{2} 0* \cos(v) + (1 + \cos(u) - 8\sin(u) 0 * \sin(v) 0 + 2 \cos(u) * (18 \pi - 5v + 10\pi * \cos(u) 0 * \sin(2v) 0 0 / (64 * \pi^{2} 0 0 $$ $$ nz = (-5 * (2 \pi - v0 * (18 \pi - 5v + 5*(2 \pi - v0 * \cos(u)0 * \sin(u)0 / (32 \pi^{2}0 $$

Add the following lines of code below the comment that says: STEP 2: add normal equations here:

curr->nx = (-5*(2*M_PI - v)*(2*(20 + 18*M_PI - 5*v)*cos(u - 2*v) + 5*(2*M_PI - v)*cos(2*(u - v)) + 20*M_PI*cos(2*v) - 10*v*cos(2*v) + 10*M_PI*cos(2*(u + v)) - 5*v*cos(2*(u + v)) - 40*cos(u + 2*v) + 36*M_PI*cos(u + 2*v) - 10*v*cos(u + 2*v) + 5*sin(u - 2*v) - 10*sin(2*v) - 5*sin(u + 2*v)))/(128.*pow(M_PI,2));

curr->ny = (-5*(2*M_PI - v)*(5*pow(cos(v),2)*(1 + cos(u) - 8*sin(u)) - 5*sin(v)*(4*v*pow(cos(u),2)*cos(v) + (1 + cos(u) - 8*sin(u))*sin(v)) + 2*cos(u)*(18*M_PI - 5*v + 10*M_PI*cos(u))*sin(2*v)))/(64.*pow(M_PI,2));

curr->nz = (-5*(2*M_PI - v)*(18*M_PI - 5*v + 5*(2*M_PI - v)*cos(u))*sin(u))/(32.*pow(M_PI,2));

// Normalize the results
float invLen = 1.0f / sqrtf(curr->nx * curr->nx + curr->ny * curr->ny + curr->nz * curr->nz);
curr->nx *= invLen;
curr->ny *= invLen;
curr->nz *= invLen;

Lastly, we need to combine triplets of vertices into triangles. We can do that by running through our vertices in squares and creating 2 triangles for each one.

// Generate indices
unsigned short *idx = indices;
unsigned short stripStart = 0;

for (NSInteger i = 0; i < SUBDIVISIONS; ++i, stripStart += (SUBDIVISIONS + 1)) {
    for (NSInteger j = 0; j < SUBDIVISIONS; ++j) {
		unsigned short v1 = stripStart + j;
		unsigned short v2 = stripStart + j + 1;
		unsigned short v3 = stripStart + (SUBDIVISIONS+1) + j;
		unsigned short v4 = stripStart + (SUBDIVISIONS+1) + j + 1;
		
		*idx++ = v4;
		*idx++ = v2;
		*idx++ = v3;
		*idx++ = v1;
		*idx++ = v3;
		*idx++ = v2;
    }
}

Finally, we create and return the geometry.

NSData *data = [NSData dataWithBytes:vertices length:vertexCount * sizeof(Vertex)];
free(vertices);

SCNGeometrySource *source = [SCNGeometrySource geometrySourceWithData:data
                                                             semantic:SCNGeometrySourceSemanticVertex
                                                          vectorCount:vertexCount
                                                      floatComponents:YES
                                                  componentsPerVector:3
                                                    bytesPerComponent:sizeof(float)
                                                           dataOffset:0
                                                           dataStride:sizeof(Vertex)];

SCNGeometrySource *normalSource = [SCNGeometrySource geometrySourceWithData:data
                                                                   semantic:SCNGeometrySourceSemanticNormal
                                                                vectorCount:vertexCount
                                                            floatComponents:YES
                                                        componentsPerVector:3
                                                          bytesPerComponent:sizeof(float)
                                                                 dataOffset:offsetof(Vertex, nx)
                                                                 dataStride:sizeof(Vertex)];

SCNGeometryElement *element = [SCNGeometryElement geometryElementWithData:[NSData dataWithBytes:indices length:indexCount * sizeof(unsigned short)]
                                                            primitiveType:SCNGeometryPrimitiveTypeTriangles
                                                           primitiveCount:indexCount/3
                                                            bytesPerIndex:sizeof(unsigned short)];

free(indices);

return [SCNGeometry geometryWithSources:@[source, normalSource] elements:@[element]];

You can find a sample program that renders the geometry on GitHub.