Normal projection:
authorAndre Susano Pinto <andresusanopinto@gmail.com>
Wed, 7 May 2008 12:45:02 +0000 (12:45 +0000)
committerAndre Susano Pinto <andresusanopinto@gmail.com>
Wed, 7 May 2008 12:45:02 +0000 (12:45 +0000)
+added option to remove faces where all vertices got unprojected

Nearest surface point
+15% faster closest point on point-tri function
(archived by projecting the point on tri-plane and solving the problem on 2D)
(its still using bruteforce on triangles.. I'll add the right data structure later)

source/blender/blenkernel/BKE_shrinkwrap.h
source/blender/blenkernel/intern/shrinkwrap.c
source/blender/makesdna/DNA_modifier_types.h
source/blender/src/buttons_editing.c

index 03c3f897dd4b322e589c72a549025d50e12f0c0b..babdcd78261b3d5442ed5800c4e228e9e4de8e18 100644 (file)
@@ -1,5 +1,5 @@
 /**
- * shrinkwrap.c
+ * BKE_shrinkwrap.h
  *
  * ***** BEGIN GPL LICENSE BLOCK *****
  *
 #ifndef BKE_SHRINKWRAP_H
 #define BKE_SHRINKWRAP_H
 
+/* bitset stuff */
+//TODO: should move this to other generic lib files?
+typedef char* BitSet;
+#define bitset_memsize(size)           (sizeof(char)*((size+7)>>3))
+
+#define bitset_new(size,name)          ((BitSet)MEM_callocN( bitset_memsize(size) , name))
+#define bitset_free(set)                       (MEM_freeN((void*)set))
+
+#define bitset_get(set,index)  ((set)[(index)>>3] & (1 << ((index)&0x7)))
+#define bitset_set(set,index)  ((set)[(index)>>3] |= (1 << ((index)&0x7)))
+
+
 struct Object;
 struct DerivedMesh;
 struct ShrinkwrapModifierData;
 
 
+
 typedef struct ShrinkwrapCalcData
 {
        ShrinkwrapModifierData *smd;    //shrinkwrap modifier data
@@ -50,7 +63,7 @@ typedef struct ShrinkwrapCalcData
 
        float keptDist;                                 //Distance to kept from target (units are in local space)
        //float *weights;                               //weights of vertexs
-       unsigned char *moved;                   //boolean indicating if vertex has moved (TODO use bitmaps)
+       BitSet moved;                                   //BitSet indicating if vertex has moved
 
 } ShrinkwrapCalcData;
 
index 122256aef6f9a8b4d3ea6460556654602d759137..5337a932e792a2547bcb045b647bd5c15f8690a0 100644 (file)
@@ -22,7 +22,7 @@
  *
  * The Original Code is: all of this file.
  *
- * Contributor(s): none yet.
+ * Contributor(s): AndrĂ© Pinto
  *
  * ***** END GPL LICENSE BLOCK *****
  */
 #include "BLI_kdtree.h"
 
 #include "RE_raytrace.h"
+#include "MEM_guardedalloc.h"
 
+
+/* Util macros */
 #define TO_STR(a)      #a
 #define JOIN(a,b)      a##b
 
@@ -82,7 +85,6 @@
 
 #endif
 
-#define CONST
 typedef void ( *Shrinkwrap_ForeachVertexCallback) (DerivedMesh *target, float *co, float *normal);
 
 
@@ -189,7 +191,7 @@ static RayTree* raytree_create_from_mesh(DerivedMesh *mesh)
                //Theres some nasty thing with non-coplanar quads (that I can't find the issue)
                //so we split quads (an odd numbered face represents the second triangle of the quad)
                if(face[i-1].v4)
-                       RE_ray_tree_add_face(tree, 0, i*2+1);
+                       RE_ray_tree_add_face(tree, 0, (RayFace*)(i*2+1));
        }
 
        RE_ray_tree_done(tree);
@@ -228,7 +230,7 @@ static float raytree_cast_ray(RayTree *tree, const float *coord, const float *di
 
        isec.labda = ABS(isec.labda);
        VECADDFAC(isec.end, isec.start, isec.vec, isec.labda);
-       return VecLenf(coord, isec.end);
+       return VecLenf((float*)coord, (float*)isec.end);
 }
 
 /*
@@ -239,7 +241,7 @@ static float raytree_cast_ray(RayTree *tree, const float *coord, const float *di
  *
  * Returns FLT_MIN in parallel case
  */
