svn merge ^/trunk/blender -r42333:42361
authorCampbell Barton <ideasman42@gmail.com>
Fri, 2 Dec 2011 23:02:21 +0000 (23:02 +0000)
committerCampbell Barton <ideasman42@gmail.com>
Fri, 2 Dec 2011 23:02:21 +0000 (23:02 +0000)
1  2 
source/blender/blenlib/BLI_math_geom.h
source/blender/blenlib/intern/math_geom.c
source/blender/blenlib/intern/scanfill.c
source/blender/bmesh/intern/bmesh_interp.c
source/blender/bmesh/intern/bmesh_polygon.c
source/blender/editors/armature/meshlaplacian.c
source/blender/editors/object/object_relations.c
source/blender/editors/uvedit/uvedit_unwrap_ops.c
source/blender/render/intern/source/render_texture.c

index 65aaa22,0000000..00e51fb
mode 100644,000000..100644
--- /dev/null
@@@ -1,959 -1,0 +1,955 @@@
-       int i, xn, yn, zn, ax, ay;
 +/*
 + * ***** BEGIN GPL LICENSE BLOCK *****
 + *
 + * This program is free software; you can redistribute it and/or
 + * modify it under the terms of the GNU General Public License
 + * as published by the Free Software Foundation; either version 2
 + * of the License, or (at your option) any later version.
 + *
 + * This program is distributed in the hope that it will be useful,
 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 + * GNU General Public License for more details.
 + *
 + * You should have received a copy of the GNU General Public License
 + * along with this program; if not, write to the Free Software Foundation,
 + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 + *
 + * The Original Code is Copyright (C) 2007 Blender Foundation.
 + * All rights reserved.
 + *
 + * The Original Code is: all of this file.
 + *
 + * Contributor(s): Geoffrey Bantle.
 + *
 + * ***** END GPL LICENSE BLOCK *****
 + */
 +
 +/** \file blender/bmesh/intern/bmesh_interp.c
 + *  \ingroup bmesh
 + *
 + * Functions for interpolating data across the surface of a mesh.
 + */
 +
 +#include "MEM_guardedalloc.h"
 +
 +#include "DNA_mesh_types.h"
 +#include "DNA_meshdata_types.h"
 +
 +#include "BKE_customdata.h" 
 +#include "BKE_utildefines.h"
 +#include "BKE_multires.h"
 +
 +#include "BLI_array.h"
 +#include "BLI_math.h"
 +#include "BLI_cellalloc.h"
 +
 +#include "bmesh.h"
 +#include "bmesh_private.h"
 +
 +/**
 + *                    bmesh_data_interp_from_verts
 + *
 + *  Interpolates per-vertex data from two sources to a target.
 + * 
 + *  Returns -
 + *    Nothing
 + */
 +void BM_Data_Interp_From_Verts(BMesh *bm, BMVert *v1, BMVert *v2, BMVert *v, float fac)
 +{
 +      void *src[2];
 +      float w[2];
 +      if (v1->head.data && v2->head.data) {
 +              src[0]= v1->head.data;
 +              src[1]= v2->head.data;
 +              w[0] = 1.0f-fac;
 +              w[1] = fac;
 +              CustomData_bmesh_interp(&bm->vdata, src, w, NULL, 2, v->head.data);
 +      }
 +}
 +
 +/*
 +    BM Data Vert Average
 +
 +    Sets all the customdata (e.g. vert, loop) associated with a vert
 +    to the average of the face regions surrounding it.
 +*/
 +
 +
 +static void UNUSED_FUNCTION(BM_Data_Vert_Average)(BMesh *UNUSED(bm), BMFace *UNUSED(f))
 +{
 +      // BMIter iter;
 +}
 +
 +/**
 + *                    bmesh_data_facevert_edgeinterp
 + *
 + *  Walks around the faces of an edge and interpolates the per-face-edge
 + *  data between two sources to a target.
 + * 
 + *  Returns -
 + *    Nothing
 +*/
 + 
 +void BM_Data_Facevert_Edgeinterp(BMesh *bm, BMVert *v1, BMVert *UNUSED(v2), BMVert *v, BMEdge *e1, float fac)
 +{
 +      void *src[2];
 +      float w[2];
 +      BMLoop *l=NULL, *v1loop = NULL, *vloop = NULL, *v2loop = NULL;
 +      
 +      w[1] = 1.0f - fac;
 +      w[0] = fac;
 +
 +      if(!e1->l) return;
 +      l = e1->l;
 +      do{
 +              if(l->v == v1){ 
 +                      v1loop = l;
 +                      vloop = v1loop->next;
 +                      v2loop = vloop->next;
 +              }else if(l->v == v){
 +                      v1loop = l->next;
 +                      vloop = l;
 +                      v2loop = l->prev;
 +              }
 +              
 +              if (!v1loop || !v2loop)
 +                      return;
 +              
 +              src[0] = v1loop->head.data;
 +              src[1] = v2loop->head.data;                                     
 +
 +              CustomData_bmesh_interp(&bm->ldata, src,w, NULL, 2, vloop->head.data);                          
 +              l = l->radial_next;
 +      }while(l!=e1->l);
 +}
 +
 +void BM_loops_to_corners(BMesh *bm, Mesh *me, int findex,
 +                         BMFace *f, int numTex, int numCol) 
 +{
 +      BMLoop *l;
 +      BMIter iter;
 +      MTFace *texface;
 +      MTexPoly *texpoly;
 +      MCol *mcol;
 +      MLoopCol *mloopcol;
 +      MLoopUV *mloopuv;
 +      int i, j;
 +
 +      for(i=0; i < numTex; i++){
 +              texface = CustomData_get_n(&me->fdata, CD_MTFACE, findex, i);
 +              texpoly = CustomData_bmesh_get_n(&bm->pdata, f->head.data, CD_MTEXPOLY, i);
 +              
 +              texface->tpage = texpoly->tpage;
 +              texface->flag = texpoly->flag;
 +              texface->transp = texpoly->transp;
 +              texface->mode = texpoly->mode;
 +              texface->tile = texpoly->tile;
 +              texface->unwrap = texpoly->unwrap;
 +
 +              j = 0;
 +              BM_ITER(l, &iter, bm, BM_LOOPS_OF_FACE, f) {
 +                      mloopuv = CustomData_bmesh_get_n(&bm->ldata, l->head.data, CD_MLOOPUV, i);
 +                      texface->uv[j][0] = mloopuv->uv[0];
 +                      texface->uv[j][1] = mloopuv->uv[1];
 +
 +                      j++;
 +              }
 +
 +      }
 +
 +      for(i=0; i < numCol; i++){
 +              mcol = CustomData_get_n(&me->fdata, CD_MCOL, findex, i);
 +
 +              j = 0;
 +              BM_ITER(l, &iter, bm, BM_LOOPS_OF_FACE, f) {
 +                      mloopcol = CustomData_bmesh_get_n(&bm->ldata, l->head.data, CD_MLOOPCOL, i);
 +                      mcol[j].r = mloopcol->r;
 +                      mcol[j].g = mloopcol->g;
 +                      mcol[j].b = mloopcol->b;
 +                      mcol[j].a = mloopcol->a;
 +
 +                      j++;
 +              }
 +      }
 +}
 +
 +/**
 + *                    BM_data_interp_from_face
 + *
 + *  projects target onto source, and pulls interpolated customdata from
 + *  source.
 + * 
 + *  Returns -
 + *    Nothing
 +*/
 +void BM_face_interp_from_face(BMesh *bm, BMFace *target, BMFace *source)
 +{
 +      BMLoop *l1, *l2;
 +      void **blocks=NULL;
 +      float (*cos)[3]=NULL, *w=NULL;
 +      BLI_array_staticdeclare(cos,     BM_NGON_STACK_SIZE);
 +      BLI_array_staticdeclare(w,       BM_NGON_STACK_SIZE);
 +      BLI_array_staticdeclare(blocks,  BM_NGON_STACK_SIZE);
 +      
 +      BM_Copy_Attributes(bm, bm, source, target);
 +      
 +      l2 = bm_firstfaceloop(source);
 +      do {
 +              BLI_array_growone(cos);
 +              copy_v3_v3(cos[BLI_array_count(cos)-1], l2->v->co);
 +              BLI_array_growone(w);
 +              BLI_array_append(blocks, l2->head.data);
 +              l2 = l2->next;
 +      } while (l2 != bm_firstfaceloop(source));
 +
 +      l1 = bm_firstfaceloop(target);
 +      do {
 +              interp_weights_poly_v3(w, cos, source->len, l1->v->co);
 +              CustomData_bmesh_interp(&bm->ldata, blocks, w, NULL, BLI_array_count(blocks), l1->head.data);
 +              l1 = l1->next;
 +      } while (l1 != bm_firstfaceloop(target));
 +
 +      BLI_array_free(cos);
 +      BLI_array_free(w);
 +      BLI_array_free(blocks);
 +}
 +
 +/****some math stuff for dealing with doubles, put here to
 +  avoid merge errors - joeedh ****/
 +
 +#define VECMUL(a, b) (((a)[0] = (a)[0] * (b)), ((a)[1] = (a)[1] * (b)), ((a)[2] = (a)[2] * (b)))
 +#define VECADD2(a, b) (((a)[0] = (a)[0] + (b)[0]), ((a)[1] = (a)[1] + (b)[1]), ((a)[2] = (a)[2] + (b)[2]))
 +#define VECSUB2(a, b) (((a)[0] = (a)[0] - (b)[0]), ((a)[1] = (a)[1] - (b)[1]), ((a)[2] = (a)[2] - (b)[2]))
 +
 +/* find closest point to p on line through l1,l2 and return lambda,
 + * where (0 <= lambda <= 1) when cp is in the line segement l1,l2
 + */
 +static double closest_to_line_v3_d(double cp[3], const double p[3], const double l1[3], const double l2[3])
 +{
 +      double h[3],u[3],lambda;
 +      VECSUB(u, l2, l1);
 +      VECSUB(h, p, l1);
 +      lambda =INPR(u,h)/INPR(u,u);
 +      cp[0] = l1[0] + u[0] * lambda;
 +      cp[1] = l1[1] + u[1] * lambda;
 +      cp[2] = l1[2] + u[2] * lambda;
 +      return lambda;
 +}
 +
 +/* point closest to v1 on line v2-v3 in 3D */
 +static void UNUSED_FUNCTION(closest_to_line_segment_v3_d)(double *closest, double v1[3], double v2[3], double v3[3])
 +{
 +      double lambda, cp[3];
 +
 +      lambda= closest_to_line_v3_d(cp,v1, v2, v3);
 +
 +      if(lambda <= 0.0) {
 +              VECCOPY(closest, v2);
 +      } else if(lambda >= 1.0) {
 +              VECCOPY(closest, v3);
 +      } else {
 +              VECCOPY(closest, cp);
 +      }
 +}
 +
 +static double UNUSED_FUNCTION(len_v3_d)(const double a[3])
 +{
 +      return sqrt(INPR(a, a));
 +}
 +
 +static double UNUSED_FUNCTION(len_v3v3_d)(const double a[3], const double b[3])
 +{
 +      double d[3];
 +
 +      VECSUB(d, b, a);
 +      return sqrt(INPR(d, d));
 +}
 +
 +static void cent_quad_v3_d(double *cent, double *v1, double *v2, double *v3, double *v4)
 +{
 +      cent[0]= 0.25*(v1[0]+v2[0]+v3[0]+v4[0]);
 +      cent[1]= 0.25*(v1[1]+v2[1]+v3[1]+v4[1]);
 +      cent[2]= 0.25*(v1[2]+v2[2]+v3[2]+v4[2]);
 +}
 +
 +static void UNUSED_FUNCTION(cent_tri_v3_d)(double *cent, double *v1, double *v2, double *v3)
 +{
 +      cent[0]= 0.33333*(v1[0]+v2[0]+v3[0]);
 +      cent[1]= 0.33333*(v1[1]+v2[1]+v3[1]);
 +      cent[2]= 0.33333*(v1[2]+v2[2]+v3[2]);
 +}
 +
 +static void UNUSED_FUNCTION(cross_v3_v3v3_d)(double r[3], const double a[3], const double b[3])
 +{
 +      r[0]= a[1]*b[2] - a[2]*b[1];
 +      r[1]= a[2]*b[0] - a[0]*b[2];
 +      r[2]= a[0]*b[1] - a[1]*b[0];
 +}
 +
 +/* distance v1 to line-piece v2-v3 */
 +static double UNUSED_FUNCTION(dist_to_line_segment_v2_d)(double v1[3], double v2[3], double v3[3])
 +{
 +      double labda, rc[2], pt[2], len;
 +      
 +      rc[0]= v3[0]-v2[0];
 +      rc[1]= v3[1]-v2[1];
 +      len= rc[0]*rc[0]+ rc[1]*rc[1];
 +      if(len==0.0) {
 +              rc[0]= v1[0]-v2[0];
 +              rc[1]= v1[1]-v2[1];
 +              return sqrt(rc[0]*rc[0]+ rc[1]*rc[1]);
 +      }
 +      
 +      labda= (rc[0]*(v1[0]-v2[0]) + rc[1]*(v1[1]-v2[1]))/len;
 +      if(labda<=0.0) {
 +              pt[0]= v2[0];
 +              pt[1]= v2[1];
 +      }
 +      else if(labda>=1.0) {
 +              pt[0]= v3[0];
 +              pt[1]= v3[1];
 +      }
 +      else {
 +              pt[0]= labda*rc[0]+v2[0];
 +              pt[1]= labda*rc[1]+v2[1];
 +      }
 +
 +      rc[0]= pt[0]-v1[0];
 +      rc[1]= pt[1]-v1[1];
 +      return sqrt(rc[0]*rc[0]+ rc[1]*rc[1]);
 +}
 +
 +
 +MINLINE double line_point_side_v2_d(const double *l1, const double *l2, const double *pt)
 +{
 +      return  ((l1[0]-pt[0]) * (l2[1]-pt[1])) -
 +                      ((l2[0]-pt[0]) * (l1[1]-pt[1]));
 +}
 +
 +/* point in quad - only convex quads */
 +static int isect_point_quad_v2_d(double pt[2], double v1[2], double v2[2], double v3[2], double v4[2])
 +{
 +      if (line_point_side_v2_d(v1,v2,pt)>=0.0) {
 +              if (line_point_side_v2_d(v2,v3,pt)>=0.0) {
 +                      if (line_point_side_v2_d(v3,v4,pt)>=0.0) {
 +                              if (line_point_side_v2_d(v4,v1,pt)>=0.0) {
 +                                      return 1;
 +                              }
 +                      }
 +              }
 +      } else {
 +              if (! (line_point_side_v2_d(v2,v3,pt)>=0.0)) {
 +                      if (! (line_point_side_v2_d(v3,v4,pt)>=0.0)) {
 +                              if (! (line_point_side_v2_d(v4,v1,pt)>=0.0)) {
 +                                      return -1;
 +                              }
 +                      }
 +              }
 +      }
 +      
 +      return 0;
 +}
 +
 +/***** multires interpolation*****
 +
 +mdisps is a grid of displacements, ordered thus:
 +
 + v1/center----v4/next -> x
 +     |           |
 +     |           |
 +  v2/prev------v3/cur
 +     |
 +     V
 +     y
 +*/
 +
 +static int compute_mdisp_quad(BMLoop *l, double v1[3], double v2[3], double v3[3], double v4[3], double e1[3], double e2[3])
 +{
 +      double cent[3] = {0.0, 0.0, 0.0}, n[3], p[3];
 +      BMLoop *l2;
 +      
 +      /*computer center*/
 +      l2 = bm_firstfaceloop(l->f);
 +      do {
 +              cent[0] += (double)l2->v->co[0];
 +              cent[1] += (double)l2->v->co[1];
 +              cent[2] += (double)l2->v->co[2];
 +              l2 = l2->next;
 +      } while (l2 != bm_firstfaceloop(l->f));
 +      
 +      VECMUL(cent, (1.0/(double)l->f->len));
 +      
 +      VECADD(p, l->prev->v->co, l->v->co);
 +      VECMUL(p, 0.5);
 +      VECADD(n, l->next->v->co, l->v->co);
 +      VECMUL(n, 0.5);
 +      
 +      VECCOPY(v1, cent);
 +      VECCOPY(v2, p);
 +      VECCOPY(v3, l->v->co);
 +      VECCOPY(v4, n);
 +      
 +      VECSUB(e1, v2, v1);
 +      VECSUB(e2, v3, v4);
 +      
 +      return 1;
 +}
 +
 +/*funnily enough, I think this is identical to face_to_crn_interp, heh*/
 +static double quad_coord(double aa[3], double bb[3], double cc[3], double dd[3], int a1, int a2)
 +{
 +      double x, y, z, f1;
 +      
 +      x = aa[a1]*cc[a2]-cc[a1]*aa[a2];
 +      y = aa[a1]*dd[a2]+bb[a1]*cc[a2]-cc[a1]*bb[a2]-dd[a1]*aa[a2];
 +      z = bb[a1]*dd[a2]-dd[a1]*bb[a2];
 +      
 +      if (fabs(2*(x-y+z)) > DBL_EPSILON*10.0) {
 +              double f2;
 +
 +              f1 = (sqrt(y*y-4.0*x*z) - y + 2.0*z) / (2.0*(x-y+z));
 +              f2 = (-sqrt(y*y-4.0*x*z) - y + 2.0*z) / (2.0*(x-y+z));
 +
 +              f1= fabs(f1);
 +              f2= fabs(f2);
 +              f1 = MIN2(f1, f2);
 +              CLAMP(f1, 0.0, 1.0+DBL_EPSILON);
 +      } else {
 +              f1 = -z/(y - 2*z);
 +              CLAMP(f1, 0.0, 1.0+DBL_EPSILON);
 +              
 +              if (isnan(f1) || f1 > 1.0 || f1 < 0.0) {
 +                      int i;
 +                      
 +                      for (i=0; i<2; i++) {
 +                              if (fabsf(aa[i]) < FLT_EPSILON*100)
 +                                      return aa[(i+1)%2] / fabs(bb[(i+1)%2] - aa[(i+1)%2]);
 +                              if (fabsf(cc[i]) < FLT_EPSILON*100)
 +                                      return cc[(i+1)%2] / fabs(dd[(i+1)%2] - cc[(i+1)%2]);
 +                      }
 +              }
 +      }
 +
 +      return f1;
 +}
 +
 +static int quad_co(double *x, double *y, double v1[3], double v2[3], double v3[3], double v4[3], double p[3], float n[3])
 +{
 +      float projverts[5][3], n2[3];
 +      double dprojverts[4][3], origin[3]={0.0f, 0.0f, 0.0f};
 +      int i;
 +
 +      /*project points into 2d along normal*/
 +      VECCOPY(projverts[0], v1);
 +      VECCOPY(projverts[1], v2);
 +      VECCOPY(projverts[2], v3);
 +      VECCOPY(projverts[3], v4);
 +      VECCOPY(projverts[4], p);
 +
 +      normal_quad_v3(n2, projverts[0], projverts[1], projverts[2], projverts[3]);
 +
 +      if (INPR(n, n2) < -FLT_EPSILON)
 +              return 0;
 +
 +      /*rotate*/      
 +      poly_rotate_plane(n, projverts, 5);
 +      
 +      /*flatten*/
 +      for (i=0; i<5; i++) projverts[i][2] = 0.0f;
 +      
 +      /*subtract origin*/
 +      for (i=0; i<4; i++) {
 +              VECSUB2(projverts[i], projverts[4]);
 +      }
 +      
 +      VECCOPY(dprojverts[0], projverts[0]);
 +      VECCOPY(dprojverts[1], projverts[1]);
 +      VECCOPY(dprojverts[2], projverts[2]);
 +      VECCOPY(dprojverts[3], projverts[3]);
 +
 +      if (!isect_point_quad_v2_d(origin, dprojverts[0], dprojverts[1], dprojverts[2], dprojverts[3]))
 +              return 0;
 +      
 +      *y = quad_coord(dprojverts[1], dprojverts[0], dprojverts[2], dprojverts[3], 0, 1);
 +      *x = quad_coord(dprojverts[2], dprojverts[1], dprojverts[3], dprojverts[0], 0, 1);
 +
 +      return 1;
 +}
 +
 +
 +/*tl is loop to project onto, l is loop whose internal displacement, co, is being
 +  projected.  x and y are location in loop's mdisps grid of point co.*/
 +static int mdisp_in_mdispquad(BMesh *bm, BMLoop *l, BMLoop *tl, double p[3], double *x, double *y, int res)
 +{
 +      double v1[3], v2[3], c[3], v3[3], v4[3], e1[3], e2[3];
 +      double eps = FLT_EPSILON*4000;
 +      
 +      if (len_v3(l->v->no) == 0.0f)
 +              BM_Vert_UpdateAllNormals(bm, l->v);
 +      if (len_v3(tl->v->no) == 0.0f)
 +              BM_Vert_UpdateAllNormals(bm, tl->v);
 +              
 +      compute_mdisp_quad(tl, v1, v2, v3, v4, e1, e2);
 +
 +      /*expand quad a bit*/
 +      cent_quad_v3_d(c, v1, v2, v3, v4);
 +      
 +      VECSUB2(v1, c); VECSUB2(v2, c);
 +      VECSUB2(v3, c); VECSUB2(v4, c);
 +      VECMUL(v1, 1.0+eps); VECMUL(v2, 1.0+eps);
 +      VECMUL(v3, 1.0+eps); VECMUL(v4, 1.0+eps);
 +      VECADD2(v1, c); VECADD2(v2, c);
 +      VECADD2(v3, c); VECADD2(v4, c);
 +      
 +      if (!quad_co(x, y, v1, v2, v3, v4, p, l->v->no))
 +              return 0;
 +      
 +      *x *= res-1;
 +      *y *= res-1;
 +                
 +      return 1;
 +}
 +
 +static void bmesh_loop_interp_mdisps(BMesh *bm, BMLoop *target, BMFace *source)
 +{
 +      MDisps *mdisps;
 +      BMLoop *l2;
 +      double x, y, d, v1[3], v2[3], v3[3], v4[3] = {0.0f, 0.0f, 0.0f}, e1[3], e2[3];
 +      int ix, iy, res;
 +      
 +      /*ignore 2-edged faces*/
 +      if (target->f->len < 3)
 +              return;
 +      
 +      if (!CustomData_has_layer(&bm->ldata, CD_MDISPS))
 +              return;
 +      
 +      mdisps = CustomData_bmesh_get(&bm->ldata, target->head.data, CD_MDISPS);
 +      compute_mdisp_quad(target, v1, v2, v3, v4, e1, e2);
 +      
 +      /*if no disps data allocate a new grid, the size of the first grid in source.*/
 +      if (!mdisps->totdisp) {
 +              MDisps *md2 = CustomData_bmesh_get(&bm->ldata, bm_firstfaceloop(source)->head.data, CD_MDISPS);
 +              
 +              mdisps->totdisp = md2->totdisp;
 +              if (mdisps->totdisp)
 +                      mdisps->disps = BLI_cellalloc_calloc(sizeof(float)*3*mdisps->totdisp, "mdisp->disps in bmesh_loop_intern_mdisps");
 +              else 
 +                      return;
 +      }
 +      
 +      res = (int)sqrt(mdisps->totdisp);
 +      d = 1.0/(double)(res-1);
 +      for (x=0.0f, ix=0; ix<res; x += d, ix++) {
 +              for (y=0.0f, iy=0; iy<res; y+= d, iy++) {
 +                      double co1[3], co2[3], co[3];
 +                      /* double xx, yy; */ /* UNUSED */
 +                      
 +                      VECCOPY(co1, e1);
 +                      
 +                      /* if (!iy) yy = y + (double)FLT_EPSILON*20; */
 +                      /* else yy = y - (double)FLT_EPSILON*20; */
 +                      
 +                      VECMUL(co1, y);
 +                      VECADD2(co1, v1);
 +                      
 +                      VECCOPY(co2, e2);
 +                      VECMUL(co2, y);
 +                      VECADD2(co2, v4);
 +                      
 +                      /* if (!ix) xx = x + (double)FLT_EPSILON*20; */ /* UNUSED */
 +                      /* else xx = x - (double)FLT_EPSILON*20; */ /* UNUSED */
 +                      
 +                      VECSUB(co, co2, co1);
 +                      VECMUL(co, x);
 +                      VECADD2(co, co1);
 +                      
 +                      l2 = bm_firstfaceloop(source);
 +                      do {
 +                              double x2, y2;
 +                              MDisps *md1, *md2;
 +
 +                              md1 = CustomData_bmesh_get(&bm->ldata, target->head.data, CD_MDISPS);
 +                              md2 = CustomData_bmesh_get(&bm->ldata, l2->head.data, CD_MDISPS);
 +                              
 +                              if (mdisp_in_mdispquad(bm, target, l2, co, &x2, &y2, res)) {
 +                                      /* int ix2 = (int)x2; */ /* UNUSED */
 +                                      /* int iy2 = (int)y2; */ /* UNUSED */
 +                                      
 +                                      old_mdisps_bilinear(md1->disps[iy*res+ix], md2->disps, res, (float)x2, (float)y2);
 +                              }
 +                              l2 = l2->next;
 +                      } while (l2 != bm_firstfaceloop(source));
 +              }
 +      }
 +}
 +
 +void BM_multires_smooth_bounds(BMesh *bm, BMFace *f)
 +{
 +      BMLoop *l;
 +      BMIter liter;
 +      
 +      if (!CustomData_has_layer(&bm->ldata, CD_MDISPS))
 +              return;
 +      
 +      BM_ITER(l, &liter, bm, BM_LOOPS_OF_FACE, f) {
 +              MDisps *mdp = CustomData_bmesh_get(&bm->ldata, l->prev->head.data, CD_MDISPS);
 +              MDisps *mdl = CustomData_bmesh_get(&bm->ldata, l->head.data, CD_MDISPS);
 +              MDisps *mdn = CustomData_bmesh_get(&bm->ldata, l->next->head.data, CD_MDISPS);
 +              float co1[3];
 +              int sides;
 +              int y;
 +              
 +              /*****
 +              mdisps is a grid of displacements, ordered thus:
 +              
 +                                 v4/next
 +                                   |                          
 +               |      v1/cent-----mid2 ---> x
 +               |         |         | 
 +               |         |         |
 +              v2/prev---mid1-----v3/cur
 +                         |
 +                         V
 +                         y
 +              *****/
 +                
 +              sides = (int)sqrt(mdp->totdisp);
 +              for (y=0; y<sides; y++) {
 +                      add_v3_v3v3(co1, mdn->disps[y*sides], mdl->disps[y]);
 +                      mul_v3_fl(co1, 0.5);
 +
 +                      copy_v3_v3(mdn->disps[y*sides], co1);
 +                      copy_v3_v3(mdl->disps[y], co1);
 +              }
 +      }
 +      
 +      BM_ITER(l, &liter, bm, BM_LOOPS_OF_FACE, f) {
 +              MDisps *mdl1 = CustomData_bmesh_get(&bm->ldata, l->head.data, CD_MDISPS);
 +              MDisps *mdl2;
 +              float co1[3], co2[3], co[3];
 +              int sides;
 +              int y;
 +              
 +              /*****
 +              mdisps is a grid of displacements, ordered thus:
 +              
 +                                 v4/next
 +                                   |                          
 +               |      v1/cent-----mid2 ---> x
 +               |         |         | 
 +               |         |         |
 +              v2/prev---mid1-----v3/cur
 +                         |
 +                         V
 +                         y
 +              *****/
 +               
 +              if (l->radial_next == l)
 +                      continue;
 +
 +              if (l->radial_next->v == l->v)
 +                      mdl2 = CustomData_bmesh_get(&bm->ldata, l->radial_next->head.data, CD_MDISPS);
 +              else
 +                      mdl2 = CustomData_bmesh_get(&bm->ldata, l->radial_next->next->head.data, CD_MDISPS);
 +                      
 +              sides = (int)sqrt(mdl1->totdisp);
 +              for (y=0; y<sides; y++) {
 +                      int a1, a2, o1, o2;
 +                      
 +                      if (l->v != l->radial_next->v) {
 +                              a1 = sides*y + sides-2;
 +                              a2 = (sides-2)*sides + y;
 +                              
 +                              o1 = sides*y + sides-1;
 +                              o2 = (sides-1)*sides + y;
 +                      } else {
 +                              a1 = sides*y + sides-2;
 +                              a2 = sides*y + sides-2;
 +                              o1 = sides*y + sides-1;
 +                              o2 = sides*y + sides-1;
 +                      }
 +                      
 +                      /*magic blending numbers, hardcoded!*/
 +                      add_v3_v3v3(co1, mdl1->disps[a1], mdl2->disps[a2]);
 +                      mul_v3_fl(co1, 0.18);
 +                      
 +                      add_v3_v3v3(co2, mdl1->disps[o1], mdl2->disps[o2]);
 +                      mul_v3_fl(co2, 0.32);
 +                      
 +                      add_v3_v3v3(co, co1, co2);
 +                      
 +                      copy_v3_v3(mdl1->disps[o1], co);
 +                      copy_v3_v3(mdl2->disps[o2], co); //*/
 +              }
 +      }
 +}
 +
 +void BM_loop_interp_multires(BMesh *bm, BMLoop *target, BMFace *source)
 +{
 +      bmesh_loop_interp_mdisps(bm, target, source);
 +}
 +
 +void BM_loop_interp_from_face(BMesh *bm, BMLoop *target, BMFace *source, 
 +                              int do_vertex, int do_multires)
 +{
 +      BMLoop *l;
 +      void **blocks=NULL;
 +      void **vblocks=NULL;
 +      float (*cos)[3]=NULL, co[3], *w=NULL, cent[3] = {0.0f, 0.0f, 0.0f};
 +      BLI_array_staticdeclare(cos,      BM_NGON_STACK_SIZE);
 +      BLI_array_staticdeclare(w,        BM_NGON_STACK_SIZE);
 +      BLI_array_staticdeclare(blocks,   BM_NGON_STACK_SIZE);
 +      BLI_array_staticdeclare(vblocks,  BM_NGON_STACK_SIZE);
-          the 2d projected coords are the same and faster to compute */
-       xn= fabsf(source->no[0]);
-       yn= fabsf(source->no[1]);
-       zn= fabsf(source->no[2]);
-       if (zn >= xn && zn >= yn)       { ax= 0; ay= 1; }
-       else if (yn >= xn && yn >= zn)  { ax= 0; ay= 2; }
-       else                            { ax= 1; ay= 2; }
-       
++      int i, ax, ay;
 +      
 +      BM_Copy_Attributes(bm, bm, source, target->f);
 +      
 +      l = bm_firstfaceloop(source);
 +      do {
 +              BLI_array_growone(cos);
 +              copy_v3_v3(cos[BLI_array_count(cos)-1], l->v->co);
 +              add_v3_v3(cent, cos[BLI_array_count(cos)-1]);
 +              
 +              BLI_array_append(w, 0.0f);
 +              BLI_array_append(blocks, l->head.data);
 +      
 +              if (do_vertex)
 +                      BLI_array_append(vblocks, l->v->head.data);
 +      
 +              l = l->next;
 +      } while (l != bm_firstfaceloop(source));
 +
 +      /* find best projection of face XY, XZ or YZ: barycentric weights of
++       * the 2d projected coords are the same and faster to compute */
++
++      axis_dominant_v3(&ax, &ay, source->no);
++
 +      /* scale source face coordinates a bit, so points sitting directonly on an
 +         edge will work.*/
 +      mul_v3_fl(cent, 1.0f/(float)source->len);
 +      for (i=0; i<source->len; i++) {
 +              float vec[3], tmp[3];
 +              sub_v3_v3v3(vec, cent, cos[i]);
 +              mul_v3_fl(vec, 0.001);
 +              add_v3_v3(cos[i], vec);
 +              
 +              copy_v3_v3(tmp, cos[i]);
 +              cos[i][0] = tmp[ax];
 +              cos[i][1] = tmp[ay];
 +              cos[i][2] = 0.0;
 +      }
 +      
 +      
 +      /*interpolate*/
 +      co[0] = target->v->co[ax];
 +      co[1] = target->v->co[ay];
 +      co[2] = 0.0f;
 +      
 +      interp_weights_poly_v3(w, cos, source->len, co);
 +      CustomData_bmesh_interp(&bm->ldata, blocks, w, NULL, source->len, target->head.data);
 +      if (do_vertex) {
 +              CustomData_bmesh_interp(&bm->vdata, vblocks, w, NULL, source->len, target->v->head.data);
 +              BLI_array_free(vblocks);
 +      }
 +      
 +      BLI_array_free(cos);
 +      BLI_array_free(w);
 +      BLI_array_free(blocks);
 +      
 +      if (do_multires) {
 +              if (CustomData_has_layer(&bm->ldata, CD_MDISPS)) {
 +                      bmesh_loop_interp_mdisps(bm, target, source);
 +              }
 +      }
 +}
 +
 +
 +void BM_vert_interp_from_face(BMesh *bm, BMVert *v, BMFace *source)
 +{
 +      BMLoop *l;
 +      void **blocks=NULL;
 +      float (*cos)[3]=NULL, *w=NULL, cent[3] = {0.0f, 0.0f, 0.0f};
 +      BLI_array_staticdeclare(cos,     BM_NGON_STACK_SIZE);
 +      BLI_array_staticdeclare(w,       BM_NGON_STACK_SIZE);
 +      BLI_array_staticdeclare(blocks,  BM_NGON_STACK_SIZE);
 +      int i;
 +      
 +      l = bm_firstfaceloop(source);
 +      do {
 +              BLI_array_growone(cos);
 +              copy_v3_v3(cos[BLI_array_count(cos)-1], l->v->co);
 +              add_v3_v3(cent, cos[BLI_array_count(cos)-1]);
 +              
 +              BLI_array_append(w, 0.0f);
 +              BLI_array_append(blocks, l->v->head.data);
 +              l = l->next;
 +      } while (l != bm_firstfaceloop(source));
 +
 +      /*scale source face coordinates a bit, so points sitting directonly on an
 +      edge will work.*/
 +      mul_v3_fl(cent, 1.0f/(float)source->len);
 +      for (i=0; i<source->len; i++) {
 +              float vec[3];
 +              sub_v3_v3v3(vec, cent, cos[i]);
 +              mul_v3_fl(vec, 0.01);
 +              add_v3_v3(cos[i], vec);
 +      }
 +      
 +      /*interpolate*/
 +      interp_weights_poly_v3(w, cos, source->len, v->co);
 +      CustomData_bmesh_interp(&bm->vdata, blocks, w, NULL, source->len, v->head.data);
 +      
 +      BLI_array_free(cos);
 +      BLI_array_free(w);
 +      BLI_array_free(blocks);
 +}
 +
 +static void update_data_blocks(BMesh *bm, CustomData *olddata, CustomData *data)
 +{
 +      BMIter iter;
 +      BLI_mempool *oldpool = olddata->pool;
 +      void *block;
 +
 +      CustomData_bmesh_init_pool(data, data==&bm->ldata ? 2048 : 512);
 +
 +      if (data == &bm->vdata) {
 +              BMVert *eve;
 +              
 +              BM_ITER(eve, &iter, bm, BM_VERTS_OF_MESH, NULL) {
 +                      block = NULL;
 +                      CustomData_bmesh_set_default(data, &block);
 +                      CustomData_bmesh_copy_data(olddata, data, eve->head.data, &block);
 +                      CustomData_bmesh_free_block(olddata, &eve->head.data);
 +                      eve->head.data= block;
 +              }
 +      }
 +      else if (data == &bm->edata) {
 +              BMEdge *eed;
 +
 +              BM_ITER(eed, &iter, bm, BM_EDGES_OF_MESH, NULL) {
 +                      block = NULL;
 +                      CustomData_bmesh_set_default(data, &block);
 +                      CustomData_bmesh_copy_data(olddata, data, eed->head.data, &block);
 +                      CustomData_bmesh_free_block(olddata, &eed->head.data);
 +                      eed->head.data= block;
 +              }
 +      }
 +      else if (data == &bm->pdata || data == &bm->ldata) {
 +              BMIter liter;
 +              BMFace *efa;
 +              BMLoop *l;
 +
 +              BM_ITER(efa, &iter, bm, BM_FACES_OF_MESH, NULL) {
 +                      if (data == &bm->pdata) {
 +                              block = NULL;
 +                              CustomData_bmesh_set_default(data, &block);
 +                              CustomData_bmesh_copy_data(olddata, data, efa->head.data, &block);
 +                              CustomData_bmesh_free_block(olddata, &efa->head.data);
 +                              efa->head.data= block;
 +                      }
 +
 +                      if (data == &bm->ldata) {
 +                              BM_ITER(l, &liter, bm, BM_LOOPS_OF_FACE, efa) {
 +                                      block = NULL;
 +                                      CustomData_bmesh_set_default(data, &block);
 +                                      CustomData_bmesh_copy_data(olddata, data, l->head.data, &block);
 +                                      CustomData_bmesh_free_block(olddata, &l->head.data);
 +                                      l->head.data= block;
 +                              }
 +                      }
 +              }
 +      }
 +
 +      if (oldpool) {
 +              /* this should never happen but can when dissolve fails - [#28960] */
 +              BLI_assert(data->pool != oldpool);
 +
 +              BLI_mempool_destroy(oldpool);
 +      }
 +}
 +
 +void BM_add_data_layer(BMesh *bm, CustomData *data, int type)
 +{
 +      CustomData olddata;
 +
 +      olddata= *data;
 +      olddata.layers= (olddata.layers)? MEM_dupallocN(olddata.layers): NULL;
 +
 +      /* the pool is now owned by olddata and must not be shared */
 +      data->pool = NULL;
 +
 +      CustomData_add_layer(data, type, CD_DEFAULT, NULL, 0);
 +
 +      update_data_blocks(bm, &olddata, data);
 +      if (olddata.layers) MEM_freeN(olddata.layers);
 +}
 +
 +void BM_add_data_layer_named(BMesh *bm, CustomData *data, int type, const char *name)
 +{
 +      CustomData olddata;
 +
 +      olddata= *data;
 +      olddata.layers= (olddata.layers)? MEM_dupallocN(olddata.layers): NULL;
 +
 +      /* the pool is now owned by olddata and must not be shared */
 +      data->pool = NULL;
 +
 +      CustomData_add_layer_named(data, type, CD_DEFAULT, NULL, 0, name);
 +
 +      update_data_blocks(bm, &olddata, data);
 +      if (olddata.layers) MEM_freeN(olddata.layers);
 +}
 +
 +void BM_free_data_layer(BMesh *bm, CustomData *data, int type)
 +{
 +      CustomData olddata;
 +
 +      olddata= *data;
 +      olddata.layers= (olddata.layers)? MEM_dupallocN(olddata.layers): NULL;
 +
 +      /* the pool is now owned by olddata and must not be shared */
 +      data->pool = NULL;
 +
 +      CustomData_free_layer_active(data, type, 0);
 +
 +      update_data_blocks(bm, &olddata, data);
 +      if (olddata.layers) MEM_freeN(olddata.layers);
 +}
 +
 +void BM_free_data_layer_n(BMesh *bm, CustomData *data, int type, int n)
 +{
 +      CustomData olddata;
 +
 +      olddata= *data;
 +      olddata.layers= (olddata.layers)? MEM_dupallocN(olddata.layers): NULL;
 +
 +      /* the pool is now owned by olddata and must not be shared */
 +      data->pool = NULL;
 +
 +      CustomData_free_layer(data, type, 0, CustomData_get_layer_index_n(data, type, n));
 +      
 +      update_data_blocks(bm, &olddata, data);
 +      if (olddata.layers) MEM_freeN(olddata.layers);
 +}
 +
 +float BM_GetCDf(CustomData *cd, void *element, int type)
 +{
 +      if (CustomData_has_layer(cd, type)) {
 +              float *f = CustomData_bmesh_get(cd, ((BMHeader*)element)->data, type);
 +              return *f;
 +      }
 +
 +      return 0.0;
 +}
 +
 +void BM_SetCDf(CustomData *cd, void *element, int type, float val)
 +{
 +      if (CustomData_has_layer(cd, type)) {
 +              float *f = CustomData_bmesh_get(cd, ((BMHeader*)element)->data, type);
 +              *f = val;
 +      }
 +
 +      return;
 +}
