Cycles: Stitching of subdivided and displaced meshes

This patch stitches the vertices along patch edges so that cracks can
no longer form when applying subdivision or displacement a mesh.

Subpatches are now formed in a way that ensures vertex indices along
subpatch edges are equal for adjacent subpatches. A mapping of vertices
along patch edges is built to preform stitching. Overall performance is
roughly the same, some gains were made in splitting, but some was lost
in stitching.

This fixes:
- T49049 (cracks between patches from material and uv seams)
- T49048 (discontinuous normals with true displacement)

Reviewers: sergey, brecht

Differential Revision: https://developer.blender.org/D3692

7 files changed, 958 insertions, 364 deletions

diff --git a/intern/cycles/subd/CMakeLists.txt b/intern/cycles/subd/CMakeLists.txt
index f5ceaa0436d..b874c5c3c2d 100644
--- a/intern/cycles/subd/CMakeLists.txt
+++ b/intern/cycles/subd/CMakeLists.txt
@@ -19,6 +19,7 @@ set(SRC_HEADERS
   subd_patch.h
   subd_patch_table.h
   subd_split.h
+  subd_subpatch.h
 )
 
 set(LIB
diff --git a/intern/cycles/subd/subd_dice.cpp b/intern/cycles/subd/subd_dice.cpp
index 914b408911e..714ff6ed629 100644
--- a/intern/cycles/subd/subd_dice.cpp
+++ b/intern/cycles/subd/subd_dice.cpp
@@ -38,107 +38,91 @@ EdgeDice::EdgeDice(const SubdParams &params_) : params(params_)
   }
 }
 
-void EdgeDice::reserve(int num_verts)
+void EdgeDice::reserve(int num_verts, int num_triangles)
 {
   Mesh *mesh = params.mesh;
 
   vert_offset = mesh->verts.size();
   tri_offset = mesh->num_triangles();
 
-  /* todo: optimize so we can reserve in advance, this is like push_back_slow() */
-  if (vert_offset + num_verts > mesh->verts.capacity()) {
-    mesh->reserve_mesh(size_t((vert_offset + num_verts) * 1.2), mesh->num_triangles());
-  }
-
-  mesh->resize_mesh(vert_offset + num_verts, tri_offset);
+  mesh->resize_mesh(mesh->verts.size() + num_verts, mesh->num_triangles());
+  mesh->reserve_mesh(mesh->verts.size() + num_verts, mesh->num_triangles() + num_triangles);
 
   Attribute *attr_vN = mesh->attributes.add(ATTR_STD_VERTEX_NORMAL);
 
-  mesh_P = mesh->verts.data();
-  mesh_N = attr_vN->data_float3();
+  mesh_P = mesh->verts.data() + vert_offset;
+  mesh_N = attr_vN->data_float3() + vert_offset;
+
+  params.mesh->num_subd_verts += num_verts;
 }
 
-int EdgeDice::add_vert(Patch *patch, float2 uv)
+void EdgeDice::set_vert(Patch *patch, int index, float2 uv)
 {
   float3 P, N;
 
   patch->eval(&P, NULL, NULL, &N, uv.x, uv.y);
 
-  assert(vert_offset < params.mesh->verts.size());
-
-  mesh_P[vert_offset] = P;
-  mesh_N[vert_offset] = N;
-  params.mesh->vert_patch_uv[vert_offset] = make_float2(uv.x, uv.y);
-
-  if (params.ptex) {
-    Attribute *attr_ptex_uv = params.mesh->attributes.add(ATTR_STD_PTEX_UV);
-    params.mesh->attributes.resize();
+  assert(index < params.mesh->verts.size());
 
-    float3 *ptex_uv = attr_ptex_uv->data_float3();
-    ptex_uv[vert_offset] = make_float3(uv.x, uv.y, 0.0f);
-  }
-
-  params.mesh->num_subd_verts++;
-
-  return vert_offset++;
+  mesh_P[index] = P;
+  mesh_N[index] = N;
+  params.mesh->vert_patch_uv[index + vert_offset] = make_float2(uv.x, uv.y);
 }
 
 void EdgeDice::add_triangle(Patch *patch, int v0, int v1, int v2)
 {
   Mesh *mesh = params.mesh;
 
-  /* todo: optimize so we can reserve in advance, this is like push_back_slow() */
-  if (mesh->triangles.size() == mesh->triangles.capacity())
-    mesh->reserve_mesh(mesh->verts.size(), size_t(max(mesh->num_triangles() + 1, 1) * 1.2));
-
-  mesh->add_triangle(v0, v1, v2, patch->shader, true);
+  mesh->add_triangle(v0 + vert_offset, v1 + vert_offset, v2 + vert_offset, patch->shader, true);
   params.mesh->triangle_patch[params.mesh->num_triangles() - 1] = patch->patch_index;
 
-  if (params.ptex) {
-    Attribute *attr_ptex_face_id = params.mesh->attributes.add(ATTR_STD_PTEX_FACE_ID);
-    params.mesh->attributes.resize();
-
-    float *ptex_face_id = attr_ptex_face_id->data_float();
-    ptex_face_id[tri_offset] = (float)patch->ptex_face_id();
-  }
-
   tri_offset++;
 }
 
-void EdgeDice::stitch_triangles(Patch *patch, vector<int> &outer, vector<int> &inner)
+void EdgeDice::stitch_triangles(Subpatch &sub, int edge)
 {
-  if (inner.size() == 0 || outer.size() == 0)
+  int Mu = max(sub.edge_u0.T, sub.edge_u1.T);
+  int Mv = max(sub.edge_v0.T, sub.edge_v1.T);
+  Mu = max(Mu, 2);
+  Mv = max(Mv, 2);
+
+  int outer_T = sub.edges[edge].T;
+  int inner_T = ((edge % 2) == 0) ? Mv - 2 : Mu - 2;
+
+  if (inner_T < 0 || outer_T < 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();) {
+  for (size_t i = 0, j = 0; i < inner_T || j < outer_T;) {
     int v0, v1, v2;
 
-    v0 = inner[i];
-    v1 = outer[j];
+    v0 = sub.get_vert_along_grid_edge(edge, i);
+    v1 = sub.get_vert_along_edge(edge, j);
 
-    if (j + 1 == outer.size()) {
-      v2 = inner[++i];
+    if (j == outer_T) {
+      v2 = sub.get_vert_along_grid_edge(edge, ++i);
     }
-    else if (i + 1 == inner.size()) {
-      v2 = outer[++j];
+    else if (i == inner_T) {
+      v2 = sub.get_vert_along_edge(edge, ++j);
     }
     else {
       /* length of diagonals */
-      float len1 = len_squared(mesh_P[inner[i]] - mesh_P[outer[j + 1]]);
-      float len2 = len_squared(mesh_P[outer[j]] - mesh_P[inner[i + 1]]);
+      float len1 = len_squared(mesh_P[sub.get_vert_along_grid_edge(edge, i)] -
+                               mesh_P[sub.get_vert_along_edge(edge, j + 1)]);
+      float len2 = len_squared(mesh_P[sub.get_vert_along_edge(edge, j)] -
+                               mesh_P[sub.get_vert_along_grid_edge(edge, i + 1)]);
 
       /* use smallest diagonal */
       if (len1 < len2)
-        v2 = outer[++j];
+        v2 = sub.get_vert_along_edge(edge, ++j);
       else
-        v2 = inner[++i];
+        v2 = sub.get_vert_along_grid_edge(edge, ++i);
     }
 
-    add_triangle(patch, v0, v1, v2);
+    add_triangle(sub.patch, v1, v0, v2);
   }
 }
 
@@ -148,22 +132,15 @@ QuadDice::QuadDice(const SubdParams &params_) : EdgeDice(params_)
 {
 }
 
-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);
-}
-
-float2 QuadDice::map_uv(SubPatch &sub, float u, float v)
+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);
+  float2 d0 = interp(sub.c00, sub.c01, v);
+  float2 d1 = interp(sub.c10, sub.c11, v);
   return interp(d0, d1, u);
 }
 
