Subsurf: Support subdivision of loose elements

Applies to vertices and edges. Biggest annoyance here is that OpenSubdiv's
topology converter expects that there is no loose geometry, otherwise it
is getting confused.

For now solution is to create some sort of mapping from real Mesh vertex
and edge index to a non-loose-index. Now the annoying part is that this
is an extra step to calculate before we can compare topology, meaning FPS
will not be as great as if we knew for sure that topology didn't change.

Loose edges subdivision is different from what it used to be with old
subdivision code, but probably nice feature now is that endpoints of loose
edges are stay at the coarse vertex locations. This allows to have things
like plane with hair strands, without need to duplicate edge vertices at
endpoints.

All this required some re-work of topology refiner creation, which is now
only passing edges and vertices which are adjacent to face. This is how
topology refiner is supposed to be used, and this is how its validator
also works. Vertices which are adjacent to loose edges are marked as
infinite sharp. This seems to be good-enough approximation for now. In the
future we might tweaks things a bit and push such vertices in average
direction of loose edges, to match old subdivision code closer.

6 files changed, 473 insertions, 103 deletions

diff --git a/intern/opensubdiv/internal/opensubdiv_converter_factory.cc b/intern/opensubdiv/internal/opensubdiv_converter_factory.cc
index 48516cc80b7..06707cbf418 100644
--- a/intern/opensubdiv/internal/opensubdiv_converter_factory.cc
+++ b/intern/opensubdiv/internal/opensubdiv_converter_factory.cc
@@ -346,16 +346,10 @@ inline bool TopologyRefinerFactory<TopologyRefinerData>::assignComponentTags(
   const int num_vertices = converter->getNumVertices(converter);
   for (int vertex_index = 0; vertex_index < num_vertices; ++vertex_index) {
     ConstIndexArray vertex_edges = getBaseVertexEdges(refiner, vertex_index);
-    for (int i = 0; i < vertex_edges.size(); ++i) {
-      const int edge_index = vertex_edges[i];
-      ConstIndexArray edge_faces = getBaseEdgeFaces(refiner, edge_index);
-      if (edge_faces.size() == 0) {
-        setBaseVertexSharpness(refiner, vertex_index,
-                               Crease::SHARPNESS_INFINITE);
-        break;
-      }
-    }
-    if (vertex_edges.size() == 2) {
+    if (converter->isInfiniteSharpVertex(converter, vertex_index)) {
+      setBaseVertexSharpness(
+          refiner, vertex_index, Crease::SHARPNESS_INFINITE);
+    } else if (vertex_edges.size() == 2) {
       const int edge0 = vertex_edges[0], edge1 = vertex_edges[1];
       const float sharpness0 = converter->getEdgeSharpness(converter, edge0);
       const float sharpness1 = converter->getEdgeSharpness(converter, edge1);
diff --git a/intern/opensubdiv/internal/opensubdiv_internal.h b/intern/opensubdiv/internal/opensubdiv_internal.h
index b9d196bbe9e..16365896edf 100644
--- a/intern/opensubdiv/internal/opensubdiv_internal.h
+++ b/intern/opensubdiv/internal/opensubdiv_internal.h
@@ -24,9 +24,6 @@
 // Never do for release builds.
 #  undef OPENSUBDIV_VALIDATE_TOPOLOGY
 #else
-// TODO(sergey): Always disabled for now, the check doesn't handle multiple
-// non-manifolds from the OpenSubdiv side currently.
-// #  undef OPENSUBDIV_VALIDATE_TOPOLOGY
 #  define OPENSUBDIV_VALIDATE_TOPOLOGY
 #endif
 
diff --git a/intern/opensubdiv/opensubdiv_converter_capi.h b/intern/opensubdiv/opensubdiv_converter_capi.h
index a939f1117e0..58a231deb30 100644
--- a/intern/opensubdiv/opensubdiv_converter_capi.h
+++ b/intern/opensubdiv/opensubdiv_converter_capi.h
@@ -19,6 +19,8 @@
 #ifndef OPENSUBDIV_CONVERTER_CAPI_H_
 #define OPENSUBDIV_CONVERTER_CAPI_H_
 
+#include <stdint.h>  // for bool
+
 #include "opensubdiv_capi_type.h"
 
 #ifdef __cplusplus
@@ -34,6 +36,7 @@ typedef struct OpenSubdiv_Converter {
 
   //////////////////////////////////////////////////////////////////////////////
   // Global geometry counters.
+
   // Number of faces/edges/vertices in the base mesh.
   int (*getNumFaces)(const struct OpenSubdiv_Converter* converter);
   int (*getNumEdges)(const struct OpenSubdiv_Converter* converter);
@@ -92,6 +95,11 @@ typedef struct OpenSubdiv_Converter {
                          const int vertex_index,
                          int* vertex_faces);
 
+  // Check whether vertex is to be marked as an infinite sharp.
+  // This is a way to make sharp vertices which are adjacent to a loose edges.
+  bool (*isInfiniteSharpVertex)(const struct OpenSubdiv_Converter* converter,
+                                const int vertex_index);
+
   //////////////////////////////////////////////////////////////////////////////
   // Face-varying data.
 
diff --git a/source/blender/blenkernel/intern/subdiv_converter_mesh.c b/source/blender/blenkernel/intern/subdiv_converter_mesh.c
index 1a2c26b3564..50143dd46e1 100644
--- a/source/blender/blenkernel/intern/subdiv_converter_mesh.c
+++ b/source/blender/blenkernel/intern/subdiv_converter_mesh.c
@@ -31,6 +31,7 @@
 #include "DNA_meshdata_types.h"
 
 #include "BLI_utildefines.h"
+#include "BLI_bitmap.h"
 #include "BLI_math_vector.h"
 
 #include "BKE_customdata.h"
@@ -68,6 +69,26 @@ typedef struct ConverterStorage {
 	 */
 	int *loop_uv_indices;
 	int num_uv_coordinates;
+
+	/* Indexed by coarse mesh elements, gives index of corresponding element
+	 * with ignoring all non-manifold entities.
+	 *
+	 * NOTE: This isn't strictly speaking manifold, this is more like non-loose
+	 * geometry index. As in, index of element as if there were no loose edges
+	 * or vertices in the mesh.
+	 */
+	int *manifold_vertex_index;
+	int *manifold_edge_index;
+	/* Indexed by vertex index from mesh, corresponds to whether this vertex has
+	 * infinite sharpness due to non-manifol topology.
+	 */
+	BLI_bitmap *infinite_sharp_vertices_map;
+	/* Reverse mapping to above. */
+	int *manifold_vertex_index_reverse;
+	int *manifold_edge_index_reverse;
+	/* Number of non-loose elements. */
+	int num_manifold_vertices;
+	int num_manifold_edges;
 } ConverterStorage;
 
 static OpenSubdiv_SchemeType get_scheme_type(
@@ -98,70 +119,81 @@ static int get_num_faces(const OpenSubdiv_Converter *converter)
 static int get_num_edges(const OpenSubdiv_Converter *converter)
 {
 	ConverterStorage *storage = converter->user_data;
-	return storage->mesh->totedge;
+	return storage->num_manifold_edges;
 }
 
-static int get_num_verts(const OpenSubdiv_Converter *converter)
+static int get_num_vertices(const OpenSubdiv_Converter *converter)
 {
 	ConverterStorage *storage = converter->user_data;
-	return storage->mesh->totvert;
+	return storage->num_manifold_vertices;
 }
 
-static int get_num_face_verts(const OpenSubdiv_Converter *converter, int face)
+static int get_num_face_vertices(const OpenSubdiv_Converter *converter,
+                                 int manifold_face_index)
 {
 	ConverterStorage *storage = converter->user_data;
-	return storage->mesh->mpoly[face].totloop;
+	return storage->mesh->mpoly[manifold_face_index].totloop;
 }
 
-static void get_face_verts(const OpenSubdiv_Converter *converter,
-                           int face,
-                           int *face_verts)
+static void get_face_vertices(const OpenSubdiv_Converter *converter,
+                              int manifold_face_index,
+                              int *manifold_face_vertices)
 {
 	ConverterStorage *storage = converter->user_data;
-	const MPoly *mp = &storage->mesh->mpoly[face];
+	const MPoly *poly = &storage->mesh->mpoly[manifold_face_index];
 	const MLoop *mloop = storage->mesh->mloop;
-	for (int loop = 0; loop < mp->totloop; loop++) {
-		face_verts[loop] = mloop[mp->loopstart + loop].v;
+	for (int corner = 0; corner < poly->totloop; corner++) {
+		manifold_face_vertices[corner] = storage->manifold_vertex_index[
+		        mloop[poly->loopstart + corner].v];
 	}
 }
 
 static void get_face_edges(const OpenSubdiv_Converter *converter,
-                           int face,
-                           int *face_edges)
+                           int manifold_face_index,
+                           int *manifold_face_edges)
 {
 	ConverterStorage *storage = converter->user_data;
-	const MPoly *mp = &storage->mesh->mpoly[face];
+	const MPoly *poly = &storage->mesh->mpoly[manifold_face_index];
 	const MLoop *mloop = storage->mesh->mloop;
-	for (int loop = 0; loop < mp->totloop; loop++) {
-		face_edges[loop] = mloop[mp->loopstart + loop].e;
+	for (int corner = 0; corner < poly->totloop; corner++) {
+		manifold_face_edges[corner] =
+		        storage->manifold_edge_index[mloop[poly->loopstart + corner].e];
 	}
 }
 
-static void get_edge_verts(const OpenSubdiv_Converter *converter,
-                           int edge,
-                           int *edge_verts)
+static void get_edge_vertices(const OpenSubdiv_Converter *converter,
+                              int manifold_edge_index,
+                              int *manifold_edge_vertices)
 {
 	ConverterStorage *storage = converter->user_data;
-	const MEdge *me = &storage->mesh->medge[edge];
-	edge_verts[0] = me->v1;
-	edge_verts[1] = me->v2;
+	const int edge_index =
+	        storage->manifold_edge_index_reverse[manifold_edge_index];
+	const MEdge *edge = &storage->mesh->medge[edge_index];
+	manifold_edge_vertices[0] = storage->manifold_vertex_index[edge->v1];
+	manifold_edge_vertices[1] = storage->manifold_vertex_index[edge->v2];
 }
 
-static int get_num_edge_faces(const OpenSubdiv_Converter *converter, int edge)
+static int get_num_edge_faces(const OpenSubdiv_Converter *converter,
+                              int manifold_edge_index)
 {
 	ConverterStorage *storage = converter->user_data;
+	const int edge_index =
+	        storage->manifold_edge_index_reverse[manifold_edge_index];
 #ifdef USE_MESH_ELEMENT_MAPPING
-	return storage->edge_poly_map[edge].count;
+	return storage->edge_poly_map[edge_index].count;
 #else
 	const Mesh *mesh = storage->mesh;
 	const MPoly *mpoly = mesh->mpoly;
 	const MLoop *mloop = mesh->mloop;
 	int num = 0;
-	for (int poly = 0; poly < mesh->totpoly; poly++) {
-		const MPoly *mp = &mpoly[poly];
-		for (int loop = 0; loop < mp->totloop; loop++) {
-			const MLoop *ml = &mloop[mp->loopstart + loop];
-			if (ml->e == edge) {
+	for (int poly_index = 0; poly_index < mesh->totpoly; poly_index++) {
+		const MPoly *poly = &mpoly[poly_index];
+		for (int corner = 0; corner < poly->totloop; corner++) {
+			const MLoop *loop = &mloop[poly->loopstart + corner];
+			if (storage->manifold_edge_index[loop->e] == -1) {
+				continue;
+			}
+			if (loop->e == edge_index) {
 				++num;
 				break;
 			}
@@ -172,25 +204,30 @@ static int get_num_edge_faces(const OpenSubdiv_Converter *converter, int edge)
 }
 
 static void get_edge_faces(const OpenSubdiv_Converter *converter,
-                           int edge,
-                           int *edge_faces)
+                           int manifold_edge_index,
+                           int *manifold_edge_faces)
 {
 	ConverterStorage *storage = converter->user_data;
+	const int edge_index =
+	        storage->manifold_edge_index_reverse[manifold_edge_index];
 #ifdef USE_MESH_ELEMENT_MAPPING
-	memcpy(edge_faces,
-	       storage->edge_poly_map[edge].indices,
-	       sizeof(int) * storage->edge_poly_map[edge].count);
+	memcpy(manifold_edge_faces,
+	       storage->edge_poly_map[edge_index].indices,
+	       sizeof(int) * storage->edge_poly_map[edge_index].count);
 #else
 	const Mesh *mesh = storage->mesh;
 	const MPoly *mpoly = mesh->mpoly;
 	const MLoop *mloop = mesh->mloop;
 	int num = 0;
-	for (int poly = 0; poly < mesh->totpoly; poly++) {
-		const MPoly *mp = &mpoly[poly];
-		for (int loop = 0; loop < mpoly->totloop; loop++) {
-			const MLoop *ml = &mloop[mp->loopstart + loop];
-			if (ml->e == edge) {
-				edge_faces[num++] = poly;
+	for (int poly_index = 0; poly_index < mesh->totpoly; poly_index++) {
+		const MPoly *poly = &mpoly[poly_index];
+		for (int corner = 0; corner < mpoly->totloop; corner++) {
+			const MLoop *loop = &mloop[poly->loopstart + corner];
+			if (storage->manifold_edge_index[loop->e] == -1) {
+				continue;
+			}
+			if (loop->e == edge_index) {
+				manifold_edge_faces[num++] = poly_index;
 				break;
 			}
 		}
@@ -198,26 +235,46 @@ static void get_edge_faces(const OpenSubdiv_Converter *converter,
 #endif
 }
 
-static float get_edge_sharpness(const OpenSubdiv_Converter *converter, int edge)
+static float get_edge_sharpness(const OpenSubdiv_Converter *converter,
+                                int manifold_edge_index)
 {
 	ConverterStorage *storage = converter->user_data;
+	const int edge_index =
+	        storage->manifold_edge_index_reverse[manifold_edge_index];
 	const MEdge *medge = storage->mesh->medge;
-	const float edge_crease =  (float)medge[edge].crease / 255.0f;
+	const float edge_crease =  (float)medge[edge_index].crease / 255.0f;
 	return edge_crease * storage->settings.level;
 }
 
-static int get_num_vert_edges(const OpenSubdiv_Converter *converter, int vert)
+static int get_num_vertex_edges(const OpenSubdiv_Converter *converter,
+                                int manifold_vertex_index)
 {
 	ConverterStorage *storage = converter->user_data;
+	const int vertex_index =
+	        storage->manifold_vertex_index_reverse[manifold_vertex_index];
 #ifdef USE_MESH_ELEMENT_MAPPING
-	return storage->vert_edge_map[vert].count;
+	const int num_vertex_edges = storage->vert_edge_map[vertex_index].count;
+	int num_manifold_vertex_edges = 0;
+	for (int i = 0; i < num_vertex_edges; i++) {
+		const int edge_index = storage->vert_edge_map[vertex_index].indices[i];
+		const int manifold_edge_index =
+		        storage->manifold_edge_index[edge_index];
+		if (manifold_edge_index == -1) {
+			continue;
+		}
+		num_manifold_vertex_edges++;
+	}
+	return num_manifold_vertex_edges;
 #else
 	const Mesh *mesh = storage->mesh;
 	const MEdge *medge = mesh->medge;
 	int num = 0;
-	for (int edge = 0; edge < mesh->totedge; edge++) {
-		const MEdge *me = &medge[edge];
-		if (me->v1 == vert || me->v2 == vert) {
+	for (int edge_index = 0; edge_index < mesh->totedge; edge_index++) {
+		const MEdge *edge = &medge[edge_index];
+		if (storage->manifold_edge_index[edge_index] == -1) {
+			continue;
+		}
+		if (edge->v1 == vertex_index || edge->v2 == vertex_index) {
 			++num;
 		}
 	}
@@ -225,43 +282,61 @@ static int get_num_vert_edges(const OpenSubdiv_Converter *converter, int vert)
 #endif
 }
 
-static void get_vert_edges(const OpenSubdiv_Converter *converter,
-                           int vert,
-                           int *vert_edges)
+static void get_vertex_edges(const OpenSubdiv_Converter *converter,
+                             int manifold_vertex_index,
+                             int *manifold_vertex_edges)
 {
 	ConverterStorage *storage = converter->user_data;
+	const int vertex_index =
+	        storage->manifold_vertex_index_reverse[manifold_vertex_index];
 #ifdef USE_MESH_ELEMENT_MAPPING
-	memcpy(vert_edges,
-	       storage->vert_edge_map[vert].indices,
-	       sizeof(int) * storage->vert_edge_map[vert].count);
+	const int num_vertex_edges = storage->vert_edge_map[vertex_index].count;
+	int num_manifold_vertex_edges = 0;
+	for (int i = 0; i < num_vertex_edges; i++) {
+		const int edge_index = storage->vert_edge_map[vertex_index].indices[i];
+		const int manifold_edge_index =
+		        storage->manifold_edge_index[edge_index];
+		if (manifold_edge_index == -1) {
+			continue;
+		}
+		manifold_vertex_edges[num_manifold_vertex_edges] = manifold_edge_index;
+		num_manifold_vertex_edges++;
+	}
 #else
 	const Mesh *mesh = storage->mesh;
 	const MEdge *medge = mesh->medge;
 	int num = 0;
-	for (int edge = 0; edge < mesh->totedge; edge++) {
-		const MEdge *me = &medge[edge];
-		if (me->v1 == vert || me->v2 == vert) {
-			vert_edges[num++] = edge;
+	for (int edge_index = 0; edge_index < mesh->totedge; edge_index++) {
+		const MEdge *edge = &medge[edge_index];
+		if (storage->manifold_edge_index[edge_index] == -1) {
+			continue;
+		}
+		if (edge->v1 == vertex_index || edge->v2 == vertex_index) {
+			manifold_vertex_edges[num++] =
+			        storage->manifold_edge_index[edge_index];
 		}
 	}
 #endif
 }
 
-static int get_num_vert_faces(const OpenSubdiv_Converter *converter, int vert)
+static int get_num_vertex_faces(const OpenSubdiv_Converter *converter,
+                                int manifold_vertex_index)
 {
 	ConverterStorage *storage = converter->user_data;
+	const int vertex_index =
+	        storage->manifold_vertex_index_reverse[manifold_vertex_index];
 #ifdef USE_MESH_ELEMENT_MAPPING
-	return storage->vert_poly_map[vert].count;
+	return storage->vert_poly_map[vertex_index].count;
 #else
 	const Mesh *mesh = storage->mesh;
 	const MPoly *mpoly = mesh->mpoly;
 	const MLoop *mloop = mesh->mloop;
 	int num = 0;
-	for (int poly = 0; poly < mesh->totpoly; poly++) {
-		const MPoly *mp = &mpoly[poly];
-		for (int loop = 0; loop < mpoly->totloop; loop++) {
-			const MLoop *ml = &mloop[mp->loopstart + loop];
-			if (ml->v == vert) {
+	for (int poly_index = 0; poly_index < mesh->totpoly; poly_index++) {
+		const MPoly *poly = &mpoly[poly_index];
+		for (int corner = 0; corner < mpoly->totloop; corner++) {
+			const MLoop *loop = &mloop[poly->loopstart + corner];
+			if (loop->v == vertex_index) {
 				++num;
 				break;
 			}
@@ -271,26 +346,28 @@ static int get_num_vert_faces(const OpenSubdiv_Converter *converter, int vert)
 #endif
 }
 
-static void get_vert_faces(const OpenSubdiv_Converter *converter,
-                           int vert,
-                           int *vert_faces)
+static void get_vertex_faces(const OpenSubdiv_Converter *converter,
+                            int manifold_vertex_index,
+                            int *manifold_vertex_faces)
 {
 	ConverterStorage *storage = converter->user_data;
+	const int vertex_index =
+	        storage->manifold_vertex_index_reverse[manifold_vertex_index];
 #ifdef USE_MESH_ELEMENT_MAPPING
-	memcpy(vert_faces,
-	       storage->vert_poly_map[vert].indices,
-	       sizeof(int) * storage->vert_poly_map[vert].count);
+	memcpy(manifold_vertex_faces,
+	       storage->vert_poly_map[vertex_index].indices,
+	       sizeof(int) * storage->vert_poly_map[vertex_index].count);
 #else
 	const Mesh *mesh = storage->mesh;
 	const MPoly *mpoly = mesh->mpoly;
 	const MLoop *mloop = mesh->mloop;
 	int num = 0;
-	for (int poly = 0; poly < mesh->totpoly; poly++) {
-		const MPoly *mp = &mpoly[poly];
-		for (int loop = 0; loop < mpoly->totloop; loop++) {
-			const MLoop *ml = &mloop[mp->loopstart + loop];
-			if (ml->v == vert) {
-				vert_faces[num++] = poly;
+	for (int poly_index = 0; poly_index < mesh->totpoly; poly_index++) {
+		const MPoly *poly = &mpoly[poly_index];
+		for (int corner = 0; corner < mpoly->totloop; corner++) {
+			const MLoop *loop = &mloop[poly->loopstart + corner];
+			if (loop->v == vertex_index) {
+				manifold_vertex_faces[num++] = poly_index;
 				break;
 			}
 		}
@@ -298,6 +375,16 @@ static void get_vert_faces(const OpenSubdiv_Converter *converter,
 #endif
 }
 
+static bool is_infinite_sharp_vertex(const OpenSubdiv_Converter *converter,
+                                     int manifold_vertex_index)
+{
+	ConverterStorage *storage = converter->user_data;
+	const int vertex_index =
+	        storage->manifold_vertex_index_reverse[manifold_vertex_index];
+	return BLI_BITMAP_TEST_BOOL(storage->infinite_sharp_vertices_map,
+	                            vertex_index);
+}
+
 static int get_num_uv_layers(const OpenSubdiv_Converter *converter)
 {
 	ConverterStorage *storage = converter->user_data;
@@ -382,6 +469,11 @@ static void free_user_data(const OpenSubdiv_Converter *converter)
 	MEM_freeN(user_data->edge_poly_map);
 	MEM_freeN(user_data->edge_poly_mem);
 #endif
+	MEM_freeN(user_data->manifold_vertex_index);
+	MEM_freeN(user_data->manifold_edge_index);
+	MEM_freeN(user_data->infinite_sharp_vertices_map);
+	MEM_freeN(user_data->manifold_vertex_index_reverse);
+	MEM_freeN(user_data->manifold_edge_index_reverse);
 	MEM_freeN(user_data);
 }
 
@@ -393,21 +485,22 @@ static void init_functions(OpenSubdiv_Converter *converter)
 
 	converter->getNumFaces = get_num_faces;
 	converter->getNumEdges = get_num_edges;
-	converter->getNumVertices = get_num_verts;
+	converter->getNumVertices = get_num_vertices;
 
-	converter->getNumFaceVertices = get_num_face_verts;
-	converter->getFaceVertices = get_face_verts;
+	converter->getNumFaceVertices = get_num_face_vertices;
+	converter->getFaceVertices = get_face_vertices;
 	converter->getFaceEdges = get_face_edges;
 
-	converter->getEdgeVertices = get_edge_verts;
+	converter->getEdgeVertices = get_edge_vertices;
 	converter->getNumEdgeFaces = get_num_edge_faces;
 	converter->getEdgeFaces = get_edge_faces;
 	converter->getEdgeSharpness = get_edge_sharpness;
 
-	converter->getNumVertexEdges = get_num_vert_edges;
-	converter->getVertexEdges = get_vert_edges;
-	converter->getNumVertexFaces = get_num_vert_faces;
-	converter->getVertexFaces = get_vert_faces;
+	converter->getNumVertexEdges = get_num_vertex_edges;
+	converter->getVertexEdges = get_vertex_edges;
+	converter->getNumVertexFaces = get_num_vertex_faces;
+	converter->getVertexFaces = get_vertex_faces;
+	converter->isInfiniteSharpVertex = is_infinite_sharp_vertex;
 
 	converter->getNumUVLayers = get_num_uv_layers;
 	converter->precalcUVLayer = precalc_uv_layer;
@@ -444,6 +537,74 @@ static void create_element_maps_if_needed(ConverterStorage *storage)
 #endif
 }
 
+static void initialize_manifold_index_array(const BLI_bitmap *used_map,
+                                            const int num_elements,
+                                            int **indices_r,
+                                            int **indices_reverse_r,
+                                            int *num_manifold_elements_r)
+{
+	int *indices = MEM_malloc_arrayN(
+	        num_elements, sizeof(int), "manifold indices");
+	int *indices_reverse = MEM_malloc_arrayN(
+	        num_elements, sizeof(int), "manifold indices reverse");
+	int offset = 0;
+	for (int i = 0; i < num_elements; i++) {
+		if (BLI_BITMAP_TEST_BOOL(used_map, i)) {
+			indices[i] = i - offset;
+			indices_reverse[i - offset] = i;
+		}
+		else {
+			indices[i] = -1;
+			offset++;
+		}
+	}
+	*indices_r = indices;
+	*indices_reverse_r = indices_reverse;
+	*num_manifold_elements_r = num_elements - offset;
+}
+
+static void initialize_manifold_indices(ConverterStorage *storage)
+{
+	const Mesh *mesh = storage->mesh;
+	const MEdge *medge = mesh->medge;
+	const MLoop *mloop = mesh->mloop;
+	const MPoly *mpoly = mesh->mpoly;
+	/* Set bits of elements which are not loose. */
+	BLI_bitmap *vert_used_map = BLI_BITMAP_NEW(mesh->totvert, "vert used map");
+	BLI_bitmap *edge_used_map = BLI_BITMAP_NEW(mesh->totedge, "edge used map");
+	for (int poly_index = 0; poly_index < mesh->totpoly; poly_index++) {
+		const MPoly *poly = &mpoly[poly_index];
+		for (int corner = 0; corner < poly->totloop; corner++) {
+			const MLoop *loop = &mloop[poly->loopstart + corner];
+			BLI_BITMAP_ENABLE(vert_used_map, loop->v);
+			BLI_BITMAP_ENABLE(edge_used_map, loop->e);
+		}
+	}
+	initialize_manifold_index_array(vert_used_map,
+	                                mesh->totvert,
+	                                &storage->manifold_vertex_index,
+	                                &storage->manifold_vertex_index_reverse,
+	                                &storage->num_manifold_vertices);
+	initialize_manifold_index_array(edge_used_map,
+	                                mesh->totedge,
+	                                &storage->manifold_edge_index,
+	                                &storage->manifold_edge_index_reverse,
+	                                &storage->num_manifold_edges);
+	/* Initialize infinite sharp mapping. */
+	storage->infinite_sharp_vertices_map =
+	        BLI_BITMAP_NEW(mesh->totvert, "vert used map");
+	for (int edge_index = 0; edge_index < mesh->totedge; edge_index++) {
+		if (!BLI_BITMAP_TEST_BOOL(edge_used_map, edge_index)) {
+			const MEdge *edge = &medge[edge_index];
+			BLI_BITMAP_ENABLE(storage->infinite_sharp_vertices_map, edge->v1);
+			BLI_BITMAP_ENABLE(storage->infinite_sharp_vertices_map, edge->v2);
+		}
+	}
+	/* Free working variables. */
+	MEM_freeN(vert_used_map);
+	MEM_freeN(edge_used_map);
+}
+
 static void init_user_data(OpenSubdiv_Converter *converter,
                            const SubdivSettings *settings,
                            const Mesh *mesh)
@@ -454,6 +615,7 @@ static void init_user_data(OpenSubdiv_Converter *converter,
 	user_data->mesh = mesh;
 	user_data->loop_uv_indices = NULL;
 	create_element_maps_if_needed(user_data);
+	initialize_manifold_indices(user_data);
 	converter->user_data = user_data;
 }
 #endif
diff --git a/source/blender/blenkernel/intern/subdiv_eval.c b/source/blender/blenkernel/intern/subdiv_eval.c
index 4c2ed07cdfa..7c90a2a25cd 100644
--- a/source/blender/blenkernel/intern/subdiv_eval.c
+++ b/source/blender/blenkernel/intern/subdiv_eval.c
@@ -33,10 +33,13 @@
 #include "DNA_meshdata_types.h"
 
 #include "BLI_utildefines.h"
+#include "BLI_bitmap.h"
 #include "BLI_math_vector.h"
 
 #include "BKE_customdata.h"
 
+#include "MEM_guardedalloc.h"
+
 #ifdef WITH_OPENSUBDIV
 #  include "opensubdiv_evaluator_capi.h"
 #  include "opensubdiv_topology_refiner_capi.h"
@@ -60,6 +63,42 @@ void BKE_subdiv_eval_begin(Subdiv *subdiv)
 }
 
 #ifdef WITH_OPENSUBDIV
+static void set_coarse_positions(Subdiv *subdiv, const Mesh *mesh)
+{
+	const MVert *mvert = mesh->mvert;
+	const MLoop *mloop = mesh->mloop;
+	const MPoly *mpoly = mesh->mpoly;
+	/* Mark vertices which needs new coordinates. */
+	/* TODO(sergey): This is annoying to calculate this on every update,
+	 * maybe it's better to cache this mapping. Or make it possible to have
+	 * OpenSubdiv's vertices match mesh ones?
+	 */
+	BLI_bitmap *vertex_used_map =
+	        BLI_BITMAP_NEW(mesh->totvert, "vert used map");
+	for (int poly_index = 0; poly_index < mesh->totpoly; poly_index++) {
+		const MPoly *poly = &mpoly[poly_index];
+		for (int corner = 0; corner < poly->totloop; corner++) {
+			const MLoop *loop = &mloop[poly->loopstart + corner];
+			BLI_BITMAP_ENABLE(vertex_used_map, loop->v);
+		}
+	}
+	for (int vertex_index = 0, manifold_veretx_index = 0;
+	     vertex_index < mesh->totvert;
+	     vertex_index++)
+	{
+		if (!BLI_BITMAP_TEST_BOOL(vertex_used_map, vertex_index)) {
+			continue;
+		}
+		const MVert *vertex = &mvert[vertex_index];
+		subdiv->evaluator->setCoarsePositions(
+		        subdiv->evaluator,
+		        vertex->co,
+		        manifold_veretx_index, 1);
+		manifold_veretx_index++;
+	}
+	MEM_freeN(vertex_used_map);
+}
+
 static void set_face_varying_data_from_uv(Subdiv *subdiv,
                                           const MLoopUV *mloopuv,
                                           const int layer_index)
@@ -94,12 +133,7 @@ void BKE_subdiv_eval_update_from_mesh(Subdiv *subdiv, const Mesh *mesh)
 #ifdef WITH_OPENSUBDIV
 	BKE_subdiv_eval_begin(subdiv);
 	/* Set coordinates of base mesh vertices. */
-	subdiv->evaluator->setCoarsePositionsFromBuffer(
-	        subdiv->evaluator,
-	        mesh->mvert,
-	        offsetof(MVert, co),
-	        sizeof(MVert),
-	        0, mesh->totvert);
+	set_coarse_positions(subdiv, mesh);
 	/* Set face-varyign data to UV maps. */
 	const int num_uv_layers =
 	        CustomData_number_of_layers(&mesh->ldata, CD_MLOOPUV);
diff --git a/source/blender/blenkernel/intern/subdiv_mesh.c b/source/blender/blenkernel/intern/subdiv_mesh.c
index 33813905105..befd2b9847c 100644
--- a/source/blender/blenkernel/intern/subdiv_mesh.c
+++ b/source/blender/blenkernel/intern/subdiv_mesh.c
@@ -33,6 +33,7 @@
 
 #include "DNA_mesh_types.h"
 #include "DNA_meshdata_types.h"
+#include "DNA_key_types.h"
 
 #include "BLI_alloca.h"
 #include "BLI_bitmap.h"
@@ -40,6 +41,7 @@
 #include "BLI_task.h"
 
 #include "BKE_mesh.h"
+#include "BKE_key.h"
 
 #include "MEM_guardedalloc.h"
 
@@ -183,6 +185,7 @@ static void subdiv_mesh_ctx_count(SubdivMeshContext *ctx)
 	ctx->num_subdiv_vertices = coarse_mesh->totvert;
 	ctx->num_subdiv_edges =
 	        coarse_mesh->totedge * (num_subdiv_vertices_per_coarse_edge + 1);
+	/* Calculate extra vertices and edges createdd by non-loose geometry. */
 	for (int poly_index = 0; poly_index < coarse_mesh->totpoly; poly_index++) {
 		const MPoly *coarse_poly = &coarse_mpoly[poly_index];
 		const int num_ptex_faces_per_poly =
@@ -224,6 +227,12 @@ static void subdiv_mesh_ctx_count(SubdivMeshContext *ctx)
 			        num_polys_per_ptex_get(no_quad_patch_resolution);
 		}
 	}
+	/* Calculate extra edges createdd by loose edges. */
+	for (int edge_index = 0; edge_index < coarse_mesh->totedge; edge_index++) {
+		if (!BLI_BITMAP_TEST_BOOL(ctx->coarse_edges_used_map, edge_index)) {
+			ctx->num_subdiv_vertices += num_subdiv_vertices_per_coarse_edge;
+		}
+	}
 	ctx->num_subdiv_loops = ctx->num_subdiv_polygons * 4;
 }
 
@@ -2186,6 +2195,164 @@ static void subdiv_create_polys(SubdivMeshContext *ctx, int poly_index)
 }
 
 /* =============================================================================
+ * Loose elements subdivision process.
+ */
+
+static void subdiv_create_loose_vertices_task(
+        void *__restrict userdata,
+        const int vertex_index,
+        const ParallelRangeTLS *__restrict UNUSED(tls))
+{
+	SubdivMeshContext *ctx = userdata;
+	if (BLI_BITMAP_TEST_BOOL(ctx->coarse_vertices_used_map, vertex_index)) {
+		/* Vertex is not loose, was handled when handling polygons. */
+		return;
+	}
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MVert *coarse_mvert = coarse_mesh->mvert;
+	const MVert *coarse_vertex = &coarse_mvert[vertex_index];
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	MVert *subdiv_mvert = subdiv_mesh->mvert;
+	MVert *subdiv_vertex = &subdiv_mvert[
+	        ctx->vertices_corner_offset + vertex_index];
+	subdiv_vertex_data_copy(ctx, coarse_vertex, subdiv_vertex);
+}
+
+/* Get neighbor edges of the given one.
+ * - neighbors[0] is an edge adjacent to edge->v1.
+ * - neighbors[1] is an edge adjacent to edge->v1.
+ */
+static void find_edge_neighbors(const SubdivMeshContext *ctx,
+                                const MEdge *edge,
+                                const MEdge *neighbors[2])
+{
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MEdge *coarse_medge = coarse_mesh->medge;
+	neighbors[0] = NULL;
+	neighbors[1] = NULL;
+	for (int edge_index = 0; edge_index < coarse_mesh->totedge; edge_index++) {
+		if (BLI_BITMAP_TEST_BOOL(ctx->coarse_edges_used_map, edge_index)) {
+			continue;
+		}
+		const MEdge *current_edge = &coarse_medge[edge_index];
+		if (current_edge == edge) {
+			continue;
+		}
+		if (ELEM(edge->v1, current_edge->v1, current_edge->v2)) {
+			neighbors[0] = current_edge;
+		}
+		if (ELEM(edge->v2, current_edge->v1, current_edge->v2)) {
+			neighbors[1] = current_edge;
+		}
+	}
+}
+
+static void points_for_loose_edges_interpolation_get(
+        SubdivMeshContext *ctx,
+        const MEdge *coarse_edge,
+        const MEdge *neighbors[2],
+		float points_r[4][3])
+{
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MVert *coarse_mvert = coarse_mesh->mvert;
+	/* Middle points corresponds to the edge. */
+	copy_v3_v3(points_r[1], coarse_mvert[coarse_edge->v1].co);
+	copy_v3_v3(points_r[2], coarse_mvert[coarse_edge->v2].co);
+	/* Start point, duplicate from edge start if no neighbor. */
+	if (neighbors[0] != NULL) {
+		if (neighbors[0]->v1 == coarse_edge->v1) {
+			copy_v3_v3(points_r[0], coarse_mvert[neighbors[0]->v2].co);
+		}
+		else {
+			copy_v3_v3(points_r[0], coarse_mvert[neighbors[0]->v1].co);
+		}
+	}
+	else {
+		sub_v3_v3v3(points_r[0], points_r[1], points_r[2]);
+		add_v3_v3(points_r[0], points_r[1]);
+	}
+	/* End point, duplicate from edge end if no neighbor. */
+	if (neighbors[1] != NULL) {
+		if (neighbors[1]->v1 == coarse_edge->v2) {
+			copy_v3_v3(points_r[3], coarse_mvert[neighbors[1]->v2].co);
+		}
+		else {
+			copy_v3_v3(points_r[3], coarse_mvert[neighbors[1]->v1].co);
+		}
+	}
+	else {
+		sub_v3_v3v3(points_r[3], points_r[2], points_r[1]);
+		add_v3_v3(points_r[3], points_r[2]);
+	}
+}
+
+static void subdiv_create_vertices_of_loose_edges_task(
+        void *__restrict userdata,
+        const int edge_index,
+        const ParallelRangeTLS *__restrict UNUSED(tls))
+{
+	SubdivMeshContext *ctx = userdata;
+	if (BLI_BITMAP_TEST_BOOL(ctx->coarse_edges_used_map, edge_index)) {
+		/* Vertex is not loose, was handled when handling polygons. */
+		return;
+	}
+	const int resolution = ctx->settings->resolution;
+	const int resolution_1 = resolution - 1;
+	const float inv_resolution_1 = 1.0f / (float)resolution_1;
+	const int num_subdiv_vertices_per_coarse_edge = resolution - 2;
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MEdge *coarse_edge = &coarse_mesh->medge[edge_index];
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	MVert *subdiv_mvert = subdiv_mesh->mvert;
+	/* Find neighbors of the current loose edge. */
+	const MEdge *neighbors[2];
+	find_edge_neighbors(ctx, coarse_edge, neighbors);
+	/* Get points for b-spline interpolation. */
+	float points[4][3];
+	points_for_loose_edges_interpolation_get(
+	        ctx, coarse_edge, neighbors, points);
+	/* Subdivion verticies which corresponds to edge's v1 and v2. */
+	MVert *subdiv_v1 = &subdiv_mvert[
+	        ctx->vertices_corner_offset + coarse_edge->v1];
+	MVert *subdiv_v2 = &subdiv_mvert[
+	        ctx->vertices_corner_offset + coarse_edge->v2];
+	/* First subdivided inner vertex of the edge.  */
+	MVert *subdiv_start_vertex = &subdiv_mvert[
+	        ctx->vertices_edge_offset +
+	        edge_index * num_subdiv_vertices_per_coarse_edge];
+	/* Perform interpolation. */
+	for (int i = 0; i < resolution; i++) {
+		const float u = i * inv_resolution_1;
+		float weights[4];
+		key_curve_position_weights(u, weights, KEY_BSPLINE);
+
+		MVert *subdiv_vertex;
+		if (i == 0) {
+			subdiv_vertex = subdiv_v1;
+		}
+		else if (i == resolution - 1) {
+			subdiv_vertex = subdiv_v2;
+		}
+		else {
+			subdiv_vertex = &subdiv_start_vertex[i - 1];
+		}
+		interp_v3_v3v3v3v3(subdiv_vertex->co,
+		                   points[0],
+		                   points[1],
+		                   points[2],
+		                   points[3],
+		                   weights);
+		/* Reset flags and such. */
+		subdiv_vertex->flag = 0;
+		subdiv_vertex->bweight = 0.0f;
+		/* Reset normal. */
+		subdiv_vertex->no[0] = 0.0f;
+		subdiv_vertex->no[1] = 0.0f;
+		subdiv_vertex->no[2] = 1.0f;
+	}
+}
+
+/* =============================================================================
  * Subdivision process entry points.
  */
 
@@ -2245,6 +2412,14 @@ Mesh *BKE_subdiv_to_mesh(
 	                        &ctx,
 	                        subdiv_eval_task,
 	                        &parallel_range_settings);
+	BLI_task_parallel_range(0, coarse_mesh->totvert,
+	                        &ctx,
+	                        subdiv_create_loose_vertices_task,
+	                        &parallel_range_settings);
+	BLI_task_parallel_range(0, coarse_mesh->totedge,
+	                        &ctx,
+	                        subdiv_create_vertices_of_loose_edges_task,
+	                        &parallel_range_settings);
 	BLI_task_parallel_range(0, coarse_mesh->totedge,
 	                        &ctx,
 	                        subdiv_create_boundary_edges_task,