diff --git a/intern/cycles/render/mesh.cpp b/intern/cycles/render/mesh.cpp
index c7e0430fe7e..3d2cc98d38b 100644
--- a/intern/cycles/render/mesh.cpp
+++ b/intern/cycles/render/mesh.cpp
@@ -560,6 +560,9 @@ void Mesh::clear(bool preserve_voxel_data)
 
   used_shaders.clear();
 
+  vert_to_stitching_key_map.clear();
+  vert_stitching_map.clear();
+
   if (!preserve_voxel_data) {
     geometry_flags = GEOMETRY_NONE;
   }
diff --git a/intern/cycles/render/mesh.h b/intern/cycles/render/mesh.h
index 50b1f0e69f1..4a24a9c2656 100644
--- a/intern/cycles/render/mesh.h
+++ b/intern/cycles/render/mesh.h
@@ -283,6 +283,14 @@ class Mesh : public Node {
 
   size_t num_subd_verts;
 
+ private:
+  unordered_map<int, int> vert_to_stitching_key_map; /* real vert index -> stitching index */
+  unordered_multimap<int, int>
+      vert_stitching_map; /* stitching index -> multiple real vert indices */
+  friend class DiagSplit;
+  friend class MeshManager;
+
+ public:
   /* Functions */
   Mesh();
   ~Mesh();
diff --git a/intern/cycles/render/mesh_displace.cpp b/intern/cycles/render/mesh_displace.cpp
index 5ae9348d83e..6a6c2fbb3eb 100644
--- a/intern/cycles/render/mesh_displace.cpp
+++ b/intern/cycles/render/mesh_displace.cpp
@@ -22,7 +22,9 @@
 #include "render/shader.h"
 
 #include "util/util_foreach.h"
+#include "util/util_map.h"
 #include "util/util_progress.h"
+#include "util/util_set.h"
 
 CCL_NAMESPACE_BEGIN
 
@@ -184,6 +186,38 @@ bool MeshManager::displace(
 
   d_output.free();
 
+  /* stitch */
+  unordered_set<int> stitch_keys;
+  for (pair<int, int> i : mesh->vert_to_stitching_key_map) {
+    stitch_keys.insert(i.second); /* stitching index */
+  }
+
+  typedef unordered_multimap<int, int>::iterator map_it_t;
+
+  for (int key : stitch_keys) {
+    pair<map_it_t, map_it_t> verts = mesh->vert_stitching_map.equal_range(key);
+
+    float3 pos = make_float3(0.0f, 0.0f, 0.0f);
+    int num = 0;
+
+    for (map_it_t v = verts.first; v != verts.second; ++v) {
+      int vert = v->second;
+
+      pos += mesh->verts[vert];
+      num++;
+    }
+
+    if (num <= 1) {
+      continue;
+    }
+
+    pos *= 1.0f / num;
+
+    for (map_it_t v = verts.first; v != verts.second; ++v) {
+      mesh->verts[v->second] = pos;
+    }
+  }
+
   /* for displacement method both, we only need to recompute the face
    * normals, as bump mapping in the shader will already alter the
    * vertex normal, so we start from the non-displaced vertex normals
@@ -238,7 +272,25 @@ bool MeshManager::displace(
     for (size_t i = 0; i < num_triangles; i++) {
       if (tri_has_true_disp[i]) {
         for (size_t j = 0; j < 3; j++) {
-          vN[mesh->get_triangle(i).v[j]] += fN[i];
+          int vert = mesh->get_triangle(i).v[j];
+          vN[vert] += fN[i];
+
+          /* add face normals to stitched vertices */
+          if (stitch_keys.size()) {
+            map_it_t key = mesh->vert_to_stitching_key_map.find(vert);
+
+            if (key != mesh->vert_to_stitching_key_map.end()) {
+              pair<map_it_t, map_it_t> verts = mesh->vert_stitching_map.equal_range(key->second);
+
+              for (map_it_t v = verts.first; v != verts.second; ++v) {
+                if (v->second == vert) {
+                  continue;
+                }
+
+                vN[v->second] += fN[i];
+              }
+            }
+          }
         }
       }
     }