-float3 QuadDice::eval_projected(SubPatch &sub, float u, float v)
+float3 QuadDice::eval_projected(Subpatch &sub, float u, float v)
 {
   float2 uv = map_uv(sub, u, v);
   float3 P;
@@ -175,70 +152,40 @@ float3 QuadDice::eval_projected(SubPatch &sub, float u, float v)
   return P;
 }
 
-int QuadDice::add_vert(SubPatch &sub, float u, float v)
+void QuadDice::set_vert(Subpatch &sub, int index, float u, float v)
 {
-  return EdgeDice::add_vert(sub.patch, map_uv(sub, u, v));
+  EdgeDice::set_vert(sub.patch, index, 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)
+void QuadDice::set_side(Subpatch &sub, int edge)
 {
-  outer.clear();
-  inner.clear();
+  int t = sub.edges[edge].T;
 
   /* 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));
+  for (int i = 0; i < t; i++) {
+    float f = i / (float)t;
+
+    float u, v;
+    switch (edge) {
+      case 0:
+        u = 0;
+        v = f;
+        break;
+      case 1:
+        u = f;
+        v = 1;
+        break;
+      case 2:
+        u = 1;
+        v = 1.0f - f;
+        break;
+      case 3:
+        u = 1.0f - f;
+        v = 0;
+        break;
+    }
 
-  /* 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));
+    set_vert(sub, sub.get_vert_along_edge(edge, i), u, v);
   }
 }
 
@@ -247,7 +194,7 @@ float QuadDice::quad_area(const float3 &a, const float3 &b, const float3 &c, con
   return triangle_area(a, b, d) + triangle_area(a, d, c);
 }
 
-float QuadDice::scale_factor(SubPatch &sub, EdgeFactors &ef, int Mu, int Mv)
+float QuadDice::scale_factor(Subpatch &sub, int Mu, int Mv)
 {
   /* estimate area as 4x largest of 4 quads */
   float3 P[3][3];
@@ -269,23 +216,14 @@ float QuadDice::scale_factor(SubPatch &sub, EdgeFactors &ef, int Mu, int Mv)
   // 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 N = 0.5f * (Ntris - (sub.edge_u0.T + sub.edge_u1.T + sub.edge_v0.T + sub.edge_v1.T));
   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 */
-  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)
+void QuadDice::add_grid(Subpatch &sub, int Mu, int Mv, int offset)
 {
   /* create inner grid */
   float du = 1.0f / (float)Mu;
@@ -296,13 +234,13 @@ void QuadDice::add_grid(SubPatch &sub, int Mu, int Mv, int offset)
       float u = i * du;
       float v = j * dv;
 
-      add_vert(sub, u, v);
+      set_vert(sub, offset + (i - 1) + (j - 1) * (Mu - 1), 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);
+        int i1 = offset + (i - 1) + (j - 1) * (Mu - 1);
+        int i2 = offset + i + (j - 1) * (Mu - 1);
+        int i3 = offset + i + j * (Mu - 1);
+        int i4 = offset + (i - 1) + j * (Mu - 1);
 
         add_triangle(sub.patch, i1, i2, i3);
         add_triangle(sub.patch, i1, i3, i4);
@@ -311,11 +249,11 @@ void QuadDice::add_grid(SubPatch &sub, int Mu, int Mv, int offset)
   }
 }
 
-void QuadDice::dice(SubPatch &sub, EdgeFactors &ef)
+void QuadDice::dice(Subpatch &sub)
 {
   /* compute inner grid size with scale factor */
-  int Mu = max(ef.tu0, ef.tu1);
-  int Mv = max(ef.tv0, ef.tv1);
+  int Mu = max(sub.edge_u0.T, sub.edge_u1.T);
+  int Mv = max(sub.edge_v0.T, sub.edge_v1.T);
 
 #if 0 /* Doesn't work very well, especially at grazing angles. */
   float S = scale_factor(sub, ef, Mu, Mv);
@@ -326,33 +264,19 @@ void QuadDice::dice(SubPatch &sub, EdgeFactors &ef)
   Mu = max((int)ceilf(S * Mu), 2);  // XXX handle 0 & 1?
   Mv = max((int)ceilf(S * Mv), 2);  // XXX handle 0 & 1?
 
-  /* reserve space for new verts */
-  int offset = params.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(sub.patch, outer, inner);
-
-  /* top side */
-  add_side_u(sub, outer, inner, Mu, Mv, ef.tu1, 1, offset);
-  stitch_triangles(sub.patch, inner, outer);
-
-  /* left side */
-  add_side_v(sub, outer, inner, Mu, Mv, ef.tv0, 0, offset);
-  stitch_triangles(sub.patch, inner, outer);
+  /* inner grid */
+  add_grid(sub, Mu, Mv, sub.inner_grid_vert_offset);
 
-  /* right side */
-  add_side_v(sub, outer, inner, Mu, Mv, ef.tv1, 1, offset);
-  stitch_triangles(sub.patch, outer, inner);
+  /* sides */
+  set_side(sub, 0);
+  set_side(sub, 1);
+  set_side(sub, 2);
+  set_side(sub, 3);
 
-  assert(vert_offset == params.mesh->verts.size());
+  stitch_triangles(sub, 0);
+  stitch_triangles(sub, 1);
+  stitch_triangles(sub, 2);
+  stitch_triangles(sub, 3);
 }
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/subd/subd_dice.h b/intern/cycles/subd/subd_dice.h
index eee54e01861..ee63403d40c 100644
--- a/intern/cycles/subd/subd_dice.h
+++ b/intern/cycles/subd/subd_dice.h
@@ -25,6 +25,8 @@
 #include "util/util_types.h"
 #include "util/util_vector.h"
 
+#include "subd/subd_subpatch.h"
+
 CCL_NAMESPACE_BEGIN
 
 class Camera;
@@ -67,78 +69,33 @@ class EdgeDice {
 
   explicit EdgeDice(const SubdParams &params);
 
-  void reserve(int num_verts);
+  void reserve(int num_verts, int num_triangles);
 
-  int add_vert(Patch *patch, float2 uv);
+  void set_vert(Patch *patch, int index, float2 uv);
   void add_triangle(Patch *patch, int v0, int v1, int v2);
 
-  void stitch_triangles(Patch *patch, vector<int> &outer, vector<int> &inner);
+  void stitch_triangles(Subpatch &sub, int edge);
 };
 
-/* Quad EdgeDice
- *
- * Edge tessellation factors and subpatch coordinates are as follows:
- *
- *            tu1
- *     P01 --------- P11
- *     |               |
- * tv0 |               | tv1
- *     |               |
- *     P00 --------- P10
- *            tu0
- */
+/* Quad EdgeDice */
 
 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;
