the make ngon function's overlap test needed some work, the API function
[blender.git] / source / blender / bmesh / intern / bmesh_polygon.c
index b5a65932743a9d1f62c9c2a55d512843dafc04d6..3c4b738a84082cd33ee6bb5d0cb5239ee672e3d1 100644 (file)
@@ -1,9 +1,16 @@
-#include "bmesh.h"
-#include "bmesh_private.h"
+#include <string.h>
+#include <math.h>
+#include <stdlib.h>
 
-#include "BLI_arithb.h"
 #include "BKE_utildefines.h"
-#include <string.h>
+
+#include "BLI_arithb.h"
+#include "BLI_blenlib.h"
+
+#include "MEM_guardedalloc.h"
+
+#include "bmesh.h"
+#include "bmesh_private.h"
 
 /*
  *
  *
 */
 
-static short testedgeside(float *v1, float *v2, float *v3)
+static short testedgeside(double *v1, double *v2, double *v3)
 /* is v3 to the right of v1-v2 ? With exception: v3==v1 || v3==v2 */
 {
-       float inp;
+       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]);
@@ -45,7 +52,7 @@ static short testedgeside(float *v1, float *v2, float *v3)
        return 1;
 }
 
-static int point_in_triangle(float *v1, float *v2, float *v3, float *pt)
+static int point_in_triangle(double *v1, double *v2, double *v3, double *pt)
 {
        if(testedgeside(v1,v2,pt) && testedgeside(v2,v3,pt) && testedgeside(v3,v1,pt))
                return 1;
@@ -88,24 +95,47 @@ static int convexangle(float *v1t, float *v2t, float *v3t)
 static void compute_poly_normal(float normal[3], float (*verts)[3], int nverts)
 {
 
-       float *u,  *v;
+       float *u,  *v;/*, *w, v1[3], v2[3];*/
+       double n[3] = {0.0, 0.0, 0.0}, l;
        int i;
 
-       normal[0] = 0.0;
-       normal[1] = 0.0;
-       normal[2] = 0.0;
-       
        for(i = 0; i < nverts; i++){
                u = verts[i];
                v = verts[(i+1) % nverts];
+               /*w = verts[(i+2) % nverts];
+
+               VecSubf(v1, u, v);
+               VecSubf(v2, w, v);
+               Crossf(normal, v1, v2);
+               Normalize(normal);
                
-               normal[0] += (u[1] - v[1]) * (u[2] + v[2]);
-               normal[1] += (u[2] - v[2]) * (u[0] + v[0]);
-               normal[2] += (u[0] - v[0]) * (u[1] + v[1]);
-               i++;
-       }
+               return;*/
+               
+               /* 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]
 
-       Normalize(normal);
+               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]);
+       }
+       
+       l = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
+       n[0] /= l;
+       n[1] /= l;
+       n[2] /= l;
+       
+       normal[0] = n[0];
+       normal[1] = n[1];
+       normal[2] = n[2];
 }
 
 /*
@@ -140,11 +170,11 @@ static int compute_poly_center(float center[3], float *area, float (*verts)[3],
        }
 
        if(area)
-               *area = atmp / 2.0    
+               *area = atmp / 2.0f;    
        
        if (atmp != 0){
-               center[0] = xtmp /  (3.0 * atmp);
-               center[1] = xtmp /  (3.0 * atmp);
+               center[0] = xtmp /  (3.0f * atmp);
+               center[1] = xtmp /  (3.0f * atmp);
                return 1;
        }
        return 0;
@@ -223,14 +253,23 @@ void compute_poly_plane(float (*verts)[3], int nverts)
 void poly_rotate_plane(float normal[3], float (*verts)[3], int nverts)
 {
 
-       float up[3] = {0.0,0.0,1.0}, axis[3], angle, q[4];
+       float up[3] = {0.0f,0.0f,1.0f}, axis[3], q[4];
        float mat[3][3];
+       double angle;
        int i;
 
+       compute_poly_normal(normal, verts, nverts);
+
        Crossf(axis, up, normal);
-       angle = VecAngle2(normal, up);
-       
-       AxisAngleToQuat(q, axis, angle);
+       axis[0] *= -1;
+       axis[1] *= -1;
+       axis[2] *= -1;
+
+       angle = saacos(normal[0]*up[0]+normal[1]*up[1] + normal[2]*up[2]);
+
+       if (angle == 0.0f) return;
+
+       AxisAngleToQuatd(q, axis, angle);
        QuatToMat3(q, mat);
 
        for(i = 0;  i < nverts;  i++)
@@ -248,6 +287,56 @@ void poly_rotate_plane(float normal[3], float (*verts)[3], int nverts)
  *
 */
 
+void BM_Face_UpdateNormal(BMesh *bm, BMFace *f)
+{
+       float projverts[12][3];
+       float (*proj)[3] = f->len < 12 ? projverts : MEM_mallocN(sizeof(float)*f->len*3, "projvertsn");
+       BMLoop *l = f->loopbase;
+       int i=0;
+
+       if (f->len < 3) return;
+       
+       do {
+               VECCOPY(proj[i], l->v->co);
+               i += 1;
+       } while (l != f->loopbase);
+
+       bmesh_update_face_normal(bm, f, proj);
+
+       if (projverts != proj) MEM_freeN(proj);
+}
+
+void BM_Edge_UpdateNormals(BMesh *bm, BMEdge *e)
+{
+       BMIter *iter;
+       BMFace *f;
+       
+       f = BMIter_New(&iter, bm, BM_FACES_OF_EDGE, e);
+       for (; f; f=BMIter_Step(&iter)) {
+               BM_Face_UpdateNormal(bm, f);
+       }
+
+       BM_Vert_UpdateNormal(bm, e->v1);
+       BM_Vert_UpdateNormal(bm, e->v2);
+}
+
+void BM_Vert_UpdateNormal(BMesh *bm, BMVert *v)
+{
+       BMIter iter;
+       BMFace *f;
+       float norm[3] = {0.0f, 0.0f, 0.0f};
+       int len=0;
+
+       f = BMIter_New(&iter, bm, BM_FACES_OF_VERT, v);
+       for (; f; f=BMIter_Step(&iter), len++) {
+               VecAddf(norm, f->no, norm);
+       }
+
+       if (!len) return;
+
+       VecMulf(norm, 1.0f/(int)len);
+}
+
 void bmesh_update_face_normal(BMesh *bm, BMFace *f, float (*projectverts)[3])
 {
        BMLoop *l;
@@ -302,6 +391,79 @@ void BM_flip_normal(BMesh *bm, BMFace *f)
 
 
 
+int winding(double *a, double *b, double *c)
+{
+       double v1[3], v2[3], v[3];
+       
+       VECSUB(v1, b, a);
+       VECSUB(v2, b, c);
+       
+       v1[2] = 0;
+       v2[2] = 0;
+       
+       Normalize_d(v1);
+       Normalize_d(v2);
+       
+       Crossd(v, v1, v2);
+
+       /*!! (turns nonzero into 1) is likely not necassary, 
+         since '>' I *think* should always
+         return 0 or 1, but I'm not totally sure. . .*/
+       return !!(v[2] > 0);
+}
+
+/* detects if two line segments cross each other (intersects).
+   note, there could be more winding cases then there needs to be. */
+int linecrosses(double *v1, double *v2, double *v3, double *v4)
+{
+       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 = winding(v1, v3, v2);
+       w2 = winding(v2, v4, v1);
+       w3 = !winding(v1, v2, v3);
+       w4 = winding(v3, v2, v4);
+       w5 = !winding(v3, v1, v4);
+       return w1 == w2 && w2 == w3 && w3 == w4 && w4==w5;
+}
+
+int goodline(float (*projectverts)[3], BMFace *f, int v1i,
+            int v2i, int v3i, int nvert) {
+       BMLoop *l = f->loopbase;
+       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 0;
+       
+       //for (i=0; i<nvert; i++) {
+       do {
+               i = l->v->head.eflag2;
+               if (i == v1i || i == v2i || i == v3i) {
+                       l = l->head.next;
+                       continue;
+               }
+               
+               VECCOPY(pv1, projectverts[l->v->head.eflag2]);
+               VECCOPY(pv2, projectverts[((BMLoop*)l->head.next)->v->head.eflag2]);
+               
+               //if (linecrosses(pv1, pv2, v1, v3)) return 0;
+               if (point_in_triangle(v1, v2, v3, pv1)) return 0;
+               if (point_in_triangle(v3, v2, v1, pv1)) return 0;
+
+               l = l->head.next;
+       } while (l != f->loopbase);
+       return 1;
+}
 /*
  * FIND EAR
  *
@@ -311,12 +473,18 @@ void BM_flip_normal(BMesh *bm, BMFace *f)
  *
 */
 
-static BMLoop *find_ear(BMFace *f, float (*verts)[3])
+typedef struct quadtree {
+       int *dsdf;
+       int *dsfds;
+} quadtree;
+
+static BMLoop *find_ear(BMesh *bm, BMFace *f, float (*verts)[3], 
+                       int nvert)
 {
        BMVert *v1, *v2, *v3;
-       BMLoop *bestear = NULL, *l, *l2;
+       BMLoop *bestear = NULL, *l;
        float angle, bestangle = 180.0f;
-       int isear;
+       int isear, i=0;
        
        l = f->loopbase;
        do{
@@ -328,23 +496,20 @@ static BMLoop *find_ear(BMFace *f, float (*verts)[3])
 
                if (BM_Edge_Exist(v1, v3)) isear = 0;
 
-               if(isear && convexangle(verts[v1->head.eflag2], verts[v2->head.eflag2], verts[v3->head.eflag2])){
-                       for(l2 = ((BMLoop*)(l->head.next->next)); l2 != ((BMLoop*)(l->head.prev)); l2 = ((BMLoop*)(l2->head.next)) ){
-                               if(point_in_triangle(verts[v1->head.eflag2], verts[v2->head.eflag2],verts[v3->head.eflag2], l2->v->co)){
-                                       isear = 0;
-                                       break;
-                               }
-                       }
-               } else isear = 0;
+               if (isear && !goodline(verts, f, v1->head.eflag2, v2->head.eflag2,
+                                      v3->head.eflag2, nvert))
+                       isear = 0;
 
                if(isear){
                        angle = VecAngle3(verts[v1->head.eflag2], verts[v2->head.eflag2], verts[v3->head.eflag2]);
-                       if(!bestear || angle < bestangle){
+                       if(!bestear || ABS(angle-45.0f) < bestangle){
                                bestear = l;
-                               bestangle = angle;
+                               bestangle = ABS(45.0f-angle);
                        }
-                       if(angle < 90.0)
-                               break;
+                       
+                       if (angle > 20 && angle < 90) break;
+                       if (angle < 100 && i > 5) break;
+                       i += 1;
                }
                l = (BMLoop*)(l->head.next);
        }
@@ -369,10 +534,13 @@ static BMLoop *find_ear(BMFace *f, float (*verts)[3])
  * -Modify this to try and find ears that will not create a non-manifold face after conversion back to editmesh
  *
 */
-void BM_Triangulate_Face(BMesh *bm, BMFace *f, float (*projectverts)[3], int newedgeflag, int newfaceflag)
+void BM_Triangulate_Face(BMesh *bm, BMFace *f, float (*projectverts)[3], 
+                        int newedgeflag, int newfaceflag)
 {
-       int i, done;
+       int i, done, nvert;
+       float no[3];
        BMLoop *l, *newl, *nextloop;
+       BMVert *v;
 
        /*copy vertex coordinates to vertspace array*/
        i = 0;
@@ -384,20 +552,51 @@ void BM_Triangulate_Face(BMesh *bm, BMFace *f, float (*projectverts)[3], int new
                l = (BMLoop*)(l->head.next);
        }while(l != f->loopbase);
        
-       bmesh_update_face_normal(bm, f, projectverts);
+       ///bmesh_update_face_normal(bm, f, projectverts);
 
-       compute_poly_plane(projectverts, i);
+       /*this fixes some weird numerical error*/
+       projectverts[0][0] += 0.0001f;
+       projectverts[0][1] += 0.0001f;
+       projectverts[0][2] += 0.0001f;
+
+       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 = find_ear(f, projectverts);
+               l = find_ear(bm, f, projectverts, nvert);
                if(l) {
                        done = 0;
-                       f = bmesh_sfme(bm, f, ((BMLoop*)(l->head.prev))->v, ((BMLoop*)(l->head.next))->v, &newl);
+                       v = l->v;
+                       f = BM_Split_Face(bm, l->f, ((BMLoop*)(l->head.prev))->v, 
+                                         ((BMLoop*)(l->head.next))->v, 
+                                         &newl, NULL, 0);
+                       VECCOPY(f->no, l->f->no);
+
+                       if (!f) {
+                               printf("yeek! triangulator failed to split face!\n");
+                               break;
+                       }
+
                        BMO_SetFlag(bm, newl->e, newedgeflag);
                        BMO_SetFlag(bm, f, newfaceflag);
+
+                       /*l = f->loopbase;
+                       do {
+                               if (l->v == v) {
+                                       f->loopbase = l;
+                                       break;
+                               }
+                               l = l->head.next;
+                       } while (l != f->loopbase);*/
                }
        }
 
@@ -405,7 +604,13 @@ void BM_Triangulate_Face(BMesh *bm, BMFace *f, float (*projectverts)[3], int new
                l = f->loopbase;
                while (l->f->len > 3){
                        nextloop = ((BMLoop*)(l->head.next->next));
-                       bmesh_sfme(bm, l->f, l->v,nextloop->v, &newl);
+                       f = BM_Split_Face(bm, l->f, l->v, nextloop->v, 
+                                         &newl, NULL, 0);
+                       if (!f) {
+                               printf("triangle fan step of triangulator failed.\n");
+                               return;
+                       }
+
                        BMO_SetFlag(bm, newl->e, newedgeflag);
                        BMO_SetFlag(bm, f, newfaceflag);
                        l = nextloop;