index edbc152,0000000..9e65261
mode 100644,000000..100644
--- /dev/null
@@@ -1,991 -1,0 +1,985 @@@
-       int xn, yn, zn, ax, ay;
 +/*
 + * ***** 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 <string.h>
 +#include <math.h>
 +#include <stdlib.h>
 +
 +#include "BKE_utildefines.h"
 +
 +#include "BLI_math.h"
 +#include "BLI_blenlib.h"
 +#include "BLI_array.h"
 +#include "BLI_utildefines.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(double *v1, double *v2, double *v3)
 +/* 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 0;
 +      else if(inp==0) {
 +              if(v1[0]==v3[0] && v1[1]==v3[1]) return 0;
 +              if(v2[0]==v3[0] && v2[1]==v3[1]) return 0;
 +      }
 +      return 1;
 +}
 +
 +static short testedgesidef(float *v1, float *v2, float *v3)
 +/* 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 0;
 +      else if(inp==0) {
 +              if(v1[0]==v3[0] && v1[1]==v3[1]) return 0;
 +              if(v2[0]==v3[0] && v2[1]==v3[1]) return 0;
 +      }
 +      return 1;
 +}
 +
 +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;
 +      return 0;
 +}
 +
 +/*
 + * 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 error*/
 +      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 *area, float (*verts)[3], int nverts)
 +{
 +      int i, j;
 +      float atmp = 0.0, xtmp = 0.0, ytmp = 0.0, ai;
 +
 +      zero_v3(center);
 +
 +      if(nverts < 3) 
 +              return 0;
 +
 +      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(area)
 +              *area = atmp / 2.0f;    
 +      
 +      if (atmp != 0){
 +              center[0] = xtmp /  (3.0f * atmp);
 +              center[1] = xtmp /  (3.0f * atmp);
 +              return 1;
 +      }
 +      return 0;
 +}
 +
 +float BM_face_area(BMFace *f)
 +{
 +      BMLoop *l;
 +      BMIter iter;
 +      float (*verts)[3];
 +      float area, center[3];
 +      int i;
 +
 +      BLI_array_fixedstack_declare(verts, BM_NGON_STACK_SIZE, f->len, __func__);
 +
 +      i = 0;
 +      BM_ITER(l, &iter, NULL, 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_Compute_Face_CenterBounds(BMesh *bm, BMFace *f, float r_cent[3])
 +{
 +      BMIter iter;
 +      BMLoop *l;
 +      float min[3], max[3];
 +      int i;
 +
 +      INIT_MINMAX(min, max);
 +      l = BMIter_New(&iter, bm, BM_LOOPS_OF_FACE, f);
 +      for (i=0; l; l=BMIter_Step(&iter), i++) {
 +              DO_MINMAX(l->v->co, min, max);
 +      }
 +
 +      mid_v3_v3v3(r_cent, min, max);
 +}
 +
 +void BM_Compute_Face_CenterMean(BMesh *bm, BMFace *f, float r_cent[3])
 +{
 +      BMIter iter;
 +      BMLoop *l;
 +      int i;
 +
 +      zero_v3(r_cent);
 +
 +      l = BMIter_New(&iter, bm, BM_LOOPS_OF_FACE, f);
 +      for (i=0; l; l=BMIter_Step(&iter), i++) {
 +              add_v3_v3(r_cent, l->v->co);
 +      }
 +
 +      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, float *v2, 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_UpdateNormal(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, proj);
 +
 +              BLI_array_fixedstack_free(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 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_Compute_Normals */
 +                              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_UpdateAllNormals(BMesh *bm, BMVert *v)
 +{
 +      BMIter iter;
 +      BMFace *f;
 +      int len=0;
 +
 +      f = BMIter_New(&iter, bm, BM_FACES_OF_VERT, v);
 +      for (; f; f=BMIter_Step(&iter), len++) {
 +              BM_Face_UpdateNormal(bm, f);
 +      }
 +
 +      BM_Vert_UpdateNormal(bm, v);
 +}
 +
 +void bmesh_update_face_normal(BMesh *bm, BMFace *f, float (*projectverts)[3])
 +{
 +      BMIter iter;
 +      BMLoop *l;
 +
 +      if(f->len > 4) {
 +              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(f->no, projectverts, f->len);       
 +      }
 +      else if(f->len == 3){
 +              BMVert *v1, *v2, *v3;
 +              v1 = bm_firstfaceloop(f)->v;
 +              v2 = bm_firstfaceloop(f)->next->v;
 +              v3 = bm_firstfaceloop(f)->next->next->v;
 +              normal_tri_v3( f->no,v1->co, v2->co, v3->co);
 +      }
 +      else if(f->len == 4){
 +              BMVert *v1, *v2, *v3, *v4;
 +              v1 = bm_firstfaceloop(f)->v;
 +              v2 = bm_firstfaceloop(f)->next->v;
 +              v3 = bm_firstfaceloop(f)->next->next->v;
 +              v4 = bm_firstfaceloop(f)->prev->v;
 +              normal_quad_v3( f->no,v1->co, v2->co, v3->co, v4->co);
 +      }
 +      else{ /*horrible, two sided face!*/
 +              zero_v3(f->no);
 +      }
 +
 +}
 +
 +
 +/*
 + * BMESH FLIP NORMAL
 + * 
 + *  Reverses the winding of a face.
 + *  Note that this updates the calculated 
 + *  normal.
 +*/
 +void BM_flip_normal(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)(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 = 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(float *v1, float *v2, float *v3, float *v4)
 +{
 +      int w1, w2, w3, w4, w5 /*, ret*/;
 +      float mv1[2], mv2[2], mv3[2], mv4[2];
 +      
 +      /*now test winding*/
 +      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 1;
 +      
 +#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 axes*/
 +      /*first do x axis*/
 +      #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 axis*/
 +      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 0; 
 +}
 +
 +/*
 +   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_Point_In_Face(BMesh *bm, BMFace *f, const float co[3])
 +{
-          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)
-       */
-       xn= fabsf(f->no[0]);
-       yn= fabsf(f->no[1]);
-       zn= fabsf(f->no[2]);
-       if (zn >= xn && zn >= yn)       { ax= 0; ay= 1; }
-       else if (yn >= xn && yn >= zn)  { ax= 0; ay= 2; }
-       else                            { ax= 1; ay= 2; }
++      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;
 +      int crosses = 0;
 +      float eps = 1.0f+(float)FLT_EPSILON*150.0f;
 +      
 +      if (dot_v3v3(f->no, f->no) <= FLT_EPSILON*10)
 +              BM_Face_UpdateNormal(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 = bm_firstfaceloop(f);
 +      do {
 +              cent[0] += l->v->co[ax];
 +              cent[1] += l->v->co[ay];
 +              l = l->next;
 +      } while (l != bm_firstfaceloop(f));
 +      
 +      mul_v2_fl(cent, 1.0f/(float)f->len);
 +      
 +      l = bm_firstfaceloop(f);
 +      do {
 +              float v1[3], v2[3];
 +              
 +              v1[0] = (l->prev->v->co[ax] - cent[ax])*eps + cent[ax];
 +              v1[1] = (l->prev->v->co[ay] - cent[ay])*eps + cent[ay];
 +              v1[2] = 0.0f;
 +              
 +              v2[0] = (l->v->co[ax] - cent[ax])*eps + cent[ax];
 +              v2[1] = (l->v->co[ay] - cent[ay])*eps + cent[ay];
 +              v2[2] = 0.0f;
 +              
 +              crosses += linecrossesf(v1, v2, co2, out) != 0;
 +              
 +              l = l->next;
 +      } while (l != bm_firstfaceloop(f));
 +      
 +      return crosses%2 != 0;
 +}
 +
 +static int goodline(float (*projectverts)[3], BMFace *f, int v1i,
 +                    int v2i, int v3i, int UNUSED(nvert))
 +{
 +      BMLoop *l = bm_firstfaceloop(f);
 +      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 = BM_GetIndex(l->v);
 +              if (i == v1i || i == v2i || i == v3i) {
 +                      l = l->next;
 +                      continue;
 +              }
 +              
 +              VECCOPY(pv1, projectverts[BM_GetIndex(l->v)]);
 +              VECCOPY(pv2, projectverts[BM_GetIndex(l->next->v)]);
 +              
 +              //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->next;
 +      } while (l != bm_firstfaceloop(f));
 +      return 1;
 +}
 +/*
 + * 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],
 +                      int nvert)
 +{
 +      BMVert *v1, *v2, *v3;
 +      BMLoop *bestear = NULL, *l;
 +      /*float angle, bestangle = 180.0f;*/
 +      int isear /*, i=0*/;
 +      
 +      l = bm_firstfaceloop(f);
 +      do {
 +              isear = 1;
 +              
 +              v1 = l->prev->v;
 +              v2 = l->v;
 +              v3 = l->next->v;
 +
 +              if (BM_Edge_Exist(v1, v3)) isear = 0;
 +
 +              if (isear && !goodline(verts, f, BM_GetIndex(v1), BM_GetIndex(v2),
 +                                     BM_GetIndex(v3), nvert))
 +                      isear = 0;
 +              
 +              if(isear) {
 +                      /*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;*/
 +                      bestear = l;
 +                      break;
 +              }
 +              l = l->next;
 +      }
 +      while(l != bm_firstfaceloop(f));
 +
 +      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_Triangulate_Face(BMesh *bm, BMFace *f, float (*projectverts)[3], 
 +                         int newedgeflag, int newfaceflag, BMFace **newfaces)
 +{
 +      int i, done, nvert, nf_i = 0;
 +      BMLoop *l, *newl, *nextloop;
 +      /* BMVert *v; */ /* UNUSED */
 +
 +      /*copy vertex coordinates to vertspace array*/
 +      i = 0;
 +      l = bm_firstfaceloop(f);
 +      do{
 +              copy_v3_v3(projectverts[i], l->v->co);
 +              BM_SetIndex(l->v, i); /* set dirty! */
 +              i++;
 +              l = l->next;
 +      }while(l != bm_firstfaceloop(f));
 +
 +      bm->elem_index_dirty |= BM_VERT; /* see above */
 +
 +      ///bmesh_update_face_normal(bm, f, 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 = find_ear(bm, f, projectverts, nvert);
 +              if(l) {
 +                      done = 0;
 +                      /* v = l->v; */ /* UNUSED */
 +                      f = BM_Split_Face(bm, l->f, l->prev->v,
 +                                        l->next->v,
 +                                        &newl, NULL);
 +                      copy_v3_v3(f->no, l->f->no);
 +
 +                      if (!f) {
 +                              fprintf(stderr, "%s: triangulator failed to split face! (bmesh internal error)\n", __func__);
 +                              break;
 +                      }
 +
 +                      BMO_SetFlag(bm, newl->e, newedgeflag);
 +                      BMO_SetFlag(bm, f, newfaceflag);
 +                      
 +                      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 = bm_firstfaceloop(f);
 +              while (l->f->len > 3){
 +                      nextloop = l->next->next;
 +                      f = BM_Split_Face(bm, l->f, l->v, nextloop->v, 
 +                                        &newl, NULL);
 +                      if (!f) {
 +                              printf("triangle fan step of triangulator failed.\n");
 +
 +                              /*NULL-terminate*/
 +                              if (newfaces) newfaces[nf_i] = NULL;
 +                              return;
 +                      }
 +
 +                      if (newfaces) newfaces[nf_i++] = f;
 +                      
 +                      BMO_SetFlag(bm, newl->e, newedgeflag);
 +                      BMO_SetFlag(bm, f, newfaceflag);
 +                      l = nextloop;
 +              }
 +      }
 +      
 +      /*NULL-terminate*/
 +      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_LegalSplits(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 = BMIter_New(&iter, bm, BM_LOOPS_OF_FACE, f);
 +      for (; l; l=BMIter_Step(&iter)) {
 +              BM_SetIndex(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);
 +      
 +      l = bm_firstfaceloop(f);
 +      for (i=0; i<f->len; i++) {
 +              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);
 +
 +              l = l->next;
 +      }
 +      
 +      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 tests*/
 +      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);
 +}
@@@ -1339,13 -1320,10 +1339,13 @@@ static int cube_project_exec(bContext *
  {
        Scene *scene= CTX_data_scene(C);
        Object *obedit= CTX_data_edit_object(C);
 -      EditMesh *em= BKE_mesh_get_editmesh((Mesh*)obedit->data);
 -      EditFace *efa;
 -      MTFace *tf;
 -      float no[3], cube_size, *loc, dx, dy;
 +      BMEditMesh *em= ((Mesh*)obedit->data)->edit_btmesh;
 +      BMFace *efa;
 +      BMLoop *l;
 +      BMIter iter, liter;
 +      /* MTexPoly *tf; */ /* UNUSED */
 +      MLoopUV *luv;
-       float no[3], cube_size, *loc, dx, dy;
++      float cube_size, *loc, dx, dy;
        int cox, coy;
  
        /* add uvs if they don't exist yet */
        /* choose x,y,z axis for projection depending on the largest normal
         * component, but clusters all together around the center of map. */
  
 -      for(efa= em->faces.first; efa; efa= efa->next) {
 -              if(efa->f & SELECT) {
 -                      tf= CustomData_em_get(&em->fdata, efa->data, CD_MTFACE);
 -                      normal_tri_v3( no,efa->v1->co, efa->v2->co, efa->v3->co);
 -
 -                      axis_dominant_v3(&cox, &coy, no);
 -
 -                      tf->uv[0][0]= 0.5f+0.5f*cube_size*(loc[cox] + efa->v1->co[cox]);
 -                      tf->uv[0][1]= 0.5f+0.5f*cube_size*(loc[coy] + efa->v1->co[coy]);
 -                      dx = floor(tf->uv[0][0]);
 -                      dy = floor(tf->uv[0][1]);
 -                      tf->uv[0][0] -= dx;
 -                      tf->uv[0][1] -= dy;
 -                      tf->uv[1][0]= 0.5f+0.5f*cube_size*(loc[cox] + efa->v2->co[cox]);
 -                      tf->uv[1][1]= 0.5f+0.5f*cube_size*(loc[coy] + efa->v2->co[coy]);
 -                      tf->uv[1][0] -= dx;
 -                      tf->uv[1][1] -= dy;
 -                      tf->uv[2][0]= 0.5f+0.5f*cube_size*(loc[cox] + efa->v3->co[cox]);
 -                      tf->uv[2][1]= 0.5f+0.5f*cube_size*(loc[coy] + efa->v3->co[coy]);
 -                      tf->uv[2][0] -= dx;
 -                      tf->uv[2][1] -= dy;
 -
 -                      if(efa->v4) {
 -                              tf->uv[3][0]= 0.5f+0.5f*cube_size*(loc[cox] + efa->v4->co[cox]);
 -                              tf->uv[3][1]= 0.5f+0.5f*cube_size*(loc[coy] + efa->v4->co[coy]);
 -                              tf->uv[3][0] -= dx;
 -                              tf->uv[3][1] -= dy;
 +      BM_ITER(efa, &iter, em->bm, BM_FACES_OF_MESH, NULL) {
 +              int first=1;
 +
 +              /* tf = CustomData_bmesh_get(&em->bm->pdata, efa->head.data, CD_MTEXPOLY); */ /* UNUSED */
 +              if (!BM_TestHFlag(efa, BM_SELECT))
 +                      continue;
-               
-               VECCOPY(no, efa->no);
 +
-               no[0]= fabs(no[0]);
-               no[1]= fabs(no[1]);
-               no[2]= fabs(no[2]);
-               
-               cox=0; coy= 1;
-               if(no[2]>=no[0] && no[2]>=no[1]);
-               else if(no[1]>=no[0] && no[1]>=no[2]) coy= 2;
-               else { cox= 1; coy= 2; }
-               
++              axis_dominant_v3(&cox, &coy, efa->no);
++
 +              dx = dy = 0;
 +              BM_ITER(l, &liter, em->bm, BM_LOOPS_OF_FACE, efa) {
 +                      luv = CustomData_bmesh_get(&em->bm->ldata, l->head.data, CD_MLOOPUV);
 +
 +                      luv->uv[0] = 0.5f+0.5f*cube_size*(loc[cox] + l->v->co[cox]);
 +                      luv->uv[1] = 0.5f+0.5f*cube_size*(loc[coy] + l->v->co[coy]);
 +                      
 +                      if (first) {
 +                              dx = floor(luv->uv[0]);
 +                              dy = floor(luv->uv[1]);
 +                              first = 0;
                        }
 +                      
 +
 +                      luv->uv[0] -= dx;
 +                      luv->uv[1] -= dy;
                }
        }