-  };
-
   explicit QuadDice(const SubdParams &params);
 
-  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);
+  float3 eval_projected(Subpatch &sub, float u, float v);
+
+  float2 map_uv(Subpatch &sub, float u, float v);
+  void set_vert(Subpatch &sub, int index, float u, float v);
+
+  void add_grid(Subpatch &sub, int Mu, int Mv, int offset);
+
+  void set_side(Subpatch &sub, int edge);
 
   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);
+  float scale_factor(Subpatch &sub, int Mu, int Mv);
 
-  void dice(SubPatch &sub, EdgeFactors &ef);
+  void dice(Subpatch &sub);
 };
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/subd/subd_patch.h b/intern/cycles/subd/subd_patch.h
index 5209d4d0b07..2e3484f95e3 100644
--- a/intern/cycles/subd/subd_patch.h
+++ b/intern/cycles/subd/subd_patch.h
@@ -24,18 +24,13 @@ CCL_NAMESPACE_BEGIN
 
 class Patch {
  public:
-  virtual ~Patch()
-  {
-  }
+  virtual ~Patch() = default;
+
   virtual void eval(float3 *P, float3 *dPdu, float3 *dPdv, float3 *N, float u, float v) = 0;
-  virtual BoundBox bound() = 0;
-  virtual int ptex_face_id()
-  {
-    return -1;
-  }
 
   int patch_index;
   int shader;
+  bool from_ngon;
 };
 
 /* Linear Quad Patch */
diff --git a/intern/cycles/subd/subd_split.cpp b/intern/cycles/subd/subd_split.cpp
index e5b85fcfd60..593bb92b0ff 100644
--- a/intern/cycles/subd/subd_split.cpp
+++ b/intern/cycles/subd/subd_split.cpp
@@ -21,6 +21,9 @@
 #include "subd/subd_patch.h"
 #include "subd/subd_split.h"
 
+#include "util/util_algorithm.h"
+#include "util/util_foreach.h"
+#include "util/util_hash.h"
 #include "util/util_math.h"
 #include "util/util_types.h"
 
@@ -28,14 +31,12 @@ CCL_NAMESPACE_BEGIN
 
 /* DiagSplit */
 
-DiagSplit::DiagSplit(const SubdParams &params_) : params(params_)
-{
-}
+#define DSPLIT_NON_UNIFORM -1
+#define STITCH_NGON_CENTER_VERT_INDEX_OFFSET 0x60000000
+#define STITCH_NGON_SPLIT_EDGE_CENTER_VERT_TAG (0x60000000 - 1)
 
-void DiagSplit::dispatch(QuadDice::SubPatch &sub, QuadDice::EdgeFactors &ef)
+DiagSplit::DiagSplit(const SubdParams &params_) : params(params_)
 {
-  subpatches_quad.push_back(sub);
-  edgefactors_quad.push_back(ef);
 }
 
 float3 DiagSplit::to_world(Patch *patch, float2 uv)
@@ -49,45 +50,62 @@ float3 DiagSplit::to_world(Patch *patch, float2 uv)
   return P;
 }
 
