Some comments.

This reverts commit f0672fb586.

CMakeLists.txt

@@ -273,6 +273,7 @@ option(WITH_MOD_REMESH          "Enable Remesh Modifier" ON)
 # option(WITH_MOD_CLOTH_ELTOPO    "Enable Experimental cloth solver" OFF)  # this is now only available in a branch
 # mark_as_advanced(WITH_MOD_CLOTH_ELTOPO)
 option(WITH_MOD_OCEANSIM        "Enable Ocean Modifier" OFF)
+option(WITH_MOD_QUADREMESH      "Enable Quad Remesh Modifier" ON)
 
 # Image format support
 option(WITH_OPENIMAGEIO         "Enable OpenImageIO Support (http://www.openimageio.org)" ON)
@@ -2847,6 +2848,7 @@ if(FIRST_RUN)
 	info_cfg_option(WITH_MOD_REMESH)
 	info_cfg_option(WITH_MOD_FLUID)
 	info_cfg_option(WITH_MOD_OCEANSIM)
+	info_cfg_option(WITH_MOD_QUADREMESH)
 
 	info_cfg_text("OpenGL:")
 	info_cfg_option(WITH_GLEW_ES)

build_files/cmake/macros.cmake

@@ -655,6 +655,10 @@ macro(SETUP_BLENDER_SORTED_LIBS)
 	if(WITH_BULLET AND NOT WITH_SYSTEM_BULLET)
 		list_insert_after(BLENDER_SORTED_LIBS "ge_logic_ngnetwork" "extern_bullet")
 	endif()
+	
+	if(WITH_MOD_QUADREMESH)
+		list(APPEND BLENDER_SORTED_LIBS bf_intern_quadremesh)
+	endif()
 
 	foreach(SORTLIB ${BLENDER_SORTED_LIBS})
 		set(REMLIB ${SORTLIB})

intern/CMakeLists.txt

@@ -77,6 +77,10 @@ if(WITH_OPENNL)
 	add_subdirectory(opennl)
 endif()
 
+if(WITH_MOD_QUADREMESH)
+	add_subdirectory(quadremesh)
+endif()
+
 # only windows needs utf16 converter
 if(WIN32)
 	add_subdirectory(utfconv)

intern/quadremesh/CMakeLists.txt

@@ -0,0 +1,50 @@
+# ***** 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.
+#
+# ***** END GPL LICENSE BLOCK *****
+
+set(INC
+	.
+	../../source/blender/blenkernel
+	../../source/blender/blenlib
+	../../source/blender/makesdna
+	../../source/blender/modifiers/intern
+	../../source/blender/modifiers
+	../../source/blender/depsgraph
+	../opennl/extern
+	../guardedalloc
+)
+
+set(INC_SYS
+	${EIGEN3_INCLUDE_DIRS}
+)
+
+set(SRC
+	QRM_quadremesh.c
+	quadremesh_util.c
+	quadremesh_line.c
+	quadremesh_link.c
+	quadremesh_edge.c
+	quadremesh_input.c
+	quadremesh_harmonic.c
+	quadremesh_curvature.cpp
+	
+	QRM_quadremesh.h
+	quadremesh_util.h
+)
+
+blender_add_lib(bf_intern_quadremesh "${SRC}" "${INC}" "${INC_SYS}")
+

intern/quadremesh/QRM_quadremesh.c

@@ -0,0 +1,613 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Alexander Pinzon Fernandez
+ *                 Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+/** \file intern/quadremesh/MOD_quadremesh.c
+ *  \ingroup quadremesh
+ */
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_modifier_types.h"
+#include "DNA_object_types.h"
+
+#include "BLI_math.h"
+#include "BLI_rand.h"
+#include "BLI_linklist.h"
+#include "BLI_memarena.h"
+
+#include "BKE_cdderivedmesh.h"
+#include "BKE_mesh_mapping.h"
+
+#include "PIL_time.h"
+
+#include "MOD_util.h"
+#include "quadremesh_util.h"
+
+#define QR_MAKEPOLYS
+//#define QR_SHOWORIENTATIONS
+
+#if 0 /* UNUSED ROUTINES */
+/**
+* Random point, P, uniformly from within triangle ABC, method given by
+* Robert Osada, Thomas Funkhouser, Bernard Chazelle, and David Dobkin. 2002. Shape distributions. ACM Trans. Graph. 21,
+* 4 (October 2002), 807-832. DOI=10.1145/571647.571648 http://doi.acm.org/10.1145/571647.571648
+* a,b,c are the triangle points
+* r1, r2, are the randon numbers betwen [0, 1]
+* P = (1 − sqrt(r1)) A + sqrt(r1)(1 − r2) B + sqrt(r1) * r2 C
+*/
+static void uniformRandomPointWithinTriangle(float r[3], float a[3], float b[3], float c[3])
+{
+	float va, vb, vc;
+	float pa[3], pb[3], pc[3];
+	float r1, r2;
+	r1 = BLI_frand();
+	r2 = BLI_frand();
+	va = 1.0f - sqrtf(r1);
+	vb = sqrtf(r1) * (1.0f - r2);
+	vc = sqrtf(r1) * r2;
+	mul_v3_v3fl(pa, a, va);
+	mul_v3_v3fl(pb, b, vb);
+	mul_v3_v3fl(pc, c, vc);
+}
+
+static void uniformRandomPointWithinFace(float r[3], LaplacianSystem *sys, int indexf){
+	int *vin;
+	vin = sys->faces[indexf];
+	uniformRandomPointWithinTriangle(r, sys->co[vin[0]], sys->co[vin[1]], sys->co[vin[2]]);
+}
+
+/*
+* Compute the normal curvature
+* k = dot(2*no, (pi - pj)) / (|pi - pj|)^2
+* no = normal on vertex pi
+* pi - pj is a vector direction on this case the gradient field direction
+* the gradient field direction on some vertex is computed how the average of the faces around vertex
+*/
+static void computeSampleDistanceFunctions(LaplacianSystem *sys, float user_h, float user_alpha) {
+	int i, j, *fin, lin;
+	float avg1[3], avg2[3], no[3], k1, k2, h1, h2;
+
+	sys->h1 = MEM_mallocN(sizeof(float) * im->num_verts, "QuadRemeshH1");
+	sys->h2 = MEM_mallocN(sizeof(float) * im->num_verts, "QuadRemeshH2");
+
+	for (i = 0; i < sys->total_verts; i++) {
+		zero_v3(avg1);
+		zero_v3(avg2);
+
+		fin = sys->ringf_map[i].indices;
+		lin = sys->ringf_map[i].count;
+
+		for (j = 0; j < lin; j++) {
+			add_v3_v3(avg1, sys->gf1[j]);
+			add_v3_v3(avg2, sys->gf2[j]);
+		}
+
+		mul_v3_fl(avg1, 1.0f / ((float)lin));
+		mul_v3_fl(avg2, 1.0f / ((float)lin));
+
+		copy_v3_v3(no, sys->no[i]);
+		mul_v3_fl(no, 2.0f);
+
+		k1 = dot_v3v3(no, avg1) / dot_v3v3(avg1, avg1);
+		k2 = dot_v3v3(no, avg2) / dot_v3v3(avg2, avg2);
+
+		h1 = user_h / (1.0f + user_alpha * (logf(1.0f + k1*k1)));
+		h2 = user_h / (1.0f + user_alpha * (logf(1.0f + k2*k2)));
+
+		sys->h1[i] = h1;
+		sys->h2[i] = h2;
+	}
+}
+bool nextPointFromVertex(float r_co[3], int *r_face, int *r_edge, GradientFlowSystem *gfsys, int in_v)
+{
+	int i, f, vf1, vf2;
+	float gf[3], a[3], b[3], alfa, beta, gamma, dummy[3];
+	QuadRemeshSystem *sys = gfsys->sys;
+
+	for (i = 0; i < sys->ringf_map[in_v].count; i++) {
+		f = sys->ringf_map[in_v].indices[i];
+		if (!getSecondAndThirdVert(&vf1, &vf2, sys, f, in_v)) continue;
+
+		sub_v3_v3v3(a, sys->co[vf1], sys->co[in_v]);
+		sub_v3_v3v3(b, sys->co[vf2], sys->co[in_v]);
+		copy_v3_v3(gf, gfsys->gfield[f]);
+
+		alfa = angle_v3v3(a, b);
+		beta = angle_v3v3(a, gf);
+		gamma = angle_v3v3(gf, b);
+
+		if (beta + gamma >= alfa - 0.001f && beta + gamma <= alfa + 0.001f) {
+			/* vertex on this face */
+			add_v3_v3(gf, sys->co[in_v]);
+			isect_line_line_v3(sys->co[vf1], sys->co[vf2], sys->co[in_v], gf, r_co, dummy);
+			*r_edge = getEdgeFromVerts(sys, vf1, vf2);
+			*r_face = f;
+
+			return true;
+		}
+	}
+
+	return false;
+}
+
+/*
+* alpha is degree of anisotropic curvature sensitivity
+* h is the desired distance
+* return ve[0] number of vertices
+* return ve[1] number of edges
+*/
+void estimateNumberGFVerticesEdges(int ve[2], QuadRemeshSystem *sys, float h)
+{
+	int i, totalv, totale;
+	float area = 0.0f;
+	float sqrtarea;
+	for (i = 0; i < sys->total_faces; i++) {
+		area += area_tri_v3(sys->co[sys->faces[i][0]], sys->co[sys->faces[i][1]], sys->co[sys->faces[i][2]]);
+	}
+	sqrtarea = sqrtf(area);
+	if (h > 0.0f) {
+		totalv = ((sqrtarea / h) + 1.0f);
+		totale = totalv * sqrtarea * 2.0f;
+		totalv = totalv * totalv;
+	}
+	else{
+		totalv = sqrtarea + 1.0f;
+		totale = totalv * sqrtarea * 2.0f;
+		totalv = totalv * totalv;
+	}
+	ve[0] = totalv;
+	ve[1] = totale;
+}
+
+static bool isPointOnSegment(float in_v[3], float in_l1[3], float in_l2[3])
+{
+	float a[3], b[3];
+	sub_v3_v3v3(a, in_v, in_l1);
+	sub_v3_v3v3(b, in_l2, in_v);
+	return dot_v3v3(a, b) >= 0.0f;
+}
+
+
+static int getFaceFromTwoEdges(InputMesh *im, int in_e1, int in_e2)
+{
+	if (im->faces_edge[in_e1][0] == im->faces_edge[in_e2][0])
+		return im->faces_edge[in_e2][0];
+
+	return im->faces_edge[in_e2][1];
+}
+
+static int getEdgeOppositeToVertex(InputMesh *im, int in_v, int in_f)
+{
+	int i;
+
+	for (i = 0; i < 3; i++)
+		if (im->faces[in_f][i] == in_v)
+			break;
+
+	return getEdgeFromVerts(im, im->faces[in_f][(i + 1) % 3], im->faces[in_f][(i + 2) % 3]);
+}
+
+static float getSamplingDistanceFunctionOnFace(GradientFlowSystem *gfsys, int in_f, float in_co[3])
+{
+	float h1, h2, h3, uv[2];
+	InputMesh *im = &gfsys->sys->input_mesh;
+
+	resolve_tri_uv_v3(uv, in_co, im->co[im->faces[in_f][0]], im->co[im->faces[in_f][1]], im->co[im->faces[in_f][2]]);
+
+	h1 = gfsys->hfunction[im->faces[in_f][0]];
+	h2 = gfsys->hfunction[im->faces[in_f][1]];
+	h3 = gfsys->hfunction[im->faces[in_f][2]];
+
+	return uv[0] * h1 + uv[1] * h2 + (1.0f - uv[0] - uv[1]) * h3;
+}
+#endif // UNUSED ROUTINES
+
+/* MESH GENERATION */
+
+static void makeFeatureEdges(OutputMesh *om, InputMesh *im)
+{
+	int i, j, f1, f2;
+	float angle;
+	GFPoint p1, p2;
+
+	for (i = 0; i < im->num_edges; i++) {
+		f1 = im->faces_edge[i][0];
+		f2 = im->faces_edge[i][1];
+		if (f1 == QR_NO_FACE || f2 == QR_NO_FACE)
+			angle = M_PI * 2.0f;
+		else
+			angle = angle_normalized_v3v3(im->no[f1], im->no[f2]);
+
+		if (angle >= M_PI_2 * 0.7f) {
+			p1.type = eVert;
+			p1.v = im->edges[i][0];
+			copy_v3_v3(p1.co, im->co[im->edges[i][0]]);
+			addGFPoint(im, om, &p1);
+
+			for (j = 0; j < im->ringe_map[p1.v].count; j++)
+				insertOnQREdge(om, &om->ringe[im->ringe_map[p1.v].indices[j]], p1.id);
+
+			p2.type = eVert;
+			p2.v = im->edges[i][1];
+			copy_v3_v3(p2.co, im->co[im->edges[i][1]]);
+			addGFPoint(im, om, &p2);
+
+			for (j = 0; j < im->ringe_map[p2.v].count; j++)
+				insertOnQREdge(om, &om->ringe[im->ringe_map[p2.v].indices[j]], p2.id);
+
+			linkOnQREdge(om, GFSYSNONE, &om->ringe[i], p1.id, p2.id);
+		}
+	}
+}
+
+static void deleteDegenerateVerts(OutputMesh *om)
+{
+	int i, n;
+	bool *donevert, doverts = true;
+	QRVertID *vertmap;
+	QRVert *newverts;
+	QRDiskLink *it;
+
+	donevert = MEM_callocN(sizeof(QRVertID) * om->num_verts, __func__);
+
+	while (doverts) {
+		doverts = false;
+		for (i = 0; i < om->num_verts; i++) {
+			if (donevert[i])
+				continue;
+
+			donevert[i] = true;
+
+			if (om->verts[i].num_links == 1) {
+				doverts = true;
+				donevert[om->verts[i].link->v] = false;
+				unlinkVerts(om, om->verts[i].link);
+			}
+			else if (om->verts[i].num_links == 2) {
+				doverts = true;
+				donevert[om->verts[i].link->v] = false;
+				donevert[om->verts[i].link->next->v] = false;
+				dissolveVert(om, i);
+			}
+		}
+	}
+
+	MEM_SAFE_FREE(donevert);
+	vertmap = MEM_mallocN(sizeof(QRVertID) * om->num_verts, __func__);
+
+	for (i = 0, n = 0; i < om->num_verts; i++) {
+		if (om->verts[i].num_links == 0)
+			vertmap[i] = -1;
+		else
+			vertmap[i] = n++;
+	}
+
+	newverts = MEM_callocN(sizeof(QRVert) * n, "Newlinks");
+
+	for (i = 0, n = 0; i < om->num_verts; i++) {
+		if (vertmap[i] == -1) continue;
+
+		memcpy(&newverts[n], &om->verts[i], sizeof(QRVert));
+
+		it = newverts[n].link;
+		do {
+			it->v = vertmap[it->v];
+			it = it->next;
+		} while (it != newverts[n].link);
+		n++;
+	}
+
+	MEM_SAFE_FREE(om->verts);
+	MEM_SAFE_FREE(vertmap);
+	om->num_verts = n;
+	om->verts = newverts;
+}
+
+static void hideEdgesOnFaces(OutputMesh *om, InputMesh *im)
+{
+	int f;
+	QREdge *e;
+	QREdgeLink *itel;
+	LinkNode *it;
+
+	for (f = 0; f < im->num_faces; f++) {
+		for (it = om->ringf[f]; it; it = it->next) {
+			e = (QREdge*)it->link;
+			for (itel = e->v1; itel; itel = itel->next)
+				if (itel->elink) {
+					itel->elink->e = 0;
+					itel->elink->brother->e = 0;
+				}
+		}
+	}
+}
+
+static void makeEdges(OutputMesh *om, MEdge *r_edges)
+{
+	int i, j, e;
+	QRDiskLink *it;
+
+	for (i = 0, e = 0; i < om->num_verts; i++) {
+		for (j = 0, it = om->verts[i].link;
+			 j < om->verts[i].num_links;
+			 j++, it = it->next)
+		{
+			if (it->e < 0) {
+				it->brother->e = it->e = e;
+
+				r_edges[e].v1 = i;
+				r_edges[e].v2 = it->v;
+				r_edges[e].bweight = 0;
+				r_edges[e].crease = 0;
+				r_edges[e].flag = ME_EDGEDRAW;
+				e++;
+			}
+		}
+	}
+}
+
+static void makePolys(OutputMesh *om, MPoly *r_polys, MLoop *r_loops)
+{
+	int i, j, s, num_loops, num_polys, max_num_polys;
+	QRVertID v;
+	QRDiskLink *it, *lit;
+
+	num_loops = 0;
+	num_polys = 0;
+	max_num_polys = 2 + om->num_edges - om->num_verts;
+
+	for (i = 0; i < om->num_verts; i++) {
+		for (j = 0, it = om->verts[i].link;
+			 j < om->verts[i].num_links;
+			 j++, it = it->next)
+		{
+			if (it->poly_on_right) continue;
+
+			s = num_loops;
+			v = i;
+			lit = it;
+			do {
+				r_loops[num_loops].v = v;
+				r_loops[num_loops].e = lit->e;
+				num_loops++;
+
+				v = lit->v;
+				lit->poly_on_right = true;
+				lit = lit->brother->next;
+			} while (v != i);
+			
+			r_polys[num_polys].loopstart = s;
+			r_polys[num_polys].totloop = num_loops - s;
+			r_polys[num_polys].mat_nr = 0;
+			r_polys[num_polys].flag = 0;
+			if (++num_polys == max_num_polys) {
+				printf("Num polys: %d / %d  ", num_polys, max_num_polys);
+				return;
+			}
+		}
+	}
+
+	printf("Num polys: %d / %d  ", num_polys, max_num_polys);
+}
+
+#if 0
+static void makeNormals(OutputMesh *om)
+{
+	int i, tv = om->vert;
+	float npos[3];
+
+	for (i = 0; i < tv; i++) {
+		normal_short_to_float_v3(npos, om->verts[i].no);
+		mul_v3_fl(npos, 0.1f);
+		add_v3_v3(npos, om->verts[i].co);
+		addEdge(om, i, addVert(om, npos, NULL));
+	}
+}
+#endif // 0
+
+static void initOutputMesh(OutputMesh *om, InputMesh *im)
+{
+	om->num_verts = om->alloc_verts = om->num_edges = 0;
+
+	om->memarena = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "Output Mesh");
+	om->vonvs = MEM_mallocN(sizeof(QRVertID) * im->num_verts, __func__);
+	copy_vn_i(om->vonvs, im->num_verts, -1);
+	om->ringe = MEM_callocN(sizeof(QREdge) * im->num_edges, "GFListEdges");
+	om->ringf = MEM_callocN(sizeof(LinkNode *) * im->num_faces, "GFListFaces");
+}
+
+static void freeOutputMesh(OutputMesh *om)
+{
+	om->num_edges = om->num_verts = om->alloc_verts = 0;
+	MEM_SAFE_FREE(om->verts);
+	MEM_SAFE_FREE(om->vonvs);
+	MEM_SAFE_FREE(om->ringe);
+	MEM_SAFE_FREE(om->ringf);
+	if (om->memarena) {
+		BLI_memarena_free(om->memarena);
+		om->memarena = NULL;
+	}
+}
+
+static DerivedMesh *makeOrientationsMesh(QuadRemeshSystem *sys)
+{
+	int i;
+	float mid[3], gra1[3], gra2[3];
+	DerivedMesh *ret;
+	InputMesh *im = &sys->input_mesh;
+	MVert *verts;
+	MEdge *edges;
+	
+	ret = CDDM_new(im->num_faces * 3, im->num_faces * 2, 0, 0, 0);
+
+	verts = ret->getVertArray(ret);
+	edges = ret->getEdgeArray(ret);
+
+	for (i = 0; i < im->num_faces; i++) {
+		mid_v3_v3v3v3(mid, im->co[im->faces[i][0]], im->co[im->faces[i][1]], im->co[im->faces[i][2]]);
+		madd_v3_v3v3fl(gra1, mid, sys->gfsys[0]->gf[i], 0.04f);
+		madd_v3_v3v3fl(gra2, mid, sys->gfsys[1]->gf[i], 0.04f);
+
+		copy_v3_v3(verts[i * 3].co, mid);
+		copy_v3_v3(verts[i * 3 + 1].co, gra1);
+		copy_v3_v3(verts[i * 3 + 2].co, gra2);
+
+		edges[i * 2].v1 = i * 3;
+		edges[i * 2].v2 = i * 3 + 1;
+		edges[i * 2].flag = ME_EDGEDRAW;
+		edges[i * 2 + 1].v1 = i * 3;
+		edges[i * 2 + 1].v2 = i * 3 + 2;
+		edges[i * 2 + 1].flag = ME_EDGEDRAW;
+	}
+
+	return ret;
+}
+
+static DerivedMesh *remesh(QuadRemeshSystem *sys)
+{
+	int i, num_loops, num_polys;
+	double start_time;
+	MVert *verts;
+	OutputMesh *om = &sys->output_mesh;
+	DerivedMesh *ret;
+
+	start_time = PIL_check_seconds_timer();
+
+	sys->rng = BLI_rng_new(sys->qmd->rng_seed);
+
+	initOutputMesh(om, &sys->input_mesh);
+	sys->gfsys[0]->id = GFSYS1;
+	sys->gfsys[1]->id = GFSYS2;
+
+	makeFeatureEdges(om, &sys->input_mesh);
+	computeFlowLines(sys);
+	//hideEdgesOnFaces(om, &sys->input_mesh);
+	//makeNormals(om);
+
+	if (!om->num_verts) {
+		freeOutputMesh(om);
+		BLI_rng_free(sys->rng);
+
+		return NULL;
+	}
+
+#ifdef QR_MAKEPOLYS
+	deleteDegenerateVerts(om);
+
+	num_loops = om->num_edges * 2;
+	num_polys = 2 + om->num_edges - om->num_verts;
+#else
+	num_loops = num_polys = 0;
+#endif
+
+	ret = CDDM_new(om->num_verts, om->num_edges, 0, num_loops, num_polys);
+
+	verts = ret->getVertArray(ret);
+
+	for (i = 0; i < om->num_verts; i++) {
+		copy_v3_v3(verts[i].co, om->verts[i].co);
+		normal_float_to_short_v3(verts[i].no, om->verts[i].no);
+		verts[i].bweight = 0;
+		verts[i].flag = 0;
+	}
+
+	makeEdges(om, ret->getEdgeArray(ret));
+
+#ifdef QR_MAKEPOLYS
+	makePolys(om, ret->getPolyArray(ret), ret->getLoopArray(ret));
+
+	CDDM_recalc_tessellation(ret);
+	//CDDM_calc_edges_tessface(ret);
+#endif
+
+	freeOutputMesh(om);
+	BLI_rng_free(sys->rng);
+
+	printf("Mesh generation time: %f ms\n", (PIL_check_seconds_timer() - start_time) * 1000.0f);
+
+	return ret;
+}
+
+DerivedMesh *makeResultMesh(QuadRemeshSystem *sys, Object *ob, DerivedMesh *in)
+{
+	/* Get input mesh */
+	if (!sys->input_mesh.is_alloc) {
+		sys->qmd->flag |= MOD_QUADREMESH_INPUT_DIRTY;
+	}
+
+	if (sys->qmd->flag & MOD_QUADREMESH_INPUT_DIRTY) {
+		if (sys->input_mesh.is_alloc) {
+			freeInputMesh(&sys->input_mesh);
+		}
+
+		getInput(sys, ob, in);
+		sys->qmd->flag |= MOD_QUADREMESH_FIELD_DIRTY;
+	}
+
+	/* Compute harmonic field */
+	if (!sys->is_alloc) {
+		sys->qmd->flag |= MOD_QUADREMESH_FIELD_DIRTY;
+	}
+
+	if (sys->qmd->flag & MOD_QUADREMESH_FIELD_DIRTY) {
+		if (sys->is_alloc) {
+			MEM_SAFE_FREE(sys->U_field);
+			freeGradientFlowSystem(sys->gfsys[0]);
+			freeGradientFlowSystem(sys->gfsys[1]);
+
+			sys->is_alloc = false;
+		}
+
+		getHarmonicGradients(sys);
+		//getPrincipalCurvatures(sys);
+
+		if (!sys->has_solution) {
+			modifier_setError((ModifierData*)sys->qmd, "No soultion found given these constraints");
+			return in;
+		}
+
+		sys->qmd->flag |= MOD_QUADREMESH_REMESH;
+	}
+
+	/* Remesh */
+	if (sys->qmd->flag & MOD_QUADREMESH_REMESH) {
+		if (sys->cache_mesh) {
+			sys->cache_mesh->release(sys->cache_mesh);
+		}
+
+		sys->cache_mesh = remesh(sys);
+	}
+
+	sys->qmd->flag = 0;
+
+#ifndef QR_SHOWORIENTATIONS
+	if (sys->cache_mesh) {
+		return CDDM_copy(sys->cache_mesh);
+	}
+#else
+	return makeOrientationsMesh(sys);
+#endif
+
+	return NULL;
+}

