Merge from trunk
authorAndre Susano Pinto <andresusanopinto@gmail.com>
Fri, 30 May 2008 20:33:33 +0000 (20:33 +0000)
committerAndre Susano Pinto <andresusanopinto@gmail.com>
Fri, 30 May 2008 20:33:33 +0000 (20:33 +0000)
svn merge  -r 14954:15064 https://svn.blender.org/svnroot/bf-blender/trunk/blender

release/scripts/scripttemplate_ipo_gen.py [deleted file]
source/blender/blenkernel/BKE_shrinkwrap.h [new file with mode: 0644]
source/blender/blenkernel/intern/modifier.c
source/blender/blenkernel/intern/shrinkwrap.c [new file with mode: 0644]
source/blender/blenlib/BLI_kdopbvh.h [new file with mode: 0644]
source/blender/blenlib/intern/BLI_kdopbvh.c [new file with mode: 0644]
source/blender/makesdna/DNA_modifier_types.h
source/blender/src/buttons_editing.c

diff --git a/release/scripts/scripttemplate_ipo_gen.py b/release/scripts/scripttemplate_ipo_gen.py
deleted file mode 100644 (file)
index 791eaed..0000000
+++ /dev/null
@@ -1,92 +0,0 @@
-#!BPY
-"""
-Name: 'Camera/Object Example'
-Blender: 245
-Group: 'ScriptTemplate'
-Tooltip: 'Script template for setting the camera direction'
-"""
-
-from Blender import Window
-import bpy
-
-script_data = \
-'''#!BPY
-"""
-Name: 'My Ipo Script'
-Blender: 245
-Group: 'Animation'
-Tooltip: 'Put some useful info here'
-"""
-
-# Add a licence here if you wish to re-distribute, we recommend the GPL
-
-from Blender import Ipo, Mathutils, Window
-import bpy, BPyMessages
-
-def makeRandomIpo(object, firstFrame, numberOfFrames, frameStep):
-       # Create an new Ipo Curve of name myIpo and type Object
-       myIpo = bpy.data.ipos.new('myIpo', 'Object')
-       
-       # Create LocX, LocY, and LocZ Ipo curves in our new Curve Object
-       # and store them so we can access them later
-       myIpo_x = myIpo.addCurve('LocX')
-       myIpo_y = myIpo.addCurve('LocY')
-       myIpo_z = myIpo.addCurve('LocZ')
-       
-       # What value we want to scale our random value by
-       ipoScale = 4
-       
-       # This Calculates the End Frame for use in an xrange() expression
-       endFrame = firstFrame + (numberOfFrames * frameStep) + frameStep
-       
-       for frame in xrange(firstFrame, endFrame, frameStep):
-               
-               # Use the Mathutils Rand() function to get random numbers
-               ipoValue_x = Mathutils.Rand(-1, 1) * ipoScale
-               ipoValue_y = Mathutils.Rand(-1, 1) * ipoScale
-               ipoValue_z = Mathutils.Rand(-1, 1) * ipoScale
-               
-               # Append to the Ipo curve at location frame, with the value ipoValue_x
-               # Note that we should pass the append function a tuple or a BezTriple
-               myIpo_x.append((frame, ipoValue_x))
-       
-               # Similar to above
-               myIpo_y.append((frame, ipoValue_y))
-               myIpo_z.append((frame, ipoValue_z))
-       
-       # Link our new Ipo Curve to the passed object
-       object.setIpo(myIpo)
-       print object
-       
-       
-def main():
-       
-       # Get the active scene, since there can be multiple ones
-       sce = bpy.data.scenes.active
-       
-       # Get the active object
-       object = sce.objects.active
-       
-       # If there is no active object, pop up an error message
-       if not object:
-               BPyMessages.Error_NoActive()
-               
-       Window.WaitCursor(1)
-       
-       # Call our makeRandomIpo function
-       # Pass it our object, Tell it to keys from the start frame until the end frame, at a step of 10 frames
-       # between them
-       
-       makeRandomIpo(object, sce.render.sFrame, sce.render.eFrame, 10)
-       
-       Window.WaitCursor(0)
-
-if __name__ == '__main__':
-       main()
-
-'''
-
-new_text = bpy.data.texts.new('ipo_template.py')
-new_text.write(script_data)
-bpy.data.texts.active = new_text
-Window.RedrawAll()
diff --git a/source/blender/blenkernel/BKE_shrinkwrap.h b/source/blender/blenkernel/BKE_shrinkwrap.h
new file mode 100644 (file)
index 0000000..babdcd7
--- /dev/null
@@ -0,0 +1,78 @@
+/**
+ * BKE_shrinkwrap.h
+ *
+ * ***** 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) Blender Foundation.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): none yet.
+ *
+ * ***** END 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
+
+       struct Object *ob;                              //object we are applying shrinkwrap to
+       struct DerivedMesh *original;   //mesh before shrinkwrap
+       struct DerivedMesh *final;              //initially a copy of original mesh.. mesh thats going to be shrinkwrapped
+
+       struct DerivedMesh *target;             //mesh we are shrinking to
+       
+       //matrixs for local<->target space transform
+       float local2target[4][4];               
+       float target2local[4][4];
+
+       float keptDist;                                 //Distance to kept from target (units are in local space)
+       //float *weights;                               //weights of vertexs
+       BitSet moved;                                   //BitSet indicating if vertex has moved
+
+} ShrinkwrapCalcData;
+
+void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *data);
+void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *data);
+void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *data);
+
+struct DerivedMesh *shrinkwrapModifier_do(struct ShrinkwrapModifierData *smd, struct Object *ob, struct DerivedMesh *dm, int useRenderParams, int isFinalCalc);
+
+#endif
+
+
index f9f17e7762da2b9409b8dd3e4857be8577478909..73005f3708e3235f2fca4ac09b031929242216ee 100644 (file)
@@ -97,6 +97,7 @@
 #include "BKE_utildefines.h"
 #include "depsgraph_private.h"
 #include "BKE_bmesh.h"
+#include "BKE_shrinkwrap.h"
 
 #include "LOD_DependKludge.h"
 #include "LOD_decimation.h"
@@ -7011,6 +7012,63 @@ static void meshdeformModifier_deformVertsEM(
                dm->release(dm);
 }
 
+
+/* Shrinkwrap */
+
+static void shrinkwrapModifier_initData(ModifierData *md)
+{
+       ShrinkwrapModifierData *smd = (ShrinkwrapModifierData*) md;
+       smd->shrinkType = MOD_SHRINKWRAP_NEAREST_SURFACE;
+       smd->shrinkOpts = MOD_SHRINKWRAP_ALLOW_DEFAULT_NORMAL;
+       smd->keptDist   = 0.0f;
+}
+
+static void shrinkwrapModifier_copyData(ModifierData *md, ModifierData *target)
+{
+       memcpy(target, md, sizeof(MeshDeformModifierData));
+}
+
+CustomDataMask shrinkwrapModifier_requiredDataMask(ModifierData *md)
+{
+       ShrinkwrapModifierData *smd = (ShrinkwrapModifierData *)md;
+       CustomDataMask dataMask = 0;
+
+       /* ask for vertexgroups if we need them */
+       if(smd->vgroup_name[0]) dataMask |= (1 << CD_MDEFORMVERT);
+
+       return dataMask;
+}
+
+static void shrinkwrapModifier_foreachObjectLink(ModifierData *md, Object *ob, ObjectWalkFunc walk, void *userData)
+{
+       ShrinkwrapModifierData *smd = (ShrinkwrapModifierData*) md;
+
+       walk(userData, ob, &smd->target);
+}
+
+static DerivedMesh *shrinkwrapModifier_applyModifier(ModifierData *md, Object *ob, DerivedMesh *derivedData, int useRenderParams, int isFinalCalc)
+{
+       return shrinkwrapModifier_do((ShrinkwrapModifierData*)md,ob,derivedData,useRenderParams,isFinalCalc);
+}
+
+static DerivedMesh *shrinkwrapModifier_applyModifierEM(ModifierData *md, Object *ob, EditMesh *editData, DerivedMesh *derivedData)
+{
+       return shrinkwrapModifier_do((ShrinkwrapModifierData*)md,ob,derivedData,0,0);
+}
+
+static void shrinkwrapModifier_updateDepgraph(ModifierData *md, DagForest *forest, Object *ob, DagNode *obNode)
+{
+       ShrinkwrapModifierData *smd = (ShrinkwrapModifierData*) md;
+
+       if (smd->target) {
+               DagNode *curNode = dag_get_node(forest, smd->target);
+
+               dag_add_relation(forest, curNode, obNode, DAG_RL_OB_DATA | DAG_RL_DATA_DATA,
+                       "Shrinkwrap Modifier");
+       }
+}
+
+
 /***/
 
 static ModifierTypeInfo typeArr[NUM_MODIFIER_TYPES];
@@ -7331,6 +7389,20 @@ ModifierTypeInfo *modifierType_getInfo(ModifierType type)
                mti->requiredDataMask = explodeModifier_requiredDataMask;
                mti->applyModifier = explodeModifier_applyModifier;
 
+               mti = INIT_TYPE(Shrinkwrap);
+               mti->type = eModifierTypeType_Nonconstructive;
+               mti->flags = eModifierTypeFlag_AcceptsMesh
+                               | eModifierTypeFlag_SupportsEditmode
+                               | eModifierTypeFlag_EnableInEditmode;
+
+               mti->initData = shrinkwrapModifier_initData;
+               mti->copyData = shrinkwrapModifier_copyData;
+               mti->requiredDataMask = shrinkwrapModifier_requiredDataMask;
+               mti->foreachObjectLink = shrinkwrapModifier_foreachObjectLink;
+               mti->applyModifier = shrinkwrapModifier_applyModifier;
+               mti->applyModifierEM = shrinkwrapModifier_applyModifierEM;
+               mti->updateDepgraph = shrinkwrapModifier_updateDepgraph;
+
                typeArrInit = 0;
 #undef INIT_TYPE
        }