@@ -289,8 +341,27 @@ bool MeshManager::displace(
         for (size_t i = 0; i < num_triangles; i++) {
           if (tri_has_true_disp[i]) {
             for (size_t j = 0; j < 3; j++) {
+              int vert = mesh->get_triangle(i).v[j];
               float3 fN = compute_face_normal(mesh->get_triangle(i), mP);
-              mN[mesh->get_triangle(i).v[j]] += fN;
+              mN[vert] += fN;
+
+              /* add face normals to stitched vertices */
+              if (stitch_keys.size()) {
+                map_it_t key = mesh->vert_to_stitching_key_map.find(vert);
+
+                if (key != mesh->vert_to_stitching_key_map.end()) {
+                  pair<map_it_t, map_it_t> verts = mesh->vert_stitching_map.equal_range(
+                      key->second);
+
+                  for (map_it_t v = verts.first; v != verts.second; ++v) {
+                    if (v->second == vert) {
+                      continue;
+                    }
+
+                    mN[v->second] += fN;
+                  }
+                }
+              }
             }
           }
         }
diff --git a/intern/cycles/render/mesh_subdivision.cpp b/intern/cycles/render/mesh_subdivision.cpp
index 46c8240fb71..a5a35fc049e 100644
--- a/intern/cycles/render/mesh_subdivision.cpp
+++ b/intern/cycles/render/mesh_subdivision.cpp
@@ -24,6 +24,7 @@
 
 #include "util/util_foreach.h"
 #include "util/util_algorithm.h"
+#include "util/util_hash.h"
 
 CCL_NAMESPACE_BEGIN
 
@@ -318,6 +319,9 @@ class OsdData {
 struct OsdPatch : Patch {
   OsdData *osd_data;
 
+  OsdPatch()
+  {
+  }
   OsdPatch(OsdData *data) : osd_data(data)
   {
   }
@@ -358,11 +362,6 @@ struct OsdPatch : Patch {
       *N = (t != 0.0f) ? *N / t : make_float3(0.0f, 0.0f, 1.0f);
     }
   }
-
-  BoundBox bound()
-  {
-    return BoundBox::empty;
-  }
 };
 
 #endif
@@ -397,29 +396,62 @@ void Mesh::tessellate(DiagSplit *split)
   Attribute *attr_vN = subd_attributes.find(ATTR_STD_VERTEX_NORMAL);
   float3 *vN = attr_vN->data_float3();
 
+  /* count patches */
+  int num_patches = 0;
   for (int f = 0; f < num_faces; f++) {
     SubdFace &face = subd_faces[f];
 
     if (face.is_quad()) {
-      /* quad */
-      QuadDice::SubPatch subpatch;
+      num_patches++;
+    }
+    else {
+      num_patches += face.num_corners;
+    }
+  }
 
-      LinearQuadPatch quad_patch;
+  /* build patches from faces */
 #ifdef WITH_OPENSUBDIV
-      OsdPatch osd_patch(&osd_data);
+  if (subdivision_type == SUBDIVISION_CATMULL_CLARK) {
+    vector<OsdPatch> osd_patches(num_patches, &osd_data);
+    OsdPatch *patch = osd_patches.data();
 
-      if (subdivision_type == SUBDIVISION_CATMULL_CLARK) {
-        osd_patch.patch_index = face.ptex_offset;
+    for (int f = 0; f < num_faces; f++) {
+      SubdFace &face = subd_faces[f];
 
-        subpatch.patch = &osd_patch;
+      if (face.is_quad()) {
+        patch->patch_index = face.ptex_offset;
+        patch->from_ngon = false;
+        patch->shader = face.shader;
+        patch++;
       }
-      else
-#endif
-      {
-        float3 *hull = quad_patch.hull;
-        float3 *normals = quad_patch.normals;
+      else {
+        for (int corner = 0; corner < face.num_corners; corner++) {
+          patch->patch_index = face.ptex_offset + corner;
+          patch->from_ngon = true;
+          patch->shader = face.shader;
+          patch++;
+        }
+      }
+    }
 
-        quad_patch.patch_index = face.ptex_offset;
+    /* split patches */
+    split->split_patches(osd_patches.data(), sizeof(OsdPatch));
+  }
+  else
+#endif
+  {
+    vector<LinearQuadPatch> linear_patches(num_patches);
+    LinearQuadPatch *patch = linear_patches.data();
+
+    for (int f = 0; f < num_faces; f++) {
+      SubdFace &face = subd_faces[f];
+
+      if (face.is_quad()) {
+        float3 *hull = patch->hull;
+        float3 *normals = patch->normals;
+
+        patch->patch_index = face.ptex_offset;
+        patch->from_ngon = false;
 
         for (int i = 0; i < 4; i++) {
           hull[i] = verts[subd_face_corners[face.start_corner + i]];
@@ -440,55 +472,11 @@ void Mesh::tessellate(DiagSplit *split)
         swap(hull[2], hull[3]);
         swap(normals[2], normals[3]);
 
-        subpatch.patch = &quad_patch;
+        patch->shader = face.shader;
+        patch++;
       }
-
-      subpatch.patch->shader = face.shader;
-
-      /* Quad faces need to be split at least once to line up with split ngons, we do this
-       * here in this manner because if we do it later edge factors may end up slightly off.
-       */
-      subpatch.P00 = make_float2(0.0f, 0.0f);
-      subpatch.P10 = make_float2(0.5f, 0.0f);
-      subpatch.P01 = make_float2(0.0f, 0.5f);
-      subpatch.P11 = make_float2(0.5f, 0.5f);
-      split->split_quad(subpatch.patch, &subpatch);
-
-      subpatch.P00 = make_float2(0.5f, 0.0f);
-      subpatch.P10 = make_float2(1.0f, 0.0f);
-      subpatch.P01 = make_float2(0.5f, 0.5f);
-      subpatch.P11 = make_float2(1.0f, 0.5f);
-      split->split_quad(subpatch.patch, &subpatch);
-
-      subpatch.P00 = make_float2(0.0f, 0.5f);
-      subpatch.P10 = make_float2(0.5f, 0.5f);
-      subpatch.P01 = make_float2(0.0f, 1.0f);
-      subpatch.P11 = make_float2(0.5f, 1.0f);
-      split->split_quad(subpatch.patch, &subpatch);
-
-      subpatch.P00 = make_float2(0.5f, 0.5f);
-      subpatch.P10 = make_float2(1.0f, 0.5f);
-      subpatch.P01 = make_float2(0.5f, 1.0f);
-      subpatch.P11 = make_float2(1.0f, 1.0f);
-      split->split_quad(subpatch.patch, &subpatch);
-    }
-    else {
-      /* ngon */
-#ifdef WITH_OPENSUBDIV
-      if (subdivision_type == SUBDIVISION_CATMULL_CLARK) {
-        OsdPatch patch(&osd_data);
-
-        patch.shader = face.shader;
-
-        for (int corner = 0; corner < face.num_corners; corner++) {
-          patch.patch_index = face.ptex_offset + corner;
-
-          split->split_quad(&patch);
-        }
-      }
-      else
-#endif
-      {
+      else {
+        /* ngon */
         float3 center_vert = make_float3(0.0f, 0.0f, 0.0f);
         float3 center_normal = make_float3(0.0f, 0.0f, 0.0f);
 
@@ -499,13 +487,13 @@ void Mesh::tessellate(DiagSplit *split)
         }
 
         for (int corner = 0; corner < face.num_corners; corner++) {
-          LinearQuadPatch patch;
-          float3 *hull = patch.hull;
-          float3 *normals = patch.normals;
+          float3 *hull = patch->hull;
+          float3 *normals = patch->normals;
 
-          patch.patch_index = face.ptex_offset + corner;
+          patch->patch_index = face.ptex_offset + corner;
+          patch->from_ngon = true;
 
-          patch.shader = face.shader;
+          patch->shader = face.shader;
 
           hull[0] =
               verts[subd_face_corners[face.start_corner + mod(corner + 0, face.num_corners)]];
@@ -537,10 +525,13 @@ void Mesh::tessellate(DiagSplit *split)
             }
           }
 
-          split->split_quad(&patch);
+          patch++;
         }
       }
     }
+
+    /* split patches */
+    split->split_patches(linear_patches.data(), sizeof(LinearQuadPatch));
   }
 
   /* interpolate center points for attributes */
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());
+  assert(index < 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();
-
-    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 = 0; i < t; i++) {
+    float f = i / (float)t;
 
-  for (int i = 1; i < tu; i++) {
-    float u = i / (float)tu;
-    float v = (side) ? 1.0f : 0.0f;
+    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;
+    }
 
-    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));
+    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);
+  /* inner grid */
+  add_grid(sub, Mu, Mv, sub.inner_grid_vert_offset);
 