intern/quadremesh/QRM_quadremesh.h

@@ -0,0 +1,52 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Alexander Pinzon Fernandez
+ *                 Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+/**
+ * Known issues:
+ *  - caching the u-field
+ *  - semi automatic extrema placement
+ *  - mesh anisotropy
+ */
+
+/** \file intern/quadremesh/MOD_quadremesh.h
+ *  \ingroup quadremesh
+ */
+
+#ifndef __MOD_QUADREMESH_H__
+#define __MOD_QUADREMESH_H__
+
+/* Default values for the modifier */
+#define QR_MINDIST 0.04f
+#define QR_SEED 1
+#define QR_AUTOUPDATES 1
+
+typedef struct QuadRemeshSystem QuadRemeshSystem;
+
+void initQuadRemeshSystem(struct QuadRemeshModifierData *qmd);
+struct DerivedMesh *makeResultMesh(QuadRemeshSystem *sys, struct Object *ob, struct DerivedMesh *dm);
+void freeQuadRemeshSystem(QuadRemeshSystem *sys);
+
+void getUField(QuadRemeshSystem *sys, struct MDeformVert *dvert, int defgrp_index);
+
+#endif /*__MOD_QUADREMESH_H__*/

intern/quadremesh/quadremesh_curvature.cpp

@@ -0,0 +1,191 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+extern "C" {
+
+#include "MEM_guardedalloc.h"
+
+#include "BLI_math.h"
+#include "BLI_utildefines.h"
+
+#include "BKE_mesh_mapping.h"
+
+#include "quadremesh_util.h"
+
+}
+
+#include <iostream>
+
+#include "Eigen/Core"
+#include "Eigen/Eigenvalues"
+
+using namespace Eigen;
+
+static int getVertexOppositeToEdge(InputMesh *im, int in_e, int in_f)
+{
+	int i;
+	for (i = 0; i < 3; i++) {
+		if (im->faces[in_f][i] != im->edges[in_e][0] &&
+			im->faces[in_f][i] != im->edges[in_e][1])
+		{
+			return im->faces[in_f][i];
+		}
+	}
+	return 0;
+}
+
+static void alignAtEdge(QuadRemeshSystem *sys, int in_e)
+{
+	int f1, f2;
+	float ev[3], vec[3], tmp[3], best[3], best_dot = -1.0f;
+	float ang1, ang2, ang_diff, power_sum, ang_mid;
+	InputMesh *im = &sys->input_mesh;
+
+	f1 = im->faces_edge[in_e][0];
+	f2 = im->faces_edge[in_e][1];
+
+	if (f1 == QR_NO_FACE || f2 == QR_NO_FACE)
+		return;
+
+	//if (sys->gfsys[0]->h[f1] < FLT_EPSILON && sys->gfsys[0]->h[f2] < FLT_EPSILON)
+	//	return;
+
+	sub_v3_v3v3(ev, im->co[im->edges[in_e][0]],
+					im->co[im->edges[in_e][1]]);
+	normalize_v3(ev);
+
+#if 0
+	ang1 = angle_signed_on_axis_normalized_fast_v3v3_v3(ev, sys->cf[f1][0], im->no[f1]);
+	while (ang1 > M_PI_4)
+		ang1 -= M_PI_2;
+
+	ang2 = angle_signed_on_axis_normalized_fast_v3v3_v3(ev, sys->cf[f2][0], im->no[f2]);
+	while (ang2 > M_PI_4)
+		ang2 -= M_PI_2;
+#endif
+
+	ang_diff = ang1 - ang2;
+
+	/*power_sum = sys->gfsys[0]->h[f1] + sys->gfsys[0]->h[f2];
+	ang_mid = ang_diff * (sys->gfsys[0]->h[f2] / power_sum);
+
+	rotate_normalized_v3_v3v3fl(tmp, sys->gfsys[0]->gf[f1], im->no[f1], ang_mid);
+	copy_v3_v3(sys->gfsys[0]->gf[f1], tmp);
+	rotate_normalized_v3_v3v3fl(tmp, sys->gfsys[0]->gf[f2], im->no[f2], -(ang_diff - ang_mid));
+	copy_v3_v3(sys->gfsys[0]->gf[f2], tmp);
+
+
+	if (sys->gfsys[0]->h[f2] > sys->gfsys[0]->h[f1]) {
+		sys->gfsys[0]->h[f1] = power_sum * 0.05f;
+	}
+	else {
+		sys->gfsys[0]->h[f2] = power_sum * 0.05f;
+	}*/
+}
+
+static void getCurvatures(QuadRemeshSystem *sys)
+{
+	int f, e;
+	float ev[3], angle, vec[3];
+	InputMesh *im = &sys->input_mesh;
+	Vector3f edge_vec;
+	Matrix3f edge_c, *face_curvatures = new Matrix3f[im->num_faces];
+
+	for (f = 0; f < im->num_faces; f++) {
+		face_curvatures[f].setZero();
+	}
+
+	for (e = 0; e < im->num_edges; e++) {
+		edge_c.setZero();
+
+		sub_v3_v3v3(ev, im->co[im->edges[e][0]],
+						im->co[im->edges[e][1]]);
+		
+		normalize_v3(ev);
+		angle = angle_normalized_v3v3(im->no[im->faces_edge[e][0]],
+									  im->no[im->faces_edge[e][1]]);
+		sub_v3_v3v3(vec, im->co[getVertexOppositeToEdge(&sys->input_mesh, e, im->faces_edge[e][0])], im->co[im->edges[e][0]]);
+		if (dot_v3v3(vec, im->no[im->faces_edge[e][1]]) > 0.0f) {
+			angle = -angle;
+		}
+
+		edge_vec << ev[0], ev[1], ev[2];
+		edge_c = angle * (edge_vec * edge_vec.transpose());
+
+		if (im->faces_edge[e][0] != QR_NO_FACE)
+			face_curvatures[im->faces_edge[e][0]] += edge_c;
+		if (im->faces_edge[e][1] != QR_NO_FACE)
+			face_curvatures[im->faces_edge[e][1]] += edge_c;
+	}
+
+	for (f = 0; f < im->num_faces; f++) {
+		//std::cout << "vert_c: " << std::endl << vert_c << std::endl;
+		EigenSolver<Matrix3f> es(face_curvatures[f]);
+		Matrix3f eigenv = es.pseudoEigenvectors();
+		Vector3cf eigenvalues = es.eigenvalues();
+		//std::cout << "Eigenvalues: " << std::endl << es.eigenvalues() << std::endl;
+		//std::cout << "Pseudo-eigenvectors: " << std::endl << es.pseudoEigenvectors() << std::endl;
+#if 0
+		int max = MAX3_PAIR(eigenvalues(0).real(), eigenvalues(1).real(), eigenvalues(2).real(), 0, 1, 2);
+		if (eigenvalues(max).real() > FLT_EPSILON) {
+			copy_v3_v3(sys->cf[f][0], eigenv.col(max).data());
+			normalize_v3(sys->cf[f][0]);
+			//sys->gfsys[0]->h[f] = SQUARE(eigenvalues(max).real());
+		}
+		else {
+			sub_v3_v3v3(sys->cf[f][0], im->co[im->faces[f][0]], im->co[im->faces[f][1]]);
+			normalize_v3(sys->cf[f][0]);
+			//sys->gfsys[0]->h[f] = 0.0f;
+		}
+#endif
+	}
+
+	delete[] face_curvatures;
+
+	for (int i = 0; i < 20; i++)
+		for (int e = 0; e < im->num_edges; e++) {
+			alignAtEdge(sys, e);
+		}
+}
+
+void getPrincipalCurvatures(QuadRemeshSystem *sys)
+{
+	int i;
+	InputMesh *im = &sys->input_mesh;
+
+	sys->has_solution = true;
+
+	sys->U_field = static_cast<float*>(MEM_mallocN(sizeof(float) * im->num_verts, "QuadRemeshUField"));
+
+	for (i = 0; i < im->num_verts; i++) {
+		sys->U_field[i] = 0.0f;
+	}
+
+	//sys->cf = static_cast<float(*)[4][3]>(MEM_mallocN(sizeof(float[4][3]) * im->num_faces, __func__));
+	//sys->gfsys[0]->h = static_cast<float*>(MEM_mallocN(sizeof(float) * im->num_faces, __func__));
+	//sys->gfsys[1]->h = static_cast<float*>(MEM_mallocN(sizeof(float) * im->num_faces, __func__));
+
+	getCurvatures(sys);
+
+	sys->is_alloc = true;
+}

