Subsurf: Begin new subdivision surface module

The idea is to use this as a replacement of old CCG, now it is
based on OpenSubdiv. The goal is to reduce any possible overhead
which was happening with OpenSubdiv used by CCG.

Currently implemented/supported:

- Creation from mesh, including topology on OpenSubdiv side,
  its refinement.

- Evaluation of limit point, first order derivatives, normal,
  and face-varying data for individual coarse position.

- Evaluation of whole patches.

  Currently not optimized, uses evaluation of individual coarse
  positions.

- Creation of Mesh from subdiv, with all geometry being real:
  all mvert, medge, mloop, and mpoly.

  This includes custom data interpolation, but all faces currently
  are getting separated (they are converted to ptex patches, which
  we need to weld back together).

Still need to support lighter weights grids and such, but this
is already a required part to have subsurf working in the middle
of modifier stack.

Annoying part is ifdef all over the place, to keep it compilable
when OpenSubdiv is disabled. More cleaner approach would be to
have stub API for OpenSubdiv, so everything gets ifdef-ed in a
much fewer places.

1 files changed, 910 insertions, 0 deletions

diff --git a/source/blender/blenkernel/intern/subdiv_mesh.c b/source/blender/blenkernel/intern/subdiv_mesh.c
new file mode 100644
index 00000000000..8a58605ba19
--- /dev/null
+++ b/source/blender/blenkernel/intern/subdiv_mesh.c
@@ -0,0 +1,910 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2018 by Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Sergey Sharybin.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+/** \file blender/blenkernel/intern/subdiv_mesh.c
+ *  \ingroup bke
+ */
+
+#include "BKE_subdiv.h"
+
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+
+#include "BLI_alloca.h"
+#include "BLI_math_vector.h"
+#include "BLI_task.h"
+
+#include "BKE_mesh.h"
+
+/* TODO(sergey): Somehow move this to subdiv code? */
+static int mpoly_ptex_faces_count_get(const MPoly *mp)
+{
+	if (mp->totloop == 4) {
+		return 1;
+	}
+	else {
+		return mp->totloop;
+	}
+}
+
+static int num_edges_per_ptex_get(const int resolution)
+{
+	return 2 * (resolution - 1) * resolution;
+}
+
+static int num_polys_per_ptex_get(const int resolution)
+{
+	return (resolution - 1) * (resolution - 1);
+}
+
+typedef struct SubdivMeshContext {
+	const Mesh *coarse_mesh;
+	Subdiv *subdiv;
+	Mesh *subdiv_mesh;
+	const SubdivToMeshSettings *settings;
+	/* Cached custom data arrays for fastter access. */
+	int *vert_origindex;
+	int *edge_origindex;
+	int *loop_origindex;
+	int *poly_origindex;
+	/* UV layers interpolation. */
+	int num_uv_layers;
+	MLoopUV *uv_layers[MAX_MTFACE];
+} SubdivMeshContext;
+
+typedef struct LoopsOfPtex {
+	/* First loop of the ptex, starts at ptex (0, 0) and goes in u direction. */
+	const MLoop *first_loop;
+	/* Last loop of the ptex, starts at ptex (0, 0) and goes in v direction. */
+	const MLoop *last_loop;
+	/* For quad coarse faces only. */
+	const MLoop *second_loop;
+	const MLoop *third_loop;
+} LoopsOfPtex;
+
+static void loops_of_ptex_get(
+        const SubdivMeshContext *ctx,
+        LoopsOfPtex *loops_of_ptex,
+        const MPoly *coarse_poly,
+        const int ptex_face_index)
+{
+	const MLoop *coarse_mloop = ctx->coarse_mesh->mloop;
+	const int first_ptex_loop_index = coarse_poly->loopstart + ptex_face_index;
+	/* Loop which look in the (opposite) V direction of the current
+	 * ptex face.
+	 *
+	 * TOOD(sergey): Get rid of using module on every iteration.
+	 */
+	const int last_ptex_loop_index =
+	        coarse_poly->loopstart +
+	        (ptex_face_index + coarse_poly->totloop - 1) % coarse_poly->totloop;
+	loops_of_ptex->first_loop = &coarse_mloop[first_ptex_loop_index];
+	loops_of_ptex->last_loop = &coarse_mloop[last_ptex_loop_index];
+	if (coarse_poly->totloop == 4) {
+		loops_of_ptex->second_loop = loops_of_ptex->first_loop + 1;
+		loops_of_ptex->third_loop = loops_of_ptex->first_loop + 2;
+	}
+	else {
+		loops_of_ptex->second_loop = NULL;
+		loops_of_ptex->third_loop = NULL;
+	}
+}
+
+typedef struct EdgesOfPtex {
+	/* First edge of the ptex, starts at ptex (0, 0) and goes in u direction. */
+	const MEdge *first_edge;
+	/* Last edge of the ptex, starts at ptex (0, 0) and goes in v direction. */
+	const MEdge *last_edge;
+	/* For quad coarse faces only. */
+	const MEdge *second_edge;
+	const MEdge *third_edge;
+} EdgesOfPtex;
+
+static void edges_of_ptex_get(
+        const SubdivMeshContext *ctx,
+        EdgesOfPtex *edges_of_ptex,
+        const MPoly *coarse_poly,
+        const int ptex_face_index)
+{
+	const MEdge *coarse_medge = ctx->coarse_mesh->medge;
+	LoopsOfPtex loops_of_ptex;
+	loops_of_ptex_get(ctx, &loops_of_ptex, coarse_poly, ptex_face_index);
+	edges_of_ptex->first_edge = &coarse_medge[loops_of_ptex.first_loop->e];
+	edges_of_ptex->last_edge = &coarse_medge[loops_of_ptex.last_loop->e];
+	if (coarse_poly->totloop == 4) {
+		edges_of_ptex->second_edge =
+		        &coarse_medge[loops_of_ptex.second_loop->e];
+		edges_of_ptex->third_edge =
+		        &coarse_medge[loops_of_ptex.third_loop->e];
+	}
+	else {
+		edges_of_ptex->second_edge = NULL;
+		edges_of_ptex->third_edge = NULL;
+	}
+}
+
+/* TODO(sergey): Somehow de-duplicate with loops storage, without too much
+ * exception cases all over the code.
+ */
+
+typedef struct VerticesForInterpolation {
+	/* This field points to a vertex data which is to be used for interpolation.
+	 * The idea is to avoid unnecessary allocations for regular faces, where
+	 * we can simply 
+	 */
+	const CustomData *vertex_data;
+	/* Vertices data calculated for ptex corners. There are always 4 elements
+	 * in this custom data, aligned the following way:
+	 *
+	 *   index 0 -> uv (0, 0)
+	 *   index 1 -> uv (0, 1)
+	 *   index 2 -> uv (1, 1)
+	 *   index 3 -> uv (1, 0)
+	 *
+	 * Is allocated for non-regular faces (triangles and n-gons).
+	 */
+	CustomData vertex_data_storage;
+	bool vertex_data_storage_allocated;
+	/* Infices within vertex_data to interpolate for. The indices are aligned
+	 * with uv coordinates in a similar way as indices in loop_data_storage.
+	 */
+	int vertex_indices[4];
+} VerticesForInterpolation;
+
+static void vertex_interpolation_init(
+        const SubdivMeshContext *ctx,
+        VerticesForInterpolation *vertex_interpolation,
+		const MPoly *coarse_poly)
+{
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MLoop *coarse_mloop = coarse_mesh->mloop;
+	if (coarse_poly->totloop == 4) {
+		vertex_interpolation->vertex_data = &coarse_mesh->vdata;
+		vertex_interpolation->vertex_indices[0] =
+		        coarse_mloop[coarse_poly->loopstart + 0].v;
+		vertex_interpolation->vertex_indices[1] =
+		        coarse_mloop[coarse_poly->loopstart + 1].v;
+		vertex_interpolation->vertex_indices[2] =
+		        coarse_mloop[coarse_poly->loopstart + 2].v;
+		vertex_interpolation->vertex_indices[3] =
+		        coarse_mloop[coarse_poly->loopstart + 3].v;
+		vertex_interpolation->vertex_data_storage_allocated = false;
+	}
+	else {
+		vertex_interpolation->vertex_data =
+		        &vertex_interpolation->vertex_data_storage;
+		/* Allocate storage for loops corresponding to ptex corners. */
+		CustomData_copy(&ctx->coarse_mesh->vdata,
+		                &vertex_interpolation->vertex_data_storage,
+		                CD_MASK_EVERYTHING,
+		                CD_CALLOC,
+		                4);
+		/* Initialize indices. */
+		vertex_interpolation->vertex_indices[0] = 0;
+		vertex_interpolation->vertex_indices[1] = 1;
+		vertex_interpolation->vertex_indices[2] = 2;
+		vertex_interpolation->vertex_indices[3] = 3;
+		vertex_interpolation->vertex_data_storage_allocated = true;
+		/* Interpolate center of poly right away, it stays unchanged for all
+		 * ptex faces.
+		 */
+		const float weight = 1.0f / (float)coarse_poly->totloop;
+		float *weights = BLI_array_alloca(weights, coarse_poly->totloop);
+		int *indices = BLI_array_alloca(indices, coarse_poly->totloop);
+		for (int i = 0; i < coarse_poly->totloop; ++i) {
+			weights[i] = weight;
+			indices[i] = coarse_poly->loopstart + i;
+		}
+		CustomData_interp(&coarse_mesh->vdata,
+		                  &vertex_interpolation->vertex_data_storage,
+		                  indices,
+		                  weights, NULL,
+		                  coarse_poly->totloop,
+		                  2);
+	}
+}
+
+static void vertex_interpolation_from_ptex(
+        const SubdivMeshContext *ctx,
+        VerticesForInterpolation *vertex_interpolation,
+        const MPoly *coarse_poly,
+        const int ptex_face_index)
+{
+	if (coarse_poly->totloop == 4) {
+		/* Nothing to do, all indices and data is already assigned. */
+	} else {
+		const CustomData *vertex_data = &ctx->coarse_mesh->vdata;
+		const Mesh *coarse_mesh = ctx->coarse_mesh;
+		const MLoop *coarse_mloop = coarse_mesh->mloop;
+		LoopsOfPtex loops_of_ptex;
+		loops_of_ptex_get(ctx, &loops_of_ptex, coarse_poly, ptex_face_index);
+		/* Ptex face corner corresponds to a poly loop with same index. */
+		CustomData_copy_data(
+		        vertex_data,
+		        &vertex_interpolation->vertex_data_storage,
+		        coarse_mloop[coarse_poly->loopstart + ptex_face_index].v,
+		        0,
+		        1);
+		/* Interpolate remaining ptex face corners, which hits loops
+		 * middle points.
+		 *
+		 * TODO(sergey): Re-use one of interpolation results from previous
+		 * iteration.
+		 */
+		const float weights[2] = {0.5f, 0.5f};
+		const int first_indices[2] = {
+		        coarse_mloop[loops_of_ptex.first_loop - coarse_mloop].v,
+		        coarse_mloop[(loops_of_ptex.first_loop + 1 - coarse_mloop) %
+		                coarse_poly->totloop].v};
+		const int last_indices[2] = {
+		        coarse_mloop[loops_of_ptex.last_loop - coarse_mloop].v,
+		        coarse_mloop[loops_of_ptex.first_loop - coarse_mloop].v};
+		CustomData_interp(vertex_data,
+		                  &vertex_interpolation->vertex_data_storage,
+		                  first_indices,
+		                  weights, NULL,
+		                  2,
+		                  1);
+		CustomData_interp(vertex_data,
+		                  &vertex_interpolation->vertex_data_storage,
+		                  last_indices,
+		                  weights, NULL,
+		                  2,
+		                  3);
+	}
+}
+
+static void vertex_interpolation_end(
+        VerticesForInterpolation *vertex_interpolation)
+{
+	if (vertex_interpolation->vertex_data_storage_allocated) {
+		CustomData_free(&vertex_interpolation->vertex_data_storage, 4);
+	}
+}
+
+typedef struct LoopsForInterpolation {
+ /* This field points to a loop data which is to be used for interpolation.
+	 * The idea is to avoid unnecessary allocations for regular faces, where
+	 * we can simply 
+	 */
+	const CustomData *loop_data;
+	/* Loops data calculated for ptex corners. There are always 4 elements
+	 * in this custom data, aligned the following way:
+	 *
+	 *   index 0 -> uv (0, 0)
+	 *   index 1 -> uv (0, 1)
+	 *   index 2 -> uv (1, 1)
+	 *   index 3 -> uv (1, 0)
+	 *
+	 * Is allocated for non-regular faces (triangles and n-gons).
+	 */
+	CustomData loop_data_storage;
+	bool loop_data_storage_allocated;
+	/* Infices within loop_data to interpolate for. The indices are aligned with
+	 * uv coordinates in a similar way as indices in loop_data_storage.
+	 */
+	int loop_indices[4];
+} LoopsForInterpolation;
+
+static void loop_interpolation_init(
+        const SubdivMeshContext *ctx,
+        LoopsForInterpolation *loop_interpolation,
+		const MPoly *coarse_poly)
+{
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	if (coarse_poly->totloop == 4) {
+		loop_interpolation->loop_data = &coarse_mesh->ldata;
+		loop_interpolation->loop_indices[0] = coarse_poly->loopstart + 0;
+		loop_interpolation->loop_indices[1] = coarse_poly->loopstart + 1;
+		loop_interpolation->loop_indices[2] = coarse_poly->loopstart + 2;
+		loop_interpolation->loop_indices[3] = coarse_poly->loopstart + 3;
+		loop_interpolation->loop_data_storage_allocated = false;
+	}
+	else {
+		loop_interpolation->loop_data = &loop_interpolation->loop_data_storage;
+		/* Allocate storage for loops corresponding to ptex corners. */
+		CustomData_copy(&ctx->coarse_mesh->ldata,
+		                &loop_interpolation->loop_data_storage,
+		                CD_MASK_EVERYTHING,
+		                CD_CALLOC,
+		                4);
+		/* Initialize indices. */
+		loop_interpolation->loop_indices[0] = 0;
+		loop_interpolation->loop_indices[1] = 1;
+		loop_interpolation->loop_indices[2] = 2;
+		loop_interpolation->loop_indices[3] = 3;
+		loop_interpolation->loop_data_storage_allocated = true;
+		/* Interpolate center of poly right away, it stays unchanged for all
+		 * ptex faces.
+		 */
+		const float weight = 1.0f / (float)coarse_poly->totloop;
+		float *weights = BLI_array_alloca(weights, coarse_poly->totloop);
+		int *indices = BLI_array_alloca(indices, coarse_poly->totloop);
+		for (int i = 0; i < coarse_poly->totloop; ++i) {
+			weights[i] = weight;
+			indices[i] = coarse_poly->loopstart + i;
+		}
+		CustomData_interp(&coarse_mesh->ldata,
+		                  &loop_interpolation->loop_data_storage,
+		                  indices,
+		                  weights, NULL,
+		                  coarse_poly->totloop,
+		                  2);
+	}
+}
+
+static void loop_interpolation_from_ptex(
+        const SubdivMeshContext *ctx,
+        LoopsForInterpolation *loop_interpolation,
+        const MPoly *coarse_poly,
+        const int ptex_face_index)
+{
+	if (coarse_poly->totloop == 4) {
+		/* Nothing to do, all indices and data is already assigned. */
+	} else {
+		const CustomData *loop_data = &ctx->coarse_mesh->ldata;
+		const Mesh *coarse_mesh = ctx->coarse_mesh;
+		const MLoop *coarse_mloop = coarse_mesh->mloop;
+		LoopsOfPtex loops_of_ptex;
+		loops_of_ptex_get(ctx, &loops_of_ptex, coarse_poly, ptex_face_index);
+		/* Ptex face corner corresponds to a poly loop with same index. */
+		CustomData_copy_data(loop_data,
+		                     &loop_interpolation->loop_data_storage,
+		                     coarse_poly->loopstart + ptex_face_index,
+		                     0,
+		                     1);
+		/* Interpolate remaining ptex face corners, which hits loops
+		 * middle points.
+		 *
+		 * TODO(sergey): Re-use one of interpolation results from previous
+		 * iteration.
+		 */
+		const float weights[2] = {0.5f, 0.5f};
+		const int first_indices[2] = {
+		        loops_of_ptex.first_loop - coarse_mloop,
+		        (loops_of_ptex.first_loop + 1 - coarse_mloop) %
+		                coarse_poly->totloop};
+		const int last_indices[2] = {
+		        loops_of_ptex.last_loop - coarse_mloop,
+		        loops_of_ptex.first_loop - coarse_mloop};
+		CustomData_interp(loop_data,
+		                  &loop_interpolation->loop_data_storage,
+		                  first_indices,
+		                  weights, NULL,
+		                  2,
+		                  1);
+		CustomData_interp(loop_data,
+		                  &loop_interpolation->loop_data_storage,
+		                  last_indices,
+		                  weights, NULL,
+		                  2,
+		                  3);
+	}
+}
+
+static void loop_interpolation_end(LoopsForInterpolation *loop_interpolation)
+{
+	if (loop_interpolation->loop_data_storage_allocated) {
+		CustomData_free(&loop_interpolation->loop_data_storage, 4);
+	}
+}
+
+static void subdiv_copy_vertex_data(
+        const SubdivMeshContext *ctx,
+        MVert *subdiv_vertex,
+        const VerticesForInterpolation *vertex_interpolation,
+        const float u, const float v)
+{
+	const int subdiv_vertex_index = subdiv_vertex - ctx->subdiv_mesh->mvert;
+	const float weights[4] = {(1.0f - u) * (1.0f - v),
+	                          u * (1.0f - v),
+	                          u * v,
+	                          (1.0f - u) * v};
+	CustomData_interp(vertex_interpolation->vertex_data,
+	                  &ctx->subdiv_mesh->vdata,
+	                  vertex_interpolation->vertex_indices,
+	                  weights, NULL,
+	                  4,
+	                  subdiv_vertex_index);
+	/* TODO(sergey): Set ORIGINDEX. */
+}
+
+static void subdiv_evaluate_vertices(SubdivMeshContext *ctx,
+                                     const int poly_index)
+{
+	Subdiv *subdiv = ctx->subdiv;
+	const int resolution = ctx->settings->resolution;
+	const int resolution2 = resolution * resolution;
+	const float inv_resolution_1 = 1.0f / (float)(resolution - 1);
+	/* Base/coarse mesh information. */
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MPoly *coarse_polyoly = coarse_mesh->mpoly;
+	const MPoly *coarse_poly = &coarse_polyoly[poly_index];
+	const int num_poly_ptex_faces = mpoly_ptex_faces_count_get(coarse_poly);
+	/* Hi-poly subdivided mesh. */
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	MVert *subdiv_vertert = subdiv_mesh->mvert;
+	const int ptex_face_index = subdiv->face_ptex_offset[poly_index];
+	/* Actual evaluation. */
+	VerticesForInterpolation vertex_interpolation;
+	vertex_interpolation_init(ctx, &vertex_interpolation, coarse_poly);
+	MVert *subdiv_vert = &subdiv_vertert[ptex_face_index * resolution2];
+	for (int ptex_of_poly_index = 0;
+	     ptex_of_poly_index < num_poly_ptex_faces;
+	     ptex_of_poly_index++)
+	{
+		vertex_interpolation_from_ptex(ctx,
+		                               &vertex_interpolation,
+		                               coarse_poly,
+		                               ptex_of_poly_index);
+		const int current_ptex_face_index =
+		        ptex_face_index + ptex_of_poly_index;
+		BKE_subdiv_eval_limit_patch_resolution_point_and_short_normal(
+		        subdiv,
+		        current_ptex_face_index,
+		        resolution,
+		        subdiv_vert, offsetof(MVert, co), sizeof(MVert),
+		        subdiv_vert, offsetof(MVert, no), sizeof(MVert));
+		for (int y = 0; y < resolution; y++) {
+			const float v = y * inv_resolution_1;
+			for (int x = 0; x < resolution; x++, subdiv_vert++) {
+				const float u = x * inv_resolution_1;
+				subdiv_copy_vertex_data(ctx,
+				                        subdiv_vert,
+				                        &vertex_interpolation,
+				                        u, v);
+			}
+		}
+	}
+	vertex_interpolation_end(&vertex_interpolation);
+}
+
+static void subdiv_copy_edge_data(SubdivMeshContext *ctx,
+	                              MEdge *subdiv_edge,
+                                  const MEdge *coarse_edge)
+{
+	if (coarse_edge == NULL) {
+		subdiv_edge->crease = 0;
+		subdiv_edge->bweight = 0;
+		subdiv_edge->flag = 0;
+		return;
+	}
+	const int coarse_edge_index = coarse_edge - ctx->coarse_mesh->medge;
+	const int subdiv_edge_index = subdiv_edge - ctx->subdiv_mesh->medge;
+	CustomData_copy_data(&ctx->coarse_mesh->edata,
+	                     &ctx->subdiv_mesh->edata,
+	                     coarse_edge_index,
+	                     subdiv_edge_index,
+	                     1);
+	if (ctx->edge_origindex != NULL) {
+		ctx->edge_origindex[subdiv_edge_index] = coarse_edge_index;
+	}
+}
+
+static MEdge *subdiv_create_edges_row(SubdivMeshContext *ctx,
+                                      MEdge *subdiv_edge,
+                                      const MEdge *coarse_edge,
+                                      const int start_vertex_index,
+                                      const int resolution)
+{
+	int vertex_index = start_vertex_index;
+	for (int edge_index = 0;
+	     edge_index < resolution - 1;
+	     edge_index++, subdiv_edge++)
+	{
+		subdiv_copy_edge_data(ctx, subdiv_edge, coarse_edge);
+		subdiv_edge->v1 = vertex_index;
+		subdiv_edge->v2 = vertex_index + 1;
+		vertex_index += 1;
+	}
+	return subdiv_edge;
+}
+
+static MEdge *subdiv_create_edges_column(SubdivMeshContext *ctx,
+                                         MEdge *subdiv_edge,
+                                         const MEdge *coarse_start_edge,
+                                         const MEdge *coarse_end_edge,
+                                         const int start_vertex_index,
+                                         const int resolution)
+{
+	int vertex_index = start_vertex_index;
+	for (int edge_index = 0;
+	     edge_index < resolution;
+	     edge_index++, subdiv_edge++)
+	{
+		const MEdge *coarse_edge = NULL;
+		if (edge_index == 0) {
+			coarse_edge = coarse_start_edge;
+		}
+		else if (edge_index == resolution - 1) {
+			coarse_edge = coarse_end_edge;
+		}
+		subdiv_copy_edge_data(ctx, subdiv_edge, coarse_edge);
+		subdiv_edge->v1 = vertex_index;
+		subdiv_edge->v2 = vertex_index + resolution;
+		vertex_index += 1;
+	}
+	return subdiv_edge;
+}
+
+static void subdiv_create_edges(SubdivMeshContext *ctx, int poly_index)
+{
+	Subdiv *subdiv = ctx->subdiv;
+	const int resolution = ctx->settings->resolution;
+	const int resolution2 = resolution * resolution;
+	const int ptex_face_index = subdiv->face_ptex_offset[poly_index];
+	const int num_edges_per_ptex = num_edges_per_ptex_get(resolution);
+	const int start_edge_index = ptex_face_index * num_edges_per_ptex;
+	/* Base/coarse mesh information. */
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MPoly *coarse_polyoly = coarse_mesh->mpoly;
+	const MPoly *coarse_poly = &coarse_polyoly[poly_index];
+	const int num_poly_ptex_faces = mpoly_ptex_faces_count_get(coarse_poly);
+	/* Hi-poly subdivided mesh. */
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	MEdge *subdiv_medge = subdiv_mesh->medge;
+	MEdge *subdiv_edge = &subdiv_medge[start_edge_index];
+	const int start_poly_vertex_index = ptex_face_index * resolution2;
+	/* Consider a subdivision of base face at level 1:
+	 *
+	 *  y
+	 *  ^
+	 *  |   (6) ---- (7) ---- (8)
+	 *  |    |        |        |
+	 *  |   (3) ---- (4) ---- (5)
+	 *  |    |        |        |
+	 *  |   (0) ---- (1) ---- (2)
+	 *  o---------------------------> x
+	 *
+	 * This is illustrate which parts of geometry is created by code below.
+	 */
+	for (int i = 0; i < num_poly_ptex_faces; i++) {
+		const int start_ptex_face_vertex_index =
+		        start_poly_vertex_index + i * resolution2;
+		EdgesOfPtex edges_of_ptex;
+		edges_of_ptex_get(ctx, &edges_of_ptex, coarse_poly, i);
+		/* Create bottom row of edges (0-1, 1-2). */
+		subdiv_edge = subdiv_create_edges_row(ctx,
+		                                      subdiv_edge,
+		                                      edges_of_ptex.first_edge,
+	                                          start_ptex_face_vertex_index,
+		                                      resolution);
+		/* Create remaining edges. */
+		for (int row = 0; row < resolution - 1; row++) {
+			const int start_row_vertex_index =
+			    start_ptex_face_vertex_index + row * resolution;
+			/* Create vertical columns.
+			 *
+			 * At first iteration it will be edges (0-3. 1-4, 2-5), then it
+			 * will be (3-6, 4-7, 5-8) and so on.
+			 */
+			subdiv_edge = subdiv_create_edges_column(
+			        ctx,
+			        subdiv_edge,
+			        edges_of_ptex.last_edge,
+			        edges_of_ptex.second_edge,
+			        start_row_vertex_index,
+			        resolution);
+			/* Create horizontal edge row.
+			 *
+			 * At first iteration it will be edges (3-4, 4-5), then it will be
+			 * (6-7, 7-8) and so on.
+			 */
+			subdiv_edge = subdiv_create_edges_row(
+			        ctx,
+			        subdiv_edge,
+			        (row == resolution - 2) ? edges_of_ptex.third_edge
+			                                : NULL,
+			        start_row_vertex_index + resolution,
+			        resolution);
+		}
+	}
+}
+
+static void subdiv_copy_loop_data(
+        const SubdivMeshContext *ctx,
+        MLoop *subdiv_loop,
+        const LoopsForInterpolation *loop_interpolation,
+        const float u, const float v)
+{
+	const int subdiv_loop_index = subdiv_loop - ctx->subdiv_mesh->mloop;
+	const float weights[4] = {(1.0f - u) * (1.0f - v),
+	                          u * (1.0f - v),
+	                          u * v,
+	                          (1.0f - u) * v};
+	CustomData_interp(loop_interpolation->loop_data,
+	                  &ctx->subdiv_mesh->ldata,
+	                  loop_interpolation->loop_indices,
+	                  weights, NULL,
+	                  4,
+	                  subdiv_loop_index);
+	/* TODO(sergey): Set ORIGINDEX. */
+}
+
+static void subdiv_eval_uv_layer(SubdivMeshContext *ctx,
+                                 MLoop *subdiv_loop,
+                                 const int ptex_face_index,
+                                 const float u, const float v,
+                                 const float inv_resolution_1)
+{
+	if (ctx->num_uv_layers == 0) {
+		return;
+	}
+	Subdiv *subdiv = ctx->subdiv;
+	const int mloop_index = subdiv_loop - ctx->subdiv_mesh->mloop;
+	const float du = inv_resolution_1;
+	const float dv = inv_resolution_1;
+	for (int layer_index = 0; layer_index < ctx->num_uv_layers; layer_index++) {
+		MLoopUV *subdiv_loopuv = &ctx->uv_layers[layer_index][mloop_index];
+		BKE_subdiv_eval_face_varying(subdiv,
+		                             ptex_face_index,
+		                             u, v,
+		                             subdiv_loopuv[0].uv);
+		BKE_subdiv_eval_face_varying(subdiv,
+		                             ptex_face_index,
+		                             u + du, v,
+		                             subdiv_loopuv[1].uv);
+		BKE_subdiv_eval_face_varying(subdiv,
+		                             ptex_face_index,
+		                             u + du, v + dv,
+		                             subdiv_loopuv[2].uv);
+		BKE_subdiv_eval_face_varying(subdiv,
+		                             ptex_face_index,
+		                             u, v + dv,
+		                             subdiv_loopuv[3].uv);
+		/* TODO(sergey): Currently evaluator only has single UV layer, so can
+		 * not evaluate more than that. Need to be solved.
+		 */
+		break;
+	}
+}
+
+static void subdiv_create_loops(SubdivMeshContext *ctx, int poly_index)
+{
+	Subdiv *subdiv = ctx->subdiv;
+	const int resolution = ctx->settings->resolution;
+	const int resolution2 = resolution * resolution;
+	const float inv_resolution_1 = 1.0f / (float)(resolution - 1);
+	const int ptex_face_index = subdiv->face_ptex_offset[poly_index];
+	const int num_edges_per_ptex = num_edges_per_ptex_get(resolution);
+	const int start_edge_index = ptex_face_index * num_edges_per_ptex;
+	const int num_polys_per_ptex = num_polys_per_ptex_get(resolution);
+	const int start_poly_index = ptex_face_index * num_polys_per_ptex;
+	const int start_loop_index = 4 * start_poly_index;
+	const int start_vert_index = ptex_face_index * resolution2;
+	const float du = inv_resolution_1;
+	const float dv = inv_resolution_1;
+	/* Base/coarse mesh information. */
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MPoly *coarse_polyoly = coarse_mesh->mpoly;
+	const MPoly *coarse_poly = &coarse_polyoly[poly_index];
+	const int num_poly_ptex_faces = mpoly_ptex_faces_count_get(coarse_poly);
+	/* Hi-poly subdivided mesh. */
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	MLoop *subdiv_loopoop = subdiv_mesh->mloop;
+	MLoop *subdiv_loop = &subdiv_loopoop[start_loop_index];
+	LoopsForInterpolation loop_interpolation;
+	loop_interpolation_init(ctx, &loop_interpolation, coarse_poly);
+	for (int ptex_of_poly_index = 0;
+	     ptex_of_poly_index < num_poly_ptex_faces;
+	     ptex_of_poly_index++)
+	{
+		loop_interpolation_from_ptex(ctx,
+		                             &loop_interpolation,
+		                             coarse_poly,
+		                             ptex_of_poly_index);
+		const int current_ptex_face_index =
+		        ptex_face_index + ptex_of_poly_index;
+		for (int y = 0; y < resolution - 1; y++) {
+			const float v = y * inv_resolution_1;
+			for (int x = 0; x < resolution - 1; x++, subdiv_loop += 4) {
+				const float u = x * inv_resolution_1;
+				/* Vertex indicies ordered counter-clockwise. */
+				const int v0 = start_vert_index +
+				               (ptex_of_poly_index * resolution2) +
+				               (y * resolution + x);
+				const int v1 = v0 + 1;
+				const int v2 = v0 + resolution + 1;
+				const int v3 = v0 + resolution;
+				/* Edge indicies ordered counter-clockwise. */
+				const int e0 = start_edge_index +
+				               (ptex_of_poly_index * num_edges_per_ptex) +
+				               (y * (2 * resolution - 1) + x);
+				const int e1 = e0 + resolution;
+				const int e2 = e0 + (2 * resolution - 1);
+				const int e3 = e0 + resolution - 1;
+				/* Initialize 4 loops of corresponding hi-poly poly. */
+				/* TODO(sergey): For ptex boundaries we should use loops from
+				 * coarse mesh.
+				 */
+				subdiv_copy_loop_data(ctx,
+				                      &subdiv_loop[0],
+				                      &loop_interpolation,
+				                      u, v);
+				subdiv_loop[0].v = v0;
+				subdiv_loop[0].e = e0;
+				subdiv_copy_loop_data(ctx,
+				                      &subdiv_loop[1],
+				                      &loop_interpolation,
+				                      u + du, v);
+				subdiv_loop[1].v = v1;
+				subdiv_loop[1].e = e1;
+				subdiv_copy_loop_data(ctx,
+				                      &subdiv_loop[2],
+				                      &loop_interpolation,
+				                      u + du, v + dv);
+				subdiv_loop[2].v = v2;
+				subdiv_loop[2].e = e2;
+				subdiv_copy_loop_data(ctx,
+				                      &subdiv_loop[3],
+				                      &loop_interpolation,
+				                      u, v + dv);
+				subdiv_loop[3].v = v3;
+				subdiv_loop[3].e = e3;
+				/* Interpolate UV layers using OpenSubdiv. */
+				subdiv_eval_uv_layer(ctx,
+				                     subdiv_loop,
+				                     current_ptex_face_index,
+				                     u, v,
+				                     inv_resolution_1);
+			}
+		}
+	}
+	loop_interpolation_end(&loop_interpolation);
+}
+
+static void subdiv_copy_poly_data(const SubdivMeshContext *ctx,
+                                  MPoly *subdiv_poly,
+                                  const MPoly *coarse_poly)
+{
+	const int coarse_poly_index = coarse_poly - ctx->coarse_mesh->mpoly;
+	const int subdiv_poly_index = subdiv_poly - ctx->subdiv_mesh->mpoly;
+	CustomData_copy_data(&ctx->coarse_mesh->pdata,
+	                     &ctx->subdiv_mesh->pdata,
+	                     coarse_poly_index,
+	                     subdiv_poly_index,
+	                     1);
+	if (ctx->poly_origindex != NULL) {
+		ctx->poly_origindex[subdiv_poly_index] = coarse_poly_index;
+	}
+}
+
+static void subdiv_create_polys(SubdivMeshContext *ctx, int poly_index)
+{
+	Subdiv *subdiv = ctx->subdiv;
+	const int resolution = ctx->settings->resolution;
+	const int ptex_face_index = subdiv->face_ptex_offset[poly_index];
+	const int num_polys_per_ptex = num_polys_per_ptex_get(resolution);
+	const int num_loops_per_ptex = 4 * num_polys_per_ptex;
+	const int start_poly_index = ptex_face_index * num_polys_per_ptex;
+	const int start_loop_index = 4 * start_poly_index;
+	/* Base/coarse mesh information. */
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MPoly *coarse_polyoly = coarse_mesh->mpoly;
+	const MPoly *coarse_poly = &coarse_polyoly[poly_index];
+	const int num_poly_ptex_faces = mpoly_ptex_faces_count_get(coarse_poly);
+	/* Hi-poly subdivided mesh. */
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	MPoly *subdiv_mpoly = subdiv_mesh->mpoly;
+	MPoly *subdiv_mp = &subdiv_mpoly[start_poly_index];
+	for (int ptex_of_poly_index = 0;
+	     ptex_of_poly_index < num_poly_ptex_faces;
+	     ptex_of_poly_index++)
+	{
+		for (int subdiv_poly_index = 0;
+		     subdiv_poly_index < num_polys_per_ptex;
+		     subdiv_poly_index++, subdiv_mp++)
+		{
+			subdiv_copy_poly_data(ctx, subdiv_mp, coarse_poly);
+			subdiv_mp->loopstart = start_loop_index +
+			                       (ptex_of_poly_index * num_loops_per_ptex) +
+			                       (subdiv_poly_index * 4);
+			subdiv_mp->totloop = 4;
+		}
+	}
+}
+
+static void subdiv_eval_task(
+        void *__restrict userdata,
+        const int poly_index,
+        const ParallelRangeTLS *__restrict UNUSED(tls))
+{
+	SubdivMeshContext *data = userdata;
+	/* Evaluate hi-poly vertex coordinates and normals. */
+	subdiv_evaluate_vertices(data, poly_index);
+	/* Create mesh geometry for the given base poly index. */
+	subdiv_create_edges(data, poly_index);
+	subdiv_create_loops(data, poly_index);
+	subdiv_create_polys(data, poly_index);
+}
+
+static void cache_uv_layers(SubdivMeshContext *ctx)
+{
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	ctx->num_uv_layers =
+	        CustomData_number_of_layers(&subdiv_mesh->ldata, CD_MLOOPUV);
+	for (int layer_index = 0; layer_index < ctx->num_uv_layers; ++layer_index) {
+		ctx->uv_layers[layer_index] = CustomData_get_layer_n(
+		        &subdiv_mesh->ldata, CD_MLOOPUV, layer_index);
+	}
+}
+
+static void cache_custom_data_layers(SubdivMeshContext *ctx)
+{
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	/* Pointers to original indices layers. */
+	ctx->vert_origindex = CustomData_get_layer(
+	        &subdiv_mesh->vdata, CD_ORIGINDEX);
+	ctx->edge_origindex = CustomData_get_layer(
+	        &subdiv_mesh->edata, CD_ORIGINDEX);
+	ctx->loop_origindex = CustomData_get_layer(
+	        &subdiv_mesh->ldata, CD_ORIGINDEX);
+	ctx->poly_origindex = CustomData_get_layer(
+	        &subdiv_mesh->pdata, CD_ORIGINDEX);
+	/* UV layers interpolation. */
+	cache_uv_layers(ctx);
+}
+
+Mesh *BKE_subdiv_to_mesh(
+        Subdiv *subdiv,
+        const SubdivToMeshSettings *settings,
+        const Mesh *coarse_mesh)
+{
+	/* Make sure evaluator is up to date with possible new topology, and that
+	 * is is refined for the new positions of coarse vertices.
+	 */
+	BKE_subdiv_eval_update_from_mesh(subdiv, coarse_mesh);
+	const int resolution = settings->resolution;
+	const int resolution2 = resolution * resolution;
+	const int num_result_verts = subdiv->num_ptex_faces * resolution2;
+	const int num_result_edges =
+			subdiv->num_ptex_faces * num_edges_per_ptex_get(resolution);
+	const int num_result_polys =
+			subdiv->num_ptex_faces * num_polys_per_ptex_get(resolution);
+	const int num_result_loops = 4 * num_result_polys;
+	/* Create mesh and its arrays. */
+	Mesh *result = BKE_mesh_new_nomain_from_template(
+	        coarse_mesh,
+	        num_result_verts,
+	        num_result_edges,
+	        0,
+	        num_result_loops,
+	        num_result_polys);
+	/* Evaluate subdivisions of base faces in threads. */
+	SubdivMeshContext ctx;
+	ctx.coarse_mesh = coarse_mesh;
+	ctx.subdiv = subdiv;
+	ctx.subdiv_mesh = result;
+	ctx.settings = settings;
+	cache_custom_data_layers(&ctx);
+	/* Multi-threaded evaluation. */
+	ParallelRangeSettings parallel_range_settings;
+	BLI_parallel_range_settings_defaults(&parallel_range_settings);
+	BLI_task_parallel_range(0, coarse_mesh->totpoly,
+	                        &ctx,
+	                        subdiv_eval_task,
+	                        &parallel_range_settings);
+	return result;
+}