-  /* corners and inner grid */
-  add_corners(sub);
-  add_grid(sub, Mu, Mv, offset);
+  /* sides */
+  set_side(sub, 0);
+  set_side(sub, 1);
+  set_side(sub, 2);
+  set_side(sub, 3);
 
-  /* 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);
-
-  /* right side */
-  add_side_v(sub, outer, inner, Mu, Mv, ef.tv1, 1, offset);
-  stitch_triangles(sub.patch, outer, inner);
-
-  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);
+  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);
+  float2 map_uv(Subpatch &sub, float u, float v);
+  void set_vert(Subpatch &sub, int index, float u, float v);
 
-  void add_corners(SubPatch &sub);
-  void add_grid(SubPatch &sub, int Mu, int Mv, int offset);
+  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);
+  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,16 +31,14 @@ CCL_NAMESPACE_BEGIN
 
 /* DiagSplit */
 
+#define DSPLIT_NON_UNIFORM -1
+#define STITCH_NGON_CENTER_VERT_INDEX_OFFSET 0x60000000
+#define STITCH_NGON_SPLIT_EDGE_CENTER_VERT_TAG (0x60000000 - 1)
+
 DiagSplit::DiagSplit(const SubdParams &params_) : params(params_)
 {
 }
 
-void DiagSplit::dispatch(QuadDice::SubPatch &sub, QuadDice::EdgeFactors &ef)
-{
-  subpatches_quad.push_back(sub);
-  edgefactors_quad.push_back(ef);
-}
-
 float3 DiagSplit::to_world(Patch *patch, float2 uv)
 {
   float3 P;
@@ -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;
-      }
-
-      Lsum += L;
-      Lmax = max(L, Lmax);
+    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;
+    }
+
+    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)
 {
-  float2 P00 = sub.P00;
-  float2 P01 = sub.P01;
-  float2 P10 = sub.P10;
-  float2 P11 = sub.P11;
+  int max_t = 1 << params.max_level;
+  int max_t_for_edge = int(max_t * len(Pstart - Pend));
 
-  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));
+  if (patch->from_ngon) {
+    max_t_for_edge >>= 1; /* Initial split of ngon causes edges to extend half the distance. */
+  }
 