intern/quadremesh/quadremesh_edge.c

@@ -0,0 +1,180 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Alexander Pinzon Fernandez
+ *                 Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+#include <float.h>
+
+#include "BLI_utildefines.h"
+#include "BLI_memarena.h"
+#include "BLI_linklist.h"
+#include "BLI_math.h"
+
+#include "quadremesh_util.h"
+
+static QREdgeLink *newQREdgeLink(OutputMesh *om, QRVertID in_v, float in_dist)
+{
+	QREdgeLink *newl;
+	newl = BLI_memarena_alloc(om->memarena, sizeof(QREdgeLink));
+
+	newl->v = in_v;
+	newl->dist = in_dist;
+	newl->elink = NULL;
+	newl->next = NULL;
+	newl->gfsysid = 0;
+
+	return newl;
+}
+
+static void appendOnQREdge(OutputMesh *om, QREdge *in_e, QRVertID in_v, float in_dist)
+{
+	QREdgeLink *newl;
+	newl = newQREdgeLink(om, in_v, in_dist);
+
+	if (in_e->v2 == NULL) {
+		in_e->v1 = in_e->v2 = newl;
+	}
+	else {
+		in_e->v2->next = newl;
+		in_e->v2 = newl;
+	}
+}
+
+static void prependOnQREdge(OutputMesh *om, QREdge *in_e, QRVertID in_v, float in_dist)
+{
+	QREdgeLink *newl;
+	newl = newQREdgeLink(om, in_v, in_dist);
+
+	newl->next = in_e->v1;
+	in_e->v1 = newl;
+
+	if (in_e->v2 == NULL)
+		in_e->v2 = newl;
+}
+
+static void insertAfterOnQREdge(OutputMesh *om, QREdge *in_e, QREdgeLink *in_l, QRVertID in_v, float in_dist)
+{
+	QREdgeLink *newl;
+	newl = newQREdgeLink(om, in_v, in_dist);
+
+	newl->next = in_l->next;
+	newl->gfsysid = in_l->gfsysid;
+
+	if (in_l->elink) {
+		unlinkVerts(om, in_l->elink);
+		in_l->elink = linkVerts(om, in_l->v, in_v);
+		newl->elink = linkVerts(om, in_v, in_l->next->v);
+	}
+
+	in_l->next = newl;
+}
+
+void insertOnQREdge(OutputMesh *om, QREdge *in_e, QRVertID in_vid)
+{
+	float tmp, vec[3];
+	QREdgeLink *it;
+
+	if (in_e->num_links == 0) {
+		appendOnQREdge(om, in_e, in_vid, 0.0f);
+		copy_v3_v3(in_e->orig, om->verts[in_vid].co);
+	}
+	else if (in_e->num_links == 1) {
+		sub_v3_v3v3(in_e->dir, om->verts[in_vid].co, in_e->orig);
+		tmp = len_v3(in_e->dir);
+
+		if (tmp < FLT_EPSILON)
+			return;
+
+		appendOnQREdge(om, in_e, in_vid, tmp);
+		mul_v3_fl(in_e->dir, 1.0f / tmp);
+	}
+	else {
+		sub_v3_v3v3(vec, om->verts[in_vid].co, in_e->orig);
+		tmp = dot_v3v3(vec, in_e->dir);
+
+		//if (IS_EQF(tmp, 0.0f))
+			//return;
+
+		if (tmp > in_e->v2->dist)
+			appendOnQREdge(om, in_e, in_vid, tmp);
+		else if (tmp < in_e->v1->dist)
+			prependOnQREdge(om, in_e, in_vid, tmp);
+		else {
+			for (it = in_e->v1; it->next && it->next->dist < tmp; it = it->next);
+			if (in_vid == it->v || in_vid == it->next->v) return;
+			//if (IS_EQF(it->dist, tmp) || IS_EQF(it->next->dist, tmp)) return;
+			insertAfterOnQREdge(om, in_e, it, in_vid, tmp);
+		}
+	}
+	in_e->num_links++;
+}
+
+void linkOnQREdge(OutputMesh *om, GFSysID sys_id, QREdge *in_e, QRVertID in_v1, QRVertID in_v2)
+{
+	QREdgeLink *it;
+
+	for (it = in_e->v1; it && it->next && it->v != in_v1; it = it->next)
+		if (it->v == in_v2) {
+			in_v2 = in_v1;
+			break;
+		}
+
+	for (; it && it->next && it->v != in_v2; it = it->next) {
+		if (!it->elink)
+			it->elink = linkVerts(om, it->v, it->next->v);
+
+		it->gfsysid |= sys_id;
+	}
+}
+
+QREdge *addQREdgeToFace(OutputMesh *om, InputMesh *im, GFSysID sys_id, int in_f, QRVertID in_v1, QRVertID in_v2)
+{
+	float isection[3];
+	QREdge *e, *newe;
+	QRVertID newv;
+	LinkNode *it;
+
+	newe = BLI_memarena_calloc(om->memarena, sizeof(QREdge));
+
+	insertOnQREdge(om, newe, in_v1);
+	insertOnQREdge(om, newe, in_v2);
+	linkOnQREdge(om, sys_id, newe, in_v1, in_v2);
+	BLI_linklist_prepend_arena(&om->ringf[in_f], (void*)newe, om->memarena);
+
+	for (it = om->ringf[in_f]->next; it; it = it->next) {
+		e = (QREdge*)it->link;
+
+		if ((e->v1->gfsysid & sys_id) != 0)
+			continue;
+
+		if (isect_seg_seg_unsafe_v3(isection,
+			                        om->verts[e->v1->v].co, om->verts[e->v2->v].co,
+					                om->verts[in_v1].co, om->verts[in_v2].co))
+		{
+			newv = addVert(om, isection, im->no[in_f]);
+			insertOnQREdge(om, e, newv);
+			insertOnQREdge(om, newe, newv);
+		}
+	}
+	
+	return newe;
+}

intern/quadremesh/quadremesh_harmonic.c

@@ -0,0 +1,225 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Alexander Pinzon Fernandez
+ *                 Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_object_types.h"
+#include "DNA_meshdata_types.h"
+
+#include "BLI_utildefines.h"
+#include "BLI_math.h"
+
+#include "BKE_deform.h"
+
+#include "ONL_opennl.h"
+#include "MOD_util.h"
+
+#include "quadremesh_util.h"
+
+static void initLaplacianMatrix(InputMesh *im)
+{
+	float v1[3], v2[3], v3[3];
+	float w2, w3;
+	int j, fi;
+	unsigned int idv1, idv2, idv3;
+
+	for (fi = 0; fi < im->num_faces; fi++) {
+		const unsigned int *vidf = im->faces[fi];
+
+		idv1 = vidf[0];
+		idv2 = vidf[1];
+		idv3 = vidf[2];
+
+		for (j = 0; j < 3; j++) {
+			idv1 = vidf[j];
+			idv2 = vidf[(j + 1) % 3];
+			idv3 = vidf[(j + 2) % 3];
+
+			copy_v3_v3(v1, im->co[idv1]);
+			copy_v3_v3(v2, im->co[idv2]);
+			copy_v3_v3(v3, im->co[idv3]);
+
+			w2 = cotangent_tri_weight_v3(v3, v1, v2);
+			w3 = cotangent_tri_weight_v3(v2, v3, v1);
+
+			if (im->constraints[idv1] == 1) {
+				nlMatrixAdd(idv1, idv1, w2 + w3);
+			}
+			else  {
+				nlMatrixAdd(idv1, idv2, -w2);
+				nlMatrixAdd(idv1, idv3, -w3);
+				nlMatrixAdd(idv1, idv1, w2 + w3);
+			}
+
+		}
+	}
+	
+}
+
+static bool computeScalarField(QuadRemeshSystem *sys)
+{
+	int i;
+	InputMesh *im = &sys->input_mesh;
+	
+	nlNewContext();
+
+	nlSolverParameteri(NL_NB_VARIABLES, im->num_verts);
+	nlSolverParameteri(NL_SYMMETRIC, NL_FALSE);
+	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
+	nlSolverParameteri(NL_NB_ROWS, im->num_verts);
+	nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 1);
+	nlBegin(NL_SYSTEM);
+	
+	for (i = 0; i < im->num_verts; i++) {
+		nlSetVariable(0, i, 0);
+	}
+		
+	nlBegin(NL_MATRIX);
+	
+	initLaplacianMatrix(im);
+
+	for (i = 0; i < im->num_verts; i++) {
+		if (im->constraints[i] == 1) {
+			nlRightHandSideSet(0, i, im->weights[i]);
+		}
+		else {
+			nlRightHandSideSet(0, i, 0);
+		}
+	}
+	nlEnd(NL_MATRIX);
+	nlEnd(NL_SYSTEM);
+	
+	return nlSolveAdvanced(NULL, NL_TRUE);
+}
+
+/**
+* Compute the gradient fields
+*
+* xi, xj, xk, are the vertices of the face
+* ui, uj, uk, are the values of scalar fields for every vertex of the face
+* n is the normal of the face.
+* gf1 is the unknown field gradient 1.
+* gf2 is the unknown field gradient 2.
+*
+* |xj - xi|         |uj - ui|
+* |xk - xj| * gf1 = |uk - uj|
+* |   nf  |         |   0   |
+*
+* gf2 = cross(n, gf1)
+*/
+static void computeGradientFields(QuadRemeshSystem * sys)
+{
+	int fi, i, j, k;
+	float val, a[3][3], u[3], inv_a[3][3], gf1[3], g[3], w[3];
+	InputMesh *im = &sys->input_mesh;
+
+	sys->gfsys[0]->gf = MEM_mallocN(sizeof(float[3]) * im->num_faces, "QuadRemeshGradientField1");
+	sys->gfsys[1]->gf = MEM_mallocN(sizeof(float[3]) * im->num_faces, "QuadRemeshGradientField2");
+
+	for (fi = 0; fi < im->num_faces; fi++) {
+		const unsigned int *vidf = im->faces[fi];
+		i = vidf[0];
+		j = vidf[1];
+		k = vidf[2];
+		sub_v3_v3v3(a[0], im->co[j], im->co[i]);
+		sub_v3_v3v3(a[1], im->co[k], im->co[j]);
+		copy_v3_v3(a[2], im->no[fi]);
+
+		/* Correct way*/
+		transpose_m3(a);
+		u[0] = sys->U_field[j] - sys->U_field[i];
+		u[1] = sys->U_field[k] - sys->U_field[j];
+		u[2] = 0;
+		invert_m3_m3(inv_a, a);
+		//mul_v3_m3v3(sys->gf1[fi], inv_a, u);
+		mul_v3_m3v3(gf1, inv_a, u);
+
+		/* Project Gradient fields on face*/
+		normalize_v3_v3(g, gf1);
+		val = dot_v3v3(g, im->no[fi]);
+		mul_v3_v3fl(u, im->no[fi], val);
+		sub_v3_v3v3(w, g, u);
+		normalize_v3_v3(sys->gfsys[0]->gf[fi], w);
+
+		cross_v3_v3v3(g, im->no[fi], sys->gfsys[0]->gf[fi]);
+		normalize_v3_v3(sys->gfsys[1]->gf[fi], g);
+		//cross_v3_v3v3(sys->gf2[fi], im->no[fi], sys->gf1[fi]);
+	}
+}
+
+void getHarmonicGradients(QuadRemeshSystem *sys)
+{
+	int i;
+	InputMesh *im = &sys->input_mesh;
+
+#ifdef OPENNL_THREADING_HACK
+	modifier_opennl_lock();
+#endif
+
+	sys->has_solution = computeScalarField(sys);
+
+	if (sys->has_solution) {
+		sys->U_field = MEM_mallocN(sizeof(float) * im->num_verts, "QuadRemeshUField");
+
+		for (i = 0; i < im->num_verts; i++) {
+			sys->U_field[i] = nlGetVariable(0, i);
+		}
+
+		sys->gfsys[0] = newGradientFlowSystem(sys);
+		sys->gfsys[1] = newGradientFlowSystem(sys);
+		computeGradientFields(sys);
+
+		sys->is_alloc = true;
+	}
+
+#ifdef OPENNL_THREADING_HACK
+	modifier_opennl_unlock();
+#endif
+}
+
+void getUField(QuadRemeshSystem *sys, MDeformVert *dvert, int defgrp_index)
+{
+	int i;
+	float y, mmin = 1000.0f, mmax = -1000.0f;
+	MDeformVert *dv = NULL;
+
+	BLI_assert(dvert != NULL);
+
+	if (sys->has_solution) {
+		dv = dvert;
+
+		for (i = 0; i < sys->input_mesh.num_verts; i++) {
+			mmin = min_ff(mmin, sys->U_field[i]);
+			mmax = max_ff(mmax, sys->U_field[i]);
+		}
+
+		for (i = 0; i < sys->input_mesh.num_verts; i++) {
+			y = (sys->U_field[i] - mmin) / (mmax - mmin);
+			//x = y * 60;
+			//y = (x % 2 == 0 ? 0.1 : 0.9);
+			defvert_add_index_notest(dv, defgrp_index, y);
+			dv++;
+		}
+	}
+}

