copying bmesh dir on its own from bmesh branch

1 files changed, 1069 insertions, 0 deletions

diff --git a/source/blender/bmesh/intern/bmesh_polygon.c b/source/blender/bmesh/intern/bmesh_polygon.c
new file mode 100644
index 00000000000..07945466a0c
--- /dev/null
+++ b/source/blender/bmesh/intern/bmesh_polygon.c
@@ -0,0 +1,1069 @@
+/*
+ * ***** 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): Joseph Eagar, Geoffrey Bantle, Campbell Barton
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+/** \file blender/bmesh/intern/bmesh_polygon.c
+ *  \ingroup bmesh
+ *
+ * This file contains code for dealing
+ * with polygons (normal/area calculation,
+ * tesselation, etc)
+ *
+ * BMESH_TODO:
+ *  - Add in Tesselator frontend that creates
+ *    BMTriangles from copied faces
+ *
+ *  - Add in Function that checks for and flags
+ *    degenerate faces.
+ */
+
+#include "BLI_math.h"
+#include "BLI_array.h"
+
+#include "MEM_guardedalloc.h"
+
+#include "bmesh.h"
+#include "bmesh_private.h"
+
+/*
+ * TEST EDGE SIDE and POINT IN TRIANGLE
+ *
+ * Point in triangle tests stolen from scanfill.c.
+ * Used for tesselator
+ *
+ */
+
+static short testedgeside(const double v1[2], const double v2[2], const double v3[2])
+{
+	/* is v3 to the right of v1 - v2 ? With exception: v3 == v1 || v3 == v2 */
+	double inp;
+
+	//inp = (v2[cox] - v1[cox]) * (v1[coy] - v3[coy]) + (v1[coy] - v2[coy]) * (v1[cox] - v3[cox]);
+	inp = (v2[0] - v1[0]) * (v1[1] - v3[1]) + (v1[1] - v2[1]) * (v1[0] - v3[0]);
+
+	if (inp < 0.0) {
+		return FALSE;
+	}
+	else if (inp == 0) {
+		if (v1[0] == v3[0] && v1[1] == v3[1]) return FALSE;
+		if (v2[0] == v3[0] && v2[1] == v3[1]) return FALSE;
+	}
+	return TRUE;
+}
+
+static short testedgesidef(const float v1[2], const float v2[2], const float v3[2])
+{
+	/* is v3 to the right of v1 - v2 ? With exception: v3 == v1 || v3 == v2 */
+	double inp;
+
+	//inp = (v2[cox] - v1[cox]) * (v1[coy] - v3[coy]) + (v1[coy] - v2[coy]) * (v1[cox] - v3[cox]);
+	inp = (v2[0] - v1[0]) * (v1[1] - v3[1]) + (v1[1] - v2[1]) * (v1[0] - v3[0]);
+
+	if (inp < 0.0) {
+		return FALSE;
+	}
+	else if (inp == 0) {
+		if (v1[0] == v3[0] && v1[1] == v3[1]) return FALSE;
+		if (v2[0] == v3[0] && v2[1] == v3[1]) return FALSE;
+	}
+	return TRUE;
+}
+
+static int point_in_triangle(const double v1[2], const double v2[2], const double v3[2], const double pt[2])
+{
+	if (testedgeside(v1, v2, pt) && testedgeside(v2, v3, pt) && testedgeside(v3, v1, pt)) {
+		return TRUE;
+	}
+	return FALSE;
+}
+
+/*
+ * COMPUTE POLY NORMAL
+ *
+ * Computes the normal of a planar
+ * polygon See Graphics Gems for
+ * computing newell normal.
+ *
+ */
+
+static void compute_poly_normal(float normal[3], float (*verts)[3], int nverts)
+{
+
+	float u[3], v[3], w[3]; /*, *w, v1[3], v2[3]; */
+	float n[3] = {0.0f, 0.0f, 0.0f} /*, l, v1[3], v2[3] */;
+	/* double l2; */
+	int i /*, s = 0 */;
+
+	/* this fixes some weird numerical erro */
+	add_v3_fl(verts[0], 0.0001f);
+
+	for (i = 0; i < nverts; i++) {
+		copy_v3_v3(u, verts[i]);
+		copy_v3_v3(v, verts[(i + 1) % nverts]);
+		copy_v3_v3(w, verts[(i + 2) % nverts]);
+		
+#if 0
+		sub_v3_v3v3(v1, w, v);
+		sub_v3_v3v3(v2, v, u);
+		normalize_v3(v1);
+		normalize_v3(v2);
+
+		l = dot_v3v3(v1, v2);
+		if (fabsf(l - 1.0) < 0.1f) {
+			continue;
+		}
+#endif
+		/* newell's method
+
+		so thats?:
+		(a[1] - b[1]) * (a[2] + b[2]);
+		a[1]*b[2] - b[1]*a[2] - b[1]*b[2] + a[1]*a[2]
+
+		odd.  half of that is the cross product. . .what's the
+		other half?
+
+		also could be like a[1]*(b[2] + a[2]) - b[1]*(a[2] - b[2])
+		*/
+
+		n[0] += (u[1] - v[1]) * (u[2] + v[2]);
+		n[1] += (u[2] - v[2]) * (u[0] + v[0]);
+		n[2] += (u[0] - v[0]) * (u[1] + v[1]);
+	}
+
+	if (normalize_v3_v3(normal, n) == 0.0f) {
+		normal[2] = 1.0f; /* other axis set to 0.0 */
+	}
+
+#if 0
+	l = len_v3(n);
+	/* fast square root, newton/babylonian method:
+	l2 = l * 0.1;
+
+	l2 = (l / l2 + l2) * 0.5;
+	l2 = (l / l2 + l2) * 0.5;
+	l2 = (l / l2 + l2) * 0.5;
+	*/
+
+	if (l == 0.0) {
+		normal[0] = 0.0f;
+		normal[1] = 0.0f;
+		normal[2] = 1.0f;
+
+		return;
+	}
+	else {
+		l = 1.0f / l;
+	}
+
+	mul_v3_fl(n, l);
+
+	copy_v3_v3(normal, n);
+#endif
+}
+
+/*
+ * COMPUTE POLY CENTER
+ *
+ * Computes the centroid and
+ * area of a polygon in the X/Y
+ * plane.
+ *
+ */
+
+static int compute_poly_center(float center[3], float *r_area, float (*verts)[3], int nverts)
+{
+	int i, j;
+	float atmp = 0.0f, xtmp = 0.0f, ytmp = 0.0f, ai;
+
+	zero_v3(center);
+
+	if (nverts < 3)
+		return FALSE;
+
+	i = nverts - 1;
+	j = 0;
+	
+	while (j < nverts) {
+		ai = verts[i][0] * verts[j][1] - verts[j][0] * verts[i][1];
+		atmp += ai;
+		xtmp += (verts[j][0] + verts[i][0]) * ai;
+		ytmp += (verts[j][1] + verts[i][1]) * ai;
+		i = j;
+		j += 1;
+	}
+
+	if (r_area)
+		*r_area = atmp / 2.0f;
+	
+	if (atmp != 0) {
+		center[0] = xtmp /  (3.0f * atmp);
+		center[1] = xtmp /  (3.0f * atmp);
+		return TRUE;
+	}
+	return FALSE;
+}
+
+float BM_face_area_calc(BMesh *bm, BMFace *f)
+{
+	BMLoop *l;
+	BMIter iter;
+	float (*verts)[3];
+	float center[3];
+	float area = 0.0f;
+	int i;
+
+	BLI_array_fixedstack_declare(verts, BM_NGON_STACK_SIZE, f->len, __func__);
+
+	i = 0;
+	BM_ITER(l, &iter, bm, BM_LOOPS_OF_FACE, f) {
+		copy_v3_v3(verts[i], l->v->co);
+		i++;
+	}
+
+	compute_poly_center(center, &area, verts, f->len);
+
+	BLI_array_fixedstack_free(verts);
+
+	return area;
+}
+
+/*
+ * computes center of face in 3d.  uses center of bounding box.
+ */
+void BM_face_center_bounds_calc(BMesh *bm, BMFace *f, float r_cent[3])
+{
+	BMIter iter;
+	BMLoop *l;
+	float min[3], max[3];
+	int i;
+
+	INIT_MINMAX(min, max);
+	l = BM_iter_new(&iter, bm, BM_LOOPS_OF_FACE, f);
+	for (i = 0; l; l = BM_iter_step(&iter), i++) {
+		DO_MINMAX(l->v->co, min, max);
+	}
+
+	mid_v3_v3v3(r_cent, min, max);
+}
+
+void BM_face_center_mean_calc(BMesh *bm, BMFace *f, float r_cent[3])
+{
+	BMIter iter;
+	BMLoop *l;
+	int i;
+
+	zero_v3(r_cent);
+
+	l = BM_iter_new(&iter, bm, BM_LOOPS_OF_FACE, f);
+	for (i = 0; l; l = BM_iter_step(&iter), i++) {
+		add_v3_v3(r_cent, l->v->co);
+	}
+
+	if (f->len) mul_v3_fl(r_cent, 1.0f / (float)f->len);
+}
+
+/*
+ * COMPUTE POLY PLANE
+ *
+ * Projects a set polygon's vertices to
+ * a plane defined by the average
+ * of its edges cross products
+ *
+ */
+
+void compute_poly_plane(float (*verts)[3], int nverts)
+{
+	
+	float avgc[3], norm[3], mag, avgn[3];
+	float *v1, *v2, *v3;
+	int i;
+	
+	if (nverts < 3)
+		return;
+
+	zero_v3(avgn);
+	zero_v3(avgc);
+
+	for (i = 0; i < nverts; i++) {
+		v1 = verts[i];
+		v2 = verts[(i + 1) % nverts];
+		v3 = verts[(i + 2) % nverts];
+		normal_tri_v3(norm, v1, v2, v3);
+
+		add_v3_v3(avgn, norm);
+	}
+
+	/* what was this bit for */
+	if (is_zero_v3(avgn)) {
+		avgn[0] = 0.0f;
+		avgn[1] = 0.0f;
+		avgn[2] = 1.0f;
+	}
+	else {
+		/* XXX, why is this being divided and _then_ normalized
+		 * division could be removed? - campbell */
+		avgn[0] /= nverts;
+		avgn[1] /= nverts;
+		avgn[2] /= nverts;
+		normalize_v3(avgn);
+	}
+	
+	for (i = 0; i < nverts; i++) {
+		v1 = verts[i];
+		mag = dot_v3v3(v1, avgn);
+		madd_v3_v3fl(v1, avgn, -mag);
+	}
+}
+
+/*
+ * BM LEGAL EDGES
+ *
+ * takes in a face and a list of edges, and sets to NULL any edge in
+ * the list that bridges a concave region of the face or intersects
+ * any of the faces's edges.
+ */
+#if 0 /* needs BLI math double versions of these functions */
+static void shrink_edged(double *v1, double *v2, double fac)
+{
+	double mid[3];
+
+	mid_v3_v3v3(mid, v1, v2);
+
+	sub_v3_v3v3(v1, v1, mid);
+	sub_v3_v3v3(v2, v2, mid);
+
+	mul_v3_fl(v1, fac);
+	mul_v3_fl(v2, fac);
+
+	add_v3_v3v3(v1, v1, mid);
+	add_v3_v3v3(v2, v2, mid);
+}
+#endif
+
+static void shrink_edgef(float v1[3], float v2[3], const float fac)
+{
+	float mid[3];
+
+	mid_v3_v3v3(mid, v1, v2);
+
+	sub_v3_v3v3(v1, v1, mid);
+	sub_v3_v3v3(v2, v2, mid);
+
+	mul_v3_fl(v1, fac);
+	mul_v3_fl(v2, fac);
+
+	add_v3_v3v3(v1, v1, mid);
+	add_v3_v3v3(v2, v2, mid);
+}
+
+
+/*
+ * POLY ROTATE PLANE
+ *
+ * Rotates a polygon so that it's
+ * normal is pointing towards the mesh Z axis
+ *
+ */
+
+void poly_rotate_plane(const float normal[3], float (*verts)[3], const int nverts)
+{
+
+	float up[3] = {0.0f, 0.0f, 1.0f}, axis[3], q[4];
+	float mat[3][3];
+	double angle;
+	int i;
+
+	cross_v3_v3v3(axis, normal, up);
+
+	angle = saacos(dot_v3v3(normal, up));
+
+	if (angle == 0.0) return;
+
+	axis_angle_to_quat(q, axis, (float)angle);
+	quat_to_mat3(mat, q);
+
+	for (i = 0;  i < nverts;  i++)
+		mul_m3_v3(mat, verts[i]);
+}
+
+/*
+ * BMESH UPDATE FACE NORMAL
+ *
+ * Updates the stored normal for the
+ * given face. Requires that a buffer
+ * of sufficient length to store projected
+ * coordinates for all of the face's vertices
+ * is passed in as well.
+ *
+ */
+
+void BM_face_normal_update(BMesh *bm, BMFace *f)
+{
+	if (f->len >= 3) {
+		float (*proj)[3];
+
+		BLI_array_fixedstack_declare(proj, BM_NGON_STACK_SIZE, f->len, __func__);
+
+		bmesh_update_face_normal(bm, f, f->no, proj);
+
+		BLI_array_fixedstack_free(proj);
+	}
+}
+/* same as BM_face_normal_update but takes vertex coords */
+void BM_face_normal_update_vcos(BMesh *bm, BMFace *f, float no[3], float (*vertexCos)[3])
+{
+	if (f->len >= 3) {
+		float (*proj)[3];
+
+		BLI_array_fixedstack_declare(proj, BM_NGON_STACK_SIZE, f->len, __func__);
+
+		bmesh_update_face_normal_vertex_cos(bm, f, no, proj, vertexCos);
+
+		BLI_array_fixedstack_free(proj);
+	}
+}
+
+void BM_edge_normals_update(BMesh *bm, BMEdge *e)
+{
+	BMIter iter;
+	BMFace *f;
+	
+	f = BM_iter_new(&iter, bm, BM_FACES_OF_EDGE, e);
+	for ( ; f; f = BM_iter_step(&iter)) {
+		BM_face_normal_update(bm, f);
+	}
+
+	BM_vert_normal_update(bm, e->v1);
+	BM_vert_normal_update(bm, e->v2);
+}
+
+void BM_vert_normal_update(BMesh *bm, BMVert *v)
+{
+	/* TODO, we can normalize each edge only once, then compare with previous edge */
+
+	BMIter eiter, liter;
+	BMEdge *e;
+	BMLoop *l;
+	float vec1[3], vec2[3], fac;
+	int len = 0;
+
+	zero_v3(v->no);
+
+	BM_ITER(e, &eiter, bm, BM_EDGES_OF_VERT, v) {
+		BM_ITER(l, &liter, bm, BM_LOOPS_OF_EDGE, e) {
+			if (l->v == v) {
+				/* Same calculation used in BM_mesh_normals_update */
+				sub_v3_v3v3(vec1, l->v->co, l->prev->v->co);
+				sub_v3_v3v3(vec2, l->next->v->co, l->v->co);
+				normalize_v3(vec1);
+				normalize_v3(vec2);
+
+				fac = saacos(-dot_v3v3(vec1, vec2));
+				
+				madd_v3_v3fl(v->no, l->f->no, fac);
+
+				len++;
+			}
+		}
+	}
+
+	if (len) {
+		normalize_v3(v->no);
+	}
+}
+
+void BM_vert_normal_update_all(BMesh *bm, BMVert *v)
+{
+	BMIter iter;
+	BMFace *f;
+	int len = 0;
+
+	f = BM_iter_new(&iter, bm, BM_FACES_OF_VERT, v);
+	for ( ; f; f = BM_iter_step(&iter), len++) {
+		BM_face_normal_update(bm, f);
+	}
+
+	BM_vert_normal_update(bm, v);
+}
+
+void bmesh_update_face_normal(BMesh *bm, BMFace *f, float no[3],
+                              float (*projectverts)[3])
+{
+	BMLoop *l;
+
+	/* common cases first */
+	switch (f->len) {
+		case 4:
+		{
+			BMVert *v1 = (l = BM_FACE_FIRST_LOOP(f))->v;
+			BMVert *v2 = (l = l->next)->v;
+			BMVert *v3 = (l = l->next)->v;
+			BMVert *v4 = (l->next)->v;
+			normal_quad_v3(no, v1->co, v2->co, v3->co, v4->co);
+			break;
+		}
+		case 3:
+		{
+			BMVert *v1 = (l = BM_FACE_FIRST_LOOP(f))->v;
+			BMVert *v2 = (l = l->next)->v;
+			BMVert *v3 = (l->next)->v;
+			normal_tri_v3(no, v1->co, v2->co, v3->co);
+			break;
+		}
+		case 0:
+		{
+			zero_v3(no);
+			break;
+		}
+		default:
+		{
+			BMIter iter;
+			int i = 0;
+			BM_ITER(l, &iter, bm, BM_LOOPS_OF_FACE, f) {
+				copy_v3_v3(projectverts[i], l->v->co);
+				i += 1;
+			}
+			compute_poly_normal(no, projectverts, f->len);
+			break;
+		}
+	}
+}
+/* exact same as 'bmesh_update_face_normal' but accepts vertex coords */
+void bmesh_update_face_normal_vertex_cos(BMesh *bm, BMFace *f, float no[3],
+                                         float (*projectverts)[3], float (*vertexCos)[3])
+{
+	BMLoop *l;
+
+	/* must have valid index data */
+	BLI_assert((bm->elem_index_dirty & BM_VERT) == 0);
+
+	/* common cases first */
+	switch (f->len) {
+		case 4:
+		{
+			BMVert *v1 = (l = BM_FACE_FIRST_LOOP(f))->v;
+			BMVert *v2 = (l = l->next)->v;
+			BMVert *v3 = (l = l->next)->v;
+			BMVert *v4 = (l->next)->v;
+			normal_quad_v3(no,
+			               vertexCos[BM_elem_index_get(v1)],
+			               vertexCos[BM_elem_index_get(v2)],
+			               vertexCos[BM_elem_index_get(v3)],
+			               vertexCos[BM_elem_index_get(v4)]);
+			break;
+		}
+		case 3:
+		{
+			BMVert *v1 = (l = BM_FACE_FIRST_LOOP(f))->v;
+			BMVert *v2 = (l = l->next)->v;
+			BMVert *v3 = (l->next)->v;
+			normal_tri_v3(no,
+			              vertexCos[BM_elem_index_get(v1)],
+			              vertexCos[BM_elem_index_get(v2)],
+			              vertexCos[BM_elem_index_get(v3)]);
+			break;
+		}
+		case 0:
+		{
+			zero_v3(no);
+			break;
+		}
+		default:
+		{
+			BMIter iter;
+			int i = 0;
+			BM_ITER(l, &iter, bm, BM_LOOPS_OF_FACE, f) {
+				copy_v3_v3(projectverts[i], vertexCos[BM_elem_index_get(l->v)]);
+				i += 1;
+			}
+			compute_poly_normal(no, projectverts, f->len);
+			break;
+		}
+	}
+}
+
+/*
+ * BMESH FLIP NORMAL
+ *
+ *  Reverses the winding of a face.
+ *  Note that this updates the calculated
+ *  normal.
+ */
+void BM_face_normal_flip(BMesh *bm, BMFace *f)
+{
+	bmesh_loop_reverse(bm, f);
+	negate_v3(f->no);
+}
+
+/* detects if two line segments cross each other (intersects).
+ * note, there could be more winding cases then there needs to be. */
+static int UNUSED_FUNCTION(linecrosses)(const double v1[2], const double v2[2], const double v3[2], const double v4[2])
+{
+	int w1, w2, w3, w4, w5;
+	
+	/* w1 = winding(v1, v3, v4);
+	w2 = winding(v2, v3, v4);
+	w3 = winding(v3, v1, v2);
+	w4 = winding(v4, v1, v2);
+	
+	return (w1 == w2) && (w3 == w4); */
+
+	w1 = testedgeside(v1, v3, v2);
+	w2 = testedgeside(v2, v4, v1);
+	w3 = !testedgeside(v1, v2, v3);
+	w4 = testedgeside(v3, v2, v4);
+	w5 = !testedgeside(v3, v1, v4);
+	return w1 == w2 && w2 == w3 && w3 == w4 && w4 == w5;
+}
+
+/* detects if two line segments cross each other (intersects).
+ * note, there could be more winding cases then there needs to be. */
+static int linecrossesf(const float v1[2], const float v2[2], const float v3[2], const float v4[2])
+{
+	int w1, w2, w3, w4, w5 /*, re */;
+	float mv1[2], mv2[2], mv3[2], mv4[2];
+	
+	/* now test windin */
+	w1 = testedgesidef(v1, v3, v2);
+	w2 = testedgesidef(v2, v4, v1);
+	w3 = !testedgesidef(v1, v2, v3);
+	w4 = testedgesidef(v3, v2, v4);
+	w5 = !testedgesidef(v3, v1, v4);
+	
+	if (w1 == w2 && w2 == w3 && w3 == w4 && w4 == w5) {
+		return TRUE;
+	}
+	
+#define GETMIN2_AXIS(a, b, ma, mb, axis) ma[axis] = MIN2(a[axis], b[axis]), mb[axis] = MAX2(a[axis], b[axis])
+#define GETMIN2(a, b, ma, mb) GETMIN2_AXIS(a, b, ma, mb, 0); GETMIN2_AXIS(a, b, ma, mb, 1);
+	
+	GETMIN2(v1, v2, mv1, mv2);
+	GETMIN2(v3, v4, mv3, mv4);
+	
+	/* do an interval test on the x and y axe */
+	/* first do x axi */
+#define T FLT_EPSILON * 15
+	if ( ABS(v1[1] - v2[1]) < T &&
+	     ABS(v3[1] - v4[1]) < T &&
+	     ABS(v1[1] - v3[1]) < T)
+	{
+		return (mv4[0] >= mv1[0] && mv3[0] <= mv2[0]);
+	}
+
+	/* now do y axi */
+	if ( ABS(v1[0] - v2[0]) < T &&
+	     ABS(v3[0] - v4[0]) < T &&
+	     ABS(v1[0] - v3[0]) < T)
+	{
+		return (mv4[1] >= mv1[1] && mv3[1] <= mv2[1]);
+	}
+
+	return FALSE;
+}
+
+/*
+ *  BM POINT IN FACE
+ *
+ * Projects co onto face f, and returns true if it is inside
+ * the face bounds.  Note that this uses a best-axis projection
+ * test, instead of projecting co directly into f's orientation
+ * space, so there might be accuracy issues.
+ */
+int BM_face_point_inside_test(BMesh *bm, BMFace *f, const float co[3])
+{
+	int ax, ay;
+	float co2[3], cent[3] = {0.0f, 0.0f, 0.0f}, out[3] = {FLT_MAX * 0.5f, FLT_MAX * 0.5f, 0};
+	BMLoop *l_iter;
+	BMLoop *l_first;
+	int crosses = 0;
+	float eps = 1.0f + (float)FLT_EPSILON * 150.0f;
+	
+	if (dot_v3v3(f->no, f->no) <= FLT_EPSILON * 10)
+		BM_face_normal_update(bm, f);
+	
+	/* find best projection of face XY, XZ or YZ: barycentric weights of
+	 * the 2d projected coords are the same and faster to compute
+	 *
+	 * this probably isn't all that accurate, but it has the advantage of
+	 * being fast (especially compared to projecting into the face orientation)
+	 */
+	axis_dominant_v3(&ax, &ay, f->no);
+
+	co2[0] = co[ax];
+	co2[1] = co[ay];
+	co2[2] = 0;
+	
+	l_iter = l_first = BM_FACE_FIRST_LOOP(f);
+	do {
+		cent[0] += l_iter->v->co[ax];
+		cent[1] += l_iter->v->co[ay];
+	} while ((l_iter = l_iter->next) != l_first);
+	
+	mul_v2_fl(cent, 1.0f / (float)f->len);
+	
+	l_iter = l_first = BM_FACE_FIRST_LOOP(f);
+	do {
+		float v1[3], v2[3];
+		
+		v1[0] = (l_iter->prev->v->co[ax] - cent[ax]) * eps + cent[ax];
+		v1[1] = (l_iter->prev->v->co[ay] - cent[ay]) * eps + cent[ay];
+		v1[2] = 0.0f;
+		
+		v2[0] = (l_iter->v->co[ax] - cent[ax]) * eps + cent[ax];
+		v2[1] = (l_iter->v->co[ay] - cent[ay]) * eps + cent[ay];
+		v2[2] = 0.0f;
+		
+		crosses += linecrossesf(v1, v2, co2, out) != 0;
+	} while ((l_iter = l_iter->next) != l_first);
+	
+	return crosses % 2 != 0;
+}
+
+static int goodline(float (*projectverts)[3], BMFace *f, int v1i,
+                    int v2i, int v3i, int UNUSED(nvert))
+{
+	BMLoop *l_iter;
+	BMLoop *l_first;
+	double v1[3], v2[3], v3[3], pv1[3], pv2[3];
+	int i;
+
+	VECCOPY(v1, projectverts[v1i]);
+	VECCOPY(v2, projectverts[v2i]);
+	VECCOPY(v3, projectverts[v3i]);
+	
+	if (testedgeside(v1, v2, v3)) {
+		return FALSE;
+	}
+
+	//for (i = 0; i < nvert; i++) {
+	l_iter = l_first = BM_FACE_FIRST_LOOP(f);
+	do {
+		i = BM_elem_index_get(l_iter->v);
+		if (i == v1i || i == v2i || i == v3i) {
+			continue;
+		}
+		
+		VECCOPY(pv1, projectverts[BM_elem_index_get(l_iter->v)]);
+		VECCOPY(pv2, projectverts[BM_elem_index_get(l_iter->next->v)]);
+		
+		//if (linecrosses(pv1, pv2, v1, v3)) return FALSE;
+
+		if ( point_in_triangle(v1, v2, v3, pv1) ||
+		     point_in_triangle(v3, v2, v1, pv1))
+		{
+			return FALSE;
+		}
+	} while ((l_iter = l_iter->next) != l_first);
+	return TRUE;
+}
+/*
+ * FIND EAR
+ *
+ * Used by tesselator to find
+ * the next triangle to 'clip off'
+ * of a polygon while tesselating.
+ *
+ */
+
+static BMLoop *find_ear(BMesh *UNUSED(bm), BMFace *f, float (*verts)[3], const int nvert)
+{
+	BMVert *v1, *v2, *v3;
+	BMLoop *bestear = NULL;
+
+	BMLoop *l_iter;
+	BMLoop *l_first;
+	/* float angle, bestangle = 180.0f; */
+	int isear /*, i = 0 */;
+	
+	l_iter = l_first = BM_FACE_FIRST_LOOP(f);
+	do {
+		isear = 1;
+		
+		v1 = l_iter->prev->v;
+		v2 = l_iter->v;
+		v3 = l_iter->next->v;
+
+		if (BM_edge_exists(v1, v3)) {
+			isear = 0;
+		}
+		else if (!goodline(verts, f, BM_elem_index_get(v1), BM_elem_index_get(v2), BM_elem_index_get(v3), nvert)) {
+			isear = 0;
+		}
+
+		if (isear) {
+#if 0
+			angle = angle_v3v3v3(verts[v1->head.eflag2], verts[v2->head.eflag2], verts[v3->head.eflag2]);
+			if (!bestear || ABS(angle-45.0f) < bestangle) {
+				bestear = l;
+				bestangle = ABS(45.0f - angle);
+			}
+			
+			if (angle > 20 && angle < 90) break;
+			if (angle < 100 && i > 5) break;
+			i += 1;
+#endif
+
+			bestear = l_iter;
+			break;
+		}
+	} while ((l_iter = l_iter->next) != l_first);
+
+	return bestear;
+}
+
+/*
+ * BMESH TRIANGULATE FACE
+ *
+ * Triangulates a face using a
+ * simple 'ear clipping' algorithm
+ * that tries to favor non-skinny
+ * triangles (angles less than
+ * 90 degrees). If the triangulator
+ * has bits left over (or cannot
+ * triangulate at all) it uses a
+ * simple fan triangulation
+ *
+ * newfaces, if non-null, must be an array of BMFace pointers,
+ * with a length equal to f->len.  it will be filled with the new
+ * triangles, and will be NULL-terminated.
+ */
+void BM_face_triangulate(BMesh *bm, BMFace *f, float (*projectverts)[3],
+                         const short newedge_oflag, const short newface_oflag, BMFace **newfaces)
+{
+	int i, done, nvert, nf_i = 0;
+	BMLoop *newl, *nextloop;
+	BMLoop *l_iter;
+	BMLoop *l_first;
+	/* BMVert *v; */ /* UNUSED */
+
+	/* copy vertex coordinates to vertspace arra */
+	i = 0;
+	l_iter = l_first = BM_FACE_FIRST_LOOP(f);
+	do {
+		copy_v3_v3(projectverts[i], l_iter->v->co);
+		BM_elem_index_set(l_iter->v, i); /* set dirty! */
+		i++;
+	} while ((l_iter = l_iter->next) != l_first);
+
+	bm->elem_index_dirty |= BM_VERT; /* see above */
+
+	///bmesh_update_face_normal(bm, f, f->no, projectverts);
+
+	compute_poly_normal(f->no, projectverts, f->len);
+	poly_rotate_plane(f->no, projectverts, i);
+
+	nvert = f->len;
+
+	//compute_poly_plane(projectverts, i);
+	for (i = 0; i < nvert; i++) {
+		projectverts[i][2] = 0.0f;
+	}
+
+	done = 0;
+	while (!done && f->len > 3) {
+		done = 1;
+		l_iter = find_ear(bm, f, projectverts, nvert);
+		if (l_iter) {
+			done = 0;
+			/* v = l->v; */ /* UNUSED */
+			f = BM_face_split(bm, l_iter->f, l_iter->prev->v,
+			                  l_iter->next->v,
+			                  &newl, NULL);
+			copy_v3_v3(f->no, l_iter->f->no);
+
+			if (!f) {
+				fprintf(stderr, "%s: triangulator failed to split face! (bmesh internal error)\n", __func__);
+				break;
+			}
+
+			BMO_elem_flag_enable(bm, newl->e, newedge_oflag);
+			BMO_elem_flag_enable(bm, f, newface_oflag);
+			
+			if (newfaces) newfaces[nf_i++] = f;
+
+			/* l = f->loopbase;
+			do {
+				if (l->v == v) {
+					f->loopbase = l;
+					break;
+				}
+				l = l->next;
+			} while (l != f->loopbase); */
+		}
+	}
+
+	if (f->len > 3) {
+		l_iter = BM_FACE_FIRST_LOOP(f);
+		while (l_iter->f->len > 3) {
+			nextloop = l_iter->next->next;
+			f = BM_face_split(bm, l_iter->f, l_iter->v, nextloop->v,
+			                  &newl, NULL);
+			if (!f) {
+				printf("triangle fan step of triangulator failed.\n");
+
+				/* NULL-terminat */
+				if (newfaces) newfaces[nf_i] = NULL;
+				return;
+			}
+
+			if (newfaces) newfaces[nf_i++] = f;
+			
+			BMO_elem_flag_enable(bm, newl->e, newedge_oflag);
+			BMO_elem_flag_enable(bm, f, newface_oflag);
+			l_iter = nextloop;
+		}
+	}
+	
+	/* NULL-terminat */
+	if (newfaces) newfaces[nf_i] = NULL;
+}
+
+/* each pair of loops defines a new edge, a split.  this function goes
+ * through and sets pairs that are geometrically invalid to null.  a
+ * split is invalid, if it forms a concave angle or it intersects other
+ * edges in the face, or it intersects another split.  in the case of
+ * intersecting splits, only the first of the set of intersecting
+ * splits survives */
+void BM_face_legal_splits(BMesh *bm, BMFace *f, BMLoop *(*loops)[2], int len)
+{
+	BMIter iter;
+	BMLoop *l;
+	float v1[3], v2[3], v3[3]/*, v4[3 */, no[3], mid[3], *p1, *p2, *p3, *p4;
+	float out[3] = {-234324.0f, -234324.0f, 0.0f};
+	float (*projverts)[3];
+	float (*edgeverts)[3];
+	float fac1 = 1.0000001f, fac2 = 0.9f; //9999f; //0.999f;
+	int i, j, a = 0, clen;
+
+	BLI_array_fixedstack_declare(projverts, BM_NGON_STACK_SIZE, f->len,      "projvertsb");
+	BLI_array_fixedstack_declare(edgeverts, BM_NGON_STACK_SIZE * 2, len * 2, "edgevertsb");
+	
+	i = 0;
+	l = BM_iter_new(&iter, bm, BM_LOOPS_OF_FACE, f);
+	for ( ; l; l = BM_iter_step(&iter)) {
+		BM_elem_index_set(l, i); /* set_loop */
+		copy_v3_v3(projverts[i], l->v->co);
+		i++;
+	}
+	
+	for (i = 0; i < len; i++) {
+		copy_v3_v3(v1, loops[i][0]->v->co);
+		copy_v3_v3(v2, loops[i][1]->v->co);
+
+		shrink_edgef(v1, v2, fac2);
+		
+		copy_v3_v3(edgeverts[a], v1);
+		a++;
+		copy_v3_v3(edgeverts[a], v2);
+		a++;
+	}
+	
+	compute_poly_normal(no, projverts, f->len);
+	poly_rotate_plane(no, projverts, f->len);
+	poly_rotate_plane(no, edgeverts, len * 2);
+
+	for (i = 0, l = BM_FACE_FIRST_LOOP(f); i < f->len; i++, l = l->next) {
+		p1 = projverts[i];
+		out[0] = MAX2(out[0], p1[0]) + 0.01f;
+		out[1] = MAX2(out[1], p1[1]) + 0.01f;
+		out[2] = 0.0f;
+		p1[2] = 0.0f;
+
+		//copy_v3_v3(l->v->co, p1);
+	}
+	
+	for (i = 0; i < len; i++) {
+		edgeverts[i * 2][2] = 0.0f;
+		edgeverts[i * 2 + 1][2] = 0.0f;
+	}
+
+	/* do convexity test */
+	for (i = 0; i < len; i++) {
+		copy_v3_v3(v2, edgeverts[i * 2]);
+		copy_v3_v3(v3, edgeverts[i * 2 + 1]);
+
+		mid_v3_v3v3(mid, v2, v3);
+		
+		clen = 0;
+		for (j = 0; j < f->len; j++) {
+			p1 = projverts[j];
+			p2 = projverts[(j + 1) % f->len];
+			
+			copy_v3_v3(v1, p1);
+			copy_v3_v3(v2, p2);
+
+			shrink_edgef(v1, v2, fac1);
+
+			if (linecrossesf(p1, p2, mid, out)) clen++;
+		}
+		
+		if (clen % 2 == 0) {
+			loops[i][0] = NULL;
+		}
+	}
+	
+	/* do line crossing test */
+	for (i = 0; i < f->len; i++) {
+		p1 = projverts[i];
+		p2 = projverts[(i + 1) % f->len];
+		
+		copy_v3_v3(v1, p1);
+		copy_v3_v3(v2, p2);
+
+		shrink_edgef(v1, v2, fac1);
+
+		for (j = 0; j < len; j++) {
+			if (!loops[j][0]) continue;
+
+			p3 = edgeverts[j * 2];
+			p4 = edgeverts[j * 2 + 1];
+
+			if (linecrossesf(v1, v2, p3, p4)) {
+				loops[j][0] = NULL;
+			}
+		}
+	}
+
+	for (i = 0; i < len; i++) {
+		for (j = 0; j < len; j++) {
+			if (j == i) continue;
+			if (!loops[i][0]) continue;
+			if (!loops[j][0]) continue;
+
+			p1 = edgeverts[i * 2];
+			p2 = edgeverts[i * 2 + 1];
+			p3 = edgeverts[j * 2];
+			p4 = edgeverts[j * 2 + 1];
+
+			copy_v3_v3(v1, p1);
+			copy_v3_v3(v2, p2);
+
+			shrink_edgef(v1, v2, fac1);
+
+			if (linecrossesf(v1, v2, p3, p4)) {
+				loops[i][0] = NULL;
+			}
+		}
+	}
+
+	BLI_array_fixedstack_free(projverts);
+	BLI_array_fixedstack_free(edgeverts);
+}