-  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);
+  T = (max_t_for_edge <= 1) ? 1 : min(T, max_t_for_edge);
+
+  assert(T >= 1 || T == DSPLIT_NON_UNIFORM);
 }
 
-void DiagSplit::split(QuadDice::SubPatch &sub, QuadDice::EdgeFactors &ef, int depth)
+void DiagSplit::resolve_edge_factors(Subpatch &sub)
+{
+  /* 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(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");
 
-    dispatch(sub, ef);
+    sub.edge_u0.T = 1;
+    sub.edge_u1.T = 1;
+    sub.edge_v0.T = 1;
+    sub.edge_v1.T = 1;
+
+    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];
 
-    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);
+    /* Update T values and offsets. */
+    for (int i = 0; i < 4; i++) {
+      Subpatch::edge_t &edge = subpatch.edges[i];
 
-    /* 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};
-
-    limit_edge_factors(sub0, ef0, 1 << params.max_level);
-    limit_edge_factors(sub1, ef1, 1 << params.max_level);
-
-    split(sub0, ef0, depth + 1);
-    split(sub1, ef1, depth + 1);
-  }
-  else if (split_v) {
-    /* 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};
-
-    limit_edge_factors(sub0, ef0, 1 << params.max_level);
-    limit_edge_factors(sub1, ef1, 1 << params.max_level);
-
-    split(sub0, ef0, depth + 1);
-    split(sub1, ef1, depth + 1);
+      edge.offset = edge.edge->T;
+      edge.edge->T += edge.T;
+    }
   }
   else {
-    dispatch(sub, ef);
+    /* 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 edges */
+    float2 P0, P1;
+
+    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);
+
+    /* Split */
+    *Pa = P1;
+    *Pb = P0;
+
+    *Pc = P1;
+    *Pd = P0;
+
+    int tsplit = T(sub.patch, P0, P1);
+
+    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;
   }
 }
 