-static float ray_intersect_plane(CONST float *point, CONST float *dir, CONST float *plane_point, CONST float *plane_normal)
+static float ray_intersect_plane(const float *point, const float *dir, const float *plane_point, const float *plane_normal)
 {
                float pp[3];
                float a, pp_dist;
@@ -248,61 +250,170 @@ static float ray_intersect_plane(CONST float *point, CONST float *dir, CONST flo
 
                if(fabs(a) < 1e-5f) return FLT_MIN;
 
-               VecSubf(pp, point, plane_point);
+               VECSUB(pp, point, plane_point);
                pp_dist = INPR(pp, plane_normal);
 
                return -pp_dist/a;
 }
 
 /*
- * Returns the minimum distance between the point and a triangle surface
- * Writes the nearest surface point in the given nearest
+ * This calculates the distance from point to the plane
+ * Distance is negative if point is on the back side of plane
  */
-static float nearest_point_in_tri_surface(CONST float *co, CONST float *v0, CONST float *v1, CONST float *v2, float *nearest)
+static float point_plane_distance(const float *point, const float *plane_point, const float *plane_normal)
 {
-       //TODO: make this efficient (probably this can be made with something like 3 point_in_slice())
-       if(point_in_tri_prism(co, v0, v1, v2))
-       {
-               float normal[3];
-               float dist;
+       float pp[3];
+       VECSUB(pp, point, plane_point);
+       return INPR(pp, plane_normal);
+}
+static float choose_nearest(const float v0[2], const float v1[2], const float point[2], float closest[2])
+{
+       float d[2][2], sdist[2];
+       VECSUB2D(d[0], v0, point);
+       VECSUB2D(d[1], v1, point);
 
-               CalcNormFloat(v0, v1, v2, normal);
-               dist = ray_intersect_plane(co, normal, v0, normal);
+       sdist[0] = d[0][0]*d[0][0] + d[0][1]*d[0][1];
+       sdist[1] = d[1][0]*d[1][0] + d[1][1]*d[1][1];
 
-               VECADDFAC(nearest, co, normal, dist);
-               return fabs(dist);
+       if(sdist[0] < sdist[1])
+       {
+               if(closest)
+                       VECCOPY2D(closest, v0);
+               return sdist[0];
        }
        else
        {
-               float dist = FLT_MAX, tdist;
-               float closest[3];
+               if(closest)
+                       VECCOPY2D(closest, v1);
+               return sdist[1];
+       }
+}
+/*
+ * calculates the closest point between point-tri (2D)
+ * returns that tri must be right-handed
+ * Returns square distance
+ */
+static float closest_point_in_tri2D(const float point[2], const float tri[3][2], float closest[2])
+{
+       float edge_di[2];
+       float v_point[2];
+       float proj[2];                                  //point projected over edge-dir, edge-normal (witouth normalized edge)
+       const float *v0 = tri[2], *v1;
+       float edge_slen, d;                             //edge squared length
+       int i;
+       const float *nearest_vertex = NULL;
+
+
+       //for each edge
+       for(i=0, v0=tri[2], v1=tri[0]; i < 3; v0=tri[i++], v1=tri[i])
+       {
+               VECSUB2D(edge_di,    v1, v0);
+               VECSUB2D(v_point, point, v0);
 
-               PclosestVL3Dfl(closest, co, v0, v1);
-               tdist = VecLenf(co, closest);
-               if(tdist < dist)
-               {
-                       dist = tdist;
-                       VECCOPY(nearest, closest);
-               }
+               proj[1] =  v_point[0]*edge_di[1] - v_point[1]*edge_di[0];       //dot product with edge normal
 
-               PclosestVL3Dfl(closest, co, v1, v2);
-               tdist = VecLenf(co, closest);
-               if(tdist < dist)
+               //point inside this edge
+               if(proj[1] < 0)
+                       continue;
+
+               proj[0] = v_point[0]*edge_di[0] + v_point[1]*edge_di[1];
+
+               //closest to this edge is v0
+               if(proj[0] < 0)
                {
-                       dist = tdist;
-                       VECCOPY(nearest, closest);
+                       if(nearest_vertex == NULL || nearest_vertex == v0)
+                               nearest_vertex = v0;
+                       else
+                       {
+                               //choose nearest
+                               return choose_nearest(nearest_vertex, v0, point, closest);
+                       }
+                       i++;    //We can skip next edge
+                       continue;
                }
 
-               PclosestVL3Dfl(closest, co, v2, v0);
-               tdist = VecLenf(co, closest);
-               if(tdist < dist)
+               edge_slen = edge_di[0]*edge_di[0] + edge_di[1]*edge_di[1];      //squared edge len
+               //closest to this edge is v1
+               if(proj[0] > edge_slen)
                {
-                       dist = tdist;
-                       VECCOPY(nearest, closest);
+                       if(nearest_vertex == NULL || nearest_vertex == v1)
+                               nearest_vertex = v1;
+                       else
+                       {
+                               return choose_nearest(nearest_vertex, v1, point, closest);
+                       }
+                       continue;
                }
 
-               return dist;
+               //nearest is on this edge
+               d= proj[1] / edge_slen;
+               closest[0] = point[0] - edge_di[1] * d;
+               closest[1] = point[1] + edge_di[0] * d;
+
+               return proj[1]*proj[1]/edge_slen;
+       }
+
+       if(nearest_vertex)
+       {
+               VECSUB2D(v_point, nearest_vertex, point);
+               VECCOPY2D(closest, nearest_vertex);
+               return v_point[0]*v_point[0] + v_point[1]*v_point[1];
+       }
+       else
+       {
+               VECCOPY(closest, point);        //point is already inside
+               return 0.0f;
+       }
+}
+
+/*
+ * Returns the square of the minimum distance between the point and a triangle surface
+ * If nearest is not NULL the nearest surface point is written on it
+ */
+static float nearest_point_in_tri_surface(const float *point, const float *v0, const float *v1, const float *v2, float *nearest)
+{
+       //Lets solve the 2D problem (closest point-tri)
+       float normal_dist, plane_sdist, plane_offset;
+       float du[3], dv[3], dw[3];      //orthogonal axis (du=(v0->v1), dw=plane normal)
+
+       float p_2d[2], tri_2d[3][2], nearest_2d[2];
+
+       CalcNormFloat((float*)v0, (float*)v1, (float*)v2, dw);
+
+       //point-plane distance and calculate axis
+       normal_dist = point_plane_distance(point, v0, dw);
+
+       VECSUB(du, v1, v0);
+       Normalize(du);
+       Crossf(dv, dw, du);
+       plane_offset = INPR(v0, dw);
+
+       //project stuff to 2d
+       tri_2d[0][0] = INPR(du, v0);
+       tri_2d[0][1] = INPR(dv, v0);
+
+       tri_2d[1][0] = INPR(du, v1);
+       tri_2d[1][1] = INPR(dv, v1);
+
+       tri_2d[2][0] = INPR(du, v2);
+       tri_2d[2][1] = INPR(dv, v2);
+
+       p_2d[0] = INPR(du, point);
+       p_2d[1] = INPR(dv, point);
+
+       //we always have a right-handed tri
+       //this should always happen because of the way normal is calculated
+       plane_sdist = closest_point_in_tri2D(p_2d, tri_2d, nearest_2d);
+
+       //project back to 3d
+       if(nearest)
+       {
+               nearest[0] = du[0]*nearest_2d[0] + dv[0] * nearest_2d[1] + dw[0] * plane_offset;
+               nearest[1] = du[1]*nearest_2d[0] + dv[1] * nearest_2d[1] + dw[1] * plane_offset;
+               nearest[2] = du[2]*nearest_2d[0] + dv[2] * nearest_2d[1] + dw[2] * plane_offset;
        }
+
+       return sasqrt(plane_sdist + normal_dist*normal_dist);
 }
 
 
@@ -321,7 +432,7 @@ static void bruteforce_shrinkwrap_calc_nearest_surface_point(DerivedMesh *target
        MVert *vert = target->getVertDataArray(target, CD_MVERT);
        MFace *face = target->getFaceDataArray(target, CD_MFACE);
 
-       VECCOPY(orig_co, co);
+       VECCOPY(orig_co, co);   
 
        for (i = 0; i < numFaces; i++)
        {
@@ -476,6 +587,112 @@ static void shrinkwrap_calc_foreach_vertex(ShrinkwrapCalcData *calc, Shrinkwrap_
        }
 }
 
+
+/*
+ * This function removes Unused faces, vertexs and edges from calc->target
+ *
+ * This function may modify calc->final. As so no data retrieved from
+ * it before the call to this function  can be considered valid
+ * In case it creates a new DerivedMesh, the old calc->final is freed
+ */
+//TODO memory checks on allocs
+static void shrinkwrap_removeUnused(ShrinkwrapCalcData *calc)
+{
+       int i, t;
+
+       DerivedMesh *old = calc->final, *new = NULL;
+       MFace *new_face = NULL;
+       MVert *new_vert  = NULL;
+
+       int numVerts= old->getNumVerts(old);
+       MVert *vert = old->getVertDataArray(old, CD_MVERT);
+
+       int     numFaces= old->getNumFaces(old);
+       MFace *face = old->getFaceDataArray(old, CD_MFACE);
+
+       BitSet moved_verts = calc->moved;
+
+       //Arrays to translate to new vertexs indexs
+       int *vert_index = (int*)MEM_callocN(sizeof(int)*(numVerts), "shrinkwrap used verts");
+       BitSet used_faces = bitset_new(numFaces, "shrinkwrap used faces");
+       int numUsedFaces = 0;
+
+       //calc real number of faces, and vertices
+       //Count used faces
+       for(i=0; i<numFaces; i++)
+       {
+               char res = bitset_get(moved_verts, face[i].v1)
+                                | bitset_get(moved_verts, face[i].v2)
+                                | bitset_get(moved_verts, face[i].v3)
+                                | (face[i].v4 ? bitset_get(moved_verts, face[i].v4) : 0);
+
+               if(res)
+               {
+                       bitset_set(used_faces, i);      //Mark face to maintain
+                       numUsedFaces++;
+
+                       vert_index[face[i].v1] = 1;
+                       vert_index[face[i].v2] = 1;
+                       vert_index[face[i].v3] = 1;
+                       if(face[i].v4) vert_index[face[i].v4] = 1;
+               }
+       }
+
+       //DP: Accumulate vertexs indexs.. (will calculate the new vertex index with a 1 offset)
+       for(i=1; i<numVerts; i++)
+               vert_index[i] += vert_index[i-1];
+               
+       
+       //Start creating the clean mesh
+       new = CDDM_new(vert_index[numVerts-1], 0, numUsedFaces);
+
+       //Copy vertexs (unused are are removed)
+       new_vert  = new->getVertDataArray(new, CD_MVERT);
+       for(i=0, t=0; i<numVerts; i++)
+       {
+               if(vert_index[i] != t)
+               {
+                       t = vert_index[i];
+                       memcpy(new_vert++, vert+i, sizeof(MVert));
+               }
+       }
+
+       //Copy faces
+       new_face = new->getFaceDataArray(new, CD_MFACE);
+       for(i=0, t=0; i<numFaces; i++)
+       {
+               if(bitset_get(used_faces, i))
+               {
+                       memcpy(new_face, face+i, sizeof(MFace));
+                       //update vertices indexs
+                       new_face->v1 = vert_index[new_face->v1]-1;
+                       new_face->v2 = vert_index[new_face->v2]-1;
+                       new_face->v3 = vert_index[new_face->v3]-1;
+                       if(new_face->v4)
+                       {
+                               new_face->v4 = vert_index[new_face->v4]-1;
+
+                               //Ups translated vertex ended on 0 .. TODO fix this
+                               if(new_face->v4 == 0)
+                               {
+                               }
+                       }                       
+                       new_face++;
+               }
+       }
+
+       //Free memory
+       bitset_free(used_faces);
+       MEM_freeN(vert_index);
+       old->release(old);
+
+       //Update edges
+       CDDM_calc_edges(new);
+       CDDM_calc_normals(new);
+
+       calc->final = new;
+}
+
 /* Main shrinkwrap function */
 DerivedMesh *shrinkwrapModifier_do(ShrinkwrapModifierData *smd, Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc)
 {
@@ -514,9 +731,6 @@ DerivedMesh *shrinkwrapModifier_do(ShrinkwrapModifierData *smd, Object *ob, Deri
                calc.keptDist = smd->keptDist;  //TODO: smd->keptDist is in global units.. must change to local
        }
 
-       calc.moved = NULL;
-
-
        //Projecting target defined - lets work!
        if(calc.target)
        {
@@ -547,7 +761,8 @@ DerivedMesh *shrinkwrapModifier_do(ShrinkwrapModifierData *smd, Object *ob, Deri
        //Destroy faces, edges and stuff
        if(calc.moved)
        {
-               //TODO
+               shrinkwrap_removeUnused(&calc);
+               bitset_free(calc.moved);
        }
 
        CDDM_calc_normals(calc.final);  
@@ -621,9 +836,6 @@ void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
                        dist = VecLenf(vert[i].co, tmp_co);
                        if(dist > 1e-5) weight *= (dist - calc->keptDist)/dist;
                        VecLerpf(vert[i].co, vert[i].co, nearest.co, weight);   //linear interpolation
-
-                       if(calc->moved)
-                               calc->moved[i] = TRUE;
                }
        }
 
@@ -668,6 +880,9 @@ void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
        vert    = calc->final->getVertDataArray(calc->final, CD_MVERT); 
        dvert   = calc->final->getVertDataArray(calc->final, CD_MDEFORMVERT);
 
+       if(calc->smd->shrinkOpts & MOD_SHRINKWRAP_REMOVE_UNPROJECTED_FACES)
+               calc->moved = bitset_new(numVerts, "shrinkwrap bitset data");
+
        for(i=0; i<numVerts; i++)
        {
                float dist = FLT_MAX;
@@ -718,7 +933,7 @@ void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
                        VecLerpf(vert[i].co, vert[i].co, tmp_co, weight);       //linear interpolation
 
                        if(calc->moved)
-                               calc->moved[i] = TRUE;
+                               bitset_set(calc->moved, i);
                }
 
        }
index 2a38358c3834421b5a1f60387c9383fac9efdd51..5178022ffad97c89ade7a3d0f92a33675ea9b8f7 100644 (file)
@@ -507,5 +507,6 @@ typedef struct ShrinkwrapModifierData {
 /* Shrinkwrap->shrinkOpts */
 #define MOD_SHRINKWRAP_ALLOW_DEFAULT_NORMAL            (1<<0)
 #define MOD_SHRINKWRAP_ALLOW_INVERTED_NORMAL   (1<<1)
+#define MOD_SHRINKWRAP_REMOVE_UNPROJECTED_FACES        (1<<2)
 
 #endif
index 47a2c6d613de4972c3333f0534b8fe92a54f24cd..68712da0b6116346af6490675e1a902ab2ce437a 100644 (file)
@@ -1824,7 +1824,7 @@ static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco
                        ShrinkwrapModifierData *smd = (ShrinkwrapModifierData*) md;
                        height = 86;
                        if (smd->shrinkType == MOD_SHRINKWRAP_NORMAL)
-                               height += 19*2;
+                               height += 19*3;
                }
                                                        /* roundbox 4 free variables: corner-rounding, nop, roundbox type, shade */
                uiDefBut(block, ROUNDBOX, 0, "", x-10, y-height-2, width, height-2, NULL, 5.0, 0.0, 12, 40, ""); 
@@ -2446,6 +2446,7 @@ static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco
                        if (smd->shrinkType == MOD_SHRINKWRAP_NORMAL){
                                uiDefButBitS(block, TOG, MOD_SHRINKWRAP_ALLOW_DEFAULT_NORMAL, B_MODIFIER_RECALC, "Default normal",      lx,(cy-=19),buttonWidth,19, &smd->shrinkOpts, 0, 0, 0, 0, "Allows vertices to move in the normal direction");
                                uiDefButBitS(block, TOG, MOD_SHRINKWRAP_ALLOW_INVERTED_NORMAL, B_MODIFIER_RECALC, "Invert normal",      lx,(cy-=19),buttonWidth,19, &smd->shrinkOpts, 0, 0, 0, 0, "Allows vertices to move in the inverse direction of their normal");
+                               uiDefButBitS(block, TOG, MOD_SHRINKWRAP_REMOVE_UNPROJECTED_FACES, B_MODIFIER_RECALC, "Remove faces",    lx,(cy-=19),buttonWidth,19, &smd->shrinkOpts, 0, 0, 0, 0, "Remove faces where all vertices haven't been projected");
                        }
 
                        but=uiDefBut(block, TEX, B_MODIFIER_RECALC, "VGroup: ",         lx, (cy-=19), buttonWidth,19, &smd->vgroup_name, 0.0, 31.0, 0, 0, "Vertex Group name");