intern/quadremesh/quadremesh_input.c

@@ -0,0 +1,289 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Alexander Pinzon Fernandez
+ *                 Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_object_types.h"
+
+#include "BLI_math.h"
+
+#include "BKE_mesh_mapping.h"
+#include "BKE_DerivedMesh.h"
+#include "BKE_deform.h"
+
+#include "MOD_util.h"
+
+#include "quadremesh_util.h"
+
+static void createFaceRingMap(InputMesh *im)
+{
+	int i, j, totalr = 0;
+	int *index_iter;
+
+	im->ringf_map = MEM_callocN(sizeof(MeshElemMap) * im->num_verts, "DeformRingMap");
+
+	for (i = 0; i < im->num_faces; i++) {
+		for (j = 0; j < 3; j++) {
+			im->ringf_map[im->faces[i][j]].count++;
+			totalr++;
+		}
+	}
+
+	im->ringf_indices = MEM_callocN(sizeof(int) * totalr, "DeformRingIndex");
+	index_iter = im->ringf_indices;
+	for (i = 0; i < im->num_verts; i++) {
+		im->ringf_map[i].indices = index_iter;
+		index_iter += im->ringf_map[i].count;
+		im->ringf_map[i].count = 0;
+	}
+
+	for (i = 0; i < im->num_faces; i++) {
+		for (j = 0; j < 3; j++) {
+			MeshElemMap *map = &im->ringf_map[im->faces[i][j]];
+			map->indices[map->count++] = i;
+		}
+	}
+}
+
+static void createEdgeRingMap(InputMesh *im)
+{
+	int i, totalr = 0;
+	int *index_iter;
+
+	im->ringe_map = MEM_callocN(sizeof(MeshElemMap) * im->num_verts, "DeformNeighborsMap");
+
+	for (i = 0; i < im->num_edges; i++) {
+		im->ringe_map[im->edges[i][0]].count++;
+		im->ringe_map[im->edges[i][1]].count++;
+		totalr += 2;
+	}
+
+	im->ringe_indices = MEM_callocN(sizeof(int) * totalr, "DeformNeighborsIndex");
+	index_iter = im->ringe_indices;
+	for (i = 0; i < im->num_verts; i++) {
+		im->ringe_map[i].indices = index_iter;
+		index_iter += im->ringe_map[i].count;
+		im->ringe_map[i].count = 0;
+	}
+
+	for (i = 0; i < im->num_edges; i++) {
+		MeshElemMap *map1 = &im->ringe_map[im->edges[i][0]];
+		MeshElemMap *map2 = &im->ringe_map[im->edges[i][1]];
+		map1->indices[map1->count++] = i;
+		map2->indices[map2->count++] = i;
+	}
+}
+
+static void createVertRingMap(InputMesh *im)
+{
+	int i, totalr = 0;
+	int *index_iter;
+
+	im->ringv_map = MEM_callocN(sizeof(MeshElemMap) * im->num_verts, "DeformNeighborsMap");
+
+	for (i = 0; i < im->num_edges; i++) {
+		im->ringv_map[im->edges[i][0]].count++;
+		im->ringv_map[im->edges[i][1]].count++;
+		totalr += 2;
+	}
+
+	im->ringv_indices = MEM_callocN(sizeof(int) * totalr, "DeformNeighborsIndex");
+	index_iter = im->ringv_indices;
+	for (i = 0; i < im->num_verts; i++) {
+		im->ringv_map[i].indices = index_iter;
+		index_iter += im->ringv_map[i].count;
+		im->ringv_map[i].count = 0;
+	}
+
+	for (i = 0; i < im->num_edges; i++) {
+		MeshElemMap *map1 = &im->ringv_map[im->edges[i][0]];
+		MeshElemMap *map2 = &im->ringv_map[im->edges[i][1]];
+		map1->indices[map1->count++] = im->edges[i][1];
+		map2->indices[map2->count++] = im->edges[i][0];
+	}
+}
+
+static void createFacesByEdge(InputMesh *im){
+	int e, i, j, v1, v2;
+	unsigned int *vin;
+
+	im->faces_edge = MEM_mallocN(sizeof(int[2]) * im->num_edges, "QuadRemeshFacesEdge");
+
+	for (e = 0; e < im->num_edges; e++) {
+		v1 = im->edges[e][0];
+		v2 = im->edges[e][1];
+
+		im->faces_edge[e][0] = QR_NO_FACE;
+		im->faces_edge[e][1] = QR_NO_FACE;
+
+		for (i = 0, j = 0; i < im->ringf_map[v1].count && j < 2; i++) {
+			vin = im->faces[im->ringf_map[v1].indices[i]];
+			if (vin[0] == v2 || vin[1] == v2 || vin[2] == v2) {
+				im->faces_edge[e][j++] = im->ringf_map[v1].indices[i];
+			}
+		}
+	}
+}
+
+void freeInputMesh(InputMesh *im)
+{
+	MEM_SAFE_FREE(im->faces);
+	MEM_SAFE_FREE(im->edges);
+	MEM_SAFE_FREE(im->faces_edge);
+	MEM_SAFE_FREE(im->co);
+	MEM_SAFE_FREE(im->no);
+	MEM_SAFE_FREE(im->vno);
+	MEM_SAFE_FREE(im->ev);
+	MEM_SAFE_FREE(im->constraints);
+	MEM_SAFE_FREE(im->weights);
+
+	MEM_SAFE_FREE(im->ringf_indices);
+	MEM_SAFE_FREE(im->ringe_indices);
+	MEM_SAFE_FREE(im->ringv_indices);
+	MEM_SAFE_FREE(im->ringf_map);
+	MEM_SAFE_FREE(im->ringe_map);
+	MEM_SAFE_FREE(im->ringv_map);
+
+	im->is_alloc = false;
+}
+
+static void getInputMeshData(InputMesh *im, DerivedMesh *dm)
+{
+	int i;
+	MVert *arrayvect;
+	MEdge *arrayedge;
+	MFace *tessface;
+
+	DM_ensure_tessface(dm);
+
+	/* Get vertices */
+	im->num_verts = dm->getNumVerts(dm);
+	im->co = MEM_mallocN(sizeof(float[3]) * im->num_verts, "QuadRemeshCoordinates");
+	arrayvect = dm->getVertArray(dm);
+	for (i = 0; i < dm->getNumVerts(dm); i++) {
+		copy_v3_v3(im->co[i], arrayvect[i].co);
+	}
+
+	/* Get edges */
+	im->edges = MEM_mallocN(sizeof(int[2]) * dm->getNumEdges(dm) * 2, "QuadRemeshEdges");
+	arrayedge = dm->getEdgeArray(dm);
+	for (im->num_edges = 0, i = 0; i < dm->getNumEdges(dm); i++) {
+		im->edges[im->num_edges][0] = arrayedge[i].v1;
+		im->edges[im->num_edges][1] = arrayedge[i].v2;
+		im->num_edges++;
+	}
+
+	/* Edge vectors */
+	im->ev = MEM_mallocN(sizeof(float[3]) * im->num_edges, "QuadRemeshEdgeVectors");
+	for (i = 0; i < im->num_edges; i++) {
+		sub_v3_v3v3(im->ev[i], im->co[im->edges[i][1]], im->co[im->edges[i][0]]);
+		normalize_v3(im->ev[i]);
+	}
+
+	/* Get faces */
+	tessface = dm->getTessFaceArray(dm);
+	im->faces = MEM_mallocN(sizeof(int[3]) * dm->getNumTessFaces(dm) * 2, "QuadRemeshFaces");
+	for (im->num_faces = 0, i = 0; i < dm->getNumTessFaces(dm); i++) {
+		im->faces[im->num_faces][0] = tessface[i].v1;
+		im->faces[im->num_faces][1] = tessface[i].v2;
+		im->faces[im->num_faces][2] = tessface[i].v3;
+		im->num_faces++;
+
+		if (tessface[i].v4 != 0) {
+			im->faces[im->num_faces][0] = tessface[i].v1;
+			im->faces[im->num_faces][1] = tessface[i].v3;
+			im->faces[im->num_faces][2] = tessface[i].v4;
+			im->num_faces++;
+
+			im->edges[im->num_edges][0] = tessface[i].v1;
+			im->edges[im->num_edges][1] = tessface[i].v3;
+			im->num_edges++;
+		}
+	}
+
+	/* Compute face and vertex normals */
+	im->no = MEM_callocN(sizeof(float[3]) * im->num_faces, "QuadRemeshNormals");
+	for (i = 0; i < im->num_faces; i++) {
+		normal_tri_v3(im->no[i], im->co[im->faces[i][0]],
+					             im->co[im->faces[i][1]],
+					             im->co[im->faces[i][2]]);
+	}
+
+	im->vno = MEM_callocN(sizeof(float[3]) * im->num_verts, "QuadRemeshNormals");
+	for (i = 0; i < im->num_verts; i++) {
+		zero_v3(im->vno[i]);
+	}
+	for (i = 0; i < im->num_faces; i++) {
+		add_v3_v3(im->vno[im->faces[i][0]], im->no[i]);
+		add_v3_v3(im->vno[im->faces[i][1]], im->no[i]);
+		add_v3_v3(im->vno[im->faces[i][2]], im->no[i]);
+	}
+	for (i = 0; i < im->num_verts; i++) {
+		normalize_v3(im->vno[i]);
+	}
+}
+
+static void getFeatures(InputMesh *im, MDeformVert *dvert, int defgrp_index)
+{
+	int i;
+	float wpaint;
+	MDeformVert *dv = NULL;
+	MDeformWeight *d;
+
+	BLI_assert(dvert != NULL);
+
+	/* Get features */
+	im->constraints = MEM_callocN(sizeof(int) * im->num_verts, __func__);
+	im->weights = MEM_callocN(sizeof(float) * im->num_verts, __func__);
+
+	dv = dvert;
+	for (i = 0, im->num_features = 0; i < im->num_verts; i++, dv++) {
+		d = defvert_find_index(dv, defgrp_index);
+		wpaint = defvert_find_weight(dv, defgrp_index);
+
+		if (d && (wpaint < 0.19 || wpaint > 0.89)) {
+			im->constraints[i] = 1;
+			im->weights[i] = -1.0f + wpaint * 2.0f;
+			im->num_features++;
+		}
+	}
+}
+
+void getInput(QuadRemeshSystem *sys, Object *ob, DerivedMesh *dm)
+{
+	int defgrp_index;
+	MDeformVert *dvert;
+
+	modifier_get_vgroup(ob, dm, sys->qmd->anchor_grp_name, &dvert, &defgrp_index);
+
+	getInputMeshData(&sys->input_mesh, dm);
+	getFeatures(&sys->input_mesh, dvert, defgrp_index);
+
+	createFaceRingMap(&sys->input_mesh);
+	createEdgeRingMap(&sys->input_mesh);
+	createVertRingMap(&sys->input_mesh);
+	createFacesByEdge(&sys->input_mesh);
+
+	sys->input_mesh.is_alloc = true;
+}

intern/quadremesh/quadremesh_line.c