-int DiagSplit::T(Patch *patch, float2 Pstart, float2 Pend)
+static void order_float2(float2 &a, float2 &b)
 {
-  float3 Plast = make_float3(0.0f, 0.0f, 0.0f);
+  if (b.x < a.x || b.y < a.y) {
+    swap(a, b);
+  }
+}
+
+int DiagSplit::T(Patch *patch, float2 Pstart, float2 Pend, bool recursive_resolve)
+{
+  order_float2(Pstart, Pend); /* May not be necessary, but better to be safe. */
+
   float Lsum = 0.0f;
   float Lmax = 0.0f;
 
-  for (int i = 0; i < params.test_steps; i++) {
+  float3 Plast = to_world(patch, Pstart);
+
+  for (int i = 1; i < params.test_steps; i++) {
     float t = i / (float)(params.test_steps - 1);
 
     float3 P = to_world(patch, Pstart + t * (Pend - Pstart));
 
-    if (i > 0) {
-      float L;
+    float L;
 
-      if (!params.camera) {
-        L = len(P - Plast);
-      }
-      else {
-        Camera *cam = params.camera;
-
-        float pixel_width = cam->world_to_raster_size((P + Plast) * 0.5f);
-        L = len(P - Plast) / pixel_width;
-      }
+    if (!params.camera) {
+      L = len(P - Plast);
+    }
+    else {
+      Camera *cam = params.camera;
 
-      Lsum += L;
-      Lmax = max(L, Lmax);
+      float pixel_width = cam->world_to_raster_size((P + Plast) * 0.5f);
+      L = len(P - Plast) / pixel_width;
     }
 
+    Lsum += L;
+    Lmax = max(L, Lmax);
+
     Plast = P;
   }
 
   int tmin = (int)ceilf(Lsum / params.dicing_rate);
   int tmax = (int)ceilf((params.test_steps - 1) * Lmax /
-                        params.dicing_rate);  // XXX paper says N instead of N-1, seems wrong?
+                       params.dicing_rate);  // XXX paper says N instead of N-1, seems wrong?
+  int res = max(tmax, 1);
 
-  if (tmax - tmin > params.split_threshold)
-    return DSPLIT_NON_UNIFORM;
+  if (tmax - tmin > params.split_threshold) {
+    if (!recursive_resolve) {
+      res = DSPLIT_NON_UNIFORM;
+    }
+    else {
+      float2 P = (Pstart + Pend) * 0.5f;
+      res = T(patch, Pstart, P, true) + T(patch, P, Pend, true);
+    }
+  }
 
-  return tmax;
+  limit_edge_factor(res, patch, Pstart, Pend);
+  return res;
 }
 
 void DiagSplit::partition_edge(
@@ -99,159 +117,632 @@ void DiagSplit::partition_edge(
     *t1 = T(patch, *P, Pend);
   }
   else {
+    assert(t >= 2); /* Need at least two segments to partition into. */
+
     int I = (int)floorf((float)t * 0.5f);
-    *P = interp(Pstart, Pend, (t == 0) ? 0 : I / (float)t); /* XXX is t faces or verts */
+    *P = interp(Pstart, Pend, I / (float)t);
     *t0 = I;
     *t1 = t - I;
   }
 }
 
-static void limit_edge_factors(const QuadDice::SubPatch &sub, QuadDice::EdgeFactors &ef, int max_t)
+void DiagSplit::limit_edge_factor(int &T, Patch *patch, float2 Pstart, float2 Pend)
+{
+  int max_t = 1 << params.max_level;
+  int max_t_for_edge = int(max_t * len(Pstart - Pend));
+
+  if (patch->from_ngon) {
+    max_t_for_edge >>= 1; /* Initial split of ngon causes edges to extend half the distance. */
+  }
+
+  T = (max_t_for_edge <= 1) ? 1 : min(T, max_t_for_edge);
+
+  assert(T >= 1 || T == DSPLIT_NON_UNIFORM);
+}
+
+void DiagSplit::resolve_edge_factors(Subpatch &sub)
 {
-  float2 P00 = sub.P00;
-  float2 P01 = sub.P01;
-  float2 P10 = sub.P10;
-  float2 P11 = sub.P11;
-
-  int tu0 = int(max_t * len(P10 - P00));
-  int tu1 = int(max_t * len(P11 - P01));
-  int tv0 = int(max_t * len(P01 - P00));
-  int tv1 = int(max_t * len(P11 - P10));
-
-  ef.tu0 = tu0 <= 1 ? 1 : min(ef.tu0, tu0);
-  ef.tu1 = tu1 <= 1 ? 1 : min(ef.tu1, tu1);
-  ef.tv0 = tv0 <= 1 ? 1 : min(ef.tv0, tv0);
-  ef.tv1 = tv1 <= 1 ? 1 : min(ef.tv1, tv1);
+  /* Resolve DSPLIT_NON_UNIFORM to actual T value if splitting is no longer possible. */
+  if (sub.edge_u0.T == 1 && sub.edge_u1.T == DSPLIT_NON_UNIFORM) {
+    sub.edge_u1.T = T(sub.patch, sub.c01, sub.c11, true);
+  }
+  if (sub.edge_u1.T == 1 && sub.edge_u0.T == DSPLIT_NON_UNIFORM) {
+    sub.edge_u0.T = T(sub.patch, sub.c00, sub.c10, true);
+  }
+  if (sub.edge_v0.T == 1 && sub.edge_v1.T == DSPLIT_NON_UNIFORM) {
+    sub.edge_v1.T = T(sub.patch, sub.c11, sub.c10, true);
+  }
+  if (sub.edge_v1.T == 1 && sub.edge_v0.T == DSPLIT_NON_UNIFORM) {
+    sub.edge_v0.T = T(sub.patch, sub.c01, sub.c00, true);
+  }
 }
 
-void DiagSplit::split(QuadDice::SubPatch &sub, QuadDice::EdgeFactors &ef, int depth)
+void DiagSplit::split(Subpatch &sub, int depth)
 {
   if (depth > 32) {
     /* We should never get here, but just in case end recursion safely. */
-    ef.tu0 = 1;
-    ef.tu1 = 1;
-    ef.tv0 = 1;
-    ef.tv1 = 1;
+    assert(!"diagsplit recursion limit reached");
+
+    sub.edge_u0.T = 1;
+    sub.edge_u1.T = 1;
+    sub.edge_v0.T = 1;
+    sub.edge_v1.T = 1;
 
-    dispatch(sub, ef);
+    subpatches.push_back(sub);
     return;
   }
 
-  bool split_u = (ef.tu0 == DSPLIT_NON_UNIFORM || ef.tu1 == DSPLIT_NON_UNIFORM);
-  bool split_v = (ef.tv0 == DSPLIT_NON_UNIFORM || ef.tv1 == DSPLIT_NON_UNIFORM);
+  bool split_u = (sub.edge_u0.T == DSPLIT_NON_UNIFORM || sub.edge_u1.T == DSPLIT_NON_UNIFORM);
+  bool split_v = (sub.edge_v0.T == DSPLIT_NON_UNIFORM || sub.edge_v1.T == DSPLIT_NON_UNIFORM);
 
   /* Split subpatches such that the ratio of T for opposite edges doesn't
    * exceed 1.5, this reduces over tessellation for some patches
    */
-  bool tmp_split_v = split_v;
-  if (!split_u && min(ef.tu0, ef.tu1) > 8 && min(ef.tu0, ef.tu1) * 1.5f < max(ef.tu0, ef.tu1))
+  /* clang-format off */
+  if (min(sub.edge_u0.T, sub.edge_u1.T) > 8 && /* Must be uniform and preferably greater than 8 to split. */
+      min(sub.edge_v0.T, sub.edge_v1.T) >= 2 && /* Must be uniform and at least 2 to split. */
+      max(sub.edge_u0.T, sub.edge_u1.T) / min(sub.edge_u0.T, sub.edge_u1.T) > 1.5f)
+  {
     split_v = true;
-  if (!tmp_split_v && min(ef.tu0, ef.tu1) > 8 && min(ef.tv0, ef.tv1) * 1.5f < max(ef.tv0, ef.tv1))
+  }
+  if (min(sub.edge_v0.T, sub.edge_v1.T) > 8 &&
+      min(sub.edge_u0.T, sub.edge_u1.T) >= 2 &&
+      max(sub.edge_v0.T, sub.edge_v1.T) / min(sub.edge_v0.T, sub.edge_v1.T) > 1.5f)
+  {
     split_u = true;
+  }
+  /* clang-format on */
 
-  /* alternate axis */
+  /* Alternate axis. */
   if (split_u && split_v) {
     split_u = depth % 2;
   }
 
-  if (split_u) {
-    /* partition edges */
-    QuadDice::EdgeFactors ef0, ef1;
-    float2 Pu0, Pu1;
+  if (!split_u && !split_v) {
+    /* Add the unsplit subpatch. */
+    subpatches.push_back(sub);
+    Subpatch &subpatch = subpatches[subpatches.size() - 1];
+
+    /* Update T values and offsets. */
+    for (int i = 0; i < 4; i++) {
+      Subpatch::edge_t &edge = subpatch.edges[i];
+
+      edge.offset = edge.edge->T;
+      edge.edge->T += edge.T;
+    }
+  }
+  else {
+    /* Copy into new subpatches. */
+    Subpatch sub_a = sub;
+    Subpatch sub_b = sub;
+
+    /* Pointers to various subpatch elements. */
+    Subpatch::edge_t *sub_across_0, *sub_across_1;
+    Subpatch::edge_t *sub_a_across_0, *sub_a_across_1;
+    Subpatch::edge_t *sub_b_across_0, *sub_b_across_1;
+
+    Subpatch::edge_t *sub_a_split, *sub_b_split;
+
+    float2 *Pa, *Pb, *Pc, *Pd;
+
+    /* Set pointers based on split axis. */
+    if (split_u) {
+      sub_across_0 = &sub.edge_u0;
+      sub_across_1 = &sub.edge_u1;
+      sub_a_across_0 = &sub_a.edge_u0;
+      sub_a_across_1 = &sub_a.edge_u1;
+      sub_b_across_0 = &sub_b.edge_u0;
+      sub_b_across_1 = &sub_b.edge_u1;
+
+      sub_a_split = &sub_a.edge_v1;
+      sub_b_split = &sub_b.edge_v0;
+
+      Pa = &sub_a.c11;
+      Pb = &sub_a.c10;
+      Pc = &sub_b.c01;
+      Pd = &sub_b.c00;
+    }
+    else {
+      sub_across_0 = &sub.edge_v0;
+      sub_across_1 = &sub.edge_v1;
+      sub_a_across_0 = &sub_a.edge_v0;
+      sub_a_across_1 = &sub_a.edge_v1;
+      sub_b_across_0 = &sub_b.edge_v0;
+      sub_b_across_1 = &sub_b.edge_v1;
+
+      sub_a_split = &sub_a.edge_u0;
+      sub_b_split = &sub_b.edge_u1;
+
+      Pa = &sub_a.c10;
+      Pb = &sub_a.c00;
+      Pc = &sub_b.c11;
+      Pd = &sub_b.c01;
+    }
 
-    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);
+    /* Partition edges */
+    float2 P0, P1;
 
-    /* split */
-    int tsplit = T(sub.patch, Pu0, Pu1);
-    ef0.tv0 = ef.tv0;
-    ef0.tv1 = tsplit;
+    partition_edge(
+        sub.patch, &P0, &sub_a_across_0->T, &sub_b_across_0->T, *Pd, *Pb, sub_across_0->T);
+    partition_edge(
+        sub.patch, &P1, &sub_a_across_1->T, &sub_b_across_1->T, *Pc, *Pa, sub_across_1->T);
 
-    ef1.tv0 = tsplit;
-    ef1.tv1 = ef.tv1;
+    /* Split */
+    *Pa = P1;
+    *Pb = P0;
 
-    /* create subpatches */
-    QuadDice::SubPatch sub0 = {sub.patch, sub.P00, Pu0, sub.P01, Pu1};
-    QuadDice::SubPatch sub1 = {sub.patch, Pu0, sub.P10, Pu1, sub.P11};
+    *Pc = P1;
+    *Pd = P0;
 
-    limit_edge_factors(sub0, ef0, 1 << params.max_level);
-    limit_edge_factors(sub1, ef1, 1 << params.max_level);
+    int tsplit = T(sub.patch, P0, P1);
 
-    split(sub0, ef0, depth + 1);
-    split(sub1, ef1, depth + 1);
+    if (depth == -2 && tsplit == 1) {
+      tsplit = 2; /* Ensure we can always split at depth -1. */
+    }
+
+    sub_a_split->T = tsplit;
+    sub_b_split->T = tsplit;
+
+    resolve_edge_factors(sub_a);
+    resolve_edge_factors(sub_b);
+
+    /* Create new edge */
+    Edge &edge = *alloc_edge();
+
+    sub_a_split->edge = &edge;
+    sub_b_split->edge = &edge;
+
+    sub_a_split->offset = 0;
+    sub_b_split->offset = 0;
+
+    sub_a_split->indices_decrease_along_edge = false;
+    sub_b_split->indices_decrease_along_edge = true;
+
+    sub_a_split->sub_edges_created_in_reverse_order = !split_u;
+    sub_b_split->sub_edges_created_in_reverse_order = !split_u;
+
+    edge.top_indices_decrease = sub_across_1->sub_edges_created_in_reverse_order;
+    edge.bottom_indices_decrease = sub_across_0->sub_edges_created_in_reverse_order;
+
+    /* Recurse */
+    edge.T = 0;
+    split(sub_a, depth + 1);
+
+    int edge_t = edge.T;
+    (void)edge_t;
+
+    edge.top_offset = sub_across_1->edge->T;
+    edge.bottom_offset = sub_across_0->edge->T;
+
+    edge.T = 0; /* We calculate T twice along each edge. :/ */
+    split(sub_b, depth + 1);
+
+    assert(edge.T == edge_t); /* If this fails we will crash at some later point! */
+
+    edge.top = sub_across_1->edge;
+    edge.bottom = sub_across_0->edge;
   }
-  else if (split_v) {
-    /* partition edges */
-    QuadDice::EdgeFactors ef0, ef1;
-    float2 Pv0, Pv1;
+}
+
+int DiagSplit::alloc_verts(int n)
+{
+  int a = num_alloced_verts;
+  num_alloced_verts += n;
+  return a;
+}
 
-    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);
+Edge *DiagSplit::alloc_edge()
+{
+  edges.emplace_back();
+  return &edges.back();
+}
 
-    /* split */
-    int tsplit = T(sub.patch, Pv0, Pv1);
-    ef0.tu0 = ef.tu0;
-    ef0.tu1 = tsplit;
+void DiagSplit::split_patches(Patch *patches, size_t patches_byte_stride)
+{
+  int patch_index = 0;
 
-    ef1.tu0 = tsplit;
-    ef1.tu1 = ef.tu1;
+  for (int f = 0; f < params.mesh->subd_faces.size(); f++) {
+    Mesh::SubdFace &face = params.mesh->subd_faces[f];
 
-    /* create subpatches */
-    QuadDice::SubPatch sub0 = {sub.patch, sub.P00, sub.P10, Pv0, Pv1};
-    QuadDice::SubPatch sub1 = {sub.patch, Pv0, Pv1, sub.P01, sub.P11};
+    Patch *patch = (Patch *)(((char *)patches) + patch_index * patches_byte_stride);
 
-    limit_edge_factors(sub0, ef0, 1 << params.max_level);
-    limit_edge_factors(sub1, ef1, 1 << params.max_level);
+    if (face.is_quad()) {
+      patch_index++;
 
-    split(sub0, ef0, depth + 1);
-    split(sub1, ef1, depth + 1);
+      split_quad(face, patch);
+    }
+    else {
+      patch_index += face.num_corners;
+
+      split_ngon(face, patch, patches_byte_stride);
+    }
   }
-  else {
-    dispatch(sub, ef);
+
+  params.mesh->vert_to_stitching_key_map.clear();
+  params.mesh->vert_stitching_map.clear();
+
+  post_split();
+}
+
+static Edge *create_edge_from_corner(DiagSplit *split,
+                                     const Mesh *mesh,
+                                     const Mesh::SubdFace &face,
+                                     int corner,
+                                     bool &reversed,
+                                     int v0,
+                                     int v1)
+{
+  int a = mesh->subd_face_corners[face.start_corner + mod(corner + 0, face.num_corners)];
+  int b = mesh->subd_face_corners[face.start_corner + mod(corner + 1, face.num_corners)];
+
+  reversed = !(b < a);
+
+  if (b < a) {
+    swap(a, b);
+    swap(v0, v1);
   }
+
+  Edge *edge = split->alloc_edge();
+
+  edge->is_stitch_edge = true;
+  edge->stitch_start_vert_index = a;
+  edge->stitch_end_vert_index = b;
+
+  edge->start_vert_index = v0;
+  edge->end_vert_index = v1;
+
+  edge->stitch_edge_key = {a, b};
+
+  return edge;
 }
 
-void DiagSplit::split_quad(Patch *patch, QuadDice::SubPatch *subpatch)
+void DiagSplit::split_quad(const Mesh::SubdFace &face, Patch *patch)
 {
-  QuadDice::SubPatch sub_split;
-  QuadDice::EdgeFactors ef_split;
+  Subpatch subpatch(patch);
+
+  int v = alloc_verts(4);
+
+  bool v0_reversed, u1_reversed, v1_reversed, u0_reversed;
+  subpatch.edge_v0.edge = create_edge_from_corner(
+      this, params.mesh, face, 3, v0_reversed, v + 3, v + 0);
+  subpatch.edge_u1.edge = create_edge_from_corner(
+      this, params.mesh, face, 2, u1_reversed, v + 2, v + 3);
+  subpatch.edge_v1.edge = create_edge_from_corner(
+      this, params.mesh, face, 1, v1_reversed, v + 1, v + 2);
+  subpatch.edge_u0.edge = create_edge_from_corner(
+      this, params.mesh, face, 0, u0_reversed, v + 0, v + 1);
+
+  subpatch.edge_v0.sub_edges_created_in_reverse_order = !v0_reversed;
+  subpatch.edge_u1.sub_edges_created_in_reverse_order = u1_reversed;
+  subpatch.edge_v1.sub_edges_created_in_reverse_order = v1_reversed;
+  subpatch.edge_u0.sub_edges_created_in_reverse_order = !u0_reversed;
+
+  subpatch.edge_v0.indices_decrease_along_edge = v0_reversed;
+  subpatch.edge_u1.indices_decrease_along_edge = u1_reversed;
+  subpatch.edge_v1.indices_decrease_along_edge = v1_reversed;
+  subpatch.edge_u0.indices_decrease_along_edge = u0_reversed;
+
+  /* Forces a split in both axis for quads, needed to match split of ngons into quads. */
+  subpatch.edge_u0.T = DSPLIT_NON_UNIFORM;
+  subpatch.edge_u1.T = DSPLIT_NON_UNIFORM;
+  subpatch.edge_v0.T = DSPLIT_NON_UNIFORM;
+  subpatch.edge_v1.T = DSPLIT_NON_UNIFORM;
+
+  split(subpatch, -2);
+}
+
+static Edge *create_split_edge_from_corner(DiagSplit *split,
+                                           const Mesh *mesh,
+                                           const Mesh::SubdFace &face,
+                                           int corner,
+                                           int side,
+                                           bool &reversed,
+                                           int v0,
+                                           int v1,
+                                           int vc)
+{
+  Edge *edge = split->alloc_edge();
+
+  int a = mesh->subd_face_corners[face.start_corner + mod(corner + 0, face.num_corners)];
+  int b = mesh->subd_face_corners[face.start_corner + mod(corner + 1, face.num_corners)];
 
-  if (subpatch) {
-    sub_split = *subpatch;
+  if (b < a) {
+    edge->stitch_edge_key = {b, a};
   }
   else {
-    sub_split.patch = patch;
-    sub_split.P00 = make_float2(0.0f, 0.0f);
-    sub_split.P10 = make_float2(1.0f, 0.0f);
-    sub_split.P01 = make_float2(0.0f, 1.0f);
-    sub_split.P11 = make_float2(1.0f, 1.0f);
+    edge->stitch_edge_key = {a, b};
   }
 
-  ef_split.tu0 = T(patch, sub_split.P00, sub_split.P10);
-  ef_split.tu1 = T(patch, sub_split.P01, sub_split.P11);
-  ef_split.tv0 = T(patch, sub_split.P00, sub_split.P01);
-  ef_split.tv1 = T(patch, sub_split.P10, sub_split.P11);
+  reversed = !(b < a);
 
-  limit_edge_factors(sub_split, ef_split, 1 << params.max_level);
+  if (side == 0) {
+    a = vc;
+  }
+  else {
+    b = vc;
+  }
+
+  if (!reversed) {
+    swap(a, b);
+    swap(v0, v1);
+  }
+
+  edge->is_stitch_edge = true;
+  edge->stitch_start_vert_index = a;
+  edge->stitch_end_vert_index = b;
+
+  edge->start_vert_index = v0;
+  edge->end_vert_index = v1;
+
+  return edge;
+}
+
+void DiagSplit::split_ngon(const Mesh::SubdFace &face, Patch *patches, size_t patches_byte_stride)
+{
+  Edge *prev_edge_u0 = nullptr;
+  Edge *first_edge_v0 = nullptr;
+
+  for (int corner = 0; corner < face.num_corners; corner++) {
+    Patch *patch = (Patch *)(((char *)patches) + corner * patches_byte_stride);
+
+    Subpatch subpatch(patch);
+
+    int v = alloc_verts(4);
+
+    /* Setup edges. */
+    Edge *edge_u1 = alloc_edge();
+    Edge *edge_v1 = alloc_edge();
+
+    edge_v1->is_stitch_edge = true;
+    edge_u1->is_stitch_edge = true;
+
+    edge_u1->stitch_start_vert_index = -(face.start_corner + mod(corner + 0, face.num_corners)) -
+                                       1;
+    edge_u1->stitch_end_vert_index = STITCH_NGON_CENTER_VERT_INDEX_OFFSET + face.ptex_offset;
+
+    edge_u1->start_vert_index = v + 3;
+    edge_u1->end_vert_index = v + 2;
+
+    edge_u1->stitch_edge_key = {edge_u1->stitch_start_vert_index, edge_u1->stitch_end_vert_index};
+
+    edge_v1->stitch_start_vert_index = -(face.start_corner + mod(corner + 1, face.num_corners)) -
+                                       1;
+    edge_v1->stitch_end_vert_index = STITCH_NGON_CENTER_VERT_INDEX_OFFSET + face.ptex_offset;
+
+    edge_v1->start_vert_index = v + 1;
+    edge_v1->end_vert_index = v + 2;
+
+    edge_v1->stitch_edge_key = {edge_v1->stitch_start_vert_index, edge_v1->stitch_end_vert_index};
+
+    bool v0_reversed, u0_reversed;
+
+    subpatch.edge_v0.edge = create_split_edge_from_corner(this,
+                                                          params.mesh,
+                                                          face,
+                                                          corner - 1,
+                                                          0,
+                                                          v0_reversed,
+                                                          v + 3,
+                                                          v + 0,
+                                                          STITCH_NGON_SPLIT_EDGE_CENTER_VERT_TAG);
+
+    subpatch.edge_u1.edge = edge_u1;
+    subpatch.edge_v1.edge = edge_v1;
+
+    subpatch.edge_u0.edge = create_split_edge_from_corner(this,
+                                                          params.mesh,
+                                                          face,
+                                                          corner + 0,
+                                                          1,
+                                                          u0_reversed,
+                                                          v + 0,
+                                                          v + 1,
+                                                          STITCH_NGON_SPLIT_EDGE_CENTER_VERT_TAG);
+
+    subpatch.edge_v0.sub_edges_created_in_reverse_order = !v0_reversed;
+    subpatch.edge_u1.sub_edges_created_in_reverse_order = false;
+    subpatch.edge_v1.sub_edges_created_in_reverse_order = true;
+    subpatch.edge_u0.sub_edges_created_in_reverse_order = !u0_reversed;
+
+    subpatch.edge_v0.indices_decrease_along_edge = v0_reversed;
+    subpatch.edge_u1.indices_decrease_along_edge = false;
+    subpatch.edge_v1.indices_decrease_along_edge = true;
+    subpatch.edge_u0.indices_decrease_along_edge = u0_reversed;
+
+    /* Perfrom split. */
+    {
+      subpatch.edge_u0.T = T(subpatch.patch, subpatch.c00, subpatch.c10);
+      subpatch.edge_u1.T = T(subpatch.patch, subpatch.c01, subpatch.c11);
+      subpatch.edge_v0.T = T(subpatch.patch, subpatch.c00, subpatch.c01);
+      subpatch.edge_v1.T = T(subpatch.patch, subpatch.c10, subpatch.c11);
+
+      resolve_edge_factors(subpatch);
+
+      split(subpatch, 0);
+    }
+
+    /* Update offsets after T is known from split. */
+    edge_u1->top = subpatch.edge_v0.edge;
+    edge_u1->stitch_top_offset = edge_u1->top->T * (v0_reversed ? -1 : 1);
+    edge_v1->top = subpatch.edge_u0.edge;
+    edge_v1->stitch_top_offset = edge_v1->top->T * (!u0_reversed ? -1 : 1);
+
+    if (corner == 0) {
+      first_edge_v0 = subpatch.edge_v0.edge;
+    }
+
+    if (prev_edge_u0) {
+      if (v0_reversed) {
+        subpatch.edge_v0.edge->stitch_offset = prev_edge_u0->T;
+      }
+      else {
+        prev_edge_u0->stitch_offset = subpatch.edge_v0.edge->T;
+      }
+
+      int T = subpatch.edge_v0.edge->T + prev_edge_u0->T;
+      subpatch.edge_v0.edge->stitch_edge_T = T;
+      prev_edge_u0->stitch_edge_T = T;
+    }
 
-  split(sub_split, ef_split);
+    if (corner == face.num_corners - 1) {
+      if (v0_reversed) {
+        subpatch.edge_u0.edge->stitch_offset = first_edge_v0->T;
+      }
+      else {
+        first_edge_v0->stitch_offset = subpatch.edge_u0.edge->T;
+      }
+
+      int T = first_edge_v0->T + subpatch.edge_u0.edge->T;
+      first_edge_v0->stitch_edge_T = T;
+      subpatch.edge_u0.edge->stitch_edge_T = T;
+    }
+
+    prev_edge_u0 = subpatch.edge_u0.edge;
+  }
+}
+
+void DiagSplit::post_split()
+{
+  int num_stitch_verts = 0;
+
+  /* All patches are now split, and all T values known. */
 
+  foreach (Edge &edge, edges) {
+    if (edge.second_vert_index < 0) {
+      edge.second_vert_index = alloc_verts(edge.T - 1);
+    }
+
+    if (edge.is_stitch_edge) {
+      num_stitch_verts = max(num_stitch_verts,
+                             max(edge.stitch_start_vert_index, edge.stitch_end_vert_index));
+    }
+  }
+
+  num_stitch_verts += 1;
+
+  /* Map of edge key to edge stitching vert offset. */
+  struct pair_hasher {
+    size_t operator()(const pair<int, int> &k) const
+    {
+      return hash_uint2(k.first, k.second);
+    }
+  };
+  typedef unordered_map<pair<int, int>, int, pair_hasher> edge_stitch_verts_map_t;
+  edge_stitch_verts_map_t edge_stitch_verts_map;
+
+  foreach (Edge &edge, edges) {
+    if (edge.is_stitch_edge) {
+      if (edge.stitch_edge_T == 0) {
+        edge.stitch_edge_T = edge.T;
+      }
+
+      if (edge_stitch_verts_map.find(edge.stitch_edge_key) == edge_stitch_verts_map.end()) {
+        edge_stitch_verts_map[edge.stitch_edge_key] = num_stitch_verts;
+        num_stitch_verts += edge.stitch_edge_T - 1;
+      }
+    }
+  }
+
+  /* Set start and end indices for edges generated from a split. */
+  foreach (Edge &edge, edges) {
+    if (edge.start_vert_index < 0) {
+      /* Fixup offsets. */
+      if (edge.top_indices_decrease) {
+        edge.top_offset = edge.top->T - edge.top_offset;
+      }
+
+      edge.start_vert_index = edge.top->get_vert_along_edge(edge.top_offset);
+    }
+
+    if (edge.end_vert_index < 0) {
+      if (edge.bottom_indices_decrease) {
+        edge.bottom_offset = edge.bottom->T - edge.bottom_offset;
+      }
+
+      edge.end_vert_index = edge.bottom->get_vert_along_edge(edge.bottom_offset);
+    }
+  }
+
+  int vert_offset = params.mesh->verts.size();
+
+  /* Add verts to stitching map. */
+  foreach (const Edge &edge, edges) {
+    if (edge.is_stitch_edge) {
+      int second_stitch_vert_index = edge_stitch_verts_map[edge.stitch_edge_key];
+
+      for (int i = 0; i <= edge.T; i++) {
+        /* Get proper stitching key. */
+        int key;
+
+        if (i == 0) {
+          key = edge.stitch_start_vert_index;
+        }
+        else if (i == edge.T) {
+          key = edge.stitch_end_vert_index;
+        }
+        else {
+          key = second_stitch_vert_index + i - 1 + edge.stitch_offset;
+        }
+
+        if (key == STITCH_NGON_SPLIT_EDGE_CENTER_VERT_TAG) {
+          if (i == 0) {
+            key = second_stitch_vert_index - 1 + edge.stitch_offset;
+          }
+          else if (i == edge.T) {
+            key = second_stitch_vert_index - 1 + edge.T;
+          }
+        }
+        else if (key < 0 && edge.top) { /* ngon spoke edge */
+          int s = edge_stitch_verts_map[edge.top->stitch_edge_key];
+          if (edge.stitch_top_offset >= 0) {
+            key = s - 1 + edge.stitch_top_offset;
+          }
+          else {
+            key = s - 1 + edge.top->stitch_edge_T + edge.stitch_top_offset;
+          }
+        }
+
+        /* Get real vert index. */
+        int vert = edge.get_vert_along_edge(i) + vert_offset;
+
+        /* Add to map */
+        if (params.mesh->vert_to_stitching_key_map.find(vert) ==
+            params.mesh->vert_to_stitching_key_map.end()) {
+          params.mesh->vert_to_stitching_key_map[vert] = key;
+          params.mesh->vert_stitching_map.insert({key, vert});
+        }
+      }
+    }
+  }
+
+  /* Dice; TODO(mai): Move this out of split. */
   QuadDice dice(params);
 
-  for (size_t i = 0; i < subpatches_quad.size(); i++) {
-    QuadDice::SubPatch &sub = subpatches_quad[i];
-    QuadDice::EdgeFactors &ef = edgefactors_quad[i];
+  int num_verts = num_alloced_verts;
+  int num_triangles = 0;
+
+  for (size_t i = 0; i < subpatches.size(); i++) {
+    subpatches[i].inner_grid_vert_offset = num_verts;
+    num_verts += subpatches[i].calc_num_inner_verts();
+    num_triangles += subpatches[i].calc_num_triangles();
+  }
+
+  dice.reserve(num_verts, num_triangles);
+
+  for (size_t i = 0; i < subpatches.size(); i++) {
+    Subpatch &sub = subpatches[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);
+    sub.edge_u0.T = max(sub.edge_u0.T, 1);
+    sub.edge_u1.T = max(sub.edge_u1.T, 1);
+    sub.edge_v0.T = max(sub.edge_v0.T, 1);
+    sub.edge_v1.T = max(sub.edge_v1.T, 1);
 
-    dice.dice(sub, ef);
+    dice.dice(sub);
   }
 
-  subpatches_quad.clear();
-  edgefactors_quad.clear();
+  /* Cleanup */
+  subpatches.clear();
+  edges.clear();
 }
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/subd/subd_split.h b/intern/cycles/subd/subd_split.h
index 6e68d8ee598..0fc22872c62 100644
--- a/intern/cycles/subd/subd_split.h
+++ b/intern/cycles/subd/subd_split.h
@@ -23,35 +23,51 @@
  * for more details. */
 
 #include "subd/subd_dice.h"
+#include "subd/subd_subpatch.h"
 
+#include "util/util_deque.h"
 #include "util/util_types.h"
 #include "util/util_vector.h"
 
+#include <deque>
+
 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;
-
   SubdParams params;
 
-  explicit DiagSplit(const SubdParams &params);
+  vector<Subpatch> subpatches;
+  /* deque is used so that element pointers remain vaild when size is changed. */
+  deque<Edge> edges;
 
   float3 to_world(Patch *patch, float2 uv);
-  int T(Patch *patch, float2 Pstart, float2 Pend);
+  int T(Patch *patch, float2 Pstart, float2 Pend, bool recursive_resolve=false);
+
+  void limit_edge_factor(int &T, Patch *patch, float2 Pstart, float2 Pend);
+  void resolve_edge_factors(Subpatch &sub);
+
   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 split(Subpatch &sub, int depth = 0);
+
+  int num_alloced_verts = 0;
+  int alloc_verts(int n); /* Returns start index of new verts. */
+
+ public:
+  Edge *alloc_edge();
+
+  explicit DiagSplit(const SubdParams &params);
+
+  void split_patches(Patch *patches, size_t patches_byte_stride);
+
+  void split_quad(const Mesh::SubdFace &face, Patch *patch);
+  void split_ngon(const Mesh::SubdFace &face, Patch *patches, size_t patches_byte_stride);
 
-  void split_quad(Patch *patch, QuadDice::SubPatch *subpatch = NULL);
+  void post_split();
 };
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/subd/subd_subpatch.h b/intern/cycles/subd/subd_subpatch.h
new file mode 100644
index 00000000000..4f0c4afee81
--- /dev/null
+++ b/intern/cycles/subd/subd_subpatch.h
@@ -0,0 +1,210 @@
+/*
+ * Copyright 2011-2018 Blender Foundation
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef __SUBD_SUBPATCH_H__
+#define __SUBD_SUBPATCH_H__
+
+#include "util/util_map.h"
+#include "util/util_types.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* Subpatch */
+
+class Subpatch {
+ public:
+  class Patch *patch; /* Patch this is a subpatch of. */
+  int inner_grid_vert_offset;
+
+  struct edge_t {
+    int T;
+    int offset; /* Offset along main edge, interpretation depends on the two flags below. */
+
+    bool indices_decrease_along_edge;
+    bool sub_edges_created_in_reverse_order;
+
+    struct Edge *edge;
+
+    int get_vert_along_edge(int n) const;
+  };
+
+  /*
+   *            eu1
+   *     c01 --------- c11
+   *     |               |
+   * ev0 |               | ev1
+   *     |               |
+   *     c00 --------- c10
+   *            eu0
+   */
+
+  union {
+    float2 corners[4]; /* UV within patch, clockwise starting from uv (0, 0) towards (0, 1) etc. */
+    struct {
+      float2 c00, c01, c11, c10;
+    };
+  };
+
+  union {
+    edge_t
+        edges[4]; /* Edges of this subpatch, each edge starts at the corner of the same index. */
+    struct {
+      edge_t edge_v0, edge_u1, edge_v1, edge_u0;
+    };
+  };
+
+  explicit Subpatch(Patch *patch = nullptr)
+      : patch(patch),
+        c00(make_float2(0.0f, 0.0f)),
+        c01(make_float2(0.0f, 1.0f)),
+        c11(make_float2(1.0f, 1.0f)),
+        c10(make_float2(1.0f, 0.0f))
+  {
+  }
+
+  Subpatch(Patch *patch, float2 c00, float2 c01, float2 c11, float2 c10)
+      : patch(patch), c00(c00), c01(c01), c11(c11), c10(c10)
+  {
+  }
+
+  int calc_num_inner_verts() const
+  {
+    int Mu = max(edge_u0.T, edge_u1.T);
+    int Mv = max(edge_v0.T, edge_v1.T);
+    Mu = max(Mu, 2);
+    Mv = max(Mv, 2);
+    return (Mu - 1) * (Mv - 1);
+  }
+
+  int calc_num_triangles() const
+  {
+    int Mu = max(edge_u0.T, edge_u1.T);
+    int Mv = max(edge_v0.T, edge_v1.T);
+    Mu = max(Mu, 2);
+    Mv = max(Mv, 2);
+
+    int inner_triangles = (Mu - 2) * (Mv - 2) * 2;
+    int edge_triangles = edge_u0.T + edge_u1.T + edge_v0.T + edge_v1.T + (Mu - 2) * 2 +
+                         (Mv - 2) * 2;
+
+    return inner_triangles + edge_triangles;
+  }
+
+  int get_vert_along_edge(int e, int n) const;
+
+  int get_vert_along_grid_edge(int edge, int n) const
+  {
+    int Mu = max(edge_u0.T, edge_u1.T);
+    int Mv = max(edge_v0.T, edge_v1.T);
+    Mu = max(Mu, 2);
+    Mv = max(Mv, 2);
+
+    switch (edge) {
+      case 0:
+        return inner_grid_vert_offset + n * (Mu - 1);
+      case 1:
+        return inner_grid_vert_offset + (Mu - 1) * (Mv - 2) + n;
+      case 2:
+        return inner_grid_vert_offset + ((Mu - 1) * (Mv - 1) - 1) - n * (Mu - 1);
+      case 3:
+        return inner_grid_vert_offset + (Mu - 2) - n;
+    }
+
+    return -1;
+  }
+};
+
+struct Edge {
+  int T; /* Number of segments the edge will be diced into, see DiagSplit paper. */
+
+  Edge *top, *bottom; /* top is edge adjacent to start, bottom is adjacent to end. */
+  int top_offset, bottom_offset;
+  bool top_indices_decrease, bottom_indices_decrease;
+
+  int start_vert_index;
+  int end_vert_index;
+  int second_vert_index; /* Index of the second vert from this edges corner along the edge towards the next corner */
+
+  bool is_stitch_edge; /* Verties on edge are to be stitched. */
+  /* Key to match this edge with others to be stitched with. The ints in the pair are ordered stitching indices */
+  pair<int, int> stitch_edge_key;
+
+  int stitch_edge_T; /* Full T along edge (may be larger than T for edges split from ngon edges) */
+  int stitch_offset;
+  int stitch_top_offset;
+  int stitch_start_vert_index;
+  int stitch_end_vert_index;
+
+  Edge()
+      : T(0),
+        top(nullptr),
+        bottom(nullptr),
+        top_offset(-1),
+        bottom_offset(-1),
+        top_indices_decrease(false),
+        bottom_indices_decrease(false),
+        start_vert_index(-1),
+        end_vert_index(-1),
+        second_vert_index(-1),
+        is_stitch_edge(false),
+        stitch_edge_T(0),
+        stitch_offset(0)
+  {
+  }
+
+  int get_vert_along_edge(int n) const
+  {
+    assert(n >= 0 && n <= T);
+
+    if (n == 0) {
+      return start_vert_index;
+    }
+    else if (n == T) {
+      return end_vert_index;
+    }
+
+    return second_vert_index + n - 1;
+  }
+};
+
+inline int Subpatch::edge_t::get_vert_along_edge(int n) const
+{
+  assert(n >= 0 && n <= T);
+
+  if (!indices_decrease_along_edge && !sub_edges_created_in_reverse_order) {
+    n = offset + n;
+  }
+  else if (!indices_decrease_along_edge && sub_edges_created_in_reverse_order) {
+    n = edge->T - offset - T + n;
+  }
+  else if (indices_decrease_along_edge && !sub_edges_created_in_reverse_order) {
+    n = offset + T - n;
+  }
+  else if (indices_decrease_along_edge && sub_edges_created_in_reverse_order) {
+    n = edge->T - offset - n;
+  }
+
+  return edge->get_vert_along_edge(n);
+}
+
+inline int Subpatch::get_vert_along_edge(int edge, int n) const
+{
+  return edges[edge].get_vert_along_edge(n);
+}
+
+CCL_NAMESPACE_END
+
+#endif /* __SUBD_SUBPATCH_H__ */