Cycles render engine, initial commit. This is the engine itself, blender modifications and build instructions will follow later.

Cycles uses code from some great open source projects, many thanks them:

* BVH building and traversal code from NVidia's "Understanding the Efficiency of Ray Traversal on GPUs":
http://code.google.com/p/understanding-the-efficiency-of-ray-traversal-on-gpus/
* Open Shading Language for a large part of the shading system:
http://code.google.com/p/openshadinglanguage/
* Blender for procedural textures and a few other nodes.
* Approximate Catmull Clark subdivision from NVidia Mesh tools:
http://code.google.com/p/nvidia-mesh-tools/
* Sobol direction vectors from:
http://web.maths.unsw.edu.au/~fkuo/sobol/
* Film response functions from:
http://www.cs.columbia.edu/CAVE/software/softlib/dorf.php

18 files changed, 3356 insertions, 0 deletions

diff --git a/intern/cycles/subd/CMakeLists.txt b/intern/cycles/subd/CMakeLists.txt
new file mode 100644
index 00000000000..a14bf7896d8
--- /dev/null
+++ b/intern/cycles/subd/CMakeLists.txt
@@ -0,0 +1,26 @@
+
+INCLUDE_DIRECTORIES(. ../util ../kernel ../kernel/svm ../render)
+
+SET(sources
+	subd_build.cpp
+	subd_dice.cpp
+	subd_mesh.cpp
+	subd_patch.cpp
+	subd_ring.cpp
+	subd_split.cpp
+	subd_stencil.cpp)
+
+SET(headers
+	subd_build.h
+	subd_dice.h
+	subd_edge.h
+	subd_face.h
+	subd_mesh.h
+	subd_patch.h
+	subd_ring.h
+	subd_split.h
+	subd_stencil.h
+	subd_vert.h)
+
+ADD_LIBRARY(subd ${sources} ${headers})
+
diff --git a/intern/cycles/subd/subd_build.cpp b/intern/cycles/subd/subd_build.cpp
new file mode 100644
index 00000000000..4c2a6299014
--- /dev/null
+++ b/intern/cycles/subd/subd_build.cpp
@@ -0,0 +1,666 @@
+/*
+ * Copyright 2006, NVIDIA Corporation Ignacio Castano <icastano@nvidia.com>
+ *
+ * Modifications copyright (c) 2011, Blender Foundation. All rights reserved.
+ * 
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ * 
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ * 
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ */
+
+#include "subd_build.h"
+#include "subd_edge.h"
+#include "subd_face.h"
+#include "subd_ring.h"
+#include "subd_mesh.h"
+#include "subd_patch.h"
+#include "subd_stencil.h"
+#include "subd_vert.h"
+
+#include "util_algorithm.h"
+#include "util_debug.h"
+#include "util_math.h"
+#include "util_string.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* Subd Builder */
+
+SubdBuilder *SubdBuilder::create(bool linear)
+{
+	if(linear)
+		return new SubdLinearBuilder();
+	else
+		return new SubdAccBuilder();
+}
+
+/* Gregory ACC Stencil */
+
+class GregoryAccStencil {
+public:
+	SubdFaceRing *ring;
+	StencilMask stencil[20];
+
+	GregoryAccStencil(SubdFaceRing *ring_)
+	{
+		ring = ring_;
+
+		for(int i = 0; i < 20; i++)
+			stencil[i].resize(ring->num_verts());
+	}
+
+	StencilMask& get(int i)
+	{
+		assert(i < 20);
+		return stencil[i];
+	}
+
+	float& get(int i, SubdVert *vert)
+	{
+		assert(i < 20);
+		return stencil[i][ring->vert_index(vert)];
+	}
+};
+
+static float pseudoValence(SubdVert *vert)
+{
+	float valence = (float)vert->valence();
+
+	if(vert->is_boundary()) {
+		/* we treat boundary verts as being half a closed mesh. corners are
+		   special case. n = 4 for corners and n = 2*(n-1) for boundaries. */
+		if(valence == 2) return 4;
+		return (valence - 1)*2;
+	}
+
+	return valence;
+}
+
+/* Subd ACC Builder */
+
+SubdAccBuilder::SubdAccBuilder()
+{
+}
+
+SubdAccBuilder::~SubdAccBuilder()
+{
+}
+
+Patch *SubdAccBuilder::run(SubdFace *face)
+{
+	SubdFaceRing ring(face, face->edge);
+	GregoryAccStencil stencil(&ring);
+	float3 position[20];
+
+	computeCornerStencil(&ring, &stencil);
+	computeEdgeStencil(&ring, &stencil);
+	computeInteriorStencil(&ring, &stencil);
+
+	ring.evaluate_stencils(position, stencil.stencil, 20);
+
+	if(face->num_edges() == 3) {
+		GregoryTrianglePatch *patch = new GregoryTrianglePatch();
+		memcpy(patch->hull, position, sizeof(float3)*20);
+		return patch;
+	}
+	else if(face->num_edges() == 4)  {
+		GregoryQuadPatch *patch = new GregoryQuadPatch();
+		memcpy(patch->hull, position, sizeof(float3)*20);
+		return patch;
+	}
+
+	assert(0); /* n-gons should have been split already */
+	return NULL;
+}
+
+/* Gregory Patch */
+
+void SubdAccBuilder::computeCornerStencil(SubdFaceRing *ring, GregoryAccStencil *stencil)
+{
+	const int cornerIndices[7] = {8, 11, 19, 16,   6, 9, 12};
+	int primitiveOffset = ring->is_quad()? 0: 4;
+
+	SubdEdge *firstEdge = ring->firstEdge();
+
+	/* compute corner control points */
+	int v = 0;
+
+	for(SubdFace::EdgeIterator it(firstEdge); !it.isDone(); it.advance(), v++) {
+		SubdVert *vert = it.current()->from();
+		int valence = vert->valence();
+		int cid = cornerIndices[primitiveOffset+v];
+
+		if(vert->is_boundary()) {
+			/* compute vertex limit position */
+			SubdEdge *edge0 = vert->edge;
+			SubdEdge *edge1 = vert->edge->prev;
+
+			assert(edge0->face == NULL);
+			assert(edge0->to() != vert);
+			assert(edge1->face == NULL);
+			assert(edge1->from() != vert);
+
+			stencil->get(cid, vert) = 2.0f/3.0f;
+			stencil->get(cid, edge0->to()) = 1.0f/6.0f;
+			stencil->get(cid, edge1->from()) = 1.0f/6.0f;
+
+			assert(stencil->get(cid).is_normalized());
+		}
+		else {
+			stencil->get(cid, vert) = 3.0f*valence*valence;
+
+			for(SubdVert::EdgeIterator eit(vert->edge); !eit.isDone(); eit.advance()) {
+				SubdEdge *edge = eit.current();
+				assert(vert->co == edge->from()->co);
+
+				stencil->get(cid, edge->to()) = 12.0f;
+
+				if(SubdFaceRing::is_triangle(edge->face)) {
+					/* distribute weight to all verts */
+					stencil->get(cid, vert) += 1.0f;
+					stencil->get(cid, edge->to()) += 1.0f;
+					stencil->get(cid, edge->next->to()) += 1.0f;
+				}
+				else
+					stencil->get(cid, edge->next->to()) = 3.0f;
+			}
+
+			/* normalize stencil. */
+			stencil->get(cid).normalize();
+		}
+	}
+}
+
+void SubdAccBuilder::computeEdgeStencil(SubdFaceRing *ring, GregoryAccStencil *stencil)
+{
+	const int cornerIndices[7] = {8, 11, 19, 16,    6,  9, 12};
+	const int edge1Indices[7] = {9, 13, 18, 14,    7, 10, 13};
+	const int edge2Indices[7] = {12, 10, 15, 17,    14,  8, 11};
+	int primitiveOffset = ring->is_quad()? 0: 4;
+
+	float tangentScales[14] = {
+		0.0f, 0.0f, 0.0f, 0.667791f, 1.0f,
+		1.11268f, 1.1284f, 1.10289f, 1.06062f,
+		1.01262f, 0.963949f, 0.916926f, 0.872541f, 0.831134f
+	};
+
+	SubdEdge *firstEdge = ring->firstEdge();
+
+	/* compute corner / edge control points */
+	int v = 0;
+
+	for(SubdFace::EdgeIterator it(firstEdge); !it.isDone(); it.advance(), v++) {
+		SubdVert *vert = it.current()->from();
+		int valence = vert->valence();
+		int cid = cornerIndices[primitiveOffset+v];
+
+		int i1 = 0, i2 = 0, j = 0;
+
+		for(SubdVert::EdgeIterator eit(vert->edge); !eit.isDone(); eit.advance(), j++) {
+			SubdEdge *edge = eit.current();
+
+			/* find index of "our" edge for edge control points */
+			if(edge == it.current())
+				i1 = j;
+			if(edge == it.current()->prev->pair)
+				i2 = j;
+		}
+
+		if(vert->is_boundary()) {
+			int num_verts = ring->num_verts();
+			StencilMask	r0(num_verts);
+			StencilMask	r1(num_verts);
+
+			computeBoundaryTangentStencils(ring, vert, r0, r1);
+
+			int k = valence - 1;
+			float omega = M_PI / k;
+
+			int eid1 = edge1Indices[primitiveOffset + v];
+			int eid2 = edge2Indices[primitiveOffset + v];
+
+			if(it.current()->is_boundary()) {
+				assert(it.current()->from() == vert);
+
+				stencil->get(eid1, vert) = 2.0f / 3.0f;
+				stencil->get(eid1, it.current()->to()) = 1.0f / 3.0f;
+
+				assert(stencil->get(eid1).is_normalized());
+
+				if(valence == 2) {
+					for(int i = 0; i < num_verts; i++)
+						stencil->get(eid1)[i] += r0[i] * 0.0001f;
+				}
+			}
+			else {
+				stencil->get(eid1) = stencil->get(cid);
+
+				/* compute index of it.current() around vert */
+				int idx = 0;
+
+				for(SubdVert::EdgeIterator eit(vert->edges()); !eit.isDone(); eit.advance(), idx++)
+					if(eit.current() == it.current())
+						break;
+
+				assert(idx != valence);
+
+				float c = cosf(idx * omega);
+				float s = sinf(idx * omega);
+
+				for(int i = 0; i < num_verts; i++)
+					stencil->get(eid1)[i] += (r0[i] * s + r1[i] * c) / 3.0f;
+			}
+
+			if(it.current()->prev->is_boundary()) {
+				assert(it.current()->prev->pair->from() == vert);
+
+				stencil->get(eid2, vert) = 2.0f / 3.0f;
+				stencil->get(eid2, it.current()->prev->pair->to()) = 1.0f / 3.0f;
+
+				assert(stencil->get(eid2).is_normalized());
+
+				if(valence == 2) {
+					for(int i = 0; i < num_verts; i++)
+						stencil->get(eid2)[i] += r0[i] * 0.0001f;
+				}
+			}
+			else {
+				stencil->get(eid2) = stencil->get(cid);
+
+				/* compute index of it.current() around vert */
+				int idx = 0;
+
+				for(SubdVert::EdgeIterator eit(vert->edges()); !eit.isDone(); eit.advance(), idx++)
+					if(eit.current() == it.current()->prev->pair)
+						break;
+
+				assert(idx != valence);
+
+				float c = cosf(idx * omega);
+				float s = sinf(idx * omega);
+
+				for(int i = 0; i < num_verts; i++)
+					stencil->get(eid2)[i] += (r0[i] * s + r1[i] * c) / 3;
+			}
+		}
+		else {
+			float costerm = cosf(M_PI / valence);
+			float sqrtterm = sqrtf(4.0f + costerm*costerm);
+
+			/* float tangentScale = 1.0f; */
+			float tangentScale = tangentScales[min(valence, 13U)];
+
+			float alpha = (1.0f +  costerm / sqrtterm) / (3.0f * valence) * tangentScale;
+			float beta  = 1.0f / (3.0f * valence * sqrtterm) * tangentScale;
+
+
+			int eid1 = edge1Indices[primitiveOffset + v];
+			int eid2 = edge2Indices[primitiveOffset + v];
+
+			stencil->get(eid1) = stencil->get(cid);
+			stencil->get(eid2) = stencil->get(cid);
+
+			int j = 0;
+			for(SubdVert::EdgeIterator eit(vert->edges()); !eit.isDone(); eit.advance(), j++) {
+				SubdEdge *edge = eit.current();
+				assert(vert->co == edge->from()->co);
+
+				float costerm1_a = cosf(M_PI * 2 * (j-i1) / valence);
+				float costerm1_b = cosf(M_PI * (2 * (j-i1)-1) / valence); /* -1 instead of +1 b/c of edge->next->to() */
+
+				float costerm2_a = cosf(M_PI * 2 * (j-i2) / valence);
+				float costerm2_b = cosf(M_PI * (2 * (j-i2)-1) / valence); /* -1 instead of +1 b/c of edge->next->to() */
+
+
+				stencil->get(eid1, edge->to()) += alpha * costerm1_a;
+				stencil->get(eid2, edge->to()) += alpha * costerm2_a;
+
+				if(SubdFaceRing::is_triangle(edge->face)) {
+					/* @@ this probably does not provide watertight results!! (1/3 + 1/3 + 1/3 != 1) */
+
+					/* distribute weight to all verts */
+					stencil->get(eid1, vert) += beta * costerm1_b / 3.0f;				
+					stencil->get(eid1, edge->to()) += beta * costerm1_b / 3.0f;
+					stencil->get(eid1, edge->next->to()) += beta * costerm1_b / 3.0f;
+
+					stencil->get(eid2, vert) += beta * costerm2_b / 3.0f;
+					stencil->get(eid2, edge->to()) += beta * costerm2_b / 3.0f;
+					stencil->get(eid2, edge->next->to()) += beta * costerm2_b / 3.0f;
+				}
+				else {
+					stencil->get(eid1, edge->next->to()) += beta * costerm1_b;
+					stencil->get(eid2, edge->next->to()) += beta * costerm2_b;
+				}
+			}
+		}
+	}
+}
+
+void SubdAccBuilder::computeInteriorStencil(SubdFaceRing *ring, GregoryAccStencil *stencil)
+{
+	static int corner1Indices[7] = {8, 11, 19, 16,    6, 9, 12};
+	static int corner2Indices[7] = {11, 19, 16, 8,    9, 12, 6};
+	static int edge1Indices[7] = {9, 13, 18, 14,    7, 10, 13};
+	static int edge2Indices[7] = {10, 15, 17, 12,    8, 11, 14};
+	static int interior1Indices[7] = {1, 3, 6, 4,    1, 3, 5};
+	static int interior2Indices[7] = {2, 7, 5, 0,    2, 4, 0};
+
+	int primitiveOffset = ring->is_quad()? 0: 4;
+
+	SubdFace * face = ring->face();
+	SubdEdge *firstEdge = ring->firstEdge();
+
+	/* interior control points */
+	int v = 0;
+	for(SubdFace::EdgeIterator it(firstEdge); !it.isDone(); it.advance(), v++) {
+		SubdEdge *edge = it.current();
+
+		if(edge->is_boundary()) {
+			float valence1 = pseudoValence(edge->from());
+			float valence2 = pseudoValence(edge->to());
+
+			float weights1[4];
+			float weights2[4];
+
+			if(ring->is_quad()) {
+				weights1[0] = 3 * valence1;
+				weights1[1] = 6;
+				weights1[2] = 3;
+				weights1[3] = 6;
+
+				weights2[0] = 6;
+				weights2[1] = 3 * valence2;
+				weights2[2] = 6;
+				weights2[3] = 3;
+			}
+			else {
+				assert(ring->is_triangle());
+				weights1[0] = 3 * valence1 + 1;
+				weights1[1] = 7;
+				weights1[2] = 7;
+
+				weights2[0] = 7;
+				weights2[1] = 3 * valence2 + 1;
+				weights2[2] = 7;
+			}
+
+			int idx1 = interior1Indices[primitiveOffset+v];
+			int idx2 = interior2Indices[primitiveOffset+v];
+
+			int i = 0;
+			for(SubdFace::EdgeIterator it(face->edges(edge)); !it.isDone(); it.advance(), i++) {
+				SubdVert *vert = it.current()->from();
+				stencil->get(idx1, vert) += weights1[i];
+				stencil->get(idx2, vert) += weights2[i];
+			}
+
+			stencil->get(idx1).normalize();
+			stencil->get(idx2).normalize();
+		}
+		else {
+			SubdVert *e0 = edge->from();
+			float costerm0 = cosf(2.0f * M_PI / pseudoValence(e0));
+
+			SubdVert *f0 = edge->to();
+			float costerm1 = cosf(2.0f * M_PI / pseudoValence(f0));
+
+			/*  p0 +------+ q0
+			 *	 |	  |
+			 *  f0 +======+ e0 <=== current edge
+			 *	 |	  |
+			 *  p1 +------+ q1
+			 */
+
+			SubdVert *q0 = edge->next->to();
+			SubdVert *p0 = edge->prev->from();
+
+			SubdVert *p1 = edge->pair->next->to();
+			SubdVert *q1 = edge->pair->prev->from();
+
+
+			StencilMask	x(ring->num_verts());
+			StencilMask	y(ring->num_verts());
+
+			for(int i = 0; i < ring->num_verts(); i++) {
+				x[i] =
+					(costerm1 * stencil->get(corner1Indices[primitiveOffset+v])[i] -
+					(2*costerm0 + costerm1) * stencil->get(edge1Indices[primitiveOffset+v])[i] +
+					2*costerm0 * stencil->get(edge2Indices[primitiveOffset+v])[i]) / 3.0f;
+			}
+
+			/* y = (2*( midedgeA1 - midedgeB1) + 4*(centroidA - centroidB))/18.0f; */
+			y[ring->vert_index(p0)] = 1;
+			y[ring->vert_index(p1)] = -1;
+
+			/* add centroidA */
+			if(ring->is_triangle()) {
+				y[ring->vert_index(p0)] += 4.0f / 3.0f;
+				y[ring->vert_index(e0)] += 4.0f / 3.0f;
+				y[ring->vert_index(f0)] += 4.0f / 3.0f;
+			}
+			else {
+				y[ring->vert_index(p0)] += 1;
+				y[ring->vert_index(q0)] += 1;
+				y[ring->vert_index(e0)] += 1;
+				y[ring->vert_index(f0)] += 1;
+			}
+
+			/* sub centroidB */
+			if(SubdFaceRing::is_triangle(edge->pair->face)) {
+				y[ring->vert_index(p1)] -= 4.0f / 3.0f;
+				y[ring->vert_index(e0)] -= 4.0f / 3.0f;
+				y[ring->vert_index(f0)] -= 4.0f / 3.0f;
+
+			}
+			else {
+				y[ring->vert_index(p1)] -= 1;
+				y[ring->vert_index(q1)] -= 1;
+				y[ring->vert_index(e0)] -= 1;
+				y[ring->vert_index(f0)] -= 1;
+			}
+
+			y /= 18.0f;
+
+			if(ring->is_triangle()) {
+				x *= 3.0f / 4.0f;
+				y *= 3.0f / 4.0f;
+			}
+
+			/* this change makes the triangle boundaries smoother, but distorts the quads next to them */
+			/*if(ring->is_triangle() || SubdFaceRing::is_triangle(edge->pair->face))
+			{
+				y *= 4.0f / 3.0f;
+			}*/
+
+			stencil->get(interior1Indices[primitiveOffset+v]) = stencil->get(edge1Indices[primitiveOffset+v]);
+			stencil->get(interior1Indices[primitiveOffset+v]) += x;
+			stencil->get(interior1Indices[primitiveOffset+v]) += y;
+
+			for(int i = 0; i < ring->num_verts(); i++) {
+				x[i] =
+					(costerm0 * stencil->get(corner2Indices[primitiveOffset+v])[i] -
+					(2*costerm1 + costerm0) * stencil->get(edge2Indices[primitiveOffset+v])[i] + 
+					2*costerm1 * stencil->get(edge1Indices[primitiveOffset+v])[i]) / 3.0f;
+			}
+
+			/* y = (2*( midedgeA2 - midedgeB2) + 4*(centroidA - centroidB))/18.0f; */
+			y = 0.0f;
+
+			/* (2*( midedgeA2 - midedgeB2) */
+			y[ring->vert_index(q0)] = 1;
+			y[ring->vert_index(q1)] = -1;
+
+			/* add centroidA */
+			if(ring->is_triangle()) {
+				y[ring->vert_index(p0)] += 4.0f / 3.0f;
+				y[ring->vert_index(e0)] += 4.0f / 3.0f;
+				y[ring->vert_index(f0)] += 4.0f / 3.0f;
+			}
+			else {
+				y[ring->vert_index(p0)] += 1;
+				y[ring->vert_index(q0)] += 1;
+				y[ring->vert_index(e0)] += 1;
+				y[ring->vert_index(f0)] += 1;
+			}
+
+			/* sub centroidB */
+			if(SubdFaceRing::is_triangle(edge->pair->face)) {
+				y[ring->vert_index(p1)] -= 4.0f / 3.0f;
+				y[ring->vert_index(e0)] -= 4.0f / 3.0f;
+				y[ring->vert_index(f0)] -= 4.0f / 3.0f;
+
+			}
+			else {
+				y[ring->vert_index(p1)] -= 1;
+				y[ring->vert_index(q1)] -= 1;
+				y[ring->vert_index(e0)] -= 1;
+				y[ring->vert_index(f0)] -= 1;
+			}
+
+			y /= 18.0f;
+
+			if(ring->is_triangle()) {
+				x *= 3.0f / 4.0f;
+				y *= 3.0f / 4.0f;
+			}
+
+			/* this change makes the triangle boundaries smoother, but distorts the quads next to them. */
+			/*if(ring->is_triangle() || SubdFaceRing::is_triangle(edge->pair->face))
+				y *= 4.0f / 3.0f;*/
+
+			stencil->get(interior2Indices[primitiveOffset+v]) = stencil->get(edge2Indices[primitiveOffset+v]);
+			stencil->get(interior2Indices[primitiveOffset+v]) += x;
+			stencil->get(interior2Indices[primitiveOffset+v]) += y;
+		}
+	}
+}
+
+void SubdAccBuilder::computeBoundaryTangentStencils(SubdFaceRing *ring, SubdVert *vert, StencilMask & r0, StencilMask & r1)
+{
+	assert(vert->is_boundary());
+	assert(r0.size() == ring->num_verts());
+	assert(r1.size() == ring->num_verts());
+
+	SubdEdge *edge0 = vert->edge;
+	assert(edge0->face == NULL);
+	assert(edge0->to() != vert);
+
+	SubdEdge *edgek = vert->edge->prev;
+	assert(edgek->face == NULL);
+	assert(edgek->from() != vert);
+
+	int valence = vert->valence();
+
+	int k = valence - 1;
+	float omega = M_PI / k;
+	float s = sinf(omega);
+	float c = cosf(omega);
+
+	float factor = 1.0f / (3 * k + c);
+
+	float gamma = -4 * s * factor;
+	r0[ring->vert_index(vert)] = gamma;
+	/* r1[ring->vert_index(vert)] = 0; */
+
+	float salpha0 = -((1 + 2 * c) * sqrtf(1 + c)) * factor / sqrtf(1 - c);
+	float calpha0 = 1.0f / 2.0f;
+
+	r0[ring->vert_index(edge0->to())] = salpha0;
+	r1[ring->vert_index(edge0->to())] = calpha0;
+
+	float salphak = salpha0;
+	float calphak = -1.0f / 2.0f;
+
+	r0[ring->vert_index(edgek->from())] = salphak;
+	r1[ring->vert_index(edgek->from())] = calphak;
+
+	int j = 0;
+	for(SubdVert::EdgeIterator it(vert->edges()); !it.isDone(); it.advance(), j++) {
+		SubdEdge *edge = it.current();
+
+		if(j == k) break;
+
+		SubdVert *p = edge->to();
+		SubdVert *q = edge->pair->prev->from();
+
+		float alphaj = 4 * sinf(j * omega) * factor;
+		float betaj = (sinf(j * omega) + sinf((j + 1) * omega)) * factor;
+
+		if(j != 0)
+			r0[ring->vert_index(p)] += alphaj;
+
+		if(edge->pair->prev->prev->prev == edge->pair) {
+			r0[ring->vert_index(vert)] += betaj / 3.0f;
+			r0[ring->vert_index(edge->pair->from())] += betaj / 3.0f;
+			r0[ring->vert_index(q)] += betaj / 3.0f;
+		}
+		else
+			r0[ring->vert_index(q)] += betaj;
+	}
+
+	if(valence == 2) {
+		/* r0 perpendicular to r1 */
+		r0[ring->vert_index(vert)] = -4.0f / 3.0f;
+		r0[ring->vert_index(edgek->from())] = 1.0f / 2.0f;
+		r0[ring->vert_index(edge0->to())] = 1.0f / 2.0f;
+		r0[ring->vert_index(edge0->next->to())] = 1.0f / 3.0f;
+	}
+}
+
+/* Subd Linear Builder */
+
+SubdLinearBuilder::SubdLinearBuilder()
+{
+}
+
+SubdLinearBuilder::~SubdLinearBuilder()
+{
+}
+
+Patch *SubdLinearBuilder::run(SubdFace *face)
+{
+	Patch *patch;
+	float3 *hull;
+
+	if(face->num_edges() == 3) {
+		LinearTrianglePatch *lpatch = new LinearTrianglePatch();
+		hull = lpatch->hull;
+		patch = lpatch;
+	}
+	else if(face->num_edges() == 4)  {
+		LinearQuadPatch *lpatch = new LinearQuadPatch();
+		hull = lpatch->hull;
+		patch = lpatch;
+	}
+	else {
+		assert(0); /* n-gons should have been split already */
+		return NULL;
+	}
+
+	int i = 0;
+
+	for(SubdFace::EdgeIterator it(face->edge); !it.isDone(); it.advance())
+		hull[i++] = it.current()->from()->co;
+
+	if(face->num_edges() == 4)
+		swap(hull[2], hull[3]);
+
+	return patch;
+}
+
+CCL_NAMESPACE_END
+
diff --git a/intern/cycles/subd/subd_build.h b/intern/cycles/subd/subd_build.h
new file mode 100644
index 00000000000..e93ffc4a2c4
--- /dev/null
+++ b/intern/cycles/subd/subd_build.h
@@ -0,0 +1,75 @@
+/*
+ * Copyright 2011, Blender Foundation.
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#ifndef __SUBD_BUILD_H__
+#define __SUBD_BUILD_H__
+
+CCL_NAMESPACE_BEGIN
+
+class SubdFace;
+class SubdFaceRing;
+class SubdVert;
+class GregoryAccStencil;
+class Patch;
+class StencilMask;
+
+/* Builder */
+
+class SubdBuilder
+{
+public:
+	virtual ~SubdBuilder() {};
+	virtual Patch *run(SubdFace *face) = 0;
+	static SubdBuilder *create(bool linear);
+};
+
+/* Approximate Catmull Clark using Loop's approximation */
+
+class SubdAccBuilder : public SubdBuilder
+{
+public:
+	SubdAccBuilder();
+	~SubdAccBuilder();
+
+	Patch *run(SubdFace *face);
+
+protected:
+	/* Gregory Patch */
+	void computeCornerStencil(SubdFaceRing *ring, GregoryAccStencil *stencil);
+	void computeEdgeStencil(SubdFaceRing *ring, GregoryAccStencil *stencil);
+	void computeInteriorStencil(SubdFaceRing *ring, GregoryAccStencil *stencil);
+	void computeBoundaryTangentStencils(SubdFaceRing *ring,
+		SubdVert *vert,
+		StencilMask & r0, StencilMask & r1);
+};
+
+/* Linear Subdivision */
+
+class SubdLinearBuilder : public SubdBuilder
+{
+public:
+	SubdLinearBuilder();
+	~SubdLinearBuilder();
+
+	Patch *run(SubdFace *face);
+};
+
+CCL_NAMESPACE_END
+
+#endif /* __SUBD_BUILD_H__ */
+
diff --git a/intern/cycles/subd/subd_dice.cpp b/intern/cycles/subd/subd_dice.cpp
new file mode 100644
index 00000000000..086b7b246d3
--- /dev/null
+++ b/intern/cycles/subd/subd_dice.cpp
@@ -0,0 +1,461 @@
+/*
+ * Copyright 2011, Blender Foundation.
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include "camera.h"
+#include "mesh.h"
+
+#include "subd_dice.h"
+#include "subd_patch.h"
+
+#include "util_debug.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* EdgeDice Base */
+
+EdgeDice::EdgeDice(Mesh *mesh_, int shader_, bool smooth_, float dicing_rate_)
+{
+	mesh = mesh_;
+	mesh_P = NULL;
+	mesh_N = NULL;
+	vert_offset = 0;
+	dicing_rate = dicing_rate_;
+	shader = shader_;
+	smooth = smooth_;
+	camera = NULL;
+
+	mesh->attributes.add(Attribute::STD_VERTEX_NORMAL);
+}
+
+void EdgeDice::reserve(int num_verts, int num_tris)
+{
+	vert_offset = mesh->verts.size();
+	tri_offset = mesh->triangles.size();
+
+	mesh->reserve(vert_offset + num_verts, tri_offset + num_tris);
+
+	Attribute *attr_vN = mesh->attributes.add(Attribute::STD_VERTEX_NORMAL);
+
+	mesh_P = &mesh->verts[0];
+	mesh_N = attr_vN->data_float3();
+}
+
+int EdgeDice::add_vert(Patch *patch, float2 uv)
+{
+	float3 P, N, dPdu, dPdv;
+
+	patch->eval(&P, &dPdu, &dPdv, uv.x, uv.y);
+	N = normalize(cross(dPdu, dPdv));
+
+	assert(vert_offset < mesh->verts.size());
+
+	mesh_P[vert_offset] = P;
+	mesh_N[vert_offset] = N;
+
+	return vert_offset++;
+}
+
+void EdgeDice::add_triangle(int v0, int v1, int v2)
+{
+	mesh->add_triangle(v0, v1, v2, shader, smooth);
+}
+
+void EdgeDice::stitch_triangles(vector<int>& outer, vector<int>& inner)
+{
+	if(inner.size() == 0 || outer.size() == 0)
+		return; // XXX avoid crashes for Mu or Mv == 1, missing polygons
+
+	/* stitch together two arrays of verts with triangles. at each step,
+	   we compare using the next verts on both sides, to find the split
+	   direction with the smallest diagonal, and use that in order to keep
+	   the triangle shape reasonable. */
+	for(size_t i = 0, j = 0; i+1 < inner.size() || j+1 < outer.size();) {
+		int v0, v1, v2;
+
+		v0 = inner[i];
+		v1 = outer[j];
+
+		if(j+1 == outer.size()) {
+			v2 = inner[++i];
+		}
+		else if(i+1 == inner.size()) {
+			v2 = outer[++j];
+		}
+		else {
+			/* length of diagonals */
+			float len1 = len(mesh_P[inner[i]] - mesh_P[outer[j+1]]);
+			float len2 = len(mesh_P[outer[j]] - mesh_P[inner[i+1]]);
+
+			/* use smallest diagonal */
+			if(len1 < len2)
+				v2 = outer[++j];
+			else
+				v2 = inner[++i];
+		}
+
+		add_triangle(v0, v1, v2);
+	}
+}
+
+/* QuadDice */
+
+QuadDice::QuadDice(Mesh *mesh_, int shader_, bool smooth_, float dicing_rate_)
+: EdgeDice(mesh_, shader_, smooth_, dicing_rate_)
+{
+}
+
+void QuadDice::reserve(EdgeFactors& ef, int Mu, int Mv)
+{
+	/* XXX need to make this also work for edge factor 0 and 1 */
+	int num_verts = (ef.tu0 + ef.tu1 + ef.tv0 + ef.tv1) + (Mu - 1)*(Mv - 1);
+	EdgeDice::reserve(num_verts, 0);
+}
+
+float2 QuadDice::map_uv(SubPatch& sub, float u, float v)
+{
+	/* map UV from subpatch to patch parametric coordinates */
+	float2 d0 = interp(sub.P00, sub.P01, v);
+	float2 d1 = interp(sub.P10, sub.P11, v);
+	return interp(d0, d1, u);
+}
+
+float3 QuadDice::eval_projected(SubPatch& sub, float u, float v)
+{
+	float2 uv = map_uv(sub, u, v);
+	float3 P;
+
+	sub.patch->eval(&P, NULL, NULL, uv.x, uv.y);
+	if(camera)
+		P = transform(&camera->worldtoraster, P);
+
+	return P;
+}
+
+int QuadDice::add_vert(SubPatch& sub, float u, float v)
+{
+	return EdgeDice::add_vert(sub.patch, map_uv(sub, u, v));
+}
+
+void QuadDice::add_side_u(SubPatch& sub,
+	vector<int>& outer, vector<int>& inner,
+	int Mu, int Mv, int tu, int side, int offset)
+{
+	outer.clear();
+	inner.clear();
+
+	/* set verts on the edge of the patch */
+	outer.push_back(offset + ((side)? 2: 0));
+
+	for(int i = 1; i < tu; i++) {
+		float u = i/(float)tu;
+		float v = (side)? 1.0f: 0.0f;
+
+		outer.push_back(add_vert(sub, u, v));
+	}
+
+	outer.push_back(offset + ((side)? 3: 1));
+
+	/* set verts on the edge of the inner grid */
+	for(int i = 0; i < Mu-1; i++) {
+		int j = (side)? Mv-1-1: 0;
+		inner.push_back(offset + 4 + i + j*(Mu-1));
+	}
+}
+
+void QuadDice::add_side_v(SubPatch& sub,
+	vector<int>& outer, vector<int>& inner,
+	int Mu, int Mv, int tv, int side, int offset)
+{
+	outer.clear();
+	inner.clear();
+
+	/* set verts on the edge of the patch */
+	outer.push_back(offset + ((side)? 1: 0));
+
+	for(int j = 1; j < tv; j++) {
+		float u = (side)? 1.0f: 0.0f;
+		float v = j/(float)tv;
+
+		outer.push_back(add_vert(sub, u, v));
+	}
+
+	outer.push_back(offset + ((side)? 3: 2));
+
+	/* set verts on the edge of the inner grid */
+	for(int j = 0; j < Mv-1; j++) {
+		int i = (side)? Mu-1-1: 0;
+		inner.push_back(offset + 4 + i + j*(Mu-1));
+	}
+}
+
+float QuadDice::quad_area(const float3& a, const float3& b, const float3& c, const float3& d)
+{
+	return triangle_area(a, b, d) + triangle_area(a, d, c);
+}
+
+float QuadDice::scale_factor(SubPatch& sub, EdgeFactors& ef, int Mu, int Mv)
+{
+	/* estimate area as 4x largest of 4 quads */
+	float3 P[3][3];
+
+	for(int i = 0; i < 3; i++)
+		for(int j = 0; j < 3; j++)
+			P[i][j] = eval_projected(sub, i*0.5f, j*0.5f);
+
+	float A1 = quad_area(P[0][0], P[1][0], P[0][1], P[1][1]);
+	float A2 = quad_area(P[1][0], P[2][0], P[1][1], P[2][1]);
+	float A3 = quad_area(P[0][1], P[1][1], P[0][2], P[1][2]);
+	float A4 = quad_area(P[1][1], P[2][1], P[1][2], P[2][2]);
+	float Apatch = max(A1, max(A2, max(A3, A4)))*4.0f;
+
+	/* solve for scaling factor */
+	float Atri = dicing_rate*dicing_rate*0.5f;
+	float Ntris = Apatch/Atri;
+
+	// XXX does the -sqrt solution matter
+	// XXX max(D, 0.0) is highly suspicious, need to test cases
+	// where D goes negative
+	float N = 0.5f*(Ntris - (ef.tu0 + ef.tu1 + ef.tv0 + ef.tv1));
+	float D = 4.0f*N*Mu*Mv + (Mu + Mv)*(Mu + Mv);
+	float S = (Mu + Mv + sqrtf(max(D, 0.0f)))/(2*Mu*Mv);
+
+	return S;
+}
+
+void QuadDice::add_corners(SubPatch& sub)
+{
+	/* add verts for patch corners */
+	if(sub.patch->is_triangle()) {
+		add_vert(sub, 0.0f, 0.0f);
+		add_vert(sub, 1.0f, 0.0f);
+		add_vert(sub, 0.0f, 1.0f);
+	}
+	else {
+		add_vert(sub, 0.0f, 0.0f);
+		add_vert(sub, 1.0f, 0.0f);
+		add_vert(sub, 0.0f, 1.0f);
+		add_vert(sub, 1.0f, 1.0f);
+	}
+}
+
+void QuadDice::add_grid(SubPatch& sub, int Mu, int Mv, int offset)
+{
+	/* create inner grid */
+	float du = 1.0f/(float)Mu;
+	float dv = 1.0f/(float)Mv;
+
+	for(int j = 1; j < Mv; j++) {
+		for(int i = 1; i < Mu; i++) {
+			float u = i*du;
+			float v = j*dv;
+
+			add_vert(sub, u, v);
+
+			if(i < Mu-1 && j < Mv-1) {
+				int i1 = offset + 4 + (i-1) + (j-1)*(Mu-1);
+				int i2 = offset + 4 + i + (j-1)*(Mu-1);
+				int i3 = offset + 4 + i + j*(Mu-1);
+				int i4 = offset + 4 + (i-1) + j*(Mu-1);
+
+				add_triangle(i1, i2, i3);
+				add_triangle(i1, i3, i4);
+			}
+		}
+	}
+}
+
+void QuadDice::dice(SubPatch& sub, EdgeFactors& ef)
+{
+	/* compute inner grid size with scale factor */
+	int Mu = max(ef.tu0, ef.tu1);
+	int Mv = max(ef.tv0, ef.tv1);
+
+	float S = scale_factor(sub, ef, Mu, Mv);
+	Mu = max((int)ceil(S*Mu), 2); // XXX handle 0 & 1?
+	Mv = max((int)ceil(S*Mv), 2); // XXX handle 0 & 1?
+
+	/* reserve space for new verts */
+	int offset = mesh->verts.size();
+	reserve(ef, Mu, Mv);
+
+	/* corners and inner grid */
+	add_corners(sub);
+	add_grid(sub, Mu, Mv, offset);
+
+	/* bottom side */
+	vector<int> outer, inner;
+
+	add_side_u(sub, outer, inner, Mu, Mv, ef.tu0, 0, offset);
+	stitch_triangles(outer, inner);
+
+	/* top side */
+	add_side_u(sub, outer, inner, Mu, Mv, ef.tu1, 1, offset);
+	stitch_triangles(inner, outer);
+
+	/* left side */
+	add_side_v(sub, outer, inner, Mu, Mv, ef.tv0, 0, offset);
+	stitch_triangles(inner, outer);
+
+	/* right side */
+	add_side_v(sub, outer, inner, Mu, Mv, ef.tv1, 1, offset);
+	stitch_triangles(outer, inner);
+
+	assert(vert_offset == mesh->verts.size());
+}
+
+/* TriangleDice */
+
+TriangleDice::TriangleDice(Mesh *mesh_, int shader_, bool smooth_, float dicing_rate_)
+: EdgeDice(mesh_, shader_, smooth_, dicing_rate_)
+{
+}
+
+void TriangleDice::reserve(EdgeFactors& ef, int M)
+{
+	int num_verts = ef.tu + ef.tv + ef.tw;
+
+	for(int m = M-2; m > 0; m -= 2)
+		num_verts += 3 + (m-1)*3;
+	
+	if(!(M & 1))
+		num_verts++;
+	
+	EdgeDice::reserve(num_verts, 0);
+}
+
+float2 TriangleDice::map_uv(SubPatch& sub, float2 uv)
+{
+	/* map UV from subpatch to patch parametric coordinates */
+	return uv.x*sub.Pu + uv.y*sub.Pv + (1.0f - uv.x - uv.y)*sub.Pw;
+}
+
+int TriangleDice::add_vert(SubPatch& sub, float2 uv)
+{
+	return EdgeDice::add_vert(sub.patch, map_uv(sub, uv));
+}
+
+void TriangleDice::add_grid(SubPatch& sub, EdgeFactors& ef, int M)
+{
+	// XXX normals are flipped, why?
+
+	/* grid is constructed starting from the outside edges, and adding
+	   progressively smaller inner triangles that connected to the outer
+	   one, until M = 1 or 2, the we fill up the last part. */
+	vector<int> outer_u, outer_v, outer_w;
+	int m;
+
+	/* add outer corners vertices */
+	{
+		float2 p_u = make_float2(1.0f, 0.0f);
+		float2 p_v = make_float2(0.0f, 1.0f);
+		float2 p_w = make_float2(0.0f, 0.0f);
+
+		int corner_u = add_vert(sub, p_u);
+		int corner_v = add_vert(sub, p_v);
+		int corner_w = add_vert(sub, p_w);
+
+		outer_u.push_back(corner_v);
+		outer_v.push_back(corner_w);
+		outer_w.push_back(corner_u);
+
+		for(int i = 1; i < ef.tu; i++)
+			outer_u.push_back(add_vert(sub, interp(p_v, p_w, i/(float)ef.tu)));
+		for(int i = 1; i < ef.tv; i++)
+			outer_v.push_back(add_vert(sub, interp(p_w, p_u, i/(float)ef.tv)));
+		for(int i = 1; i < ef.tw; i++)
+			outer_w.push_back(add_vert(sub, interp(p_u, p_v, i/(float)ef.tw)));
+
+		outer_u.push_back(corner_w);
+		outer_v.push_back(corner_u);
+		outer_w.push_back(corner_v);
+	}
+
+	for(m = M-2; m > 0; m -= 2) {
+		vector<int> inner_u, inner_v, inner_w;
+
+		float t = m/(float)M;
+		float2 center = make_float2(1.0f/3.0f, 1.0f/3.0f);
+
+		/* 3 corner vertices */
+		float2 p_u = interp(center, make_float2(1.0f, 0.0f), t);
+		float2 p_v = interp(center, make_float2(0.0f, 1.0f), t);
+		float2 p_w = interp(center, make_float2(0.0f, 0.0f), t);
+
+		int corner_u = add_vert(sub, p_u);
+		int corner_v = add_vert(sub, p_v);
+		int corner_w = add_vert(sub, p_w);
+
+		/* construct array of vertex indices for each side */
+		inner_u.push_back(corner_v);
+		inner_v.push_back(corner_w);
+		inner_w.push_back(corner_u);
+
+		for(int i = 1; i < m; i++) {
+			/* add vertices between corners */
+			float t = i/(float)m;
+
+			inner_u.push_back(add_vert(sub, interp(p_v, p_w, t)));
+			inner_v.push_back(add_vert(sub, interp(p_w, p_u, t)));
+			inner_w.push_back(add_vert(sub, interp(p_u, p_v, t)));
+		}
+
+		inner_u.push_back(corner_w);
+		inner_v.push_back(corner_u);
+		inner_w.push_back(corner_v);
+
+		/* stitch together inner/outer with triangles */
+		stitch_triangles(outer_u, inner_u);
+		stitch_triangles(outer_v, inner_v);
+		stitch_triangles(outer_w, inner_w);
+
+		outer_u = inner_u;
+		outer_v = inner_v;
+		outer_w = inner_w;
+	}
+
+	/* fill up last part */
+	if(m == -1) {
+		/* single triangle */
+		add_triangle(outer_w[0], outer_u[0], outer_v[0]);
+	}
+	else {
+		/* center vertex + 6 triangles */
+		int center = add_vert(sub, make_float2(1.0f/3.0f, 1.0f/3.0f));
+
+		add_triangle(outer_w[0], outer_w[1], center);
+		add_triangle(outer_w[1], outer_w[2], center);
+		add_triangle(outer_u[0], outer_u[1], center);
+		add_triangle(outer_u[1], outer_u[2], center);
+		add_triangle(outer_v[0], outer_v[1], center);
+		add_triangle(outer_v[1], outer_v[2], center);
+	}
+}
+
+void TriangleDice::dice(SubPatch& sub, EdgeFactors& ef)
+{
+	/* todo: handle 2 1 1 resolution */
+	int M = max(ef.tu, max(ef.tv, ef.tw));
+
+	reserve(ef, M);
+	add_grid(sub, ef, M);
+
+	assert(vert_offset == mesh->verts.size());
+}
+
+CCL_NAMESPACE_END
+
diff --git a/intern/cycles/subd/subd_dice.h b/intern/cycles/subd/subd_dice.h
new file mode 100644
index 00000000000..71879f1f5f1
--- /dev/null
+++ b/intern/cycles/subd/subd_dice.h
@@ -0,0 +1,159 @@
+/*
+ * Copyright 2011, Blender Foundation.
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#ifndef __SUBD_DICE_H__
+#define __SUBD_DICE_H__
+
+/* DX11 like EdgeDice implementation, with different tesselation factors for
+ * each edge for watertight tesselation, with subpatch remapping to work with
+ * DiagSplit. For more algorithm details, see the DiagSplit paper or the
+ * ARB_tessellation_shader OpenGL extension, Section 2.X.2. */
+
+#include "util_types.h"
+#include "util_vector.h"
+
+CCL_NAMESPACE_BEGIN
+
+class Camera;
+class Mesh;
+class Patch;
+
+/* EdgeDice Base */
+
+class EdgeDice {
+public:
+	Camera *camera;
+	Mesh *mesh;
+	float3 *mesh_P;
+	float3 *mesh_N;
+	float dicing_rate;
+	size_t vert_offset;
+	size_t tri_offset;
+	int shader;
+	bool smooth;
+
+	EdgeDice(Mesh *mesh, int shader, bool smooth, float dicing_rate);
+
+	void reserve(int num_verts, int num_tris);
+
+	int add_vert(Patch *patch, float2 uv);
+	void add_triangle(int v0, int v1, int v2);
+
+	void stitch_triangles(vector<int>& outer, vector<int>& inner);
+};
+
+/* Quad EdgeDice
+ *
+ * Edge tesselation factors and subpatch coordinates are as follows:
+ *
+ *            tu1
+ *     P01 --------- P11 
+ *     |               |
+ * tv0 |               | tv1
+ *     |               |
+ *     P00 --------- P10
+ *            tu0
+ */
+
+class QuadDice : public EdgeDice {
+public:
+	struct SubPatch {
+		Patch *patch;
+
+		float2 P00;
+		float2 P10;
+		float2 P01;
+		float2 P11;
+	};
+
+	struct EdgeFactors {
+		int tu0;
+		int tu1;
+		int tv0;
+		int tv1;
+	};
+
+	QuadDice(Mesh *mesh, int shader, bool smooth, float dicing_rate);
+
+	void reserve(EdgeFactors& ef, int Mu, int Mv);
+	float3 eval_projected(SubPatch& sub, float u, float v);
+
+	float2 map_uv(SubPatch& sub, float u, float v);
+	int add_vert(SubPatch& sub, float u, float v);
+
+	void add_corners(SubPatch& sub);
+	void add_grid(SubPatch& sub, int Mu, int Mv, int offset);
+
+	void add_side_u(SubPatch& sub,
+		vector<int>& outer, vector<int>& inner,
+		int Mu, int Mv, int tu, int side, int offset);
+
+	void add_side_v(SubPatch& sub,
+		vector<int>& outer, vector<int>& inner,
+		int Mu, int Mv, int tv, int side, int offset);
+
+	float quad_area(const float3& a, const float3& b, const float3& c, const float3& d);
+	float scale_factor(SubPatch& sub, EdgeFactors& ef, int Mu, int Mv);
+
+	void dice(SubPatch& sub, EdgeFactors& ef);
+};
+
+/* Triangle EdgeDice
+ *
+ * Edge tesselation factors and subpatch coordinates are as follows:
+ *
+ *        Pw
+ *        /\
+ *    tv /  \ tu
+ *      /    \
+ *     /      \
+ *  Pu -------- Pv
+ *        tw     
+ */
+
+class TriangleDice : public EdgeDice {
+public:
+	struct SubPatch {
+		Patch *patch;
+
+		float2 Pu;
+		float2 Pv;
+		float2 Pw;
+	};
+
+	struct EdgeFactors {
+		int tu;
+		int tv;
+		int tw;
+	};
+
+	TriangleDice(Mesh *mesh, int shader, bool smooth, float dicing_rate);
+
+	void reserve(EdgeFactors& ef, int M);
+
+	float2 map_uv(SubPatch& sub, float2 uv);
+	int add_vert(SubPatch& sub, float2 uv);
+
+	void add_grid(SubPatch& sub, EdgeFactors& ef, int M);
+	void dice(SubPatch& sub, EdgeFactors& ef);
+};
+
+CCL_NAMESPACE_END
+
+#endif /* __SUBD_DICE_H__ */
+
diff --git a/intern/cycles/subd/subd_edge.h b/intern/cycles/subd/subd_edge.h
new file mode 100644
index 00000000000..edf98c8a5a0
--- /dev/null
+++ b/intern/cycles/subd/subd_edge.h
@@ -0,0 +1,70 @@
+/*
+ * Original code in the public domain -- castanyo@yahoo.es
+ * 
+ * Modifications copyright (c) 2011, Blender Foundation.
+ * All rights reserved.
+ * 
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ * * Redistributions of source code must retain the above copyright
+ *   notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above copyright
+ *   notice, this list of conditions and the following disclaimer in the
+ *   documentation and/or other materials provided with the distribution.
+ * 
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+ * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+ * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef __SUBD_EDGE_H__
+#define __SUBD_EDGE_H__
+
+#include "subd_mesh.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* Subd Half Edge */
+
+class SubdEdge
+{
+public:
+	int id;
+	
+	SubdEdge *next;
+	SubdEdge *prev;
+	SubdEdge *pair;
+	SubdVert *vert;
+	SubdFace *face;
+
+	SubdEdge(int id_)
+	{
+		id = id_;
+		next = NULL;
+		prev = NULL;
+		pair = NULL;
+		vert = NULL;
+		face = NULL;
+	}
+	
+	/* vertex queries */
+	SubdVert *from() { return vert; }
+	SubdVert *to() { return next->vert; }
+
+	/* face queries */
+	bool is_boundary() { return !(face && pair->face); }
+};
+
+CCL_NAMESPACE_END
+
+#endif /* __SUBD_EDGE_H__ */
+
diff --git a/intern/cycles/subd/subd_face.h b/intern/cycles/subd/subd_face.h
new file mode 100644
index 00000000000..8f8a9ec3f47
--- /dev/null
+++ b/intern/cycles/subd/subd_face.h
@@ -0,0 +1,109 @@
+/*
+ * Original code in the public domain -- castanyo@yahoo.es
+ * 
+ * Modifications copyright (c) 2011, Blender Foundation.
+ * All rights reserved.
+ * 
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ * * Redistributions of source code must retain the above copyright
+ *   notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above copyright
+ *   notice, this list of conditions and the following disclaimer in the
+ *   documentation and/or other materials provided with the distribution.
+ * 
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+ * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+ * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef __SUBD_FACE_H__
+#define __SUBD_FACE_H__
+
+#include "subd_edge.h"
+#include "subd_mesh.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* Subd Face */
+
+class SubdFace
+{
+public:
+	int id;
+	SubdEdge *edge;
+
+	SubdFace(int id_)
+	{
+		id = id_;
+		edge = NULL;
+	}
+
+	bool contains(SubdEdge *e)
+	{
+		for(EdgeIterator it(edges()); !it.isDone(); it.advance())
+			if(it.current() == e)
+				return true;
+
+		return false;
+	}
+
+	int num_edges()
+	{
+		int num = 0;
+
+		for(EdgeIterator it(edges()); !it.isDone(); it.advance())
+			num++;
+
+		return num;
+	}
+
+	bool is_boundary()
+	{
+		for(EdgeIterator it(edges()); !it.isDone(); it.advance()) {
+			SubdEdge *edge = it.current();
+
+			if(edge->pair->face == NULL)
+				return true;
+		}
+
+		return false;
+	}
+	
+	/* iterate over edges in clockwise order */
+	class EdgeIterator
+	{
+	public:
+		EdgeIterator(SubdEdge *e) : end(NULL), cur(e) { }
+
+		virtual void advance()
+		{
+			if (end == NULL) end = cur;
+			cur = cur->next;
+		}
+
+		virtual bool isDone() { return end == cur; }
+		virtual SubdEdge *current() { return cur; }
+
+	private:
+		SubdEdge *end;
+		SubdEdge *cur;
+	};
+
+	EdgeIterator edges() { return EdgeIterator(edge); }
+	EdgeIterator edges(SubdEdge *edge) {  return EdgeIterator(edge);  }
+};
+
+CCL_NAMESPACE_END
+
+#endif /* __SUBD_FACE_H__ */
+
diff --git a/intern/cycles/subd/subd_mesh.cpp b/intern/cycles/subd/subd_mesh.cpp
new file mode 100644
index 00000000000..8ed4854b786
--- /dev/null
+++ b/intern/cycles/subd/subd_mesh.cpp
@@ -0,0 +1,309 @@
+/*
+ * Original code in the public domain -- castanyo@yahoo.es
+ * 
+ * Modifications copyright (c) 2011, Blender Foundation.
+ * All rights reserved.
+ * 
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ * * Redistributions of source code must retain the above copyright
+ *   notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above copyright
+ *   notice, this list of conditions and the following disclaimer in the
+ *   documentation and/or other materials provided with the distribution.
+ * 
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+ * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+ * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include <stdio.h>
+
+#include "subd_build.h"
+#include "subd_edge.h"
+#include "subd_face.h"
+#include "subd_mesh.h"
+#include "subd_patch.h"
+#include "subd_split.h"
+#include "subd_vert.h"
+
+#include "util_debug.h"
+#include "util_foreach.h"
+
+CCL_NAMESPACE_BEGIN
+
+SubdMesh::SubdMesh()
+{
+}
+
+SubdMesh::~SubdMesh()
+{
+	pair<Key, SubdEdge*> em;
+
+	foreach(SubdVert *vertex, verts)
+		delete vertex;
+	foreach(em, edge_map)
+		delete em.second;
+	foreach(SubdFace *face, faces)
+		delete face;
+
+	verts.clear();
+	edges.clear();
+	edge_map.clear();
+	faces.clear();
+}
+
+SubdVert *SubdMesh::add_vert(const float3& co)
+{
+	SubdVert *v = new SubdVert(verts.size());
+	v->co = co;
+	verts.push_back(v);
+
+	return v;
+}
+
+SubdFace *SubdMesh::add_face(int v0, int v1, int v2)
+{
+	int index[3] = {v0, v1, v2};
+	return add_face(index, 3);
+}
+
+SubdFace *SubdMesh::add_face(int v0, int v1, int v2, int v3)
+{
+	int index[4] = {v0, v1, v2, v3};
+	return add_face(index, 4);
+}
+
+SubdFace *SubdMesh::add_face(int *index, int num)
+{
+	/* test non-manifold cases */
+	if(!can_add_face(index, num)) {
+		/* we could try to add face in opposite winding instead .. */
+		fprintf(stderr, "Warning: non manifold mesh, invalid face '%lu'.\n", faces.size());
+		return NULL;
+	}
+	
+	SubdFace *f = new SubdFace(faces.size());
+	
+	SubdEdge *first_edge = NULL;
+	SubdEdge *last = NULL;
+	SubdEdge *current = NULL;
+
+	/* add edges */
+	for(int i = 0; i < num-1; i++) {
+		current = add_edge(index[i], index[i+1]);
+		assert(current != NULL);
+		
+		current->face = f;
+		
+		if(last != NULL) {
+			last->next = current;
+			current->prev = last;
+		}
+		else
+			first_edge = current;
+		
+		last = current;
+	}
+
+	current = add_edge(index[num-1], index[0]);
+	assert(current != NULL);
+	
+	current->face = f;
+
+	last->next = current;
+	current->prev = last;
+
+	current->next = first_edge;
+	first_edge->prev = current;
+
+	f->edge = first_edge;
+	faces.push_back(f);
+
+	return f;
+}
+
+bool SubdMesh::can_add_face(int *index, int num)
+{
+	/* manifold check */
+	for(int i = 0; i < num-1; i++)
+		if(!can_add_edge(index[i], index[i+1]))
+			return false;
+
+	return can_add_edge(index[num-1], index[0]);
+}
+
+bool SubdMesh::can_add_edge(int i, int j)
+{
+	/* check for degenerate edge */
+	if(i == j)
+		return false;
+
+	/* make sure edge has not been added yet. */
+	return find_edge(i, j) == NULL;
+}
+
+SubdEdge *SubdMesh::add_edge(int i, int j)
+{
+	SubdEdge *edge;
+
+	/* find pair */
+	SubdEdge *pair = find_edge(j, i);
+
+	if(pair != NULL) {
+		/* create edge with same id */
+		edge = new SubdEdge(pair->id + 1);
+		
+		/* link edge pairs */
+		edge->pair = pair;
+		pair->pair = edge;
+		
+		/* not sure this is necessary? */
+		pair->vert->edge = pair;
+	}
+	else {
+		/* create edge */
+		edge = new SubdEdge(2*edges.size());
+		
+		/* add only unpaired edges */
+		edges.push_back(edge);
+	}
+	
+	/* assign vertex and put into map */
+	edge->vert = verts[i];
+	edge_map[Key(i, j)] = edge;
+	
+	/* face and next are set by add_face */
+	
+	return edge;
+}
+
+SubdEdge *SubdMesh::find_edge(int i, int j)
+{
+	map<Key, SubdEdge*>::const_iterator it = edge_map.find(Key(i, j));
+
+	return (it == edge_map.end())? NULL: it->second;
+}
+
+bool SubdMesh::link_boundary()
+{
+	/* link boundary edges once the mesh has been created */
+	int num = 0;
+	
+	/* create boundary edges */
+	int num_edges = edges.size();
+
+	for(int e = 0; e < num_edges; e++) {
+		SubdEdge *edge = edges[e];
+
+		if(edge->pair == NULL) {
+			SubdEdge *pair = new SubdEdge(edge->id + 1);
+
+			int i = edge->from()->id;
+			int j = edge->to()->id;
+
+			assert(edge_map.find(Key(j, i)) == edge_map.end());
+
+			pair->vert = verts[j];
+			edge_map[Key(j, i)] = pair;
+			
+			edge->pair = pair;
+			pair->pair = edge;
+			
+			num++;
+		}
+	}
+
+	/* link boundary edges */
+	for(int e = 0; e < num_edges; e++) {
+		SubdEdge *edge = edges[e];
+
+		if(edge->pair->face == NULL)
+			link_boundary_edge(edge->pair);
+	}
+	
+	/* detect boundary intersections */
+	int boundaryIntersections = 0;
+	int num_verts = verts.size();
+
+	for(int v = 0; v < num_verts; v++) {
+		SubdVert *vertex = verts[v];
+
+		int boundarySubdEdges = 0;
+		for(SubdVert::EdgeIterator it(vertex->edges()); !it.isDone(); it.advance())
+			if(it.current()->is_boundary())
+				boundarySubdEdges++;
+
+		if(boundarySubdEdges > 2) {
+			assert((boundarySubdEdges & 1) == 0);
+			boundaryIntersections++;
+		}
+	}
+
+	if(boundaryIntersections != 0) {
+		fprintf(stderr, "Invalid mesh, boundary intersections found!\n");
+		return false;
+	}
+
+	return true;
+}
+
+void SubdMesh::link_boundary_edge(SubdEdge *edge)
+{
+	/* link this boundary edge. */
+
+	/* make sure next pointer has not been set. */
+	assert(edge->face == NULL);
+	assert(edge->next == NULL);
+		
+	SubdEdge *next = edge;
+
+	while(next->pair->face != NULL) {
+		/* get pair prev */
+		SubdEdge *e = next->pair->next;
+
+		while(e->next != next->pair)
+			e = e->next;
+
+		next = e;
+	}
+
+	edge->next = next->pair;
+	next->pair->prev = edge;
+	
+	/* adjust vertex edge, so that it's the boundary edge. (required for is_boundary()) */
+	if(edge->vert->edge != edge)
+		edge->vert->edge = edge;
+}
+
+void SubdMesh::tesselate(DiagSplit *split, bool linear, Mesh *mesh, int shader, bool smooth)
+{
+	SubdBuilder *builder = SubdBuilder::create(linear);
+	int num_faces = faces.size();
+		        
+	for(int f = 0; f < num_faces; f++) {
+		SubdFace *face = faces[f];
+		Patch *patch = builder->run(face);
+
+		if(patch->is_triangle())
+			split->split_triangle(mesh, patch, shader, smooth);
+		else
+			split->split_quad(mesh, patch, shader, smooth);
+
+		delete patch;
+	}
+
+	delete builder;
+}
+
+CCL_NAMESPACE_END
+
diff --git a/intern/cycles/subd/subd_mesh.h b/intern/cycles/subd/subd_mesh.h
new file mode 100644
index 00000000000..999d92d6daf
--- /dev/null
+++ b/intern/cycles/subd/subd_mesh.h
@@ -0,0 +1,90 @@
+/*
+ * Original code in the public domain -- castanyo@yahoo.es
+ *
+ * Modifications copyright (c) 2011, Blender Foundation.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ * * Redistributions of source code must retain the above copyright
+ *   notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above copyright
+ *   notice, this list of conditions and the following disclaimer in the
+ *   documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+ * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+ * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef __SUBD_MESH_H__
+#define __SUBD_MESH_H__
+
+#include "util_map.h"
+#include "util_types.h"
+#include "util_vector.h"
+
+CCL_NAMESPACE_BEGIN
+
+class SubdFace;
+class SubdVert;
+class SubdEdge;
+
+class DiagSplit;
+class Mesh;
+
+/* Subd Mesh, half edge based for dynamic mesh manipulation */
+
+class SubdMesh
+{
+public:
+	vector<SubdVert*> verts;
+	vector<SubdEdge*> edges;
+	vector<SubdFace*> faces;
+
+	SubdMesh();
+	~SubdMesh();
+
+	SubdVert *add_vert(const float3& co);
+	
+	SubdFace *add_face(int v0, int v1, int v2);
+	SubdFace *add_face(int v0, int v1, int v2, int v3);
+	SubdFace *add_face(int *index, int num);
+
+	bool link_boundary();
+	void tesselate(DiagSplit *split, bool linear,
+		Mesh *mesh, int shader, bool smooth);
+
+protected:
+	bool can_add_face(int *index, int num);
+	bool can_add_edge(int i, int j);
+	SubdEdge *add_edge(int i, int j);
+	SubdEdge *find_edge(int i, int j);
+	void link_boundary_edge(SubdEdge *edge);
+	
+	struct Key {
+		Key() {}
+		Key(int v0, int v1) : p0(v0), p1(v1) {}
+
+		bool operator<(const Key& k) const
+		{ return (p0 < k.p0 || (p0 == k.p0 && p1 < k.p1)); }
+
+		int p0, p1;
+	};
+
+	map<Key, SubdEdge *> edge_map;
+};
+
+CCL_NAMESPACE_END
+
+#endif /* __SUBD_MESH_H__ */
+
diff --git a/intern/cycles/subd/subd_patch.cpp b/intern/cycles/subd/subd_patch.cpp
new file mode 100644
index 00000000000..ff477296c7e
--- /dev/null
+++ b/intern/cycles/subd/subd_patch.cpp
@@ -0,0 +1,288 @@
+/*
+ * Copyright 2011, Blender Foundation.
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/* Parts adapted from code in the public domain in NVidia Mesh Tools. */
+
+#include "mesh.h"
+
+#include "subd_patch.h"
+
+#include "util_math.h"
+#include "util_types.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* De Casteljau Evaluation */
+
+static float3 decasteljau_quadratic(float t, const float3 cp[3])
+{
+	float3 d0 = cp[0] + t*(cp[1] - cp[0]);
+	float3 d1 = cp[1] + t*(cp[2] - cp[1]);
+
+	return d0 + t*(d1 - d0);
+}
+
+static void decasteljau_cubic(float3 *P, float3 *dt, float t, const float3 cp[4])
+{
+	float3 d0 = cp[0] + t*(cp[1] - cp[0]);
+	float3 d1 = cp[1] + t*(cp[2] - cp[1]);
+	float3 d2 = cp[2] + t*(cp[3] - cp[2]);
+
+	d0 += t*(d1 - d0);
+	d1 += t*(d2 - d1);
+
+	*P = d0 + t*(d1 - d0);
+	if(dt) *dt = d1 - d0;
+}
+
+static void decasteljau_bicubic(float3 *P, float3 *du, float3 *dv, const float3 cp[16], float u, float v)
+{
+	float3 ucp[4], utn[4];
+
+	/* interpolate over u */
+	decasteljau_cubic(ucp+0, utn+0, u, cp);
+	decasteljau_cubic(ucp+1, utn+1, u, cp+4);
+	decasteljau_cubic(ucp+2, utn+2, u, cp+8);
+	decasteljau_cubic(ucp+3, utn+3, u, cp+12);
+
+	/* interpolate over v */
+	decasteljau_cubic(P, dv, v, ucp);
+	if(du) decasteljau_cubic(du, NULL, v, utn);
+}
+
+static float3 decasteljau_tangent(const float3 cp[12], float u, float v)
+{
+	float3 ucp[3];
+
+	decasteljau_cubic(ucp+0, NULL, v, cp);
+	decasteljau_cubic(ucp+1, NULL, v, cp+4);
+	decasteljau_cubic(ucp+2, NULL, v, cp+8);
+
+	return decasteljau_quadratic(u, ucp);
+}
+
+/* Linear Quad Patch */
+
+void LinearQuadPatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
+{
+	float3 d0 = interp(hull[0], hull[1], u);
+	float3 d1 = interp(hull[2], hull[3], u);
+
+	*P = interp(d0, d1, v);
+
+	if(dPdu && dPdv) {
+		*dPdu = interp(hull[1] - hull[0], hull[3] - hull[2], v);
+		*dPdv = interp(hull[2] - hull[0], hull[3] - hull[1], u);
+	}
+}
+
+BoundBox LinearQuadPatch::bound()
+{
+	BoundBox bbox;
+
+	for(int i = 0; i < 4; i++)
+		bbox.grow(hull[i]);
+	
+	return bbox;
+}
+
+/* Linear Triangle Patch */
+
+void LinearTrianglePatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
+{
+	*P = u*hull[0] + v*hull[1] + (1.0f - u - v)*hull[2];
+
+	if(dPdu && dPdv) {
+		*dPdu = hull[0] - hull[2];
+		*dPdv = hull[1] - hull[2];
+	}
+}
+
+BoundBox LinearTrianglePatch::bound()
+{
+	BoundBox bbox;
+
+	for(int i = 0; i < 3; i++)
+		bbox.grow(hull[i]);
+	
+	return bbox;
+}
+
+/* Bicubic Patch */
+
+void BicubicPatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
+{
+	decasteljau_bicubic(P, dPdu, dPdv, hull, u, v);
+}
+
+BoundBox BicubicPatch::bound()
+{
+	BoundBox bbox;
+
+	for(int i = 0; i < 16; i++)
+		bbox.grow(hull[i]);
+	
+	return bbox;
+}
+
+/* Bicubic Patch with Tangent Fields */
+
+void BicubicTangentPatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
+{
+	decasteljau_bicubic(P, NULL, NULL, hull, u, v);
+
+	if(dPdu) *dPdu = decasteljau_tangent(utan, u, v);
+	if(dPdv) *dPdv = decasteljau_tangent(vtan, v, u);
+}
+
+BoundBox BicubicTangentPatch::bound()
+{
+	BoundBox bbox;
+
+	for(int i = 0; i < 16; i++)
+		bbox.grow(hull[i]);
+	
+	return bbox;
+}
+
+/* Gregory Patch */
+
+static float no_zero_div(float f)
+{
+	if(f == 0.0f) return 1.0f;
+	return f;
+}
+
+void GregoryQuadPatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
+{
+	float3 bicubic[16];
+
+	float U = 1 - u;
+	float V = 1 - v;
+
+	/*  8     9     10     11
+	 * 12   0\1     2/3    13
+	 * 14   4/5     6\7    15
+	 * 16    17     18     19
+	 */
+
+	bicubic[5] = (u*hull[1] + v*hull[0])/no_zero_div(u + v);
+	bicubic[6] = (U*hull[2] + v*hull[3])/no_zero_div(U + v);
+	bicubic[9] = (u*hull[5] + V*hull[4])/no_zero_div(u + V);
+	bicubic[10] = (U*hull[6] + V*hull[7])/no_zero_div(U + V);
+
+	// Map gregory control points to bezier control points.
+	bicubic[0] = hull[8];
+	bicubic[1] = hull[9];
+	bicubic[2] = hull[10];
+	bicubic[3] = hull[11];
+	bicubic[4] = hull[12];
+	bicubic[7] = hull[13];
+	bicubic[8] = hull[14];
+	bicubic[11] = hull[15];
+	bicubic[12] = hull[16];
+	bicubic[13] = hull[17];
+	bicubic[14] = hull[18];
+	bicubic[15] = hull[19];
+
+	decasteljau_bicubic(P, dPdu, dPdv, bicubic, u, v);
+}
+
+BoundBox GregoryQuadPatch::bound()
+{
+	BoundBox bbox;
+
+	for(int i = 0; i < 20; i++)
+		bbox.grow(hull[i]);
+	
+	return bbox;
+}
+
+void GregoryTrianglePatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
+{
+	/*		      6
+	 *		      
+	 *		14   0/1   7
+	 *						  
+	 *    13   5/4     3\2   8
+	 *		       
+	 * 12      11       10      9
+	 */
+
+	float w = 1 - u - v;
+	float uu = u * u;
+	float vv = v * v;
+	float ww = w * w;
+	float uuu = uu * u;
+	float vvv = vv * v;
+	float www = ww * w;
+
+	float U = 1 - u;
+	float V = 1 - v;
+	float W = 1 - w;
+
+	float3 C0 = ( v*U * hull[5] + u*V * hull[4] ) / no_zero_div(v*U + u*V);
+	float3 C1 = ( w*V * hull[3] + v*W * hull[2] ) / no_zero_div(w*V + v*W);
+	float3 C2 = ( u*W * hull[1] + w*U * hull[0] ) / no_zero_div(u*W + w*U);
+
+	*P =
+		(hull[12] * www + 3*hull[11] * ww*u + 3*hull[10] * w*uu + hull[ 9]*uuu) * (w + u) +
+		(hull[ 9] * uuu + 3*hull[ 8] * uu*v + 3*hull[ 7] * u*vv + hull[ 6]*vvv) * (u + v) +
+		(hull[ 6] * vvv + 3*hull[14] * vv*w + 3*hull[13] * v*ww + hull[12]*www) * (v + w) -
+		(hull[12] * www*w + hull[ 9] * uuu*u + hull[ 6] * vvv*v) +
+		12*(C0 * u*v*ww + C1 * uu*v*w   + C2 * u*vv*w);
+
+	if(dPdu || dPdv) {
+		float3 E1 = (hull[12]*www + 3*hull[11]*ww*u + 3*hull[10]*w*uu + hull[ 9]*uuu);
+		float3 E2 = (hull[ 9]*uuu + 3*hull[ 8]*uu*v + 3*hull[ 7]*u*vv + hull[ 6]*vvv);
+		float3 E3 = (hull[ 6]*vvv + 3*hull[14]*vv*w + 3*hull[13]*v*ww + hull[12]*www);
+
+		if(dPdu) {
+			float3 E1u = 3*( - hull[12]*ww + hull[11]*(ww-2*u*w) +   hull[10]*(2*u*w-uu) + hull[ 9]*uu);
+			float3 E2u = 3*(   hull[ 9]*uu + 2*hull[ 8]*u*v      +   hull[ 7]*vv		 );
+			float3 E3u = 3*(		    - hull[14]*vv		 - 2*hull[13]*v*w		- hull[12]*ww);
+			float3 Su  = 4*( -hull[12]*www + hull[9]*uuu);
+			float3 Cu  = 12*( C0*(ww*v-2*u*v*w) + C1*(2*u*v*w-uu*v) + C2*vv*(w-u) );
+
+			*dPdu = E1u*(w+u) + (E2+E2u*(u+v)) + (E3u*(v+w)-E3) - Su + Cu;
+		}
+
+		if(dPdv) {
+			float3 E1v = 3*(-hull[12]*ww  - 2*hull[11]*w*u       -   hull[10]*uu		 );
+			float3 E2v = 3*(		      hull[ 8]*uu		 + 2*hull[ 7]*u*v		+ hull[ 6]*vv);
+			float3 E3v = 3*( hull[ 6]*vv  +  hull[14]*(2*w*v-vv) +   hull[13]*(ww-2*w*v) - hull[12]*ww);
+			float3 Sv  = 4*(-hull[12]*www +  hull[ 6]*vvv);
+			float3 Cv  = 12*(C0*(u*ww-2*u*v*w) + C1*uu*(w-v) + C2*(2*u*v*w-u*vv));
+
+			*dPdv = ((E1v*(w+u)-E1) + (E2+E2v*(u+v)) + E3v*(v+w) - Sv + Cv );
+		}
+	}
+}
+
+BoundBox GregoryTrianglePatch::bound()
+{
+	BoundBox bbox;
+
+	for(int i = 0; i < 20; i++)
+		bbox.grow(hull[i]);
+	
+	return bbox;
+}
+
+CCL_NAMESPACE_END
+
diff --git a/intern/cycles/subd/subd_patch.h b/intern/cycles/subd/subd_patch.h
new file mode 100644
index 00000000000..8d4b8a1c911
--- /dev/null
+++ b/intern/cycles/subd/subd_patch.h
@@ -0,0 +1,107 @@
+/*
+ * Copyright 2011, Blender Foundation.
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#ifndef __SUBD_PATCH_H__
+#define __SUBD_PATCH_H__
+
+#include "util_boundbox.h"
+#include "util_types.h"
+
+CCL_NAMESPACE_BEGIN
+
+class Mesh;
+
+/* Base */
+
+class Patch {
+public:
+	virtual void eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v) = 0;
+	virtual bool is_triangle() = 0;
+	virtual BoundBox bound() = 0;
+};
+
+/* Linear Quad Patch */
+
+class LinearQuadPatch : public Patch {
+public:
+	float3 hull[4];
+
+	void eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v);
+	bool is_triangle() { return false; }
+	BoundBox bound();
+};
+
+/* Linear Triangle Patch */
+
+class LinearTrianglePatch : public Patch {
+public:
+	float3 hull[3];
+
+	void eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v);
+	bool is_triangle() { return true; }
+	BoundBox bound();
+};
+
+/* Bicubic Patch */
+
+class BicubicPatch : public Patch {
+public:
+	float3 hull[16];
+
+	void eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v);
+	bool is_triangle() { return false; }
+	BoundBox bound();
+};
+
+/* Bicubic Patch with Tangent Fields */
+
+class BicubicTangentPatch : public Patch {
+public:
+	float3 hull[16];
+	float3 utan[12];
+	float3 vtan[12];
+
+	void eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v);
+	bool is_triangle() { return false; }
+	BoundBox bound();
+};
+
+/* Gregory Patches */
+
+class GregoryQuadPatch : public Patch {
+public:
+	float3 hull[20];
+
+	void eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v);
+	bool is_triangle() { return false; }
+	BoundBox bound();
+};
+
+class GregoryTrianglePatch : public Patch  {
+public:
+	float3 hull[20];
+
+	void eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v);
+	bool is_triangle() { return true; }
+	BoundBox bound();
+};
+
+CCL_NAMESPACE_END
+
+#endif /* __SUBD_PATCH_H__ */
+
diff --git a/intern/cycles/subd/subd_ring.cpp b/intern/cycles/subd/subd_ring.cpp
new file mode 100644
index 00000000000..cbd12e60da0
--- /dev/null
+++ b/intern/cycles/subd/subd_ring.cpp
@@ -0,0 +1,236 @@
+/*
+ * Copyright 2006, NVIDIA Corporation Ignacio Castano <icastano@nvidia.com>
+ * 
+ * Modifications copyright (c) 2011, Blender Foundation.
+ * All rights reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ * 
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ * 
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ */
+
+#include "subd_face.h"
+#include "subd_ring.h"
+#include "subd_stencil.h"
+#include "subd_vert.h"
+
+#include "util_algorithm.h"
+#include "util_debug.h"
+#include "util_math.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* Utility for sorting verts by index */
+
+class vertex_compare
+{
+public:
+	const vector<int>& index;
+	vertex_compare(const vector<int>& index_) : index(index_) {}
+	bool operator()(int a, int b) const { return index[a] < index[b]; }
+};
+
+/* Subd Face Ring */
+
+SubdFaceRing::SubdFaceRing(SubdFace *face, SubdEdge *firstEdge)
+{
+	m_face = face;
+	m_firstEdge = firstEdge;
+	m_num_edges = face->num_edges();
+
+	assert(m_num_edges == 3 || m_num_edges == 4);
+
+	initVerts();
+}
+
+int SubdFaceRing::num_verts()
+{
+	return m_verts.size();
+}
+
+SubdVert *SubdFaceRing::vertexAt(int i)
+{
+	return m_verts[i];
+}
+
+int SubdFaceRing::vert_index(SubdVert *vertex)
+{
+	int count = this->num_verts();
+
+	for(int i = 0; i < count; i++)
+		if(m_verts[i] == vertex)
+			return i;
+	
+	assert(0);
+	return 0;
+}
+
+void SubdFaceRing::evaluate_stencils(float3 *P, StencilMask *mask, int num)
+{
+	/* first we sort verts by id. this way verts will always be added
+	   in the same order to ensure the exact same float ops happen for control
+	   points of other patches, so we get water-tight patches */
+	int num_verts = m_verts.size();
+
+	vector<int> vmap(num_verts);
+	vector<int> vertid(num_verts);
+
+	for(int v = 0; v < num_verts; v++) {
+		vmap[v] = v;
+		vertid[v] = m_verts[v]->id;
+	}
+
+	sort(vmap.begin(), vmap.end(), vertex_compare(vertid));
+
+	/* then evaluate the stencils */
+	for(int j = 0; j < num; j++) {
+		float3 p = make_float3(0.0f, 0.0f, 0.0f);
+
+		for(int i = 0; i < num_verts; i++) {
+			int idx = vmap[i];
+			p += m_verts[idx]->co * mask[j][idx];
+		}
+
+		P[j] = p;
+	}
+}
+
+bool SubdFaceRing::is_triangle()
+{
+	return m_num_edges == 3;
+}
+
+bool SubdFaceRing::is_quad()
+{
+	return m_num_edges == 4;
+}
+
+int SubdFaceRing::num_edges()
+{
+	return m_num_edges;
+}
+
+bool SubdFaceRing::is_regular(SubdFace *face)
+{
+	if(!is_quad(face))
+		return false;
+
+	for(SubdFace::EdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+		SubdEdge *edge = it.current();
+
+		/* regular faces don't have boundary edges */
+		if(edge->is_boundary())
+			return false;
+
+		/* regular faces only have ordinary verts */
+		if(edge->vert->valence() != 4)
+			return false;
+
+		/* regular faces are only adjacent to quads */
+		for(SubdVert::EdgeIterator eit(edge); !eit.isDone(); eit.advance())
+			if(eit.current()->face == NULL || eit.current()->face->num_edges() != 4)
+				return false;
+	}
+
+	return true;
+}
+
+bool SubdFaceRing::is_triangle(SubdFace *face)
+{
+	return face->num_edges() == 3;
+}
+bool SubdFaceRing::is_quad(SubdFace *face)
+{
+	return face->num_edges() == 4;
+}
+
+bool SubdFaceRing::is_boundary(SubdFace *face)
+{
+	/* note that face->is_boundary() returns a different result. That function
+	   returns true when any of the *edges* are on the boundary. however, this
+	   function returns true if any of the face *verts* are on the boundary.  */
+
+	for(SubdFace::EdgeIterator it(face->edges()); !it.isDone(); it.advance()) {
+		SubdEdge *edge = it.current();
+
+		if(edge->vert->is_boundary())
+			return true;
+	}
+
+	return false;
+}
+
+void SubdFaceRing::initVerts()
+{
+	m_verts.reserve(16);
+
+	/* add face verts */
+	for(SubdFace::EdgeIterator it(m_firstEdge); !it.isDone(); it.advance()) {
+		SubdEdge *edge = it.current();
+		m_verts.push_back(edge->from());
+	}
+	
+	// @@ Add support for non manifold surfaces!
+	// The fix: 
+	// - not all colocals should point to the same edge.
+	// - multiple colocals could belong to different boundaries, make sure they point to the right one.
+
+	// @@ When the face neighborhood wraps that could result in overlapping stencils. 
+
+	// Add surronding verts.
+	for(SubdFace::EdgeIterator it(m_firstEdge); !it.isDone(); it.advance()) {
+		SubdEdge *firstEdge = it.current();
+	
+		int i = 0;
+
+		/* traverse edges around vertex */
+		for(SubdVert::ReverseEdgeIterator eit(firstEdge); !eit.isDone(); eit.advance(), i++) {
+			SubdEdge *edge = eit.current();
+
+			assert(edge->from()->co == firstEdge->from()->co);
+
+			add_vert(edge->to());
+
+			if(edge->face != NULL)
+				add_vert(edge->next->to());
+		}
+
+		assert(i == firstEdge->from()->valence());
+	}
+}
+
+
+void SubdFaceRing::add_vert(SubdVert *vertex)
+{
+	/* avoid having duplicate verts */
+	if(!has_vert(vertex))
+		m_verts.push_back(vertex);
+}
+
+bool SubdFaceRing::has_vert(SubdVert *vertex)
+{
+	int num_verts = m_verts.size();
+
+	for(int i = 0; i < num_verts; i++)
+		if(m_verts[i]->co == vertex->co)
+			return true;
+
+	return false;
+}
+
+CCL_NAMESPACE_END
+
diff --git a/intern/cycles/subd/subd_ring.h b/intern/cycles/subd/subd_ring.h
new file mode 100644
index 00000000000..f25342233b0
--- /dev/null
+++ b/intern/cycles/subd/subd_ring.h
@@ -0,0 +1,75 @@
+/*
+ * Copyright 2006, NVIDIA Corporation Ignacio Castano <icastano@nvidia.com>
+ *
+ * Modifications copyright (c) 2011, Blender Foundation.
+ * All rights reserved.
+ * 
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ * 
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ * 
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ */
+
+#ifndef __SUBD_FACE_RING_H__
+#define __SUBD_FACE_RING_H__
+
+CCL_NAMESPACE_BEGIN
+
+class StencilMask;
+class SubdVert;
+class SubdEdge;
+class SubdFace;
+
+class SubdFaceRing
+{
+public:
+	SubdFaceRing(SubdFace *face, SubdEdge *edge);
+
+	SubdFace *face() { return m_face; }
+	SubdEdge *firstEdge() { return m_firstEdge; }
+
+	int num_verts();
+	SubdVert *vertexAt(int i);
+	int vert_index(SubdVert *vertex);
+
+	void evaluate_stencils(float3 *P, StencilMask *mask, int num);
+
+	bool is_triangle();
+	bool is_quad();
+	int num_edges();
+
+	static bool is_regular(SubdFace *face);
+	static bool is_triangle(SubdFace *face);
+	static bool is_quad(SubdFace *face);
+	static bool is_boundary(SubdFace *face);
+
+protected:
+	void initVerts();
+	void add_vert(SubdVert *vertex);
+	bool has_vert(SubdVert *vertex);
+
+protected:
+	SubdFace *m_face;
+	SubdEdge *m_firstEdge;
+
+	int m_num_edges;
+	vector<SubdVert*> m_verts;
+};
+
+CCL_NAMESPACE_END
+
+#endif /* __SUBD_FACE_RING_H__ */
+
diff --git a/intern/cycles/subd/subd_split.cpp b/intern/cycles/subd/subd_split.cpp
new file mode 100644
index 00000000000..712bd041e64
--- /dev/null
+++ b/intern/cycles/subd/subd_split.cpp
@@ -0,0 +1,325 @@
+/*
+ * Copyright 2011, Blender Foundation.
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include "camera.h"
+#include "mesh.h"
+
+#include "subd_dice.h"
+#include "subd_patch.h"
+#include "subd_split.h"
+
+#include "util_debug.h"
+#include "util_math.h"
+#include "util_types.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* DiagSplit */
+
+DiagSplit::DiagSplit()
+{
+	test_steps = 3;
+	split_threshold = 1;
+	dicing_rate = 0.1f;
+	camera = NULL;
+}
+
+void DiagSplit::dispatch(QuadDice::SubPatch& sub, QuadDice::EdgeFactors& ef)
+{
+	subpatches_quad.push_back(sub);
+	edgefactors_quad.push_back(ef);
+}
+
+void DiagSplit::dispatch(TriangleDice::SubPatch& sub, TriangleDice::EdgeFactors& ef)
+{
+	subpatches_triangle.push_back(sub);
+	edgefactors_triangle.push_back(ef);
+}
+
+float3 DiagSplit::project(Patch *patch, float2 uv)
+{
+	float3 P;
+
+	patch->eval(&P, NULL, NULL, uv.x, uv.y);
+	if(camera)
+		P = transform(&camera->worldtoraster, P);
+
+	return P;
+}
+
+int DiagSplit::T(Patch *patch, float2 Pstart, float2 Pend)
+{
+	float3 Plast = make_float3(0.0f, 0.0f, 0.0f);
+	float Lsum = 0.0f;
+	float Lmax = 0.0f;
+
+	for(int i = 0; i < test_steps; i++) {
+		float t = i/(float)(test_steps-1);
+
+		float3 P = project(patch, Pstart + t*(Pend - Pstart));
+
+		if(i > 0) {
+			float L = len(P - Plast);
+			Lsum += L;
+			Lmax = max(L, Lmax);
+		}
+
+		Plast = P;
+	}
+
+	int tmin = ceil(Lsum/dicing_rate);
+	int tmax = ceil((test_steps-1)*Lmax/dicing_rate); // XXX paper says N instead of N-1, seems wrong?
+
+	if(tmax - tmin > split_threshold)
+		return DSPLIT_NON_UNIFORM;
+	
+	return tmax;
+}
+
+void DiagSplit::partition_edge(Patch *patch, float2 *P, int *t0, int *t1, float2 Pstart, float2 Pend, int t)
+{
+	if(t == DSPLIT_NON_UNIFORM) {
+		*P = (Pstart + Pend)*0.5f;
+		*t0 = T(patch, Pstart, *P);
+		*t1 = T(patch, *P, Pend);
+	}
+	else {
+		int I = floor(t*0.5f);
+		*P = interp(Pstart, Pend, (t == 0)? 0: I/(float)t); /* XXX is t faces or verts */
+		*t0 = I;
+		*t1 = t - I;
+	}
+}
+
+void DiagSplit::split(TriangleDice::SubPatch& sub, TriangleDice::EdgeFactors& ef, int depth)
+{
+	assert(ef.tu == T(sub.patch, sub.Pv, sub.Pw));
+	assert(ef.tv == T(sub.patch, sub.Pw, sub.Pu));
+	assert(ef.tw == T(sub.patch, sub.Pu, sub.Pv));
+
+	if(depth == 0) {
+		dispatch(sub, ef);
+		return;
+	}
+
+	if(ef.tu == DSPLIT_NON_UNIFORM) {
+		/* partition edges */
+		TriangleDice::EdgeFactors ef0, ef1;
+		float2 Psplit;
+
+		partition_edge(sub.patch,
+			&Psplit, &ef1.tu, &ef0.tu, sub.Pv, sub.Pw, ef.tu);
+
+		/* split */
+		int tsplit = T(sub.patch, sub.Pu, Psplit);
+		ef0.tv = ef.tv;
+		ef0.tw = tsplit;
+
+		ef1.tv = tsplit;
+		ef1.tw = ef.tw;
+
+		/* create subpatches */
+		TriangleDice::SubPatch sub0 = {sub.patch, sub.Pu, Psplit, sub.Pw};
+		TriangleDice::SubPatch sub1 = {sub.patch, sub.Pu, sub.Pv, Psplit};
+
+		split(sub0, ef0, depth+1);
+		split(sub1, ef1, depth+1);
+	}
+	else if(ef.tv == DSPLIT_NON_UNIFORM) {
+		/* partition edges */
+		TriangleDice::EdgeFactors ef0, ef1;
+		float2 Psplit;
+
+		partition_edge(sub.patch,
+			&Psplit, &ef1.tu, &ef0.tu, sub.Pw, sub.Pu, ef.tv);
+
+		/* split */
+		int tsplit = T(sub.patch, sub.Pv, Psplit);
+		ef0.tv = ef.tw;
+		ef0.tw = tsplit;
+
+		ef1.tv = tsplit;
+		ef1.tw = ef.tu;
+
+		/* create subpatches */
+		TriangleDice::SubPatch sub0 = {sub.patch, sub.Pv, Psplit, sub.Pu};
+		TriangleDice::SubPatch sub1 = {sub.patch, sub.Pv, sub.Pw, Psplit};
+
+		split(sub0, ef0, depth+1);
+		split(sub1, ef1, depth+1);
+	}
+	else if(ef.tw == DSPLIT_NON_UNIFORM) {
+		/* partition edges */
+		TriangleDice::EdgeFactors ef0, ef1;
+		float2 Psplit;
+
+		partition_edge(sub.patch,
+			&Psplit, &ef1.tu, &ef0.tu, sub.Pu, sub.Pv, ef.tw);
+
+		/* split */
+		int tsplit = T(sub.patch, sub.Pw, Psplit);
+		ef0.tv = ef.tu;
+		ef0.tw = tsplit;
+
+		ef1.tv = tsplit;
+		ef1.tw = ef.tv;
+
+		/* create subpatches */
+		TriangleDice::SubPatch sub0 = {sub.patch, sub.Pw, Psplit, sub.Pv};
+		TriangleDice::SubPatch sub1 = {sub.patch, sub.Pw, sub.Pu, Psplit};
+
+		split(sub0, ef0, depth+1);
+		split(sub1, ef1, depth+1);
+	}
+	else
+		dispatch(sub, ef);
+}
+
+void DiagSplit::split(QuadDice::SubPatch& sub, QuadDice::EdgeFactors& ef, int depth)
+{
+	if((ef.tu0 == DSPLIT_NON_UNIFORM || ef.tu1 == DSPLIT_NON_UNIFORM)) {
+		/* partition edges */
+		QuadDice::EdgeFactors ef0, ef1;
+		float2 Pu0, Pu1;
+
+		partition_edge(sub.patch,
+			&Pu0, &ef0.tu0, &ef1.tu0, sub.P00, sub.P10, ef.tu0);
+		partition_edge(sub.patch,
+			&Pu1, &ef0.tu1, &ef1.tu1, sub.P01, sub.P11, ef.tu1);
+
+		/* split */
+		int tsplit = T(sub.patch, Pu0, Pu1);
+		ef0.tv0 = ef.tv0;
+		ef0.tv1 = tsplit;
+
+		ef1.tv0 = tsplit;
+		ef1.tv1 = ef.tv1;
+
+		/* create subpatches */
+		QuadDice::SubPatch sub0 = {sub.patch, sub.P00, Pu0, sub.P01, Pu1};
+		QuadDice::SubPatch sub1 = {sub.patch, Pu0, sub.P10, Pu1, sub.P11};
+
+		split(sub0, ef0, depth+1);
+		split(sub1, ef1, depth+1);
+	}
+	else if(ef.tv0 == DSPLIT_NON_UNIFORM || ef.tv1 == DSPLIT_NON_UNIFORM) {
+		/* partition edges */
+		QuadDice::EdgeFactors ef0, ef1;
+		float2 Pv0, Pv1;
+
+		partition_edge(sub.patch,
+			&Pv0, &ef0.tv0, &ef1.tv0, sub.P00, sub.P01, ef.tv0);
+		partition_edge(sub.patch,
+			&Pv1, &ef0.tv1, &ef1.tv1, sub.P10, sub.P11, ef.tv1);
+
+		/* split */
+		int tsplit = T(sub.patch, Pv0, Pv1);
+		ef0.tu0 = ef.tu0;
+		ef0.tu1 = tsplit;
+
+		ef1.tu0 = tsplit;
+		ef1.tu1 = ef.tu1;
+
+		/* create subpatches */
+		QuadDice::SubPatch sub0 = {sub.patch, sub.P00, sub.P10, Pv0, Pv1};
+		QuadDice::SubPatch sub1 = {sub.patch, Pv0, Pv1, sub.P01, sub.P11};
+
+		split(sub0, ef0, depth+1);
+		split(sub1, ef1, depth+1);
+	}
+	else
+		dispatch(sub, ef);
+}
+
+void DiagSplit::split_triangle(Mesh *mesh, Patch *patch, int shader, bool smooth)
+{
+	TriangleDice::SubPatch sub;
+	TriangleDice::EdgeFactors ef;
+
+	sub.patch = patch;
+	sub.Pu = make_float2(1.0f, 0.0f);
+	sub.Pv = make_float2(0.0f, 1.0f);
+	sub.Pw = make_float2(0.0f, 0.0f);
+
+	ef.tu = T(patch, sub.Pv, sub.Pw);
+	ef.tv = T(patch, sub.Pw, sub.Pu);
+	ef.tw = T(patch, sub.Pu, sub.Pv);
+
+	split(sub, ef);
+
+	TriangleDice dice(mesh, shader, smooth, dicing_rate);
+	dice.camera = camera;
+
+	for(size_t i = 0; i < subpatches_triangle.size(); i++) {
+		TriangleDice::SubPatch& sub = subpatches_triangle[i];
+		TriangleDice::EdgeFactors& ef = edgefactors_triangle[i];
+
+		ef.tu = 4;
+		ef.tv = 4;
+		ef.tw = 4;
+
+		ef.tu = max(ef.tu, 1);
+		ef.tv = max(ef.tv, 1);
+		ef.tw = max(ef.tw, 1);
+
+		dice.dice(sub, ef);
+	}
+
+	subpatches_triangle.clear();
+	edgefactors_triangle.clear();
+}
+
+void DiagSplit::split_quad(Mesh *mesh, Patch *patch, int shader, bool smooth)
+{
+	QuadDice::SubPatch sub;
+	QuadDice::EdgeFactors ef;
+
+	sub.patch = patch;
+	sub.P00 = make_float2(0.0f, 0.0f);
+	sub.P10 = make_float2(1.0f, 0.0f);
+	sub.P01 = make_float2(0.0f, 1.0f);
+	sub.P11 = make_float2(1.0f, 1.0f);
+
+	ef.tu0 = T(patch, sub.P00, sub.P10);
+	ef.tu1 = T(patch, sub.P01, sub.P11);
+	ef.tv0 = T(patch, sub.P00, sub.P01);
+	ef.tv1 = T(patch, sub.P10, sub.P11);
+
+	split(sub, ef);
+
+	QuadDice dice(mesh, shader, smooth, dicing_rate);
+	dice.camera = camera;
+
+	for(size_t i = 0; i < subpatches_quad.size(); i++) {
+		QuadDice::SubPatch& sub = subpatches_quad[i];
+		QuadDice::EdgeFactors& ef = edgefactors_quad[i];
+
+		ef.tu0 = max(ef.tu0, 1);
+		ef.tu1 = max(ef.tu1, 1);
+		ef.tv0 = max(ef.tv0, 1);
+		ef.tv1 = max(ef.tv1, 1);
+
+		dice.dice(sub, ef);
+	}
+
+	subpatches_quad.clear();
+	edgefactors_quad.clear();
+}
+
+CCL_NAMESPACE_END
+
diff --git a/intern/cycles/subd/subd_split.h b/intern/cycles/subd/subd_split.h
new file mode 100644
index 00000000000..8ec5d24266e
--- /dev/null
+++ b/intern/cycles/subd/subd_split.h
@@ -0,0 +1,71 @@
+/*
+ * Copyright 2011, Blender Foundation.
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#ifndef __SUBD_SPLIT_H__
+#define __SUBD_SPLIT_H__
+
+/* DiagSplit: Parallel, Crack-free, Adaptive Tessellation for Micropolygon Rendering
+ * Splits up patches and determines edge tesselation factors for dicing. Patch
+ * evaluation at arbitrary points is required for this to work. See the paper
+ * for more details. */
+
+#include "subd_dice.h"
+
+#include "util_types.h"
+#include "util_vector.h"
+
+CCL_NAMESPACE_BEGIN
+
+class Mesh;
+class Patch;
+
+#define DSPLIT_NON_UNIFORM -1
+
+class DiagSplit {
+public:
+	vector<QuadDice::SubPatch> subpatches_quad;
+	vector<QuadDice::EdgeFactors> edgefactors_quad;
+	vector<TriangleDice::SubPatch> subpatches_triangle;
+	vector<TriangleDice::EdgeFactors> edgefactors_triangle;
+
+	int test_steps;
+	int split_threshold;
+	float dicing_rate;
+	Camera *camera;
+
+	DiagSplit();
+
+	float3 project(Patch *patch, float2 uv);
+	int T(Patch *patch, float2 Pstart, float2 Pend);
+	void partition_edge(Patch *patch, float2 *P, int *t0, int *t1,
+		float2 Pstart, float2 Pend, int t);
+
+	void dispatch(QuadDice::SubPatch& sub, QuadDice::EdgeFactors& ef);
+	void split(QuadDice::SubPatch& sub, QuadDice::EdgeFactors& ef, int depth=0);
+
+	void dispatch(TriangleDice::SubPatch& sub, TriangleDice::EdgeFactors& ef);
+	void split(TriangleDice::SubPatch& sub, TriangleDice::EdgeFactors& ef, int depth=0);
+
+	void split_triangle(Mesh *mesh, Patch *patch, int shader, bool smooth);
+	void split_quad(Mesh *mesh, Patch *patch, int shader, bool smooth);
+};
+
+CCL_NAMESPACE_END
+
+#endif /* __SUBD_SPLIT_H__ */
+
diff --git a/intern/cycles/subd/subd_stencil.cpp b/intern/cycles/subd/subd_stencil.cpp
new file mode 100644
index 00000000000..599a745eebb
--- /dev/null
+++ b/intern/cycles/subd/subd_stencil.cpp
@@ -0,0 +1,103 @@
+/*
+ * Copyright 2011, Blender Foundation.
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include "subd_stencil.h"
+
+#include "util_debug.h"
+#include "util_math.h"
+
+CCL_NAMESPACE_BEGIN
+
+StencilMask::StencilMask()
+{
+}
+
+StencilMask::StencilMask(int size)
+{
+	/* initialize weights to zero. */
+	weights.resize(size, 0.0f);
+}
+
+void StencilMask::resize(int size)
+{
+	weights.resize(size, 0.0f);
+}
+
+StencilMask& StencilMask::operator=(float value)
+{
+	const int size = weights.size();
+	for(int i = 0; i < size; i++)
+		weights[i] = value;
+
+	return *this;
+}
+
+void StencilMask::operator+=(const StencilMask& mask)
+{
+	assert(mask.size() == size());
+	
+	const int size = weights.size();
+	for(int i = 0; i < size; i++)
+		weights[i] += mask.weights[i];
+}
+
+void StencilMask::operator-=(const StencilMask& mask)
+{
+	assert(mask.size() == size());
+	
+	const int size = weights.size();
+	for(int i = 0; i < size; i++)
+		weights[i] -= mask.weights[i];
+}
+
+void StencilMask::operator*=(float scale)
+{
+	const int size = weights.size();
+
+	for(int i = 0; i < size; i++)
+		weights[i] *= scale;
+}
+
+void StencilMask::operator/=(float scale)
+{
+	*this *= 1.0f/scale;
+}
+
+float StencilMask::sum() const
+{
+	float total = 0.0f;
+	const int size = weights.size();
+
+	for(int i = 0; i < size; i++)
+		total += weights[i];
+	
+	return total;
+}
+
+bool StencilMask::is_normalized() const
+{
+	return fabsf(sum() - 1.0f) < 0.0001f;
+}
+
+void StencilMask::normalize()
+{
+	*this /= sum();
+}
+
+CCL_NAMESPACE_END
+
diff --git a/intern/cycles/subd/subd_stencil.h b/intern/cycles/subd/subd_stencil.h
new file mode 100644
index 00000000000..e11d8f37cd3
--- /dev/null
+++ b/intern/cycles/subd/subd_stencil.h
@@ -0,0 +1,65 @@
+/*
+ * Copyright 2006, NVIDIA Corporation Ignacio Castano <icastano@nvidia.com>
+ *
+ * Modifications copyright (c) 2011, Blender Foundation.
+ * All rights reserved.
+ * 
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ * 
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ * 
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ */
+
+#ifndef __SUBD_STENCIL__
+#define __SUBD_STENCIL__
+
+#include "util_types.h"
+#include "util_vector.h"
+
+CCL_NAMESPACE_BEGIN
+
+class StencilMask
+{
+public:
+	StencilMask();
+	StencilMask(int size);
+
+	void resize(int size);
+
+	StencilMask& operator=(float value);
+
+	void operator+=(const StencilMask& mask);
+	void operator-=(const StencilMask& mask);
+	void operator*=(float scale);
+	void operator/=(float scale);
+
+	int size() const { return weights.size(); }
+
+	float operator[](int i) const { return weights[i]; }
+	float& operator[](int i) { return weights[i]; }
+
+	float sum() const;
+	bool is_normalized() const;
+	void normalize();
+
+private:
+	vector<float> weights;
+};
+
+CCL_NAMESPACE_END
+
+#endif /* __SUBD_STENCIL__ */
+
diff --git a/intern/cycles/subd/subd_vert.h b/intern/cycles/subd/subd_vert.h
new file mode 100644
index 00000000000..638019547ae
--- /dev/null
+++ b/intern/cycles/subd/subd_vert.h
@@ -0,0 +1,121 @@
+/*
+ * Original code in the public domain -- castanyo@yahoo.es
+ * 
+ * Modifications copyright (c) 2011, Blender Foundation.
+ * All rights reserved.
+ * 
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ * * Redistributions of source code must retain the above copyright
+ *   notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above copyright
+ *   notice, this list of conditions and the following disclaimer in the
+ *   documentation and/or other materials provided with the distribution.
+ * 
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+ * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+ * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef __SUBD_VERTEX_H__
+#define __SUBD_VERTEX_H__
+
+#include "subd_edge.h"
+#include "subd_mesh.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* Subd Vertex */
+
+class SubdVert
+{
+public:
+	int id;
+	float3 co;
+	
+	SubdEdge *edge;
+	SubdVert *next;
+	SubdVert *prev;
+
+	SubdVert(int id_)
+	{
+		id = id_;
+		edge = NULL;
+		co = make_float3(0.0f, 0.0f, 0.0f);
+
+		next = this;
+		prev = this;
+	}
+	
+	/* valence */
+	int valence()
+	{
+		int num = 0;
+
+		for(EdgeIterator it(edges()); !it.isDone(); it.advance())
+			num++;
+
+		return num;
+	}
+
+	/* edge queries */
+	bool is_boundary() { return (edge && !edge->face); }
+
+	/* iterator over edges in counterclockwise order */
+	class EdgeIterator
+	{
+	public:
+		EdgeIterator(SubdEdge *e) : end(NULL), cur(e) { }
+
+		virtual void advance()
+		{
+			if(end == NULL) end = cur;
+			cur = cur->pair->next;
+			//cur = cur->prev->pair;
+		}
+
+		virtual bool isDone() { return end == cur; }
+		virtual SubdEdge *current() { return cur; }
+
+	private:
+		SubdEdge *end;
+		SubdEdge *cur;
+	};
+
+	/* iterator over edges in clockwise order */
+	class ReverseEdgeIterator
+	{
+	public:
+		ReverseEdgeIterator(SubdEdge *e) : end(NULL), cur(e) { }
+
+		virtual void advance()
+		{
+			if(end == NULL) end = cur;
+			cur = cur->prev->pair;
+		}
+
+		virtual bool isDone() { return end == cur; }
+		virtual SubdEdge *current() { return cur; }
+
+	private:
+		SubdEdge *end;
+		SubdEdge *cur;
+	};
+
+	EdgeIterator edges() { return EdgeIterator(edge); }
+	EdgeIterator edges(SubdEdge *edge) { return EdgeIterator(edge); }
+};
+
+CCL_NAMESPACE_END
+
+#endif /* __SUBD_VERTEX_H__ */
+