@@ -0,0 +1,539 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Alexander Pinzon Fernandez
+ *                 Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_modifier_types.h"
+
+#include "BLI_utildefines.h"
+#include "BLI_math.h"
+#include "BLI_heap.h"
+#include "BLI_rand.h"
+
+#include "BKE_mesh_mapping.h"
+
+#include "quadremesh_util.h"
+
+#define QR_GFLINE_QSIZE 64
+#define QR_LINELIMIT 100000
+//#define QR_SHOWQUERIES
+
+typedef struct GFLine {
+	struct GradientFlowSystem *gfsys;
+	GFPoint end, seed;
+	GFPoint *oldp, *lastchkp;
+	GFSegment lastchks;
+	int d;
+
+	float lastchklen;
+
+	float qlen;
+	int num_q;
+	GFSegment q[QR_GFLINE_QSIZE];
+} GFLine;
+
+
+static int getEdgeFromVerts(InputMesh *im, int v1, int v2)
+{
+	int *eidn, nume, i;
+
+	nume = im->ringe_map[v1].count;
+	eidn = im->ringe_map[v1].indices;
+
+	for (i = 0; i < nume; i++) {
+		if (im->edges[eidn[i]][0] == v2 || im->edges[eidn[i]][1] == v2){
+			return eidn[i];
+		}
+	}
+
+	return -1;
+}
+
+static int getOtherFaceAdjacentToEdge(InputMesh *im, int in_f, int in_e)
+{
+	if (im->faces_edge[in_e][0] == in_f) {
+		return im->faces_edge[in_e][1];
+	}
+
+	return im->faces_edge[in_e][0];
+}
+
+static void makePoint(GFPoint *r_p, GFPointType p_type, int p_vef, float p_co[3])
+{
+	r_p->type = p_type;
+	r_p->v = p_vef;
+	r_p->id = -1;
+	copy_v3_v3(r_p->co, p_co);
+}
+
+static void makeSegment(GFSegment *r, InputMesh *im, GFPoint *prev, int in_f)
+{
+	if (prev->type == eEdge && r->p.type == eEdge) {
+		if (prev->e == r->p.e) { /* eEdgeSegment */
+			r->type = eEdgeSegment;
+			r->e = r->p.e;
+		}
+		else { /* eFaceSegment */
+			r->type = eFaceSegment;
+			r->f = in_f; //getFaceFromTwoEdges(im, prev->e, p_vef);
+		}
+	}
+	else if (prev->type == eVert && r->p.type == eEdge) { /* likely eFaceSegment */
+		if (im->edges[r->p.e][0] == prev->v || im->edges[r->p.e][1] == prev->v) {
+			r->type = eEdgeSegment;
+			r->e = r->p.e;
+		}
+		else {
+			r->type = eFaceSegment;
+			r->f = in_f;
+		}
+	}
+	else if (prev->type == eEdge && r->p.type == eVert) { /* dunno */
+		if (im->edges[prev->e][0] == r->p.v || im->edges[prev->e][1] == r->p.v) {
+			r->type = eEdgeSegment;
+			r->e = prev->e;
+		}
+		else {
+			r->type = eFaceSegment;
+			r->f = in_f;
+		}
+	}
+	else if (prev->type == eVert && r->p.type == eVert) { /* eEdgeSegment */
+		r->type = eEdgeSegment;
+		r->e = getEdgeFromVerts(im, prev->v, r->p.v);
+	}
+	else { /* prev on a face or this on a face */
+		r->type = eFaceSegment;
+		r->f = in_f; //p_type == eFace ? p_vef : prev->f;
+	}
+}
+
+/**
+ * /return true - all ok
+ *         false - error
+ */
+static bool nextPointOnFace(GFPoint *r_p, InputMesh *im, GFPoint *in_p, int in_f, float in_dir[3])
+{
+	int i, pick = -1, v = -1, e;
+	bool is_on_vertex = false;
+	float a[3][3], b[3][3], c[2][3], co2[3], result[3], dummy[3];
+
+	/* check if direction is coplanar to triangle */
+	/* check if point is inside triangle */
+	/* check if triangle is degenerate */
+
+	add_v3_v3v3(co2, in_dir, in_p->co); /* second point on direction */
+
+	for (i = 0; i < 3; i++) {
+		sub_v3_v3v3(a[i], in_p->co, im->co[im->faces[in_f][i]]);
+		if (dot_v3v3(a[i], a[i]) < FLT_EPSILON)
+			v = i;
+	}
+
+	if (v != -1) {
+		is_on_vertex = true;
+		normalize_v3_v3(c[0], a[(v + 1) % 3]);
+		normalize_v3_v3(c[1], a[(v + 2) % 3]);
+		add_v3_v3(c[0], c[1]);
+		mul_v3_v3fl(a[v], c[0], -0.5f);
+	}
+
+	for (i = 0; i < 3; i++) {
+		cross_v3_v3v3(b[i], in_dir, a[i]);
+		dummy[i] = dot_v3v3(b[i], im->no[in_f]);
+	}
+
+	for (i = 0; i < 3; i++)
+		if (dummy[i] < 0.0f && dummy[(i + 1) % 3] >= 0.0f)
+			pick = i;
+
+	if (pick == -1)
+		return false;
+
+	e = getEdgeFromVerts(im, im->faces[in_f][pick], im->faces[in_f][(pick + 1) % 3]);
+
+	isect_line_line_v3(im->co[im->faces[in_f][pick]],
+					   im->co[im->faces[in_f][(pick + 1) % 3]],
+					   in_p->co, co2, result, dummy);
+
+	if (len_squared_v3v3(result, im->co[im->edges[e][0]]) < FLT_EPSILON)
+		makePoint(r_p, eVert, im->edges[e][0], im->co[im->edges[e][0]]);
+	else if (len_squared_v3v3(result, im->co[im->edges[e][1]]) < FLT_EPSILON)
+		makePoint(r_p, eVert, im->edges[e][1], im->co[im->edges[e][1]]);
+	else
+		makePoint(r_p, eEdge, e, result);
+
+	return true;
+}
+
+/**
+ * 0 - intersection found
+ * 1 - intersection not found
+ * 2 - wrong direction for this face
+ */
+static int queryDirection(GradientFlowSystem *gfsys, GFPoint *in_p, int in_f, float in_dir[3],
+						  float dist, float maxdist, bool make_seed)
+{
+	int oldf;
+	float c[3], len, actlen, prevlen, chkco[3], seedco[3], dir[3];
+	InputMesh *im = &gfsys->sys->input_mesh;
+	OutputMesh *om = &gfsys->sys->output_mesh;
+	GFPoint oldp, newp;
+
+	copy_v3_v3(dir, in_dir);
+	memcpy(&oldp, in_p, sizeof(GFPoint));
+	oldf = in_f;
+	actlen = prevlen = 0.0f;
+
+	while (in_f != QR_NO_FACE && oldp.type != eVert) {
+		project_plane_v3_v3v3(dir, dir, im->no[in_f]);
+
+		if (normalize_v3(dir) < FLT_EPSILON) 
+			break;
+
+		if (dot_v3v3(im->no[oldf], im->no[in_f]) < 0.0f)
+			mul_v3_fl(dir, -1.0f);
+
+		if (!nextPointOnFace(&newp, im, &oldp, in_f, dir))
+			break;
+
+		sub_v3_v3v3(c, newp.co, oldp.co);
+		len = len_v3(c);
+		actlen = prevlen + len;
+
+		if (actlen > maxdist && make_seed) {
+			madd_v3_v3v3fl(seedco, oldp.co, c, (maxdist - prevlen) / len);
+
+			addSeedToQueue(gfsys->seeds, seedco, eFace, in_f, -maxdist);
+		}
+		else
+			copy_v3_v3(seedco, newp.co);
+
+		if (actlen > dist && prevlen < dist)
+			madd_v3_v3v3fl(chkco, oldp.co, c, (dist - prevlen) / len);
+		else
+			copy_v3_v3(chkco, newp.co);
+
+#ifdef QR_SHOWQUERIES
+		int vf1, vf2;
+		if (make_seed) {
+			vf1 = addVert(&gfsys->sys->output_mesh, oldp.co, NULL);
+			vf2 = addVert(&gfsys->sys->output_mesh, seedco, NULL);
+		}
+		else {
+			vf1 = addVert(&gfsys->sys->output_mesh, oldp.co, NULL);
+			vf2 = addVert(&gfsys->sys->output_mesh, chkco, NULL);
+		}
+		linkVerts(&gfsys->sys->output_mesh, vf1, vf2);
+#endif
+
+		if (prevlen < dist) {
+			if (isectSegmentWithOthersOnFace(om, gfsys->id, oldp.co, chkco, in_f))
+				return 0;
+
+			if (newp.type == eEdge && actlen <= dist) {
+				if (isectPointWithQREdge(om, gfsys->id, chkco, newp.e))
+					return 0;
+			}
+		}
+		
+		if ((actlen > dist && !make_seed) || actlen > maxdist)
+			return 1;
+
+		if (newp.type == eVert)
+			break;
+
+		oldf = in_f;
+		in_f = getOtherFaceAdjacentToEdge(im, in_f, newp.e);
+		memcpy(&oldp, &newp, sizeof(GFPoint));
+		prevlen = actlen;
+	}
+
+	return 2;
+}
+
+static bool checkPoint(GradientFlowSystem *gfsys, float in_oldco[3], GFPoint *in_p, float dist, float maxdist)
+{
+	bool make_seed = BLI_rng_get_float(gfsys->sys->rng) < QR_SEEDPROB;
+
+	int d, f;
+	float seg[3], dir[3], no[3];
+	InputMesh *im = &gfsys->sys->input_mesh;
+
+	sub_v3_v3v3(seg, in_oldco, in_p->co);
+
+	if (len_squared_v3(seg) < FLT_EPSILON)
+		return true;
+
+	if (in_p->type == eEdge) {
+		getNormalAtEdge(no, im, in_p->e);
+
+		if ((f = im->faces_edge[in_p->e][0]) == QR_NO_FACE)	
+			if ((f = im->faces_edge[in_p->e][1]) == QR_NO_FACE)
+				return true;
+	}
+	else {
+		copy_v3_v3(no, im->no[in_p->f]);
+		f = in_p->f;
+	}
+
+	cross_v3_v3v3(dir, no, seg);
+	normalize_v3(dir);
+
+	for (d = 0; d < 2; d++) {
+		if (!queryDirection(gfsys, in_p, f, dir, dist, maxdist, make_seed))
+			return false;
+		mul_v3_fl(dir, -1.0f);
+	}
+
+	return true;
+}
+
+static void addSegmentToLine(GFLine *line, GFSegment *in_s)
+{
+	int i;
+	OutputMesh *om = &line->gfsys->sys->output_mesh;
+	InputMesh *im = &line->gfsys->sys->input_mesh;
+
+	addGFPoint(im, om, &in_s->p);
+
+	if (in_s->p.type == eVert) {
+		for (i = 0; i < im->ringe_map[in_s->p.v].count; i++)
+			insertOnQREdge(om, &om->ringe[im->ringe_map[in_s->p.v].indices[i]], in_s->p.id);
+	}
+	else if (in_s->p.type == eEdge)
+		insertOnQREdge(om, &om->ringe[in_s->p.e], in_s->p.id);
+
+	if (in_s->type == eEdgeSegment)
+		linkOnQREdge(om, line->gfsys->id, &om->ringe[in_s->e], line->end.id, in_s->p.id);
+	else
+		addQREdgeToFace(om, im, line->gfsys->id, in_s->f, line->end.id, in_s->p.id);
+
+	memcpy(&line->end, &in_s->p, sizeof(GFPoint));
+}
+
+static void flushGFLineQueue(GFLine *line)
+{
+	int i;
+
+	for (i = 0; i < line->num_q; i++)
+		addSegmentToLine(line, &line->q[i]);
+
+	line->num_q = 0;
+}
+
+static bool enqueueSegment(GFLine *line, GFSegment *in_s)
+{
+	if (line->num_q == QR_GFLINE_QSIZE)
+		return false;
+
+	memcpy(&line->q[line->num_q], in_s, sizeof(GFSegment));
+	line->num_q++;
+
+	return true;
+}
+
+static bool changeLineDirection(GFLine *line)
+{
+	flushGFLineQueue(line);
+
+	if (line->d == 0 && line->end.id > line->seed.id + 3 &&
+		line->end.type == eFace && line->seed.type == eFace &&
+		line->end.f == line->seed.f)
+	{
+		addQREdgeToFace(&line->gfsys->sys->output_mesh,
+						&line->gfsys->sys->input_mesh,
+						line->gfsys->id, line->end.f,
+						line->end.id, line->seed.id);
+		return false;
+	}
+
+	/* reset to original state */
+	memcpy(&line->end, &line->seed, sizeof(GFPoint));
+	line->lastchkp = &line->seed;
+	line->oldp = &line->seed;
+	line->lastchklen = line->qlen = 0.0f;
+	line->lastchks.type = eNoSegment;
+
+	return ++line->d != 2;
+}
+
+static bool initGFLine(GradientFlowSystem *gfsys, GFLine *line, GFPoint *in_seed)
+{
+	const float maxdist = gfsys->sys->qmd->max_line_dist;
+	const float seeddist = maxdist * QR_MAXDIST_TO_SEEDDIST;
+
+	int i;
+	float old[3];
+	OutputMesh *om = &gfsys->sys->output_mesh;
+	InputMesh *im = &gfsys->sys->input_mesh;
+
+	if (in_seed->type == eFace) {
+		add_v3_v3v3(old, in_seed->co, gfsys->gf[in_seed->f]);
+		if (!checkPoint(gfsys, old, in_seed, maxdist, seeddist))
+			return false;
+	}
+		
+	addGFPoint(im, om, in_seed);
+
+	memcpy(&line->seed, in_seed, sizeof(GFPoint));
+	memcpy(&line->end, in_seed, sizeof(GFPoint));
+	line->lastchkp = &line->seed;
+	line->oldp = &line->seed;
+	line->gfsys = gfsys;
+	line->lastchks.type = eNoSegment;
+
+	line->d = line->num_q = 0;
+	line->lastchklen = line->qlen = 0.0f;
+
+	if (line->seed.type == eVert) {
+		for (i = 0; i < im->ringe_map[line->seed.v].count; i++)
+			insertOnQREdge(om, &om->ringe[im->ringe_map[line->seed.v].indices[i]], line->seed.id);
+	}
+
+	return true;
+}
+
+static bool nextLineSegment(GFLine *line, GFSegment *in_s)
+{
+	const float chklen = QR_SAMPLING_RATE;
+	const float maxdist = line->gfsys->sys->qmd->max_line_dist;
+	const float seeddist = maxdist * QR_MAXDIST_TO_SEEDDIST;//(1.0f + BLI_rng_get_float(line->gfsys->sys->rng) * (QR_MAXDIST_TO_SEEDDIST - 1.0f));
+
+	float seg[3], co[3];
+	float curlen;
+	GFSegment newchks;
+	InputMesh *im = &line->gfsys->sys->input_mesh;
+
+	/* qco[0] - first point after last checked
+	 * qco[num_q - 1] - last added point */
+
+	sub_v3_v3v3(seg, in_s->p.co, line->oldp->co);
+	curlen = len_v3(seg);
+
+	while (line->qlen + curlen > line->lastchklen + chklen) {
+		mul_v3_v3fl(co, seg, (line->lastchklen + chklen - line->qlen) / curlen);
+		add_v3_v3(co, line->oldp->co);
+
+		if (in_s->type == eEdgeSegment)
+			makePoint(&newchks.p, eEdge, in_s->e, co);
+		else
+			makePoint(&newchks.p, eFace, in_s->f, co);
+
+		makeSegment(&newchks, im, line->oldp, in_s->f);
+
+		if (!checkPoint(line->gfsys, line->lastchkp->co, &newchks.p, maxdist, seeddist)) {
+			if (line->lastchks.type != eNoSegment)
+				addSegmentToLine(line, &line->lastchks);
+			line->num_q = 0;
+
+			return false;
+		}
+		
+		flushGFLineQueue(line);
+		memcpy(&line->lastchks, &newchks, sizeof(GFSegment));
+		line->lastchkp = &line->lastchks.p;
+		line->lastchklen += chklen;
+	}
+
+	if(!enqueueSegment(line, in_s))
+		return false;
+
+	line->qlen += curlen;
+	line->oldp = &line->q[line->num_q - 1].p;
+
+	return true;
+}
+
+static void computeGFLine(GFLine *line)
+{
+	int i, f = QR_NO_FACE;
+	float gf[3], dir = 1.0f;
+	QuadRemeshSystem *sys = line->gfsys->sys;
+	InputMesh *im = &sys->input_mesh;
+	GFSegment news;
+
+	do {
+		if (line->oldp->type == eFace)
+			f = line->oldp->f;
+
+		do {
+			if (line->oldp->type == eVert) {
+				for (i = 0; i < im->ringf_map[line->oldp->v].count; i++) {
+					f = im->ringf_map[line->oldp->v].indices[i];
+
+					mul_v3_v3fl(gf, line->gfsys->gf[f], dir);
+					if (!nextPointOnFace(&news.p, im, line->oldp, f, gf))
+						continue;
+
+					if (news.p.type == eEdge)
+						break;
+				}
+				if (news.p.type != eEdge)
+					break;
+			}
+			else {
+				mul_v3_v3fl(gf, line->gfsys->gf[f], dir);
+				if (!nextPointOnFace(&news.p, im, line->oldp, f, gf))
+					break;
+				
+				if (news.p.type == eEdge && line->oldp->type == eEdge && line->oldp->e == news.p.e) {
+					if (dir * sys->U_field[im->edges[news.p.e][0]] < dir * sys->U_field[im->edges[news.p.e][1]])
+						makePoint(&news.p, eVert, im->edges[news.p.e][0], im->co[im->edges[news.p.e][0]]);
+					else
+						makePoint(&news.p, eVert, im->edges[news.p.e][1], im->co[im->edges[news.p.e][1]]);
+				}
+			}
+
+			makeSegment(&news, im, line->oldp, f);
+			if (!nextLineSegment(line, &news))
+				break;
+				
+			if (line->oldp->type == eEdge)
+				f = getOtherFaceAdjacentToEdge(im, f, line->oldp->e);
+		} while (f != QR_NO_FACE);
+
+		dir = -dir;
+	} while (changeLineDirection(line));
+}
+
+void computeFlowLines(QuadRemeshSystem *sys) {
+	GFPoint *seed;
+	GFLine line;
+	int s, comp;
+	
+	for (s = 0; s < 2; s++) {
+		getInitialSeeds(sys->gfsys[s]);
+
+		comp = 0;
+		while (!BLI_heap_is_empty(sys->gfsys[s]->seeds)) {
+			seed = getTopSeedFromQueue(sys->gfsys[s]->seeds);
+
+			if (++comp < QR_LINELIMIT && initGFLine(sys->gfsys[s], &line, seed))
+				computeGFLine(&line);
+
+			MEM_SAFE_FREE(seed);
+		}
+	}
+}