-void DiagSplit::split_quad(Patch *patch, QuadDice::SubPatch *subpatch)
+int DiagSplit::alloc_verts(int n)
 {
-  QuadDice::SubPatch sub_split;
-  QuadDice::EdgeFactors ef_split;
+  int a = num_alloced_verts;
+  num_alloced_verts += n;
+  return a;
+}
 
-  if (subpatch) {
-    sub_split = *subpatch;
+Edge *DiagSplit::alloc_edge()
+{
+  edges.emplace_back();
+  return &edges.back();
+}
+
+void DiagSplit::split_patches(Patch *patches, size_t patches_byte_stride)
+{
+  int patch_index = 0;
+
+  for (int f = 0; f < params.mesh->subd_faces.size(); f++) {
+    Mesh::SubdFace &face = params.mesh->subd_faces[f];
+
+    Patch *patch = (Patch *)(((char *)patches) + patch_index * patches_byte_stride);
+
+    if (face.is_quad()) {
+      patch_index++;
+
+      split_quad(face, patch);
+    }
+    else {
+      patch_index += face.num_corners;
+
+      split_ngon(face, patch, patches_byte_stride);
+    }
+  }
+
+  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(const Mesh::SubdFace &face, Patch *patch)
+{
+  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 (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;
+  }
 
-  split(sub_split, ef_split);
+  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;
+    }
+
+    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;
 
-    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);
+  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();
   }
 
-  subpatches_quad.clear();
-  edgefactors_quad.clear();
+  dice.reserve(num_verts, num_triangles);
+
+  for (size_t i = 0; i < subpatches.size(); i++) {
+    Subpatch &sub = subpatches[i];
+
+    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);
+  }
+
+  /* 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);
 
-  void split_quad(Patch *patch, QuadDice::SubPatch *subpatch = NULL);
+  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 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__ */
diff --git a/intern/cycles/util/CMakeLists.txt b/intern/cycles/util/CMakeLists.txt
index 4f1c82897cc..a6a6b6704f0 100644
--- a/intern/cycles/util/CMakeLists.txt
+++ b/intern/cycles/util/CMakeLists.txt
@@ -55,6 +55,7 @@ set(SRC_HEADERS
   util_boundbox.h
   util_debug.h
   util_defines.h
+  util_deque.h
   util_guarded_allocator.cpp
   util_foreach.h
   util_function.h
diff --git a/intern/cycles/util/util_deque.h b/intern/cycles/util/util_deque.h
new file mode 100644
index 00000000000..ccac961aa7d
--- /dev/null
+++ b/intern/cycles/util/util_deque.h
@@ -0,0 +1,28 @@
+/*
+ * 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 __UTIL_DEQUE_H__
+#define __UTIL_DEQUE_H__
+
+#include <deque>
+
+CCL_NAMESPACE_BEGIN
+
+using std::deque;
+
+CCL_NAMESPACE_END
+
+#endif /* __UTIL_DEQUE_H__ */
diff --git a/intern/cycles/util/util_map.h b/intern/cycles/util/util_map.h
index 8385b08dd5a..f1b2522362f 100644
--- a/intern/cycles/util/util_map.h
+++ b/intern/cycles/util/util_map.h
@@ -25,6 +25,7 @@ CCL_NAMESPACE_BEGIN
 using std::map;
 using std::pair;
 using std::unordered_map;
+using std::unordered_multimap;
 
 template<typename T> static void map_free_memory(T &data)
 {