diff --git a/source/blender/blenkernel/intern/shrinkwrap.c b/source/blender/blenkernel/intern/shrinkwrap.c
new file mode 100644 (file)
index 0000000..e54b7a8
--- /dev/null
@@ -0,0 +1,1078 @@
+/**
+ * shrinkwrap.c
+ *
+ * ***** 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) Blender Foundation.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): AndrĂ© Pinto
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+#include <string.h>
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+#include <time.h>
+
+#include "DNA_object_types.h"
+#include "DNA_modifier_types.h"
+#include "DNA_meshdata_types.h"
+
+#include "BKE_shrinkwrap.h"
+#include "BKE_DerivedMesh.h"
+#include "BKE_utildefines.h"
+#include "BKE_deform.h"
+#include "BKE_cdderivedmesh.h"
+#include "BKE_global.h"
+
+#include "BLI_arithb.h"
+#include "BLI_kdtree.h"
+#include "BLI_kdopbvh.h"
+
+#include "RE_raytrace.h"
+#include "MEM_guardedalloc.h"
+
+
+/* Util macros */
+#define TO_STR(a)      #a
+#define JOIN(a,b)      a##b
+
+#define OUT_OF_MEMORY()        ((void)printf("Shrinkwrap: Out of memory\n"))
+
+/* Benchmark macros */
+#if 1
+
+#define BENCH(a)       \
+       do {                    \
+               clock_t _clock_init = clock();  \
+               (a);                                                    \
+               printf("%s: %fms\n", #a, (float)(clock()-_clock_init)*1000/CLOCKS_PER_SEC);     \
+       } while(0)
+
+#define BENCH_VAR(name)                clock_t JOIN(_bench_step,name) = 0, JOIN(_bench_total,name) = 0
+#define BENCH_BEGIN(name)      JOIN(_bench_step, name) = clock()
+#define BENCH_END(name)                JOIN(_bench_total,name) += clock() - JOIN(_bench_step,name)
+#define BENCH_RESET(name)      JOIN(_bench_total, name) = 0
+#define BENCH_REPORT(name)     printf("%s: %fms\n", TO_STR(name), JOIN(_bench_total,name)*1000.0f/CLOCKS_PER_SEC)
+
+#else
+
+#define BENCH(a)       (a)
+#define BENCH_VAR(name)
+#define BENCH_BEGIN(name)
+#define BENCH_END(name)
+#define BENCH_RESET(name)
+#define BENCH_REPORT(name)
+
+#endif
+
+typedef void ( *Shrinkwrap_ForeachVertexCallback) (DerivedMesh *target, float *co, float *normal);
+
+static float nearest_point_in_tri_surface(const float *point, const float *v0, const float *v1, const float *v2, float *nearest);
+
+static void normal_short2float(const short *ns, float *nf)
+{
+       nf[0] = ns[0] / 32767.0f;
+       nf[1] = ns[1] / 32767.0f;
+       nf[2] = ns[2] / 32767.0f;
+}
+
+static float vertexgroup_get_weight(MDeformVert *dvert, int index, int vgroup)
+{
+       if(dvert && vgroup >= 0)
+       {
+               int j;
+               for(j = 0; j < dvert[index].totweight; j++)
+                       if(dvert[index].dw[j].def_nr == vgroup)
+                               return dvert[index].dw[j].weight;
+       }
+       return 1.0;
+}
+
+/*
+ * BVH tree from mesh vertices
+ */
+static BVHTree* bvhtree_from_mesh_verts(DerivedMesh *mesh)
+{
+       int i;
+       int numVerts= mesh->getNumVerts(mesh);
+       MVert *vert     = mesh->getVertDataArray(mesh, CD_MVERT);
+
+       BVHTree *tree = BLI_bvhtree_new(numVerts, 0, 2, 6);
+       if(tree != NULL)
+       {
+               for(i = 0; i < numVerts; i++)
+                       BLI_bvhtree_insert(tree, i, vert[i].co, 1);
+
+               BLI_bvhtree_balance(tree);
+       }
+
+       return tree;
+}
+
+static BVHTree* bvhtree_from_mesh_tri(DerivedMesh *mesh)
+{
+       int i;
+       int numFaces= mesh->getNumFaces(mesh), totFaces;
+       MVert *vert     = mesh->getVertDataArray(mesh, CD_MVERT);
+       MFace *face = mesh->getFaceDataArray(mesh, CD_MFACE);
+       BVHTree *tree= NULL;
+
+       /* Count needed faces */
+       for(totFaces=numFaces, i=0; i<numFaces; i++)
+               if(face[i].v4) totFaces++;
+
+       /* Create a bvh-tree of the given target */
+       tree = BLI_bvhtree_new(totFaces, 0, 2, 6);
+       if(tree != NULL)
+       {
+               for(i = 0; i < numFaces; i++)
+               {
+                       float co[3][3];
+
+                       VECCOPY(co[0], vert[ face[i].v1 ].co);
+                       VECCOPY(co[1], vert[ face[i].v2 ].co);
+                       VECCOPY(co[2], vert[ face[i].v3 ].co);
+                       BLI_bvhtree_insert(tree, 2*i, co[0], 3);
+                       if(face[i].v4)
+                       {
+                               /* second face is v1,v3,v4 */
+                               VECCOPY(co[1], vert[ face[i].v3 ].co);
+                               VECCOPY(co[2], vert[ face[i].v4 ].co);
+                               BLI_bvhtree_insert(tree, 2*i+1, co[0], 3);
+                       }
+               }
+
+               BLI_bvhtree_balance(tree);
+       }
+
+       return tree;
+}
+
+static float mesh_tri_nearest_point(void *userdata, int index, const float *co, float *nearest)
+{
+       DerivedMesh *mesh = (DerivedMesh*)(userdata);
+       MVert *vert     = (MVert*)mesh->getVertDataArray(mesh, CD_MVERT);
+       MFace *face = (MFace*)mesh->getFaceDataArray(mesh, CD_MFACE) + index/2;
+
+       if(index & 1)
+               return nearest_point_in_tri_surface(co, vert[ face->v1 ].co, vert[ face->v3 ].co, vert[ face->v4 ].co, nearest);
+       else
+               return nearest_point_in_tri_surface(co, vert[ face->v1 ].co, vert[ face->v2 ].co, vert[ face->v3 ].co, nearest);
+}
+
+/*
+ * Raytree from mesh
+ */
+static MVert *raytree_from_mesh_verts = NULL;
+static MFace *raytree_from_mesh_faces = NULL;
+
+static int raytree_check_always(Isect *is, int ob, RayFace *face)
+{
+       return TRUE;
+}
+
+static void raytree_from_mesh_get_coords(RayFace *face, float **v1, float **v2, float **v3, float **v4)
+{
+       MFace *mface= raytree_from_mesh_faces + (int)face/2 - 1 ;
+
+       if(face == (RayFace*)(-1))
+       {
+               *v1 = NULL;
+               *v2 = NULL;
+               *v3 = NULL;
+               *v4 = NULL;
+               return;
+       }
+
+       //Nasty quad splitting
+       if(((int)face) & 1)     // we want the 2 triangle of the quad
+       {
+               *v1= raytree_from_mesh_verts[mface->v1].co;
+               *v2= raytree_from_mesh_verts[mface->v4].co;
+               *v3= raytree_from_mesh_verts[mface->v3].co;
+               *v4= NULL;
+       }
+       else
+       {
+               *v1= raytree_from_mesh_verts[mface->v1].co;
+               *v2= raytree_from_mesh_verts[mface->v2].co;
+               *v3= raytree_from_mesh_verts[mface->v3].co;
+               *v4= NULL;
+       }
+}
+
+/*
+ * Creates a raytree from the given mesh
+ * No copy of the mesh is done, so it must exist and remain
+ * imutable as long the tree is intended to be used
+ *
+ * No more than 1 raytree can exist.. since this code uses a static variable
+ * to pass data to raytree_from_mesh_get_coords
+ */
+static RayTree* raytree_create_from_mesh(DerivedMesh *mesh)
+{
+       int i;
+       float min[3], max[3];
+
+       RayTree*tree= NULL;
+
+       int numFaces= mesh->getNumFaces(mesh);
+       MFace *face = mesh->getFaceDataArray(mesh, CD_MFACE);
+       int numVerts= mesh->getNumVerts(mesh);
+
+       //Initialize static vars
+       raytree_from_mesh_verts = mesh->getVertDataArray(mesh, CD_MVERT);
+       raytree_from_mesh_faces = face;
+
+
+       //calculate bounding box
+       INIT_MINMAX(min, max);
+
+       for(i=0; i<numVerts; i++)
+               DO_MINMAX(raytree_from_mesh_verts[i].co, min, max);
+       
+       tree = RE_ray_tree_create(64, numFaces, min, max, raytree_from_mesh_get_coords, raytree_check_always, NULL, NULL);
+       if(tree == NULL)
+               return NULL;
+
+       //Add faces to the RayTree (RayTree uses face=0, with some special value to setup things)
+       for(i=1; i<=numFaces; i++)
+       {
+               RE_ray_tree_add_face(tree, 0, (RayFace*)(i*2) );
+
+                //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, (RayFace*)(i*2+1));
+       }
+
+       RE_ray_tree_done(tree);
+
+       return tree;
+}
+
+static void free_raytree_from_mesh(RayTree *tree)
+{
+       raytree_from_mesh_verts = NULL;
+       RE_ray_tree_free(tree);
+}
+
+/*
+ * Cast a ray on the specified direction
+ * Returns the distance the ray must travel until intersect something
+ * Returns FLT_MAX in case of nothing intersection
+ */
+static float raytree_cast_ray(RayTree *tree, const float *coord, const float *direction)
+{
+       Isect isec;
+
+       /* Setup intersection */
+       isec.mode               = RE_RAY_MIRROR; /* We want closest intersection */
+       isec.lay                = -1;
+       isec.face_last  = NULL;
+       isec.faceorig   = (RayFace*)(-1);
+       isec.labda              = 1e10f;
+
+       VECCOPY(isec.start, coord);
+       VECCOPY(isec.vec, direction);
+       VECADDFAC(isec.end, isec.start, isec.vec, isec.labda);
+
+       if(!RE_ray_tree_intersect(tree, &isec))
+               return FLT_MAX;
+
+       isec.labda = ABS(isec.labda);
+       VECADDFAC(isec.end, isec.start, isec.vec, isec.labda);
+       return VecLenf((float*)coord, (float*)isec.end);
+}
+
+/*
+ * Returns the squared distance between two given points
+ */
+static float squared_dist(const float *a, const float *b)
+{
+       float tmp[3];
+       VECSUB(tmp, a, b);
+       return INPR(tmp, tmp);
+}
+
+/*
+ * This calculates the distance (in dir units) that the ray must travel to intersect plane
+ * It can return negative values
+ *
+ * TODO theres probably something like this on blender code
+ *
+ * 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)
+{
+               float pp[3];
+               float a, pp_dist;
+
+               a = INPR(dir, plane_normal);
+
+               if(fabs(a) < 1e-5f) return FLT_MIN;
+
+               VECSUB(pp, point, plane_point);
+               pp_dist = INPR(pp, plane_normal);
+
+               return -pp_dist/a;
+}
+
+/*
+ * This calculates the distance from point to the plane
+ * Distance is negative if point is on the back side of plane
+ */
+static float point_plane_distance(const float *point, const float *plane_point, const float *plane_normal)
+{
+       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);
+
+       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];
+
+       if(sdist[0] < sdist[1])
+       {
+               if(closest)
+                       VECCOPY2D(closest, v0);
+               return sdist[0];
+       }
+       else
+       {
+               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);
+
+               proj[1] =  v_point[0]*edge_di[1] - v_point[1]*edge_di[0];       //dot product with edge normal
+
+               //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)
+               {
+                       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;
+               }
+
+               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)
+               {
+                       if(nearest_vertex == NULL || nearest_vertex == v1)
+                               nearest_vertex = v1;
+                       else
+                       {
+                               return choose_nearest(nearest_vertex, v1, point, closest);
+                       }
+                       continue;
+               }
+
+               //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 plane_sdist + normal_dist*normal_dist;
+}
+
+
+
+/*
+ * Shrink to nearest surface point on target mesh
+ */
+static void bruteforce_shrinkwrap_calc_nearest_surface_point(DerivedMesh *target, float *co, float *unused)
+{
+       float minDist = FLT_MAX;
+       float orig_co[3];
+
+       int i;
+       int     numFaces = target->getNumFaces(target);
+       MVert *vert = target->getVertDataArray(target, CD_MVERT);
+       MFace *face = target->getFaceDataArray(target, CD_MFACE);
+
+       VECCOPY(orig_co, co);   
+
+       for (i = 0; i < numFaces; i++)
+       {
+               float *v0, *v1, *v2, *v3;
+
+               v0 = vert[ face[i].v1 ].co;
+               v1 = vert[ face[i].v2 ].co;
+               v2 = vert[ face[i].v3 ].co;
+               v3 = face[i].v4 ? vert[ face[i].v4 ].co : 0;
+
+               while(v2)
+               {
+                       float dist;
+                       float tmp[3];
+
+                       dist = nearest_point_in_tri_surface(orig_co, v0, v1, v2, tmp);
+
+                       if(dist < minDist)
+                       {
+                               minDist = dist;
+                               VECCOPY(co, tmp);
+                       }
+
+                       v1 = v2;
+                       v2 = v3;
+                       v3 = 0;
+               }
+       }
+}
+
+/*
+ * Projects the vertex on the normal direction over the target mesh
+ */
+static void bruteforce_shrinkwrap_calc_normal_projection(DerivedMesh *target, float *co, float *vnormal)
+{
+       //TODO: this should use raycast code probably existent in blender
+       float minDist = FLT_MAX;
+       float orig_co[3];
+
+       int i;
+       int     numFaces = target->getNumFaces(target);
+       MVert *vert = target->getVertDataArray(target, CD_MVERT);
+       MFace *face = target->getFaceDataArray(target, CD_MFACE);
+
+       VECCOPY(orig_co, co);
+
+       for (i = 0; i < numFaces; i++)
+       {
+               float *v0, *v1, *v2, *v3;
+
+               v0 = vert[ face[i].v1 ].co;
+               v1 = vert[ face[i].v2 ].co;
+               v2 = vert[ face[i].v3 ].co;
+               v3 = face[i].v4 ? vert[ face[i].v4 ].co : 0;
+
+               while(v2)
+               {
+                       float dist;
+                       float pnormal[3];
+
+                       CalcNormFloat(v0, v1, v2, pnormal);
+                       dist =  ray_intersect_plane(orig_co, vnormal, v0, pnormal);
+
+                       if(fabs(dist) < minDist)
+                       {
+                               float tmp[3], nearest[3];
+                               VECADDFAC(tmp, orig_co, vnormal, dist);
+
+                               if( fabs(nearest_point_in_tri_surface(tmp, v0, v1, v2, nearest)) < 0.0001)
+                               {
+                                       minDist = fabs(dist);
+                                       VECCOPY(co, nearest);
+                               }
+                       }
+                       v1 = v2;
+                       v2 = v3;
+                       v3 = 0;
+               }
+       }
+}
+
+/*
+ * Shrink to nearest vertex on target mesh
+ */
+static void bruteforce_shrinkwrap_calc_nearest_vertex(DerivedMesh *target, float *co, float *unused)
+{
+       float minDist = FLT_MAX;
+       float orig_co[3];
+
+       int i;
+       int     numVerts = target->getNumVerts(target);
+       MVert *vert = target->getVertDataArray(target, CD_MVERT);
+
+       VECCOPY(orig_co, co);
+
+       for (i = 0; i < numVerts; i++)
+       {
+               float sdist = squared_dist( orig_co, vert[i].co);
+               
+               if(sdist < minDist)
+               {
+                       minDist = sdist;
+                       VECCOPY(co, vert[i].co);
+               }
+       }
+}
+
+
+static void shrinkwrap_calc_foreach_vertex(ShrinkwrapCalcData *calc, Shrinkwrap_ForeachVertexCallback callback)
+{
+       int i;
+       int vgroup              = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
+       int     numVerts        = 0;
+
+       MDeformVert *dvert = NULL;
+       MVert           *vert  = NULL;
+
+       numVerts = calc->final->getNumVerts(calc->final);
+       dvert = calc->final->getVertDataArray(calc->final, CD_MDEFORMVERT);
+       vert  = calc->final->getVertDataArray(calc->final, CD_MVERT);
+
+       //Shrink (calculate each vertex final position)
+       for(i = 0; i<numVerts; i++)
+       {
+               float weight = vertexgroup_get_weight(dvert, i, vgroup);
+
+               float orig[3], final[3]; //Coords relative to target
+               float normal[3];
+               float dist;
+
+               if(weight == 0.0f) continue;    //Skip vertexs where we have no influence
+
+               VecMat4MulVecfl(orig, calc->local2target, vert[i].co);
+               VECCOPY(final, orig);
+
+               //We also need to apply the rotation to normal
+               if(calc->smd->shrinkType == MOD_SHRINKWRAP_NORMAL)
+               {
+                       normal_short2float(vert[i].no, normal);
+                       Mat4Mul3Vecfl(calc->local2target, normal);
+                       Normalize(normal);      //Watch out for scaling (TODO: do we really needed a unit-len normal?)
+               }
+               (callback)(calc->target, final, normal);
+
+               VecMat4MulVecfl(final, calc->target2local, final);
+
+               dist = VecLenf(vert[i].co, final);
+               if(dist > 1e-5) weight *= (dist - calc->keptDist)/dist;
+               VecLerpf(vert[i].co, vert[i].co, final, weight);        //linear interpolation
+       }
+}
+
+
+/*
+ * 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)
+{
+
+       ShrinkwrapCalcData calc;
+       memset(&calc, 0, sizeof(calc));
+
+       //Init Shrinkwrap calc data
+       calc.smd = smd;
+
+       calc.ob = ob;
+       calc.original = dm;
+       calc.final = CDDM_copy(calc.original);
+
+       if(!calc.final)
+       {
+               OUT_OF_MEMORY();
+               return dm;
+       }
+
+       if(smd->target)
+       {
+               calc.target = (DerivedMesh *)smd->target->derivedFinal;
+
+               if(!calc.target)
+               {
+                       printf("Target derived mesh is null! :S\n");
+               }
+
+               //TODO should we reduce the number of matrix mults? by choosing applying matrixs to target or to derived mesh?
+               //Calculate matrixs for local <-> target
+               Mat4Invert (smd->target->imat, smd->target->obmat);     //inverse is outdated
+               Mat4MulSerie(calc.local2target, smd->target->imat, ob->obmat, 0, 0, 0, 0, 0, 0);
+               Mat4Invert(calc.target2local, calc.local2target);
+       
+               calc.keptDist = smd->keptDist;  //TODO: smd->keptDist is in global units.. must change to local
+       }
+
+       //Projecting target defined - lets work!
+       if(calc.target)
+       {
+               printf("Shrinkwrap (%s)%d over (%s)%d\n",
+                       calc.ob->id.name,                       calc.final->getNumVerts(calc.final),
+                       calc.smd->target->id.name,      calc.target->getNumVerts(calc.target)
+               );
+
+               switch(smd->shrinkType)
+               {
+                       case MOD_SHRINKWRAP_NEAREST_SURFACE:
+                               BENCH(shrinkwrap_calc_nearest_surface_point(&calc));
+//                             BENCH(shrinkwrap_calc_foreach_vertex(&calc, bruteforce_shrinkwrap_calc_nearest_surface_point));
+                       break;
+
+                       case MOD_SHRINKWRAP_NORMAL:
+                               BENCH(shrinkwrap_calc_normal_projection(&calc));
+//                             BENCH(shrinkwrap_calc_foreach_vertex(&calc, bruteforce_shrinkwrap_calc_normal_projection));
+                       break;
+
+                       case MOD_SHRINKWRAP_NEAREST_VERTEX:
+                               BENCH(shrinkwrap_calc_nearest_vertex(&calc));
+//                             BENCH(shrinkwrap_calc_foreach_vertex(&calc, bruteforce_shrinkwrap_calc_nearest_vertex));
+                       break;
+               }
+
+       }
+
+       //Destroy faces, edges and stuff
+       if(calc.moved)
+       {
+               shrinkwrap_removeUnused(&calc);
+               bitset_free(calc.moved);
+       }
+
+       CDDM_calc_normals(calc.final);  
+
+       return calc.final;
+}
+
+
+/*
+ * Shrinkwrap to the nearest vertex
+ *
+ * it builds a kdtree of vertexs we can attach to and then
+ * for each vertex on performs a nearest vertex search on the tree
+ */
+void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
+{
+       int i;
+       int vgroup              = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
+       float tmp_co[3];
+
+       BVHTree *tree   = NULL;
+       BVHTreeNearest nearest;
+
+       BENCH_VAR(query);
+
+       int     numVerts;
+       MVert *vert = NULL;
+       MDeformVert *dvert = NULL;
+
+
+
+       BENCH(tree = bvhtree_from_mesh_verts(calc->target));
+       if(tree == NULL) return OUT_OF_MEMORY();
+
+       //Setup nearest
+       nearest.index = -1;
+       nearest.dist = FLT_MAX;
+
+
+       //Find the nearest vertex 
+       numVerts= calc->final->getNumVerts(calc->final);
+       vert    = calc->final->getVertDataArray(calc->final, CD_MVERT); 
+       dvert   = calc->final->getVertDataArray(calc->final, CD_MDEFORMVERT);
+
+       BENCH_BEGIN(query);
+       for(i=0; i<numVerts; i++)
+       {
+               int index;
+               float weight = vertexgroup_get_weight(dvert, i, vgroup);
+               if(weight == 0.0f) continue;
+
+               VecMat4MulVecfl(tmp_co, calc->local2target, vert[i].co);
+
+               if(nearest.index != -1)
+               {
+                       nearest.dist = squared_dist(tmp_co, nearest.nearest);
+               }
+               else nearest.dist = FLT_MAX;
+
+               index = BLI_bvhtree_find_nearest(tree, tmp_co, &nearest, NULL, NULL);
+
+               if(index != -1)
+               {
+                       float dist;
+
+                       VecMat4MulVecfl(tmp_co, calc->target2local, nearest.nearest);
+                       dist = VecLenf(vert[i].co, tmp_co);
+                       if(dist > 1e-5) weight *= (dist - calc->keptDist)/dist;
+                       VecLerpf(vert[i].co, vert[i].co, tmp_co, weight);       //linear interpolation
+               }
+       }
+       BENCH_END(query);
+       BENCH_REPORT(query);
+
+       BLI_bvhtree_free(tree);
+}
+
+/*
+ * Shrinkwrap projecting vertexs allong their normals over the target
+ *
+ * it builds a RayTree from the target mesh and then performs a
+ * raycast for each vertex (ray direction = normal)
+ */
+void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
+{
+       int i;
+       int vgroup              = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
+       char use_normal = calc->smd->shrinkOpts;
+       RayTree *target = NULL;
+
+       int     numVerts;
+       MVert *vert = NULL;
+       MDeformVert *dvert = NULL;
+       float tmp_co[3], tmp_no[3];
+
+       if( (use_normal & (MOD_SHRINKWRAP_ALLOW_INVERTED_NORMAL | MOD_SHRINKWRAP_ALLOW_DEFAULT_NORMAL)) == 0)
+               return; //Nothing todo
+
+       //setup raytracing
+       target = raytree_create_from_mesh(calc->target);
+       if(target == NULL) return OUT_OF_MEMORY();
+
+
+
+       //Project each vertex along normal
+       numVerts= calc->final->getNumVerts(calc->final);
+       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;
+               float weight = vertexgroup_get_weight(dvert, i, vgroup);
+               if(weight == 0.0f) continue;
+
+               //Transform coordinates local->target
+               VecMat4MulVecfl(tmp_co, calc->local2target, vert[i].co);
+
+               normal_short2float(vert[i].no, tmp_no);
+               Mat4Mul3Vecfl(calc->local2target, tmp_no);      //Watch out for scaling on normal
+               Normalize(tmp_no);                                                      //(TODO: do we really needed a unit-len normal? and we could know the scale factor before hand?)
+
+
+               if(use_normal & MOD_SHRINKWRAP_ALLOW_DEFAULT_NORMAL)
+               {
+                       dist = raytree_cast_ray(target, tmp_co, tmp_no);
+               }
+
+               normal_short2float(vert[i].no, tmp_no);
+               Mat4Mul3Vecfl(calc->local2target, tmp_no);      //Watch out for scaling on normal
+               Normalize(tmp_no);                                                      //(TODO: do we really needed a unit-len normal? and we could know the scale factor before hand?)
+
+               if(use_normal & MOD_SHRINKWRAP_ALLOW_INVERTED_NORMAL)
+               {
+                       float inv[3]; // = {-tmp_no[0], -tmp_no[1], -tmp_no[2]};
+                       float tdist;
+
+                       inv[0] = -tmp_no[0];
+                       inv[1] = -tmp_no[1];
+                       inv[2] = -tmp_no[2];
+
+                       tdist = raytree_cast_ray(target, tmp_co, inv);
+
+                       if(ABS(tdist) < ABS(dist))
+                               dist = -tdist;
+               }
+
+               if(ABS(dist) != FLT_MAX)
+               {
+                       float dist_t;
+
+                       VECADDFAC(tmp_co, tmp_co, tmp_no, dist);
+                       VecMat4MulVecfl(tmp_co, calc->target2local, tmp_co);
+
+                       dist_t = VecLenf(vert[i].co, tmp_co);
+                       if(dist_t > 1e-5) weight *= (dist_t - calc->keptDist)/dist_t;
+                       VecLerpf(vert[i].co, vert[i].co, tmp_co, weight);       //linear interpolation
+
+                       if(calc->moved)
+                               bitset_set(calc->moved, i);
+               }
+
+       }
+
+       free_raytree_from_mesh(target);
+}
+
+/*
+ * Shrinkwrap moving vertexs to the nearest surface point on the target
+ *
+ * it builds a BVHTree from the target mesh and then performs a
+ * NN matchs for each vertex
+ */
+void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc)
+{
+       int i;
+       int vgroup              = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
+       float tmp_co[3];
+
+       BVHTree *tree   = NULL;
+       BVHTreeNearest nearest;
+
+       int     numVerts;
+       MVert *vert = NULL;
+       MDeformVert *dvert = NULL;
+
+
+       //Create a bvh-tree of the given target
+       tree = bvhtree_from_mesh_tri(calc->target);
+       if(tree == NULL) return OUT_OF_MEMORY();
+
+       //Setup nearest
+       nearest.index = -1;
+       nearest.dist = FLT_MAX;
+
+
+       //Find the nearest vertex 
+       numVerts= calc->final->getNumVerts(calc->final);
+       vert    = calc->final->getVertDataArray(calc->final, CD_MVERT); 
+       dvert   = calc->final->getVertDataArray(calc->final, CD_MDEFORMVERT);
+
+       for(i=0; i<numVerts; i++)
+       {
+               int index;
+               float weight = vertexgroup_get_weight(dvert, i, vgroup);
+               if(weight == 0.0f) continue;
+
+               VecMat4MulVecfl(tmp_co, calc->local2target, vert[i].co);
+
+               if(nearest.index != -1)
+               {
+                       nearest.dist = squared_dist(tmp_co, nearest.nearest);
+               }
+               else nearest.dist = FLT_MAX;
+
+               index = BLI_bvhtree_find_nearest(tree, tmp_co, &nearest, mesh_tri_nearest_point, calc->target);
+
+               if(index != -1)
+               {
+                       float dist;
+
+                       VecMat4MulVecfl(tmp_co, calc->target2local, nearest.nearest);
+                       dist = VecLenf(vert[i].co, tmp_co);
+                       if(dist > 1e-5) weight *= (dist - calc->keptDist)/dist;
+                       VecLerpf(vert[i].co, vert[i].co, tmp_co, weight);       //linear interpolation
+               }
+       }
+
+       BLI_bvhtree_free(tree);
+}
+
diff --git a/source/blender/blenlib/BLI_kdopbvh.h b/source/blender/blenlib/BLI_kdopbvh.h
new file mode 100644 (file)
index 0000000..055c275
--- /dev/null
@@ -0,0 +1,77 @@
+/**
+ *
+ * ***** 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2006 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): Daniel Genrich, Andre Pinto
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+
+#ifndef BLI_KDOPBVH_H
+#define BLI_KDOPBVH_H
+
+#include <float.h>
+
+struct BVHTree;
+typedef struct BVHTree BVHTree;
+
+typedef struct BVHTreeOverlap {
+       int indexA;
+       int indexB;
+} BVHTreeOverlap;
+
+typedef struct BVHTreeNearest
+{
+       int index;                      /* the index of the nearest found (untouched if none is found within a dist radius from the given coordinates) */
+       float nearest[3];       /* nearest coordinates (untouched it none is found within a dist radius from the given coordinates) */
+       float dist;                     /* squared distance to search arround */
+} BVHTreeNearest;
+
+/* returns square of the minimum distance from given co to the node, nearest point is stored on nearest */
+typedef float (*BVHTree_NearestPointCallback) (void *userdata, int index, const float *co, float *nearest);
+
+
+BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis);
+void BLI_bvhtree_free(BVHTree *tree);
+
+/* construct: first insert points, then call balance */
+int BLI_bvhtree_insert(BVHTree *tree, int index, float *co, int numpoints);
+void BLI_bvhtree_balance(BVHTree *tree);
+
+/* update: first update points/nodes, then call update_tree to refit the bounding volumes */
+int BLI_bvhtree_update_node(BVHTree *tree, int index, float *co, float *co_moving, int numpoints);
+void BLI_bvhtree_update_tree(BVHTree *tree);
+
+/* collision/overlap: check two trees if they overlap, alloc's *overlap with length of the int return value */
+BVHTreeOverlap *BLI_bvhtree_overlap(BVHTree *tree1, BVHTree *tree2, int *result);
+
+float BLI_bvhtree_getepsilon(BVHTree *tree);
+
+/* find nearest node to the given coordinates (if nearest is given it will only search nodes where square distance is smaller than nearest->dist) */
+int BLI_bvhtree_find_nearest(BVHTree *tree, float *co, BVHTreeNearest *nearest, BVHTree_NearestPointCallback callback, void *userdata);
+
+#endif // BLI_KDOPBVH_H
+
+
+
+
diff --git a/source/blender/blenlib/intern/BLI_kdopbvh.c b/source/blender/blenlib/intern/BLI_kdopbvh.c
new file mode 100644 (file)
index 0000000..ef7c13a
--- /dev/null
@@ -0,0 +1,988 @@
+/*  kdop.c      
+* 
+*
+* ***** BEGIN GPL/BL DUAL 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. The Blender
+* Foundation also sells licenses for usinle in proprietary software under
+* the Blender License.  See http://www.blender.org/BL/ for information
+* about this.
+*
+* 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*
+* The Original Code is Copyright (C) Blender Foundation
+* All rights reserved.
+*
+* The Original Code is: all of this file.
+*
+* Contributor(s): Daniel Genrich, Andre Pinto
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+
+#include "math.h"
+#include <stdio.h>
+#include <stdlib.h> 
+#include <string.h>
+
+#include "MEM_guardedalloc.h"
+
+#include "BKE_utildefines.h"
+
+#include "BLI_kdopbvh.h"
+#include "BLI_arithb.h"
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#define node_get_bv(tree, node)                ((node)->bv)
+#define node_get_child(tree,node,i)    ((node)->children[i])
+
+/* Util macros */
+#define TO_STR(a)      #a
+#define JOIN(a,b)      a##b
+
+/* Benchmark macros */
+#if 1
+
+#define BENCH(a)       \
+       do {                    \
+               clock_t _clock_init = clock();  \
+               (a);                                                    \
+               printf("%s: %fms\n", #a, (float)(clock()-_clock_init)*1000/CLOCKS_PER_SEC);     \
+} while(0)
+
+#define BENCH_VAR(name)                clock_t JOIN(_bench_step,name) = 0, JOIN(_bench_total,name) = 0
+#define BENCH_BEGIN(name)      JOIN(_bench_step, name) = clock()
+#define BENCH_END(name)                JOIN(_bench_total,name) += clock() - JOIN(_bench_step,name)
+#define BENCH_RESET(name)      JOIN(_bench_total, name) = 0
+#define BENCH_REPORT(name)     printf("%s: %fms\n", TO_STR(name), JOIN(_bench_total,name)*1000.0f/CLOCKS_PER_SEC)
+
+#else
+
+#define BENCH(a)       (a)
+#define BENCH_VAR(name)
+#define BENCH_BEGIN(name)
+#define BENCH_END(name)
+#define BENCH_RESET(name)
+#define BENCH_REPORT(name)
+
+#endif
+
+
+typedef struct BVHNode
+{
+       struct BVHNode **children;// max 8 children
+       struct BVHNode *parent; // needed for bottom - top update
+       float *bv;                              // Bounding volume of all nodes, max 13 axis
+       int index;                              // face, edge, vertex index
+       char totnode;                   // how many nodes are used, used for speedup
+       char traversed;                 // how many nodes already traversed until this level?
+       char main_axis;                 // axis that was used to split childs
+} BVHNode;
+
+struct BVHTree
+{
+       BVHNode **nodes;
+       BVHNode *nodearray;             // pre-alloc branchs
+       BVHNode **nodechild;    // pre-alloc childs for nodes
+       float   *nodebv;                // pre-alloc bounding-volumes for nodes
+       float   epsilon;                // epslion is used for inflation of the k-dop
+       int     totleaf;                // leafs
+       int     totbranch;
+       char    tree_type;              // type of tree (4 => quadtree)
+       char    axis;                   // kdop type (6 => OBB, 7 => AABB, ...)
+       char    start_axis, stop_axis; // KDOP_AXES array indices according to axis     char    start_axis, stop_axis; // KDOP_AXES array indices according to axis
+};
+
+typedef struct BVHOverlapData 
+{  
+       BVHTree *tree1, *tree2; 
+       BVHTreeOverlap *overlap; 
+       int i, max_overlap; /* i is number of overlaps */
+} BVHOverlapData;
+
+typedef struct BVHNearestData
+{
+       BVHTree *tree;
+       float   *co;
+       BVHTree_NearestPointCallback callback;
+       void    *userdata;
+       float proj[13];                 //coordinates projection over axis
+       BVHTreeNearest nearest;
+} BVHNearestData;
+////////////////////////////////////////
+
+
+////////////////////////////////////////////////////////////////////////
+// Bounding Volume Hierarchy Definition
+// 
+// Notes: From OBB until 26-DOP --> all bounding volumes possible, just choose type below
+// Notes: You have to choose the type at compile time ITM
+// Notes: You can choose the tree type --> binary, quad, octree, choose below
+////////////////////////////////////////////////////////////////////////
+
+static float KDOP_AXES[13][3] =
+{ {1.0, 0, 0}, {0, 1.0, 0}, {0, 0, 1.0}, {1.0, 1.0, 1.0}, {1.0, -1.0, 1.0}, {1.0, 1.0, -1.0},
+{1.0, -1.0, -1.0}, {1.0, 1.0, 0}, {1.0, 0, 1.0}, {0, 1.0, 1.0}, {1.0, -1.0, 0}, {1.0, 0, -1.0},
+{0, 1.0, -1.0}
+};
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////
+// Introsort 
+// with permission deriven from the following Java code:
+// http://ralphunden.net/content/tutorials/a-guide-to-introsort/
+// and he derived it from the SUN STL 
+//////////////////////////////////////////////////////////////////////////////////////////////////////
+static int size_threshold = 16;
+/*
+* Common methods for all algorithms
+*/
+static int floor_lg(int a)
+{
+       return (int)(floor(log(a)/log(2)));
+}
+
+/*
+* Insertion sort algorithm
+*/
+static void bvh_insertionsort(BVHNode **a, int lo, int hi, int axis)
+{
+       int i,j;
+       BVHNode *t;
+       for (i=lo; i < hi; i++)
+       {
+               j=i;
+               t = a[i];
+               while((j!=lo) && (t->bv[axis] < (a[j-1])->bv[axis]))
+               {
+                       a[j] = a[j-1];
+                       j--;
+               }
+               a[j] = t;
+       }
+}
+
+static int bvh_partition(BVHNode **a, int lo, int hi, BVHNode * x, int axis)
+{
+       int i=lo, j=hi;
+       while (1)
+       {
+               while ((a[i])->bv[axis] < x->bv[axis]) i++;
+               j--;
+               while (x->bv[axis] < (a[j])->bv[axis]) j--;
+               if(!(i < j))
+                       return i;
+               SWAP( BVHNode* , a[i], a[j]);
+               i++;
+       }
+}
+
+/*
+* Heapsort algorithm
+*/
+static void bvh_downheap(BVHNode **a, int i, int n, int lo, int axis)
+{
+       BVHNode * d = a[lo+i-1];
+       int child;
+       while (i<=n/2)
+       {
+               child = 2*i;
+               if ((child < n) && ((a[lo+child-1])->bv[axis] < (a[lo+child])->bv[axis]))
+               {
+                       child++;
+               }
+               if (!(d->bv[axis] < (a[lo+child-1])->bv[axis])) break;
+               a[lo+i-1] = a[lo+child-1];
+               i = child;
+       }
+       a[lo+i-1] = d;
+}
+
+static void bvh_heapsort(BVHNode **a, int lo, int hi, int axis)
+{
+       int n = hi-lo, i;
+       for (i=n/2; i>=1; i=i-1)
+       {
+               bvh_downheap(a, i,n,lo, axis);
+       }
+       for (i=n; i>1; i=i-1)
+       {
+               SWAP(BVHNode*, a[lo],a[lo+i-1]);
+               bvh_downheap(a, 1,i-1,lo, axis);
+       }
+}
+
+static BVHNode *bvh_medianof3(BVHNode **a, int lo, int mid, int hi, int axis) // returns Sortable
+{
+       if ((a[mid])->bv[axis] < (a[lo])->bv[axis])
+       {
+               if ((a[hi])->bv[axis] < (a[mid])->bv[axis])
+                       return a[mid];
+               else
+               {
+                       if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
+                               return a[hi];
+                       else
+                               return a[lo];
+               }
+       }
+       else
+       {
+               if ((a[hi])->bv[axis] < (a[mid])->bv[axis])
+               {
+                       if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
+                               return a[lo];
+                       else
+                               return a[hi];
+               }
+               else
+                       return a[mid];
+       }
+}
+/*
+* Quicksort algorithm modified for Introsort
+*/
+static void bvh_introsort_loop (BVHNode **a, int lo, int hi, int depth_limit, int axis)
+{
+       int p;
+
+       while (hi-lo > size_threshold)
+       {
+               if (depth_limit == 0)
+               {
+                       bvh_heapsort(a, lo, hi, axis);
+                       return;
+               }
+               depth_limit=depth_limit-1;
+               p=bvh_partition(a, lo, hi, bvh_medianof3(a, lo, lo+((hi-lo)/2)+1, hi-1, axis), axis);
+               bvh_introsort_loop(a, p, hi, depth_limit, axis);
+               hi=p;
+       }
+}
+
+static void sort(BVHNode **a0, int begin, int end, int axis)
+{
+       if (begin < end)
+       {
+               BVHNode **a=a0;
+               bvh_introsort_loop(a, begin, end, 2*floor_lg(end-begin), axis);
+               bvh_insertionsort(a, begin, end, axis);
+       }
+}
+void sort_along_axis(BVHTree *tree, int start, int end, int axis)
+{
+       sort(tree->nodes, start, end, axis);
+}
+
+//after a call to this function you can expect one of:
+//      every node to left of a[n] are smaller or equal to it
+//      every node to the right of a[n] are greater or equal to it
+int partition_nth_element(BVHNode **a, int _begin, int _end, int n, int axis){        
+       int begin = _begin, end = _end, cut;        
+       while(end-begin > 3)        
+       {                            
+               cut = bvh_partition(a, begin, end, bvh_medianof3(a, begin, (begin+end)/2, end-1, axis), axis );                 
+               if(cut <= n)                        
+                       begin = cut;                
+               else                        
+                       end = cut;        
+       }        
+       bvh_insertionsort(a, begin, end, axis);
+
+       return n;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void BLI_bvhtree_free(BVHTree *tree)
+{      
+       if(tree)
+       {
+               MEM_freeN(tree->nodes);
+               MEM_freeN(tree->nodearray);
+               MEM_freeN(tree->nodebv);
+               MEM_freeN(tree->nodechild);
+               MEM_freeN(tree);
+       }
+}
+
+BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis)
+{
+       BVHTree *tree;
+       int numbranches=0, i;
+       
+       // only support up to octree
+       if(tree_type > 8)
+               return NULL;
+
+       tree = (BVHTree *)MEM_callocN(sizeof(BVHTree), "BVHTree");
+       
+       if(tree)
+       {
+               tree->epsilon = epsilon;
+               tree->tree_type = tree_type; 
+               tree->axis = axis;
+               
+               if(axis == 26)
+               {
+                       tree->start_axis = 0;
+                       tree->stop_axis = 13;
+               }
+               else if(axis == 18)
+               {
+                       tree->start_axis = 7;
+                       tree->stop_axis = 13;
+               }
+               else if(axis == 14)
+               {
+                       tree->start_axis = 0;
+                       tree->stop_axis = 7;
+               }
+               else if(axis == 8)
+               {
+                       tree->start_axis = 0;
+                       tree->stop_axis = 4;
+               }
+               else if(axis == 6) // AABB
+               {
+                       tree->start_axis = 0;
+                       tree->stop_axis = 3;
+               }
+               else
+               {
+                       MEM_freeN(tree);
+                       return NULL;
+               }
+
+
+               // calculate max number of branches, our bvh kdop is "almost perfect"
+               for(i = 1; i <= (int)ceil((float)((float)log(maxsize)/(float)log(tree_type))); i++)
+                       numbranches += (pow(tree_type, i) / tree_type);
+               
+               tree->nodes = (BVHNode **)MEM_callocN(sizeof(BVHNode *)*(numbranches+maxsize + tree_type), "BVHNodes");
+               
+               if(!tree->nodes)
+               {
+                       MEM_freeN(tree);
+                       return NULL;
+               }
+               
+               tree->nodebv = (float*)MEM_callocN(sizeof(float)* axis * (numbranches+maxsize + tree_type), "BVHNodeBV");
+               if(!tree->nodebv)
+               {
+                       MEM_freeN(tree->nodes);
+                       MEM_freeN(tree);
+               }
+
+               tree->nodechild = (BVHNode**)MEM_callocN(sizeof(BVHNode*) * tree_type * (numbranches+maxsize + tree_type), "BVHNodeBV");
+               if(!tree->nodechild)
+               {
+                       MEM_freeN(tree->nodebv);
+                       MEM_freeN(tree->nodes);
+                       MEM_freeN(tree);
+               }
+
+               tree->nodearray = (BVHNode *)MEM_callocN(sizeof(BVHNode)*(numbranches+maxsize + tree_type), "BVHNodeArray");
+               
+               if(!tree->nodearray)
+               {
+                       MEM_freeN(tree->nodechild);
+                       MEM_freeN(tree->nodebv);
+                       MEM_freeN(tree->nodes);
+                       MEM_freeN(tree);
+                       return NULL;
+               }
+
+               //link the dynamic bv and child links
+               for(i=0; i< numbranches+maxsize + tree_type; i++)
+               {
+                       tree->nodearray[i].bv = tree->nodebv + i * axis;
+                       tree->nodearray[i].children = tree->nodechild + i * tree_type;
+               }
+               
+       }
+
+       return tree;
+}
+
+
+static void create_kdop_hull(BVHTree *tree, BVHNode *node, float *co, int numpoints, int moving)
+{
+       float newminmax;
+       float *bv = node_get_bv(tree, node);
+       int i, k;
+       
+       // don't init boudings for the moving case
+       if(!moving)
+       {
+               for (i = tree->start_axis; i < tree->stop_axis; i++)
+               {
+                       bv[2*i] = FLT_MAX;
+                       bv[2*i + 1] = -FLT_MAX;
+               }
+       }
+       
+       for(k = 0; k < numpoints; k++)
+       {
+               // for all Axes.
+               for (i = tree->start_axis; i < tree->stop_axis; i++)
+               {
+                       newminmax = INPR(&co[k * 3], KDOP_AXES[i]);
+                       if (newminmax < bv[2 * i])
+                               bv[2 * i] = newminmax;
+                       if (newminmax > bv[(2 * i) + 1])
+                               bv[(2 * i) + 1] = newminmax;
+               }
+       }
+}
+
+// depends on the fact that the BVH's for each face is already build
+static void refit_kdop_hull(BVHTree *tree, BVHNode *node, int start, int end)
+{
+       float newmin,newmax;
+       int i, j;
+       float *bv = node_get_bv(tree, node);
+       
+       for (i = tree->start_axis; i < tree->stop_axis; i++)
+       {
+               bv[2*i] = FLT_MAX;
+               bv[2*i + 1] = -FLT_MAX;
+       }
+
+       for (j = start; j < end; j++)
+       {
+        // for all Axes.
+               for (i = tree->start_axis; i < tree->stop_axis; i++)
+               {
+                       newmin = tree->nodes[j]->bv[(2 * i)];   
+                       if ((newmin < bv[(2 * i)]))
+                               bv[(2 * i)] = newmin;
+                       newmax = tree->nodes[j]->bv[(2 * i) + 1];
+                       if ((newmax > bv[(2 * i) + 1]))
+                               bv[(2 * i) + 1] = newmax;
+               }
+       }
+}
+
+static void inflate_kdop_hull(BVHTree *tree, BVHNode *node)
+{
+       int i;
+       float *bv = node_get_bv(tree, node);
+
+       for(i= tree->axis; i; i--)
+       {
+               *bv++ -= tree->epsilon; // minimum 
+               *bv++ += tree->epsilon; // maximum 
+       }
+}
+
+// only supports x,y,z axis in the moment
+// but we should use a plain and simple function here for speed sake
+static char kdop_hull_get_largest_axis(BVHTree *tree, BVHNode *node)
+{
+       float *bv = node_get_bv(tree, node);
+       float middle_point[3];
+
+       middle_point[0] = (bv[1]) - (bv[0]); // x axis
+       middle_point[1] = (bv[3]) - (bv[2]); // y axis
+       middle_point[2] = (bv[5]) - (bv[4]); // z axis
+       if (middle_point[0] > middle_point[1]) 
+       {
+               if (middle_point[0] > middle_point[2])
+                       return 1; // max x axis
+               else
+                       return 5; // max z axis
+       }
+       else 
+       {
+               if (middle_point[1] > middle_point[2])
+                       return 3; // max y axis
+               else
+                       return 5; // max z axis
+       }
+}
+
+
+int BLI_bvhtree_insert(BVHTree *tree, int index, float *co, int numpoints)
+{
+       BVHNode *node= NULL;
+       
+       // insert should only possible as long as tree->totbranch is 0
+       if(tree->totbranch > 0)
+               return 0;
+       
+       if(tree->totleaf+1 >= MEM_allocN_len(tree->nodes))
+               return 0;
+       
+       // TODO check if have enough nodes in array
+       
+       node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]);
+       tree->totleaf++;
+       
+       //create and inflate hull
+       create_kdop_hull(tree, node, co, numpoints, 0);
+       inflate_kdop_hull(tree, node);
+
+       node->index= index;
+       
+       return 1;
+}
+
+
+static void bvh_div_nodes(BVHTree *tree, BVHNode *node, int start, int end, char lastaxis)
+{
+       char laxis;
+       int i, tend;
+       BVHNode *tnode;
+       int slice = (end-start+tree->tree_type-1)/tree->tree_type;      //division rounded up
+       
+       // Determine which axis to split along
+       laxis = kdop_hull_get_largest_axis(tree, node);
+       
+       // split nodes along longest axis
+       for (i=0; start < end; start += slice, i++) //i counts the current child
+       {       
+               tend = start + slice;
+               
+               if(tend > end) tend = end;
+               
+               if(tend-start == 1)     // ok, we have 1 left for this node
+               {
+                       node->children[i] = tree->nodes[start];
+                       node->children[i]->parent = node;
+               }
+               else
+               {
+                       partition_nth_element(tree->nodes, start, end, tend, laxis);
+
+                       tnode = node->children[i] = tree->nodes[tree->totleaf  + tree->totbranch] = &(tree->nodearray[tree->totbranch + tree->totleaf]);
+                       tree->totbranch++;
+                       tnode->parent = node;
+                       
+                       if(tend != end)
+                               partition_nth_element(tree->nodes, start, end, tend, laxis);
+                       refit_kdop_hull(tree, tnode, start, tend);
+                       bvh_div_nodes(tree, tnode, start, tend, laxis);
+               }
+               node->totnode++;
+       }
+       
+       return;
+}
+
+static void verify_tree(BVHTree *tree)
+{
+       int i, j, check = 0;
+       
+       // check the pointer list
+       for(i = 0; i < tree->totleaf; i++)
+       {
+               if(tree->nodes[i]->parent == NULL)
+                       printf("Leaf has no parent: %d\n", i);
+               else
+               {
+                       for(j = 0; j < tree->tree_type; j++)
+                       {
+                               if(tree->nodes[i]->parent->children[j] == tree->nodes[i])
+                                       check = 1;
+                       }
+                       if(!check)
+                       {
+                               printf("Parent child relationship doesn't match: %d\n", i);
+                       }
+                       check = 0;
+               }
+       }
+       
+       // check the leaf list
+       for(i = 0; i < tree->totleaf; i++)
+       {
+               if(tree->nodearray[i].parent == NULL)
+                       printf("Leaf has no parent: %d\n", i);
+               else
+               {
+                       for(j = 0; j < tree->tree_type; j++)
+                       {
+                               if(tree->nodearray[i].parent->children[j] == &tree->nodearray[i])
+                                       check = 1;
+                       }
+                       if(!check)
+                       {
+                               printf("Parent child relationship doesn't match: %d\n", i);
+                       }
+                       check = 0;
+               }
+       }
+       
+       printf("branches: %d, leafs: %d, total: %d\n", tree->totbranch, tree->totleaf, tree->totbranch + tree->totleaf);
+}
+       
+void BLI_bvhtree_balance(BVHTree *tree)
+{
+       BVHNode *node;
+       
+       if(tree->totleaf == 0)
+               return;
+       
+       // create root node
+       node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]);
+       tree->totbranch++;
+       
+       // refit root bvh node
+       refit_kdop_hull(tree, tree->nodes[tree->totleaf], 0, tree->totleaf);
+       // create + balance tree
+       bvh_div_nodes(tree, tree->nodes[tree->totleaf], 0, tree->totleaf, 0);
+       
+       // verify_tree(tree);
+}
+
+// overlap - is it possbile for 2 bv's to collide ?
+static int tree_overlap(BVHNode *node1, BVHNode *node2, int start_axis, int stop_axis)
+{
+       float *bv1 = node1->bv;
+       float *bv2 = node2->bv;
+
+       float *bv1_end = bv1 + (stop_axis<<1);
+               
+       bv1 += start_axis<<1;
+       bv2 += start_axis<<1;
+       
+       // test all axis if min + max overlap
+       for (; bv1 != bv1_end; bv1+=2, bv2+=2)
+       {
+               if ((*(bv1) > *(bv2 + 1)) || (*(bv2) > *(bv1 + 1))) 
+                       return 0;
+       }
+       
+       return 1;
+}
+
+static void traverse(BVHOverlapData *data, BVHNode *node1, BVHNode *node2)
+{
+       int j;
+       
+       if(tree_overlap(node1, node2, MIN2(data->tree1->start_axis, data->tree2->start_axis), MIN2(data->tree1->stop_axis, data->tree2->stop_axis)))
+       {
+               // check if node1 is a leaf
+               if(node1->index)
+               {
+                       // check if node2 is a leaf
+                       if(node2->index)
+                       {
+                               
+                               if(node1 == node2)
+                               {
+                                       return;
+                               }
+                                       
+                               if(data->i >= data->max_overlap)
+                               {       
+                                       // try to make alloc'ed memory bigger
+                                       data->overlap = realloc(data->overlap, sizeof(BVHTreeOverlap)*data->max_overlap*2);
+                                       
+                                       if(!data->overlap)
+                                       {
+                                               printf("Out of Memory in traverse\n");
+                                               return;
+                                       }
+                                       data->max_overlap *= 2;
+                               }
+                               
+                               // both leafs, insert overlap!
+                               data->overlap[data->i].indexA = node1->index;
+                               data->overlap[data->i].indexB = node2->index;
+
+                               data->i++;
+                       }
+                       else
+                       {
+                               for(j = 0; j < data->tree2->tree_type; j++)
+                               {
+                                       if(node2->children[j])
+                                               traverse(data, node1, node2->children[j]);
+                               }
+                       }
+               }
+               else
+               {
+                       
+                       for(j = 0; j < data->tree2->tree_type; j++)
+                       {
+                               if(node1->children[j])
+                                       traverse(data, node1->children[j], node2);
+                       }
+               }
+       }
+       return;
+}
+
+BVHTreeOverlap *BLI_bvhtree_overlap(BVHTree *tree1, BVHTree *tree2, int *result)
+{
+       int j, total = 0;
+       BVHTreeOverlap *overlap = NULL, *to = NULL;
+       BVHOverlapData **data;
+       
+       // check for compatibility of both trees (can't compare 14-DOP with 18-DOP)
+       if((tree1->axis != tree2->axis) && ((tree1->axis == 14) || tree2->axis == 14))
+               return 0;
+       
+       // fast check root nodes for collision before doing big splitting + traversal
+       if(!tree_overlap(tree1->nodes[tree1->totleaf], tree2->nodes[tree2->totleaf], MIN2(tree1->start_axis, tree2->start_axis), MIN2(tree1->stop_axis, tree2->stop_axis)))
+               return 0;
+
+       data = MEM_callocN(sizeof(BVHOverlapData *)* tree1->tree_type, "BVHOverlapData_star");
+       
+       for(j = 0; j < tree1->tree_type; j++)
+       {
+               data[j] = (BVHOverlapData *)MEM_callocN(sizeof(BVHOverlapData), "BVHOverlapData");
+               
+               // init BVHOverlapData
+               data[j]->overlap = (BVHTreeOverlap *)malloc(sizeof(BVHTreeOverlap)*MAX2(tree1->totleaf, tree2->totleaf));
+               data[j]->tree1 = tree1;
+               data[j]->tree2 = tree2;
+               data[j]->max_overlap = MAX2(tree1->totleaf, tree2->totleaf);
+               data[j]->i = 0;
+       }
+
+#pragma omp parallel for private(j) schedule(static)
+       for(j = 0; j < tree1->tree_type; j++)
+       {
+               traverse(data[j], tree1->nodes[tree1->totleaf]->children[j], tree2->nodes[tree2->totleaf]);
+       }
+       
+       for(j = 0; j < tree1->tree_type; j++)
+               total += data[j]->i;
+       
+       to = overlap = (BVHTreeOverlap *)MEM_callocN(sizeof(BVHTreeOverlap)*total, "BVHTreeOverlap");
+       
+       for(j = 0; j < tree1->tree_type; j++)
+       {
+               memcpy(to, data[j]->overlap, data[j]->i*sizeof(BVHTreeOverlap));
+               to+=data[j]->i;
+       }
+       
+       for(j = 0; j < tree1->tree_type; j++)
+       {
+               free(data[j]->overlap);
+               MEM_freeN(data[j]);
+       }
+       MEM_freeN(data);
+       
+       (*result) = total;
+       return overlap;
+}
+
+
+// bottom up update of bvh tree:
+// join the 4 children here
+static void node_join(BVHTree *tree, BVHNode *node)
+{
+       int i, j;
+       
+       for (i = tree->start_axis; i < tree->stop_axis; i++)
+       {
+               node->bv[2*i] = FLT_MAX;
+               node->bv[2*i + 1] = -FLT_MAX;
+       }
+       
+       for (i = 0; i < tree->tree_type; i++)
+       {
+               if (node->children[i]) 
+               {
+                       for (j = tree->start_axis; j < tree->stop_axis; j++)
+                       {
+                               // update minimum 
+                               if (node->children[i]->bv[(2 * j)] < node->bv[(2 * j)]) 
+                                       node->bv[(2 * j)] = node->children[i]->bv[(2 * j)];
+                               
+                               // update maximum 
+                               if (node->children[i]->bv[(2 * j) + 1] > node->bv[(2 * j) + 1])
+                                       node->bv[(2 * j) + 1] = node->children[i]->bv[(2 * j) + 1];
+                       }
+               }
+               else
+                       break;
+       }
+}
+
+// call before BLI_bvhtree_update_tree()
+int BLI_bvhtree_update_node(BVHTree *tree, int index, float *co, float *co_moving, int numpoints)
+{
+       BVHNode *node= NULL;
+       
+       // check if index exists
+       if(index > tree->totleaf)
+               return 0;
+       
+       node = tree->nodearray + index;
+       
+       create_kdop_hull(tree, node, co, numpoints, 0);
+       
+       if(co_moving)
+               create_kdop_hull(tree, node, co_moving, numpoints, 1);
+       
+       // inflate the bv with some epsilon
+       inflate_kdop_hull(tree, node);
+       
+       return 1;
+}
+
+// call BLI_bvhtree_update_node() first for every node/point/triangle
+void BLI_bvhtree_update_tree(BVHTree *tree)
+{
+       BVHNode *leaf, *parent;
+       
+       // reset tree traversing flag
+       for (leaf = tree->nodearray + tree->totleaf; leaf != tree->nodearray + tree->totleaf + tree->totbranch; leaf++)
+               leaf->traversed = 0;
+       
+       for (leaf = tree->nodearray; leaf != tree->nodearray + tree->totleaf; leaf++)
+       {
+               for (parent = leaf->parent; parent; parent = parent->parent)
+               {
+                       parent->traversed++;    // we tried to go up in hierarchy 
+                       if (parent->traversed < parent->totnode) 
+                               break;  // we do not need to check further 
+                       else 
+                               node_join(tree, parent);
+               }
+       }
+}
+
+float BLI_bvhtree_getepsilon(BVHTree *tree)
+{
+       return tree->epsilon;
+}
+
+
+//Nearest neighbour
+static float squared_dist(const float *a, const float *b)
+{
+       float tmp[3];
+       VECSUB(tmp, a, b);
+       return INPR(tmp, tmp);
+}
+
+static float calc_nearest_point(BVHNearestData *data, BVHNode *node, float *nearest)
+{
+       int i;
+       const float *bv = node_get_bv(data->tree, node);
+
+       //nearest on AABB hull
+       for(i=0; i != 3; i++, bv += 2)
+       {
+               if(bv[0] > data->proj[i])
+                       nearest[i] = bv[0];
+               else if(bv[1] < data->proj[i])
+                       nearest[i] = bv[1];
+               else
+                       nearest[i] = data->proj[i];
+       }
+
+/*
+       //nearest on a general hull
+       VECCOPY(nearest, data->co);
+       for(i = data->tree->start_axis; i != data->tree->stop_axis; i++, bv+=2)
+       {
+               float proj = INPR( nearest, KDOP_AXES[i]);
+               float dl = bv[0] - proj;
+               float du = bv[1] - proj;
+
+               if(dl > 0)
+               {
+                       VECADDFAC(nearest, nearest, KDOP_AXES[i], dl);
+               }
+               else if(du < 0)
+               {
+                       VECADDFAC(nearest, nearest, KDOP_AXES[i], du);
+               }
+       }
+*/
+       return squared_dist(data->co, nearest);
+}
+
+
+static void dfs_find_nearest(BVHNearestData *data, BVHNode *node)
+{
+       int i;
+       float nearest[3], sdist;
+
+       sdist = calc_nearest_point(data, node, nearest);
+       if(sdist >= data->nearest.dist) return;
+
+       if(node->totnode == 0)
+       {
+               if(data->callback)
+                       sdist = data->callback(data->userdata , node->index, data->co, nearest);
+
+               if(sdist >= data->nearest.dist) return;
+
+               data->nearest.index     = node->index;
+               VECCOPY(data->nearest.nearest, nearest);
+               data->nearest.dist      = sdist;
+       }
+       else
+       {
+               if(sdist < data->nearest.dist)
+               {
+                       for(i=0; i != node->totnode; i++)
+                       {
+                               dfs_find_nearest(data, node->children[i]);
+                       }
+               }
+       }
+}
+
+int BLI_bvhtree_find_nearest(BVHTree *tree, float *co, BVHTreeNearest *nearest, BVHTree_NearestPointCallback callback, void *userdata)
+{
+       int i;
+
+       BVHNearestData data;
+
+       //init data to search
+       data.tree = tree;
+       data.co = co;
+
+       data.callback = callback;
+       data.userdata = userdata;
+
+       for(i = data.tree->start_axis; i != data.tree->stop_axis; i++)
+       {
+               data.proj[i] = INPR(data.co, KDOP_AXES[i]);
+       }
+
+       if(nearest)
+       {
+               memcpy( &data.nearest , nearest, sizeof(*nearest) );
+       }
+       else
+       {
+               data.nearest.index = -1;
+               data.nearest.dist = FLT_MAX;
+       }
+
+       //dfs search
+       dfs_find_nearest(&data, tree->nodes[tree->totleaf] );
+
+       //copy back results
+       if(nearest)
+       {
+               memcpy(nearest, &data.nearest, sizeof(*nearest));
+       }
+
+       return data.nearest.index;
+}
+
index b7b438174747b0a7d8be78c38c81e9932826438a..5178022ffad97c89ade7a3d0f92a33675ea9b8f7 100644 (file)
@@ -35,6 +35,7 @@ typedef enum ModifierType {
        eModifierType_Cloth,
        eModifierType_Collision,
        eModifierType_Bevel,
+       eModifierType_Shrinkwrap,
        NUM_MODIFIER_TYPES
 } ModifierType;
 