intern/quadremesh/quadremesh_link.c

@@ -0,0 +1,179 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Alexander Pinzon Fernandez
+ *                 Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+#include "BLI_utildefines.h"
+#include "BLI_memarena.h"
+#include "BLI_math.h"
+
+#include "quadremesh_util.h"
+
+static QRDiskLink *getLink(OutputMesh *om, QRVertID in_v1, QRVertID in_v2)
+{
+	QRDiskLink *it;
+
+	it = om->verts[in_v1].link;
+	while(it) {
+		if (it->v == in_v2)
+			return it;
+		
+		it = it->next;
+		if (it == om->verts[in_v1].link)
+			break;
+	}
+
+	return NULL;
+}
+
+static QRDiskLink *insertLink(OutputMesh *om, QRVertID in_a, QRVertID in_b)
+{
+	float vec[3];
+	QRDiskLink *it, *l;
+
+	l = BLI_memarena_alloc(om->memarena, sizeof(QRDiskLink));
+
+	l->e = -1;
+	l->v = in_b;
+	l->poly_on_right = false;
+
+	sub_v3_v3v3(vec, om->verts[in_b].co, om->verts[in_a].co);
+	project_plane_v3_v3v3(vec, vec, om->verts[in_a].no);
+	normalize_v3(vec);
+
+	if (om->verts[in_a].num_links == 0) {
+		om->verts[in_a].link = l;
+		copy_v3_v3(om->verts[in_a].vec, vec);
+
+		l->next = l;
+		l->prev = l;
+		l->ang = 0.0f;
+	}
+	else {
+		l->ang = angle_signed_on_axis_normalized_fast_v3v3_v3(om->verts[in_a].vec, vec, om->verts[in_a].no);
+
+		if (l->ang <= om->verts[in_a].link->ang) {
+			it = om->verts[in_a].link->prev;
+			om->verts[in_a].link = l;
+		}
+		else {
+			for (it = om->verts[in_a].link; it->next != om->verts[in_a].link; it = it->next) {
+				if (it->next->ang > l->ang)
+					break;
+			}
+		}
+
+		it->next->prev = l;
+		l->next = it->next;
+		l->prev = it;
+		it->next = l;
+	}
+	om->verts[in_a].num_links++;
+
+	return l;
+}
+
+QRDiskLink *linkVerts(OutputMesh *om, QRVertID in_v1, QRVertID in_v2)
+{
+	QRDiskLink *l1, *l2;
+
+	l1 = getLink(om, in_v1, in_v2);
+	if (l1)
+		return l1;
+
+	l1 = insertLink(om, in_v1, in_v2);
+	l2 = insertLink(om, in_v2, in_v1);
+
+	l1->brother = l2;
+	l2->brother = l1;
+
+	om->num_edges++;
+
+	return l1;
+}
+
+void deleteLink(OutputMesh *om, QRVert *ll, QRDiskLink *l)
+{
+	BLI_assert(ll->num_links != 0);
+
+	if (l->next == l) {
+		ll->link = NULL;
+	}
+	else {
+		if (ll->link == l) ll->link = l->next;
+		l->prev->next = l->next;
+		l->next->prev = l->prev;
+	}
+	ll->num_links--;
+}
+
+void unlinkVerts(OutputMesh *om, QRDiskLink *l)
+{
+	QRVert *ll1, *ll2;
+
+	ll2 = &om->verts[l->v];
+	ll1 = &om->verts[l->brother->v];
+
+	deleteLink(om, ll1, l);
+	deleteLink(om, ll2, l->brother);
+
+	om->num_edges--;
+}
+
+void dissolveVert(OutputMesh *om, QRVertID in_v)
+{
+	QRVertID a, b;
+	QRDiskLink *al, *bl;
+
+	a = om->verts[in_v].link->v;
+	al = om->verts[in_v].link->brother;
+	b = om->verts[in_v].link->next->v;
+	bl = om->verts[in_v].link->next->brother;
+
+	if (a != b) {
+		deleteLink(om, &om->verts[in_v], om->verts[in_v].link);
+		deleteLink(om, &om->verts[in_v], om->verts[in_v].link);
+
+		if (getLink(om, a, b)) {
+			deleteLink(om, &om->verts[a], al);
+			deleteLink(om, &om->verts[b], bl);
+
+			om->num_edges -= 2;
+		}
+		else {
+			al->v = b;
+			al->brother = bl;
+			bl->v = a;
+			bl->brother = al;
+			
+			om->num_edges--;
+		}
+	}
+	else {
+		deleteLink(om, &om->verts[in_v], om->verts[in_v].link);
+		deleteLink(om, &om->verts[in_v], om->verts[in_v].link);
+		deleteLink(om, &om->verts[a], al);
+		deleteLink(om, &om->verts[b], bl);
+
+		om->num_edges -= 2;
+	}
+}

intern/quadremesh/quadremesh_link.h

@@ -0,0 +1,31 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Alexander Pinzon Fernandez
+ *                 Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+#ifndef __QUADREMESH_LINK_H__
+#define __QUADREMESH_LINK_H__
+
+#include "BLI_utildefines.h"
+#include "quadremesh_util.h"
+
+#endif

intern/quadremesh/quadremesh_util.c

@@ -0,0 +1,277 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Alexander Pinzon Fernandez
+ *                 Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_modifier_types.h"
+
+#include "BLI_math.h"
+#include "BLI_rand.h"
+#include "BLI_heap.h"
+#include "BLI_linklist.h"
+
+#include "BKE_mesh_mapping.h"
+#include "BKE_DerivedMesh.h"
+
+#include "quadremesh_util.h"
+
+void getNormalAtEdge(float r_no[3], InputMesh *im, int in_e)
+{
+	if (im->faces_edge[in_e][0] == QR_NO_FACE)
+		copy_v3_v3(r_no, im->no[im->faces_edge[in_e][1]]);
+	else if (im->faces_edge[in_e][1] == QR_NO_FACE)
+		copy_v3_v3(r_no, im->no[im->faces_edge[in_e][0]]);
+	else {
+		add_v3_v3v3(r_no, im->no[im->faces_edge[in_e][0]],
+						  im->no[im->faces_edge[in_e][1]]);
+		mul_v3_fl(r_no, 0.5f);
+	}
+}
+
+/* ADDING STUFF TO OUTPUT MESH */
+
+QRVertID addVert(OutputMesh *om, float in_co[3], float in_no[3])
+{
+	if (om->num_verts == om->alloc_verts) {
+		om->alloc_verts = om->alloc_verts * 2 + 10;
+		om->verts = MEM_reallocN(om->verts, sizeof(QRVert) * om->alloc_verts);
+	}
+
+	copy_v3_v3(om->verts[om->num_verts].co, in_co);
+	if (in_no)
+		copy_v3_v3(om->verts[om->num_verts].no, in_no);
+
+	om->verts[om->num_verts].link = NULL;
+	om->verts[om->num_verts].num_links = 0;
+	
+	return om->num_verts++;
+}
+
+
+void addGFPoint(InputMesh *im, OutputMesh *om, GFPoint *in_p)
+{
+	int i;
+	float no[3];
+	QRVertID newv;
+
+	if (in_p->type == eVert) {
+		if (om->vonvs[in_p->v] == -1) {
+			zero_v3(no);
+			for (i = 0; i < im->ringf_map[in_p->v].count; i++)
+				add_v3_v3(no, im->no[im->ringf_map[in_p->v].indices[i]]);
+
+			mul_v3_fl(no, 1.0f / (float)im->ringf_map[in_p->v].count);
+
+			newv = addVert(om, in_p->co, no);
+			om->vonvs[in_p->v] = newv;
+		}
+		else
+			newv = om->vonvs[in_p->v];
+	}
+	else if (in_p->type == eEdge) {
+		getNormalAtEdge(no, im, in_p->e);
+		newv = addVert(om, in_p->co, no);
+	}
+	else { /* eFace */
+		copy_v3_v3(no, im->no[in_p->f]);
+		newv = addVert(om, in_p->co, no);
+	}
+
+	in_p->id = newv;
+}
+
+
+/*
+* List of vertices from original mesh with special features (edge dihedral angle less that 90) to be preserves
+* return the size of array
+*/
+static int *findFeaturesOnMesh(InputMesh *im, int size[2])
+{
+	int i, f1, f2, total, *listverts, *listdest;
+	float angle;
+
+	listverts = MEM_callocN(sizeof(int) * im->num_verts, __func__);
+	listdest = NULL;
+	total = 0;
+
+	for (i = 0; i < im->num_edges; i++) {
+		f1 = im->faces_edge[i][0];
+		f2 = im->faces_edge[i][1];
+		angle = angle_normalized_v3v3(im->no[f1], im->no[f2]);
+		if (angle >= M_PI_2) {
+			listverts[im->edges[i][0]] = 1;
+			listverts[im->edges[i][1]] = 1;
+			break;
+		}
+	}
+
+	for (i = 0; i < im->num_verts; i++) {
+		if (im->constraints[i] == 1) {
+			listverts[i] = 1;
+		}
+	}
+
+	for (i = 0; i < im->num_verts; i++) {
+		if (listverts[i] == 1) {
+			total++;
+		}
+	}
+	if (total > 0) {
+		listdest = MEM_mallocN(sizeof(int) * total, __func__);
+	}
+	total = 0;
+	for (i = 0; i < im->num_verts; i++) {
+		if (listverts[i] == 1) {
+			listdest[total++] = i;
+		}
+	}
+	MEM_SAFE_FREE(listverts);
+	size[0] = total;
+
+	return listdest;
+}
+
+void addSeedToQueue(Heap *aheap, float in_co[3], GFPointType in_type, int in_val, float weight)
+{
+	GFPoint *seed;
+
+	seed = MEM_mallocN(sizeof(GFPoint), __func__);
+
+	copy_v3_v3(seed->co, in_co);
+	seed->type = in_type;
+	seed->e = seed->f = seed->v = seed->id = -1;
+	if (in_type == eVert)
+		seed->v = in_val;
+	else
+		seed->f = in_val;
+
+	BLI_heap_insert(aheap, weight, seed);
+}
+
+void getInitialSeeds(GradientFlowSystem *gfsys)
+{
+	int i, *lverts, sizeverts[2];
+
+	lverts = findFeaturesOnMesh(&gfsys->sys->input_mesh, sizeverts);
+	
+	for (i = 0; i < sizeverts[0]; i++)
+		addSeedToQueue(gfsys->seeds, gfsys->sys->input_mesh.co[lverts[i]], eVert, lverts[i], 0.0f);
+	
+	MEM_SAFE_FREE(lverts);
+}
+
+/* SEED QUEUE */
+
+GFPoint *getTopSeedFromQueue(struct Heap *aheap)
+{
+	return (GFPoint*)BLI_heap_popmin(aheap);
+}
+
+bool isectSegmentWithOthersOnFace(OutputMesh *om, GFSysID sys_id, float in_a[3], float in_b[3], int in_f)
+{
+	QREdge *e;
+	LinkNode *iter;
+
+	for (iter = om->ringf[in_f]; iter; iter = iter->next) {
+		e = (QREdge*)iter->link;
+		
+		if ((e->v1->gfsysid & sys_id) == 0)
+			continue;
+
+		if (isect_seg_seg_unsafe_v3(NULL, in_a, in_b, om->verts[e->v1->v].co, om->verts[e->v2->v].co))
+			return true;
+	}
+
+	return false;
+}
+
+bool isectPointWithQREdge(OutputMesh *om, GFSysID sys_id, float in_co[3], int in_e)
+{
+	float vec[3], tmp;
+	QREdge *e;
+	QREdgeLink *it;
+
+	e = &om->ringe[in_e];
+
+	if (!e || e->num_links < 2)
+		return false;
+
+	sub_v3_v3v3(vec, in_co, e->orig);
+	tmp = dot_v3v3(vec, e->dir);
+
+	if (tmp <= e->v2->dist && tmp >= e->v1->dist) {
+		for (it = e->v1; it->next && it->next->dist < tmp; it = it->next);
+		
+		if (it->elink != NULL && (it->gfsysid & sys_id) != 0)
+			return true;
+	}
+
+	return false;
+}
+
+void freeGradientFlowSystem(GradientFlowSystem *gfsys)
+{
+	BLI_heap_free(gfsys->seeds, MEM_freeN);
+	MEM_SAFE_FREE(gfsys->gf);
+	MEM_SAFE_FREE(gfsys);
+}
+
+void freeQuadRemeshSystem(QuadRemeshSystem *sys)
+{
+	if (sys->input_mesh.is_alloc) {
+		freeInputMesh(&sys->input_mesh);
+	}
+
+	if (sys->is_alloc) {
+		MEM_SAFE_FREE(sys->U_field);
+		freeGradientFlowSystem(sys->gfsys[0]);
+		freeGradientFlowSystem(sys->gfsys[1]);
+	}
+
+	if (sys->cache_mesh) {
+		sys->cache_mesh->release(sys->cache_mesh);
+	}
+	
+	MEM_SAFE_FREE(sys);
+}
+
+GradientFlowSystem *newGradientFlowSystem(QuadRemeshSystem *sys)
+{
+	GradientFlowSystem *gfsys = MEM_callocN(sizeof(GradientFlowSystem), "GradientFlowSystem");
+	
+	gfsys->sys = sys;
+	gfsys->seeds = BLI_heap_new();
+	
+	return gfsys;
+}
+
+void initQuadRemeshSystem(QuadRemeshModifierData *qmd)
+{
+	QuadRemeshSystem *sys;
+
+	sys = MEM_callocN(sizeof(QuadRemeshSystem), "QuadRemeshSystem - Cache");
+
+	sys->qmd = qmd;
+	qmd->cache_system = sys;
+}

intern/quadremesh/quadremesh_util.h

