the make ngon function's overlap test needed some work, the API function
[blender.git] / source / blender / bmesh / intern / bmesh_polygon.c
index 149a356b5a6aec75c95b5114f210d46a4fcc30cf..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]
+
+               odd.  half of that is the cross product. . .what's the
+               other half?
 
-       Normalize(normal);
+               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];
 }
 
 /*
@@ -223,16 +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);
+       axis[0] *= -1;
+       axis[1] *= -1;
+       axis[2] *= -1;
 
-       if (angle == 0.0) return;
-       
-       AxisAngleToQuat(q, axis, angle);
+       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++)
@@ -250,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;
@@ -304,20 +391,20 @@ void BM_flip_normal(BMesh *bm, BMFace *f)
 
 
 
-int winding(float *a, float *b, float *c)
+int winding(double *a, double *b, double *c)
 {
-       float v1[3], v2[3], v[3];
+       double v1[3], v2[3], v[3];
        
-       VecSubf(v1, b, a);
-       VecSubf(v2, b, c);
+       VECSUB(v1, b, a);
+       VECSUB(v2, b, c);
        
        v1[2] = 0;
        v2[2] = 0;
        
-       Normalize(v1);
-       Normalize(v2);
+       Normalize_d(v1);
+       Normalize_d(v2);
        
-       Crossf(v, v1, v2);
+       Crossd(v, v1, v2);
 
        /*!! (turns nonzero into 1) is likely not necassary, 
          since '>' I *think* should always
@@ -327,7 +414,7 @@ int winding(float *a, float *b, float *c)
 
 /* detects if two line segments cross each other (intersects).
    note, there could be more winding cases then there needs to be. */
-int linecrosses(float *v1, float *v2, float *v3, float *v4)
+int linecrosses(double *v1, double *v2, double *v3, double *v4)
 {
        int w1, w2, w3, w4, w5;
        
@@ -346,104 +433,37 @@ int linecrosses(float *v1, float *v2, float *v3, float *v4)
        return w1 == w2 && w2 == w3 && w3 == w4 && w4==w5;
 }
 
-static int goodline_notworking(BMesh *bm, BMFace *f, BMVert *v1, BMVert *v2, BMVert *v3,
-            float (*p)[3], float *outv, int nvert)
-{
-       BMIter iter;
-       BMLoop *l;
-       int i, ret = 1;
-       float r = 0.00001f;
-       
-       /*cases of a vector being too close to
-         an axis can cause problems, so this is to
-         prevent that.*/
-       for (i=0; i<3; i++) {
-               p[v1->head.eflag2][i] += r;
-               p[v2->head.eflag2][i] -= r;
-               p[v3->head.eflag2][i] += r;             
-       }
-
-       //l = BMIter_New(&iter, bm, BM_LOOPS_OF_FACE, f);
-       //for (; l; l=BMIter_Step(&iter)) {
-       if (convexangle(p[v1->head.eflag2], p[v2->head.eflag2],
-           p[v3->head.eflag2])) {
-               ret = 0;
-               goto cleanup;
-       }
-
-       for (i=0; i<nvert; i++) {
-               if (i!=v1->head.eflag2 && i!=v2->head.eflag2 && 
-                   i!=v3->head.eflag2)
-               {
-                       if (point_in_triangle(p[v3->head.eflag2], 
-                           p[v2->head.eflag2], p[v1->head.eflag2], 
-                           p[i]))
-                       {
-                               ret = 0;
-                               goto cleanup;
-                       }
-               }
-       }
-
-cleanup:
-       for (i=0; i<3; i++) {
-               p[v1->head.eflag2][i] -= r;
-               p[v2->head.eflag2][i] += r;
-               p[v3->head.eflag2][i] -= r;             
-       }
-
-       return ret;
-}
+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;
 