@@ -488,4 +489,24 @@ typedef struct ExplodeModifierData {
        float protect;
 } ExplodeModifierData;
 
+typedef struct ShrinkwrapModifierData {
+       ModifierData modifier;
+
+       struct Object *target;  /* shrink target */
+       char vgroup_name[32];   /* optional vertexgroup name */
+       float keptDist;                 /* distance offset from mesh/projection point */
+       short shrinkType;               /* shrink type projection */
+       short shrinkOpts;               /* shrink options */
+} ShrinkwrapModifierData;
+
+/* Shrinkwrap->shrinkType */
+#define MOD_SHRINKWRAP_NEAREST_SURFACE 0
+#define MOD_SHRINKWRAP_NORMAL                  1
+#define MOD_SHRINKWRAP_NEAREST_VERTEX  2
+
+/* 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 1365baf075a18a2bb18c3a9371d8c5c00be85aa7..b1d2d26bc4628a78fbdfc77f2521892201e733d2 100644 (file)
@@ -1826,6 +1826,11 @@ static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco
                        height = 94;
                } else if (md->type==eModifierType_Explode) {
                        height = 94;
+               } else if (md->type==eModifierType_Shrinkwrap) {
+                       ShrinkwrapModifierData *smd = (ShrinkwrapModifierData*) md;
+                       height = 86;
+                       if (smd->shrinkType == MOD_SHRINKWRAP_NORMAL)
+                               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, ""); 
@@ -2443,6 +2448,24 @@ static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco
                        uiDefButBitS(block, TOG, eExplodeFlag_Alive, B_MODIFIER_RECALC, "Alive",        lx+buttonWidth/3, cy, buttonWidth/3,19, &emd->flag, 0, 0, 0, 0, "Show mesh when particles are alive");
                        uiDefButBitS(block, TOG, eExplodeFlag_Dead, B_MODIFIER_RECALC, "Dead",  lx+buttonWidth*2/3, cy, buttonWidth/3,19, &emd->flag, 0, 0, 0, 0, "Show mesh when particles are dead");
                        uiBlockEndAlign(block);
+               } else if (md->type==eModifierType_Shrinkwrap) {
+                       ShrinkwrapModifierData *smd = (ShrinkwrapModifierData*) md;
+
+                       char shrinktypemenu[]="Shrinkwrap type%t|nearest surface point %x0|normal projection %x1|nearest vertex %x2";
+                       uiDefButS(block, MENU, B_MODIFIER_RECALC, shrinktypemenu, lx,(cy-=19),buttonWidth,19, &smd->shrinkType, 0, 0, 0, 0, "Selects type of shrinkwrap algorithm for target position.");
+
+                       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");
+                       uiButSetCompleteFunc(but, autocomplete_vgroup, (void *)ob);
+
+                       uiDefIDPoinBut(block, modifier_testMeshObj, ID_OB, B_CHANGEDEP, "Ob: ", lx, (cy-=19), buttonWidth,19, &smd->target, "Target to shrink to");
+                       uiDefButF(block, NUM, B_MODIFIER_RECALC, "Offset:",     lx,(cy-=19),buttonWidth,19, &smd->keptDist, 0.0f, 100.0f, 1.0f, 0, "Specify distance to kept from the target");
+                       uiBlockEndAlign(block);
                }
 
                uiBlockEndAlign(block);