@@ -0,0 +1,207 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Alexander Pinzon Fernandez
+ *                 Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+#ifndef __QUADREMESH_UTIL_H__
+#define __QUADREMESH_UTIL_H__
+
+/* Probability of producing a seed while tracing
+ * a line during 'sampling' - see QR_SAMPLING_RATE.
+ *  - 1.0f meaning 100%
+ * Do not reccomend setting this much higher, every line
+ * starts producing too many seeds and the amount
+ * of active seeds quickly grows. */
+#define QR_SEEDPROB 0.2f
+
+/* Distance, after which the line will be 'probed'.
+ * Probing: lines are traced in two opposite directions,
+ * perpendicular to the line being traced. If any of them hit
+ * another line before being 'maxdist' long, the tracing stops and the line is cut.
+ * If they don't, and RNG gods say to make a seed, lines are
+ * continued and a seed is placed at QR_MAXDIST_TO_SEEDDIST * maxdist
+ * distance from mother line. */
+#define QR_SAMPLING_RATE 0.04f
+#define QR_MAXDIST_TO_SEEDDIST 2.0f
+
+#define QR_NO_FACE 0xffffffff
+
+typedef int MEdgeID;
+typedef int QRVertID;
+
+typedef enum {
+	GFSYS1 = 1 << 0,
+	GFSYS2 = 1 << 1,
+	GFSYSNONE = 1 << 2
+} GFSysID;
+
+typedef struct QRDiskLink {
+	struct QRDiskLink *next, *prev, *brother;
+	MEdgeID e;
+	QRVertID v;
+	float ang;
+	bool poly_on_right;
+} QRDiskLink;
+
+typedef struct QREdgeLink {
+	QRVertID v;
+	float dist;
+
+	struct QREdgeLink *next;
+	struct QRDiskLink *elink;
+
+	GFSysID gfsysid;
+} QREdgeLink;
+
+typedef struct QREdge {
+	QREdgeLink *v1, *v2;
+	float dir[3], orig[3];
+	int num_links;
+} QREdge;
+
+typedef struct QRVert {
+	struct QRDiskLink *link;
+	float co[3], no[3], vec[3];
+	int num_links;
+} QRVert;
+
+typedef enum {
+	eVert,
+	eEdge,
+	eFace
+} GFPointType;
+
+typedef struct GFPoint {
+	float co[3];
+	union {
+		int v, e, f;
+	};
+	QRVertID id;
+	GFPointType type;
+} GFPoint;
+
+typedef enum {
+	eEdgeSegment,
+	eFaceSegment,
+	eNoSegment
+} GFSegmentType;
+
+typedef struct GFSegment {
+	GFPoint p;
+	union {
+		int e, f;
+	};
+	GFSegmentType type;
+} GFSegment;
+
+/* GradientFlowSysten, one gfsys for every gradient field */
+typedef struct GradientFlowSystem {
+	GFSysID id;
+	struct Heap *seeds;
+	//float *h;
+	float(*gf)[3];				/* Gradient Field */
+
+	struct QuadRemeshSystem *sys;
+} GradientFlowSystem;
+
+typedef struct InputMesh {
+	bool is_alloc;
+
+	int num_verts, num_edges, num_faces, num_features;
+	float(*co)[3];					/* Original vertex coordinates */
+	float(*no)[3];					/* Original face normal */
+	float(*vno)[3];
+	float(*ev)[3]; /* Normalized edge vectors */
+
+	unsigned int(*faces)[3];		/* Copy of MFace (tessface) v1-v3, v2-v4 */
+	unsigned int(*edges)[2];		/* Copy of edges v1-v2 */
+	unsigned int(*faces_edge)[2];	/* Faces by edges  */
+
+	int *ringf_indices;				/* Indices of faces per vertex */
+	int *ringe_indices;				/* Indices of edges per vertex */
+	int *ringv_indices;				/* Indices of verts per vertex */
+	struct MeshElemMap *ringf_map;  /* Map of faces per vertex */
+	struct MeshElemMap *ringe_map;	/* Map of edges per vertex */
+	struct MeshElemMap *ringv_map;	/* Map of verts per vertex */
+
+	int *constraints;				/* Feature points constraints*/
+	float *weights;					/* Feature points weights*/
+} InputMesh;
+
+typedef struct OutputMesh {
+	struct MemArena *memarena;
+	QRVert *verts;
+	QRVertID *vonvs;
+	struct QREdge *ringe;                  /* QREdges per original edges */
+	struct LinkNode **ringf;        /* Lists of QREdge per original faces */
+	
+	int num_verts, alloc_verts;
+	int num_edges;
+} OutputMesh;
+
+typedef struct QuadRemeshSystem {
+	struct QuadRemeshModifierData *qmd;
+
+	InputMesh input_mesh;
+	
+	bool has_solution, is_alloc;
+	GradientFlowSystem *gfsys[2];
+	float *U_field;					/* Initial scalar field*/
+
+	OutputMesh output_mesh;
+	struct RNG *rng;
+
+	struct DerivedMesh *cache_mesh;
+} QuadRemeshSystem;
+
+void getNormalAtEdge(float r_no[3], InputMesh *im, int in_e);
+
+QRDiskLink *linkVerts(OutputMesh *om, QRVertID in_v1, QRVertID in_v2);
+void deleteLink(OutputMesh *om, QRVert *ll, QRDiskLink *l);
+void unlinkVerts(OutputMesh *om, QRDiskLink *l);
+void dissolveVert(OutputMesh *om, QRVertID in_v);
+
+void insertOnQREdge(OutputMesh *om, QREdge *in_e, QRVertID in_vid);
+void linkOnQREdge(OutputMesh *om, GFSysID sys_id, QREdge *in_e, QRVertID in_v1, QRVertID in_v2);
+QREdge *addQREdgeToFace(OutputMesh *om, InputMesh *im, GFSysID sys_id, int in_f, QRVertID in_v1, QRVertID in_v2);
+
+QRVertID addVert(OutputMesh *om, float in_co[3], float in_no[3]);
+void addGFPoint(InputMesh *im, OutputMesh *om, GFPoint *in_p);
+
+void addSeedToQueue(struct Heap *aheap, float in_co[3], GFPointType in_type, int in_val, float weight);
+GFPoint *getTopSeedFromQueue(struct Heap *aheap);
+
+bool isectSegmentWithOthersOnFace(OutputMesh *om, GFSysID sys_id, float in_a[3], float in_b[3], int in_f);
+bool isectPointWithQREdge(OutputMesh *om, GFSysID sys_id, float in_co[3], int in_e);
+
+void getInitialSeeds(GradientFlowSystem *gfsys);
+void computeFlowLines(QuadRemeshSystem *sys);
+
+void freeInputMesh(InputMesh *im);
+void getInput(QuadRemeshSystem *sys, struct Object *ob, struct DerivedMesh *dm);
+void getHarmonicGradients(QuadRemeshSystem *sys);
+
+GradientFlowSystem *newGradientFlowSystem(QuadRemeshSystem *sys);
+void freeGradientFlowSystem(GradientFlowSystem *gfsys);
+//void getPrincipalCurvatures(QuadRemeshSystem *sys);
+
+#endif

release/scripts/startup/bl_ui/properties_data_modifier.py

@@ -679,6 +679,32 @@ class DATA_PT_modifiers(ModifierButtonsPanel, Panel):
     def PARTICLE_SYSTEM(self, layout, ob, md):
         layout.label(text="Settings can be found inside the Particle context")
 
+    def QUADREMESH(self, layout, ob, md):
+        row = layout.row()
+        row.label(text="Features Vertex Group:")
+
+        row = layout.row()
+        row.prop_search(md, "vertex_group", ob, "vertex_groups", text="")
+        
+        layout.separator()
+        row = layout.row()
+        layout.prop(md, "max_line_dist")
+
+        row = layout.row()
+        layout.prop(md, "seed")
+        
+        layout.separator()
+        row = layout.row()
+        layout.prop(md, "is_autoupdate")
+
+        row = layout.row()
+        row.enabled = not bool(md.is_autoupdate)
+        row.operator("object.quadremesh_computeflow", text="Update all")
+
+        row = layout.row()
+        row.enabled = not bool(md.is_autoupdate)
+        row.operator("object.quadremesh_remesh", text="Update remesh")
+
     def SCREW(self, layout, ob, md):
         split = layout.split()
 

source/blender/blenlib/BLI_math_geom.h

@@ -153,6 +153,9 @@ int isect_line_line_v3(
         const float v1[3], const float v2[3],
         const float v3[3], const float v4[3],
         float i1[3], float i2[3]);
+bool isect_seg_seg_unsafe_v3(float r[3],
+                             const float v1[3], const float v2[3],
+                             const float v3[3], const float v4[3]);
 bool isect_line_line_strict_v3(const float v1[3], const float v2[3],
                                const float v3[3], const float v4[3],
                                float vi[3], float *r_lambda);

source/blender/blenlib/BLI_math_vector.h

@@ -265,6 +265,7 @@ float cos_v3v3v3(const float p1[3], const float p2[3], const float p3[3]) ATTR_W
 float angle_normalized_v3v3(const float v1[3], const float v2[3]) ATTR_WARN_UNUSED_RESULT;
 float angle_on_axis_v3v3v3_v3(const float v1[3], const float v2[3], const float v3[3], const float axis[3]) ATTR_WARN_UNUSED_RESULT;
 float angle_signed_on_axis_v3v3v3_v3(const float v1[3], const float v2[3], const float v3[3], const float axis[3]) ATTR_WARN_UNUSED_RESULT;
+float angle_signed_on_axis_normalized_fast_v3v3_v3(const float v1[3], const float v2[3], const float axis[3]) ATTR_WARN_UNUSED_RESULT;
 void angle_tri_v3(float angles[3], const float v1[3], const float v2[3], const float v3[3]);
 void angle_quad_v3(float angles[4], const float v1[3], const float v2[3], const float v3[3], const float v4[3]);
 void angle_poly_v3(float *angles, const float *verts[3], int len);

source/blender/blenlib/intern/math_geom.c

@@ -1760,6 +1760,49 @@ int isect_line_line_v3(
 	return isect_line_line_epsilon_v3(v1, v2, v3, v4, i1, i2, epsilon);
 }
 
+bool isect_seg_seg_unsafe_v3(float r[3],
+                             const float v1[3], const float v2[3],
+                             const float v3[3], const float v4[3])
+{
+	float a[3], b[3], c[3], ab[3], cb[3], ca[3];
+	float d, div;
+
+	sub_v3_v3v3(c, v3, v1);
+	sub_v3_v3v3(a, v2, v1);
+	sub_v3_v3v3(b, v4, v3);
+
+	cross_v3_v3v3(ab, a, b);
+	d = dot_v3v3(c, ab);
+	div = dot_v3v3(ab, ab);
+
+	/* test zero length line */
+	if (UNLIKELY(div == 0.0f)) {
+		return false;
+	}
+	/* test if the two lines are coplanar */
+	else if (d > -0.000001f && d < 0.000001f) {
+		float f1, f2;
+		cross_v3_v3v3(cb, c, b);
+		cross_v3_v3v3(ca, c, a);
+
+		f1 = dot_v3v3(cb, ab) / div;
+		f2 = dot_v3v3(ca, ab) / div;
+
+		if (f1 >= 0.0f && f1 <= 1.0f &&
+		    f2 >= 0.0f && f2 <= 1.0f)
+		{
+			if (r) {
+				mul_v3_fl(a, f1);
+				add_v3_v3v3(r, v1, a);
+			}
+
+			return true; /* intersection found */
+		}
+	}
+	
+	return false;
+}
+
 /** Intersection point strictly between the two lines
  * \return false when no intersection is found
  */

source/blender/blenlib/intern/math_vector.c

@@ -499,6 +499,23 @@ float angle_signed_on_axis_v3v3v3_v3(const float v1[3], const float v2[3], const
 	return angle;
 }
 
+float angle_signed_on_axis_normalized_fast_v3v3_v3(const float v1[3], const float v2[3], const float axis[3])
+{
+	float tproj[3];
+	float angle, dot;
+
+	dot = dot_v3v3(v1, v2);
+	CLAMP(dot, -1.0f, 1.0f);
+	angle = acosf(dot);
+
+	cross_v3_v3v3(tproj, v2, v1);
+	if (dot_v3v3(tproj, axis) < 0.0f) {
+		angle = ((float)(M_PI * 2.0f)) - angle;
+	}
+
+	return angle;
+}
+
 void angle_tri_v3(float angles[3], const float v1[3], const float v2[3], const float v3[3])
 {
 	float ed1[3], ed2[3], ed3[3];

source/blender/blenloader/intern/readfile.c

@@ -4986,6 +4986,11 @@ static void direct_link_modifiers(FileData *fd, ListBase *lb)
 			csmd->delta_cache = NULL;
 			csmd->delta_cache_num = 0;
 		}
+		else if (md->type == eModifierType_QuadRemesh) {
+			QuadRemeshModifierData *qmd = (QuadRemeshModifierData*)md;
+
+			qmd->cache_system = NULL;
+		}
 	}
 }
 

source/blender/editors/object/object_intern.h

@@ -184,6 +184,8 @@ void OBJECT_OT_skin_loose_mark_clear(struct wmOperatorType *ot);
 void OBJECT_OT_skin_radii_equalize(struct wmOperatorType *ot);
 void OBJECT_OT_skin_armature_create(struct wmOperatorType *ot);
 void OBJECT_OT_laplaciandeform_bind(struct wmOperatorType *ot);
+void OBJECT_OT_quadremesh_computeflow(struct wmOperatorType *ot);
+void OBJECT_OT_quadremesh_remesh(struct wmOperatorType *ot);
 
 /* object_constraint.c */
 void OBJECT_OT_constraint_add(struct wmOperatorType *ot);

source/blender/editors/object/object_modifier.c

@@ -2295,3 +2295,95 @@ void OBJECT_OT_laplaciandeform_bind(wmOperatorType *ot)
 	ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO | OPTYPE_INTERNAL;
 	edit_modifier_properties(ot);
 }
+
+
+/************************ QuadRemesh compute flow operator *********************/
+
+static int quadremesh_poll(bContext *C)
+{
+	return edit_modifier_poll_generic(C, &RNA_QuadRemeshModifier, 0);
+}
+
+static int quadremesh_computeflow_exec(bContext *C, wmOperator *op)
+{
+	Object *ob = ED_object_active_context(C);
+	QuadRemeshModifierData *lmd = (QuadRemeshModifierData *)edit_modifier_property_get(op, ob, eModifierType_QuadRemesh);
+
+	if (!lmd)
+		return OPERATOR_CANCELLED;
+
+	lmd->flag |= MOD_QUADREMESH_INPUT_DIRTY;
+
+	DAG_id_tag_update(&ob->id, OB_RECALC_DATA);
+	WM_event_add_notifier(C, NC_OBJECT | ND_MODIFIER, ob);
+
+	return OPERATOR_FINISHED;
+}
+
+static int quadremesh_computeflow_invoke(bContext *C, wmOperator *op, const wmEvent *UNUSED(event))
+{
+	if (edit_modifier_invoke_properties(C, op))
+		return quadremesh_computeflow_exec(C, op);
+	else
+		return OPERATOR_CANCELLED;
+}
+
+void OBJECT_OT_quadremesh_computeflow(wmOperatorType *ot)
+{
+	/* identifiers */
+	ot->name = "Quad Remesh Compute Flow";
+	ot->description = "Compute gradient flow in Quad Remesh modifier";
+	ot->idname = "OBJECT_OT_quadremesh_computeflow";
+
+	/* api callbacks */
+	ot->poll = quadremesh_poll;
+	ot->invoke = quadremesh_computeflow_invoke;
+	ot->exec = quadremesh_computeflow_exec;
+
+	/* flags */
+	ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO | OPTYPE_INTERNAL;
+	edit_modifier_properties(ot);
+}
+
+/************************ QuadRemesh remesh operator *********************/
+
+static int quadremesh_remesh_exec(bContext *C, wmOperator *op)
+{
+	Object *ob = ED_object_active_context(C);
+	QuadRemeshModifierData *lmd = (QuadRemeshModifierData *)edit_modifier_property_get(op, ob, eModifierType_QuadRemesh);
+
+	if (!lmd)
+		return OPERATOR_CANCELLED;
+
+	lmd->flag |= MOD_QUADREMESH_REMESH;
+
+	DAG_id_tag_update(&ob->id, OB_RECALC_DATA);
+	WM_event_add_notifier(C, NC_OBJECT | ND_MODIFIER, ob);
+
+	return OPERATOR_FINISHED;
+}
+
+static int quadremesh_remesh_invoke(bContext *C, wmOperator *op, const wmEvent *UNUSED(event))
+{
+	if (edit_modifier_invoke_properties(C, op))
+		return quadremesh_remesh_exec(C, op);
+	else
+		return OPERATOR_CANCELLED;
+}
+
+void OBJECT_OT_quadremesh_remesh(wmOperatorType *ot)
+{
+	/* identifiers */
+	ot->name = "Quad Remesh Remeshing";
+	ot->description = "Compute Remeshing modifier";
+	ot->idname = "OBJECT_OT_quadremesh_remesh";
+
+	/* api callbacks */
+	ot->poll = quadremesh_poll;
+	ot->invoke = quadremesh_remesh_invoke;
+	ot->exec = quadremesh_remesh_exec;
+
+	/* flags */
+	ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO | OPTYPE_INTERNAL;
+	edit_modifier_properties(ot);
+}