-static int goodline(float (*projectverts)[3], int v1i, int v2i, int nvert, float *outv)
-{
-       /*the hardcoded stuff here, 0.999 and 1.0001, may be problems
-         in the future, not sure. - joeedh*/
-       float v1[3], v2[3], p[3], vv1[3], vv2[3], mid[3], a[3], b[3];
-       int i = 0, lcount=0;
-       
-       VECCOPY(v1, projectverts[v1i])
-       VECCOPY(v2, projectverts[v2i])
-       
-       VecAddf(p, v1, v2);
-       VecMulf(p, 0.5f);
-       
-       VecSubf(a, v1, p);
-       VecSubf(b, v2, p);
+       VECCOPY(v1, projectverts[v1i]);
+       VECCOPY(v2, projectverts[v2i]);
+       VECCOPY(v3, projectverts[v3i]);
        
-       VecMulf(a, 0.999f);
-       VecMulf(b, 0.999f);
+       if (testedgeside(v1, v2, v3)) return 0;
        
-       VecAddf(v1, p, a);
-       VecAddf(v2, p, b);
-       
-       while (i < nvert) {
-               VECCOPY(vv1, projectverts[i]);
-               VECCOPY(vv2, projectverts[(i+1)%nvert]);
-               
-               if (linecrosses(vv1, vv2, v1, v2)) return 0;
-
-               VecAddf(mid, vv1, vv2);
-               VecMulf(mid, 0.5f);
-               
-               VecSubf(a, vv1, mid);
-               VecSubf(b, vv2, mid);
+       //for (i=0; i<nvert; i++) {
+       do {
+               i = l->v->head.eflag2;
+               if (i == v1i || i == v2i || i == v3i) {
+                       l = l->head.next;
+                       continue;
+               }
                
-               VecMulf(a, 1.0001f);
-               VecMulf(b, 1.0001f);
+               VECCOPY(pv1, projectverts[l->v->head.eflag2]);
+               VECCOPY(pv2, projectverts[((BMLoop*)l->head.next)->v->head.eflag2]);
                
-               VecAddf(vv1, mid, a);
-               VecAddf(vv2, mid, b);
-                               
-               if (linecrosses(vv1, vv2, p, outv)) lcount += 1;
-
-               i += 1;
-       }
-       if ((lcount % 2) == 0) return 0;
+               //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
  *
@@ -453,7 +473,13 @@ static int goodline(float (*projectverts)[3], int v1i, int v2i, int nvert, float
  *
 */
 
-static BMLoop *find_ear(BMesh *bm, BMFace *f, float (*verts)[3], int nvert, float *outv)
+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;
@@ -470,17 +496,19 @@ static BMLoop *find_ear(BMesh *bm, BMFace *f, float (*verts)[3], int nvert, floa
 
                if (BM_Edge_Exist(v1, v3)) isear = 0;
 
-               if (isear && !goodline(verts, v1->head.eflag2, v3->head.eflag2, nvert, outv))
+               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 || ABS(angle-40.0f) < bestangle){
+                       if(!bestear || ABS(angle-45.0f) < bestangle){
                                bestear = l;
-                               bestangle = ABS(40.0f-angle);
+                               bestangle = ABS(45.0f-angle);
                        }
                        
-                       if ((angle > 10 && angle < 140) || i > 5) break;
+                       if (angle > 20 && angle < 90) break;
+                       if (angle < 100 && i > 5) break;
                        i += 1;
                }
                l = (BMLoop*)(l->head.next);
@@ -506,11 +534,13 @@ static BMLoop *find_ear(BMesh *bm, BMFace *f, float (*verts)[3], int nvert, floa
  * -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, nvert;
+       float no[3];
        BMLoop *l, *newl, *nextloop;
-       float outv[3] = {-1.0e30f, -1.0e30f, -1.0e30f};
+       BMVert *v;
 
        /*copy vertex coordinates to vertspace array*/
        i = 0;
@@ -519,30 +549,38 @@ void BM_Triangulate_Face(BMesh *bm, BMFace *f, float (*projectverts)[3], int new
                VECCOPY(projectverts[i], l->v->co);
                l->v->head.eflag2 = i; /*warning, abuse! never duplicate in tools code! never you hear?*/ /*actually, get rid of this completely, use a new structure for this....*/
                i++;
-
-               outv[0] = MAX2(outv[0], l->v->co[0])+1.0f;
-               outv[1] = MAX2(outv[1], l->v->co[1])+1.0f;
-               outv[2] = MAX2(outv[2], l->v->co[2])+1.0f;
-
                l = (BMLoop*)(l->head.next);
        }while(l != f->loopbase);
        
-       
-       //bmesh_update_face_normal(bm, f, projectverts);
+       ///bmesh_update_face_normal(bm, f, projectverts);
+
+       /*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);       
-       compute_poly_plane(projectverts, i);
+       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(bm, f, projectverts, nvert, outv);
+               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;
@@ -550,6 +588,15 @@ void BM_Triangulate_Face(BMesh *bm, BMFace *f, float (*projectverts)[3], int new
 
                        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);*/
                }
        }
 
@@ -557,7 +604,8 @@ 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));
-                       f = 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;