source/blender/editors/object/object_ops.c

@@ -250,6 +250,8 @@ void ED_operatortypes_object(void)
 	WM_operatortype_append(OBJECT_OT_lod_remove);
 
 	WM_operatortype_append(OBJECT_OT_vertex_random);
+	WM_operatortype_append(OBJECT_OT_quadremesh_computeflow);
+	WM_operatortype_append(OBJECT_OT_quadremesh_remesh);
 
 	WM_operatortype_append(OBJECT_OT_data_transfer);
 	WM_operatortype_append(OBJECT_OT_datalayout_transfer);

source/blender/makesdna/DNA_modifier_types.h

@@ -85,6 +85,7 @@ typedef enum ModifierType {
 	eModifierType_DataTransfer      = 49,
 	eModifierType_NormalEdit        = 50,
 	eModifierType_CorrectiveSmooth  = 51,
+	eModifierType_QuadRemesh        = 52,
 	NUM_MODIFIER_TYPES
 } ModifierType;
 
@@ -1444,6 +1445,26 @@ enum {
 	MOD_WIREFRAME_CREASE        = (1 << 5),
 };
 
+typedef struct QuadRemeshModifierData {
+	ModifierData modifier;
+	char anchor_grp_name[64];  /* MAX_VGROUP_NAME */
+
+	short flag;
+	char auto_updates;
+	char pad[5];
+
+	float max_line_dist;
+	int rng_seed;
+
+	void *cache_system;  /* runtime only */
+} QuadRemeshModifierData;
+
+/* QuadRemesh modifier flags */
+enum {
+	MOD_QUADREMESH_REMESH = (1 << 1),
+	MOD_QUADREMESH_INPUT_DIRTY = (1 << 2),
+	MOD_QUADREMESH_FIELD_DIRTY = (1 << 3),
+};
 
 typedef struct DataTransferModifierData {
 	ModifierData modifier;

source/blender/makesrna/RNA_access.h

@@ -466,6 +466,7 @@ extern StructRNA RNA_PropertyGroupItem;
 extern StructRNA RNA_PropertySensor;
 extern StructRNA RNA_PythonConstraint;
 extern StructRNA RNA_PythonController;
+extern StructRNA RNA_QuadRemeshModifier;
 extern StructRNA RNA_QuickTimeSettings;
 extern StructRNA RNA_RadarSensor;
 extern StructRNA RNA_RandomSensor;

source/blender/makesrna/intern/rna_modifier.c

@@ -83,6 +83,7 @@ EnumPropertyItem modifier_type_items[] = {
 	{eModifierType_Mirror, "MIRROR", ICON_MOD_MIRROR, "Mirror", ""},
 	{eModifierType_Multires, "MULTIRES", ICON_MOD_MULTIRES, "Multiresolution", ""},
 	{eModifierType_Remesh, "REMESH", ICON_MOD_REMESH, "Remesh", ""},
+	{eModifierType_QuadRemesh, "QUADREMESH", ICON_MOD_REMESH, "Quadrilateral Remesh", "" },
 	{eModifierType_Screw, "SCREW", ICON_MOD_SCREW, "Screw", ""},
 	{eModifierType_Skin, "SKIN", ICON_MOD_SKIN, "Skin", ""},
 	{eModifierType_Solidify, "SOLIDIFY", ICON_MOD_SOLIDIFY, "Solidify", ""},
@@ -374,6 +375,8 @@ static StructRNA *rna_Modifier_refine(struct PointerRNA *ptr)
 			return &RNA_LaplacianDeformModifier;
 		case eModifierType_Wireframe:
 			return &RNA_WireframeModifier;
+		case eModifierType_QuadRemesh:
+			return &RNA_QuadRemeshModifier;
 		case eModifierType_DataTransfer:
 			return &RNA_DataTransferModifier;
 		case eModifierType_NormalEdit:
@@ -471,6 +474,7 @@ RNA_MOD_VGROUP_NAME_SET(WeightVGMix, mask_defgrp_name);
 RNA_MOD_VGROUP_NAME_SET(WeightVGProximity, defgrp_name);
 RNA_MOD_VGROUP_NAME_SET(WeightVGProximity, mask_defgrp_name);
 RNA_MOD_VGROUP_NAME_SET(Wireframe, defgrp_name);
+RNA_MOD_VGROUP_NAME_SET(QuadRemesh, anchor_grp_name);
 
 static void rna_ExplodeModifier_vgroup_get(PointerRNA *ptr, char *value)
 {
@@ -745,6 +749,12 @@ static int rna_LaplacianDeformModifier_is_bind_get(PointerRNA *ptr)
 	return ((lmd->flag & MOD_LAPLACIANDEFORM_BIND) && (lmd->cache_system != NULL));
 }
 
+static int rna_QuadRemeshModifier_is_autoupdate_get(PointerRNA *ptr)
+{
+	QuadRemeshModifierData *qmd = (QuadRemeshModifierData *)ptr->data;
+	return qmd->auto_updates;
+}
+
 /* NOTE: Curve and array modifiers requires curve path to be evaluated,
  * dependency graph will make sure that curve eval would create such a path,
  * but if curve was already evaluated we might miss path.
@@ -4236,7 +4246,45 @@ static void rna_def_modifier_wireframe(BlenderRNA *brna)
 	RNA_def_property_update(prop, 0, "rna_Modifier_update");
 }
 
+
+static void rna_def_modifier_quadremesh(BlenderRNA *brna)
+{
+	StructRNA *srna;
+	PropertyRNA *prop;
+
+	srna = RNA_def_struct(brna, "QuadRemeshModifier", "Modifier");
+	RNA_def_struct_ui_text(srna, "Quadrilateral Remesh Modifier", "Quadrilateral Remesh modifier");
+	RNA_def_struct_sdna(srna, "QuadRemeshModifierData");
+	RNA_def_struct_ui_icon(srna, ICON_MOD_REMESH);
+
+	prop = RNA_def_property(srna, "vertex_group", PROP_STRING, PROP_NONE);
+	RNA_def_property_string_sdna(prop, NULL, "anchor_grp_name");
+	RNA_def_property_ui_text(prop, "Vertex group for feature points",
+		"Name of Vertex Group which determines feature points");
+	RNA_def_property_string_funcs(prop, NULL, NULL, "rna_QuadRemeshModifier_anchor_grp_name_set");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+
+	prop = RNA_def_property(srna, "max_line_dist", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "max_line_dist");
+	RNA_def_property_range(prop, 0.001f, 10.0f);
+	RNA_def_property_ui_range(prop, 0.001f, 10.0f, 0.1, 1);
+	RNA_def_property_ui_text(prop, "Line Spacing", "Remeshing Line Spacing");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+
+	prop = RNA_def_property(srna, "seed", PROP_INT, PROP_NONE);
+	RNA_def_property_int_sdna(prop, NULL, "rng_seed");
+	RNA_def_property_range(prop, 0.0f, 1000.0f);
+	RNA_def_property_ui_text(prop, "Seed", "Seed used for remeshing");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+
+	prop = RNA_def_property(srna, "is_autoupdate", PROP_BOOLEAN, PROP_NONE);
+	RNA_def_property_boolean_sdna(prop, NULL, "auto_updates", 0);
+	RNA_def_property_ui_text(prop, "Automatic update", "Update all data automatically");
+	RNA_def_property_update(prop, 0, "rna_Modifier_update");
+}
+
 static void rna_def_modifier_datatransfer(BlenderRNA *brna)
+
 {
 	StructRNA *srna;
 	PropertyRNA *prop;
@@ -4667,6 +4715,7 @@ void RNA_def_modifier(BlenderRNA *brna)
 	rna_def_modifier_meshcache(brna);
 	rna_def_modifier_laplaciandeform(brna);
 	rna_def_modifier_wireframe(brna);
+	rna_def_modifier_quadremesh(brna);
 	rna_def_modifier_datatransfer(brna);
 	rna_def_modifier_normaledit(brna);
 }

source/blender/modifiers/CMakeLists.txt

@@ -79,6 +79,7 @@ set(SRC
 	intern/MOD_ocean.c
 	intern/MOD_particleinstance.c
 	intern/MOD_particlesystem.c
+	intern/MOD_quadremesh.c
 	intern/MOD_remesh.c
 	intern/MOD_screw.c
 	intern/MOD_shapekey.c
@@ -151,4 +152,11 @@ if(WITH_OPENNL)
 	)
 endif()
 
+if(WITH_MOD_QUADREMESH)
+	add_definitions(-DWITH_MOD_QUADREMESH)
+	list(APPEND INC
+		../../../intern/quadremesh
+	)
+endif()
+
 blender_add_lib(bf_modifiers "${SRC}" "${INC}" "${INC_SYS}")

source/blender/modifiers/MOD_modifiertypes.h

@@ -81,6 +81,7 @@ extern ModifierTypeInfo modifierType_UVWarp;
 extern ModifierTypeInfo modifierType_MeshCache;
 extern ModifierTypeInfo modifierType_LaplacianDeform;
 extern ModifierTypeInfo modifierType_Wireframe;
+extern ModifierTypeInfo modifierType_QuadRemesh;
 extern ModifierTypeInfo modifierType_DataTransfer;
 extern ModifierTypeInfo modifierType_NormalEdit;
 extern ModifierTypeInfo modifierType_CorrectiveSmooth;

source/blender/modifiers/intern/MOD_quadremesh.c

@@ -0,0 +1,156 @@
+/*
+ * ***** 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.
+ *
+ * Contributor(s): Alexander Pinzon Fernandez
+ *                 Krzysztof Recko
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ *
+ */
+
+/** \file blender/modifiers/intern/MOD_quadremesh.c
+ *  \ingroup modifiers
+ */
+
+#include "BLI_math.h"
+#include "BLI_utildefines.h"
+#include "BLI_string.h"
+#include "BLI_rand.h"
+
+#include "MEM_guardedalloc.h"
+
+#include "BKE_cdderivedmesh.h"
+#include "BKE_particle.h"
+#include "BKE_deform.h"
+
+#include "MOD_util.h"
+#include "QRM_quadremesh.h"
+
+static void initData(ModifierData *md)
+{
+	QuadRemeshModifierData *lmd = (QuadRemeshModifierData *)md;
+	lmd->anchor_grp_name[0] = '\0';
+	lmd->flag = 0;
+
+	lmd->cache_system = NULL;
+	lmd->max_line_dist = QR_MINDIST;
+	lmd->rng_seed = QR_SEED;
+	lmd->auto_updates = QR_AUTOUPDATES;
+}
+
+static void copyData(ModifierData *md, ModifierData *target)
+{
+	QuadRemeshModifierData *qmd = (QuadRemeshModifierData *)md;
+	
+	modifier_copyData_generic(md, target);
+	qmd->cache_system = NULL;
+}
+
+static bool isDisabled(ModifierData *md, int UNUSED(useRenderParams))
+{
+	QuadRemeshModifierData *lmd = (QuadRemeshModifierData *)md;
+
+	if (!lmd->anchor_grp_name[0])
+		return true;
+
+	return false;
+}
+
+static CustomDataMask requiredDataMask(Object *UNUSED(ob), ModifierData *md)
+{
+	QuadRemeshModifierData *lmd = (QuadRemeshModifierData *)md;
+	CustomDataMask dataMask = 0;
+	if (lmd->anchor_grp_name[0]) dataMask |= CD_MASK_MDEFORMVERT;
+	return dataMask;
+}
+
+#define QR_SHOW_UFIELD
+
+static DerivedMesh *applyModifier(ModifierData *md,
+	Object *ob,
+	DerivedMesh *dm,
+	ModifierApplyFlag UNUSED(flag))
+{
+	DerivedMesh *result = dm;
+	QuadRemeshModifierData *qmd = (QuadRemeshModifierData *)md;
+	QuadRemeshSystem *sys;
+
+	if (qmd->auto_updates) {
+		qmd->flag |= MOD_QUADREMESH_INPUT_DIRTY;
+	}
+
+	if (!qmd->cache_system) {
+		initQuadRemeshSystem(qmd);
+	}
+
+	sys = qmd->cache_system;
+
+	result = makeResultMesh(sys, ob, dm);
+
+#ifdef QR_SHOW_UFIELD
+	{
+		int defgrp_index;
+		MDeformVert *dvert = NULL;
+
+		if (!defgroup_find_name(ob, "QuadRemeshFlow")) {
+			BKE_defgroup_new(ob, "QuadRemeshFlow");
+		}
+
+		modifier_get_vgroup(ob, dm, "QuadRemeshFlow", &dvert, &defgrp_index);
+		getUField(sys, dvert, defgrp_index);
+	}
+#endif
+
+	return result;
+}
+
+static void freeData(ModifierData *md)
+{
+	QuadRemeshModifierData *qmd = (QuadRemeshModifierData *)md;
+	QuadRemeshSystem *sys = (QuadRemeshSystem *)qmd->cache_system;
+
+	if (sys) {
+		freeQuadRemeshSystem(sys);
+	}
+}
+
+ModifierTypeInfo modifierType_QuadRemesh = {
+	/* name */              "QuadRemesh",
+	/* structName */        "QuadRemeshModifierData",
+	/* structSize */        sizeof(QuadRemeshModifierData),
+	/* type */              eModifierTypeType_Nonconstructive,
+	/* flags */             eModifierTypeFlag_AcceptsMesh |
+							eModifierTypeFlag_AcceptsCVs,
+
+	/* copyData */          copyData,
+	/* deformVerts */       NULL,
+	/* deformMatrices */    NULL,
+	/* deformVertsEM */     NULL,
+	/* deformMatricesEM */  NULL,
+	/* applyModifier */     applyModifier,
+	/* applyModifierEM */   NULL,
+	/* initData */          initData,
+	/* requiredDataMask */  requiredDataMask,
+	/* freeData */          freeData,
+	/* isDisabled */        isDisabled,
+	/* updateDepgraph */    NULL,
+	/* dependsOnTime */     NULL,
+	/* dependsOnNormals */	NULL,
+	/* foreachObjectLink */ NULL,
+	/* foreachIDLink */     NULL,
+	/* foreachTexLink */    NULL,
+};

source/blender/modifiers/intern/MOD_util.c

@@ -304,6 +304,7 @@ void modifier_type_init(ModifierTypeInfo *types[])
 	INIT_TYPE(MeshCache);
 	INIT_TYPE(LaplacianDeform);
 	INIT_TYPE(Wireframe);
+	INIT_TYPE(QuadRemesh);
 	INIT_TYPE(DataTransfer);
 	INIT_TYPE(NormalEdit);
 	INIT_TYPE(CorrectiveSmooth);