Merged shrinkwrap modifier from soc-2008-jaguarandi

9 files changed, 1207 insertions, 36 deletions

diff --git a/source/blender/blenkernel/BKE_deform.h b/source/blender/blenkernel/BKE_deform.h
index a1975dd4265..73a9b2b5d4e 100644
--- a/source/blender/blenkernel/BKE_deform.h
+++ b/source/blender/blenkernel/BKE_deform.h
@@ -46,5 +46,8 @@ int get_defgroup_num (struct Object *ob, struct bDeformGroup *dg);
 int get_named_vertexgroup_num (Object *ob, char *name);
 void unique_vertexgroup_name (struct bDeformGroup *dg, struct Object *ob);
 
+float deformvert_get_weight(const struct MDeformVert *dvert, int group_num);
+float vertexgroup_get_vertex_weight(const struct MDeformVert *dvert, int index, int group_num);
+
 #endif
 
diff --git a/source/blender/blenkernel/BKE_shrinkwrap.h b/source/blender/blenkernel/BKE_shrinkwrap.h
new file mode 100644
index 00000000000..e8276238ff2
--- /dev/null
+++ b/source/blender/blenkernel/BKE_shrinkwrap.h
@@ -0,0 +1,146 @@
+/**
+ * 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
+
+/* mesh util */
+//TODO: move this somewhere else
+#include "BKE_customdata.h"
+struct DerivedMesh;
+struct Object;
+struct DerivedMesh *object_get_derived_final(struct Object *ob, CustomDataMask dataMask);
+
+
+/* SpaceTransform stuff */
+/*
+ * TODO: move this somewhere else
+ *
+ * this structs encapsulates all needed data to convert between 2 coordinate spaces
+ * (where conversion can be represented by a matrix multiplication)
+ *
+ * This is used to reduce the number of arguments to pass to functions that need to perform
+ * this kind of operation and make it easier for the coder, as he/she doenst needs to recode
+ * the matrix calculation.
+ *
+ * A SpaceTransform is initialized using:
+ *   space_transform_setup( &data,  ob1, ob2 )
+ *
+ * After that the following calls can be used:
+ *   space_transform_apply (&data, co); //converts a coordinate in ob1 coords space to the corresponding ob2 coords
+ *   space_transform_invert(&data, co); //converts a coordinate in ob2 coords space to the corresponding ob1 coords
+ *
+ *	//Same Concept as space_transform_apply and space_transform_invert, but no is normalized after conversion
+ *   space_transform_apply_normal (&data, &no);
+ *   space_transform_invert_normal(&data, &no);
+ *
+ */
+struct Object;
+
+typedef struct SpaceTransform
+{
+	float local2target[4][4];
+	float target2local[4][4];
+
+} SpaceTransform;
+
+void space_transform_from_matrixs(SpaceTransform *data, float local[][4], float target[][4]);
+#define space_transform_setup(data, local, target) space_transform_from_matrixs(data, (local)->obmat, (target)->obmat)
+
+void space_transform_apply (const SpaceTransform *data, float *co);
+void space_transform_invert(const SpaceTransform *data, float *co);
+
+void space_transform_apply_normal (const SpaceTransform *data, float *no);
+void space_transform_invert_normal(const SpaceTransform *data, float *no);
+
+/* Shrinkwrap stuff */
+#include "BKE_bvhutils.h"
+
+/*
+ * Shrinkwrap is composed by a set of functions and options that define the type of shrink.
+ *
+ * 3 modes are available:
+ *    - Nearest vertex
+ *	  - Nearest surface
+ *    - Normal projection
+ *
+ * ShrinkwrapCalcData encapsulates all needed data for shrinkwrap functions.
+ * (So that you dont have to pass an enormous ammount of arguments to functions)
+ */
+
+struct Object;
+struct DerivedMesh;
+struct ShrinkwrapModifierData;
+struct BVHTree;
+
+
+typedef struct ShrinkwrapCalcData
+{
+	ShrinkwrapModifierData *smd;	//shrinkwrap modifier data
+
+	struct Object *ob;				//object we are applying shrinkwrap to
+	struct DerivedMesh *original;	//mesh before shrinkwrap
+
+	float (*vertexCos)[3];			//vertexs being shrinkwraped
+	int numVerts;
+
+	struct DerivedMesh *target;		//mesh we are shrinking to	
+	SpaceTransform local2target;	//transform to move bettwem local and target space
+
+	float keepDist;					//Distance to kept from target (units are in local space)
+
+} ShrinkwrapCalcData;
+
+void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *data);
+void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *data);
+void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *data);
+
+void shrinkwrapModifier_deform(struct ShrinkwrapModifierData *smd, struct Object *ob, struct DerivedMesh *dm, float (*vertexCos)[3], int numVerts);
+
+/*
+ * This function casts a ray in the given BVHTree.. but it takes into consideration the space_transform, that is:
+ *
+ * if transf was configured with "space_transform_setup( &transf,  ob1, ob2 )"
+ * then the input (vert, dir, BVHTreeRayHit) must be defined in ob1 coordinates space
+ * and the BVHTree must be built in ob2 coordinate space.
+ *
+ * Thus it provides an easy way to cast the same ray across several trees (where each tree was built on its own coords space)
+ */
+int normal_projection_project_vertex(char options, const float *vert, const float *dir, const SpaceTransform *transf, BVHTree *tree, BVHTreeRayHit *hit, BVHTree_RayCastCallback callback, void *userdata);
+
+/*
+ * NULL initializers to local data
+ */
+#define NULL_ShrinkwrapCalcData	{NULL, }
+#define NULL_BVHTreeFromMesh	{NULL, }
+#define NULL_BVHTreeRayHit		{NULL, }
+#define NULL_BVHTreeNearest		{0, }
+
+
+#endif
+
diff --git a/source/blender/blenkernel/intern/bvhutils.c b/source/blender/blenkernel/intern/bvhutils.c
index 042e2afad53..ae449843d2a 100644
--- a/source/blender/blenkernel/intern/bvhutils.c
+++ b/source/blender/blenkernel/intern/bvhutils.c
@@ -81,7 +81,7 @@ static float sphereray_tri_intersection(const BVHTreeRay *ray, float radius, con
  * Function adapted from David Eberly's distance tools (LGPL)
  * http://www.geometrictools.com/LibFoundation/Distance/Distance.html
  */
-static float nearest_point_in_tri_surface(const float *v0,const float *v1,const float *v2,const float *p, int *v, int *e, float *d, float *nearest )
+static float nearest_point_in_tri_surface(const float *v0,const float *v1,const float *v2,const float *p, int *v, int *e, float *nearest )
 {
 	float diff[3];
 	float e0[3];
@@ -386,7 +386,7 @@ static float nearest_point_in_tri_surface(const float *v0,const float *v1,const
 		VecMulf(y, T);
 		VECADD(z, w, x);
 		VECADD(z, z, y);
-		VECSUB(d, p, z);
+		//VECSUB(d, p, z);
 		VECCOPY(nearest, z);
 		// d = p - ( v0 + S * e0 + T * e1 );
 	}
@@ -418,16 +418,16 @@ static void mesh_faces_nearest_point(void *userdata, int index, const float *co,
 	
 	do
 	{	
-		float nearest_tmp[3], col_normal[3], dist;
+		float nearest_tmp[3], dist;
 		int vertex, edge;
 		
-		dist = nearest_point_in_tri_surface(t0, t1, t2, co, &vertex, &edge, col_normal, nearest_tmp);
+		dist = nearest_point_in_tri_surface(t0, t1, t2, co, &vertex, &edge, nearest_tmp);
 		if(dist < nearest->dist)
 		{
 			nearest->index = index;
 			nearest->dist = dist;
 			VECCOPY(nearest->co, nearest_tmp);
-			VECCOPY(nearest->no, col_normal);
+			CalcNormFloat(t0, t1, t2, nearest->no);
 		}
 
 		t1 = t2;
diff --git a/source/blender/blenkernel/intern/deform.c b/source/blender/blenkernel/intern/deform.c
index ab53571b62d..3143c5e4df2 100644
--- a/source/blender/blenkernel/intern/deform.c
+++ b/source/blender/blenkernel/intern/deform.c
@@ -220,3 +220,31 @@ void unique_vertexgroup_name (bDeformGroup *dg, Object *ob)
 		}
 	}	
 }
+
+float deformvert_get_weight(const struct MDeformVert *dvert, int group_num)
+{
+	if(dvert)
+	{
+		const MDeformWeight *dw = dvert->dw;
+		int i;
+
+		for(i=dvert->totweight; i>0; i--, dw++)
+			if(dw->def_nr == group_num)
+				return dw->weight;
+	}
+
+	/* Not found */
+	return 0.0;
+}
+
+float vertexgroup_get_vertex_weight(const struct MDeformVert *dvert, int index, int group_num)
+{
+	if(group_num == -1)
+		return 1.0;
+
+	if(dvert == 0)
+		return 0.0;
+
+	return deformvert_get_weight(dvert+index, group_num);
+}
+
diff --git a/source/blender/blenkernel/intern/modifier.c b/source/blender/blenkernel/intern/modifier.c
index 777be445e63..4c74fe1cca2 100644
--- a/source/blender/blenkernel/intern/modifier.c
+++ b/source/blender/blenkernel/intern/modifier.c
@@ -101,6 +101,8 @@
 #include "BKE_utildefines.h"
 #include "depsgraph_private.h"
 #include "BKE_bmesh.h"
+#include "BKE_deform.h"
+#include "BKE_shrinkwrap.h"
 
 #include "LOD_DependKludge.h"
 #include "LOD_decimation.h"
@@ -6120,22 +6122,6 @@ CustomDataMask explodeModifier_requiredDataMask(ModifierData *md)
 	return dataMask;
 }
 
-/* this should really be put somewhere permanently */
-static float vert_weight(MDeformVert *dvert, int group)
-{
-	MDeformWeight *dw;
-	int i;
-	
-	if(dvert) {
-		dw= dvert->dw;
-		for(i= dvert->totweight; i>0; i--, dw++) {
-			if(dw->def_nr == group) return dw->weight;
-			if(i==1) break; /*otherwise dw will point to somewhere it shouldn't*/
-		}
-	}
-	return 0.0;
-}
-
 static void explodeModifier_createFacepa(ExplodeModifierData *emd,
 					 ParticleSystemModifierData *psmd,
       Object *ob, DerivedMesh *dm)
@@ -6179,7 +6165,7 @@ static void explodeModifier_createFacepa(ExplodeModifierData *emd,
 			for(i=0; i<totvert; i++){
 				val = BLI_frand();
 				val = (1.0f-emd->protect)*val + emd->protect*0.5f;
-				if(val < vert_weight(dvert+i,emd->vgroup-1))
+				if(val < deformvert_get_weight(dvert+i,emd->vgroup-1))
 					vertpa[i] = -1;
 			}
 		}
@@ -7236,6 +7222,126 @@ 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_PROJECT_ALLOW_POS_DIR;
+	smd->keepDist	= 0.0f;
+
+	smd->target		= NULL;
+	smd->auxTarget	= NULL;
+}
+
+static void shrinkwrapModifier_copyData(ModifierData *md, ModifierData *target)
+{
+	ShrinkwrapModifierData *smd  = (ShrinkwrapModifierData*)md;
+	ShrinkwrapModifierData *tsmd = (ShrinkwrapModifierData*)target;
+
+	tsmd->target	= smd->target;
+	tsmd->auxTarget = smd->auxTarget;
+
+	strcpy(tsmd->vgroup_name, smd->vgroup_name);
+
+	tsmd->keepDist	= smd->keepDist;
+	tsmd->shrinkType= smd->shrinkType;
+	tsmd->shrinkOpts= smd->shrinkOpts;
+}
+
+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);
+
+	if(smd->shrinkType == MOD_SHRINKWRAP_PROJECT
+	&& smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL)
+		dataMask |= (1 << CD_MVERT);
+		
+	return dataMask;
+}
+
+static int shrinkwrapModifier_isDisabled(ModifierData *md)
+{
+	ShrinkwrapModifierData *smd = (ShrinkwrapModifierData*) md;
+	return !smd->target;
+}
+
+
+static void shrinkwrapModifier_foreachObjectLink(ModifierData *md, Object *ob, ObjectWalkFunc walk, void *userData)
+{
+	ShrinkwrapModifierData *smd = (ShrinkwrapModifierData*) md;
+
+	walk(userData, ob, &smd->target);
+	walk(userData, ob, &smd->auxTarget);
+}
+
+static void shrinkwrapModifier_deformVerts(ModifierData *md, Object *ob, DerivedMesh *derivedData, float (*vertexCos)[3], int numVerts)
+{
+	DerivedMesh *dm = NULL;
+	CustomDataMask dataMask = shrinkwrapModifier_requiredDataMask(md);
+
+	/* We implement requiredDataMask but thats not really usefull since mesh_calc_modifiers pass a NULL derivedData or without the modified vertexs applied */
+	if(shrinkwrapModifier_requiredDataMask(md))
+	{
+		if(derivedData) dm = CDDM_copy(derivedData);
+		else if(ob->type==OB_MESH) dm = CDDM_from_mesh(ob->data, ob);
+		else return;
+
+		if(dataMask & CD_MVERT)
+		{
+			CDDM_apply_vert_coords(dm, vertexCos);
+			CDDM_calc_normals(dm);
+		}
+	}
+
+	shrinkwrapModifier_deform((ShrinkwrapModifierData*)md, ob, dm, vertexCos, numVerts);
+
+	if(dm)
+		dm->release(dm);
+}
+
+static void shrinkwrapModifier_deformVertsEM(ModifierData *md, Object *ob, EditMesh *editData, DerivedMesh *derivedData, float (*vertexCos)[3], int numVerts)
+{
+	DerivedMesh *dm = NULL;
+	CustomDataMask dataMask = shrinkwrapModifier_requiredDataMask(md);
+
+	if(dataMask)
+	{
+		if(derivedData) dm = CDDM_copy(derivedData);
+		else if(ob->type==OB_MESH) dm = CDDM_from_editmesh(editData, ob->data);
+		else return;
+
+		if(dataMask & CD_MVERT)
+		{
+			CDDM_apply_vert_coords(dm, vertexCos);
+			CDDM_calc_normals(dm);
+		}
+	}
+
+	shrinkwrapModifier_deform((ShrinkwrapModifierData*)md, ob, dm, vertexCos, numVerts);
+
+	if(dm)
+		dm->release(dm);
+}
+
+static void shrinkwrapModifier_updateDepgraph(ModifierData *md, DagForest *forest, Object *ob, DagNode *obNode)
+{
+	ShrinkwrapModifierData *smd = (ShrinkwrapModifierData*) md;
+
+	if (smd->target)
+		dag_add_relation(forest, dag_get_node(forest, smd->target),   obNode, DAG_RL_OB_DATA | DAG_RL_DATA_DATA, "Shrinkwrap Modifier");
+
+	if (smd->auxTarget)
+		dag_add_relation(forest, dag_get_node(forest, smd->auxTarget), obNode, DAG_RL_OB_DATA | DAG_RL_DATA_DATA, "Shrinkwrap Modifier");
+}
+
 /***/
 
 static ModifierTypeInfo typeArr[NUM_MODIFIER_TYPES];
@@ -7557,6 +7663,21 @@ ModifierTypeInfo *modifierType_getInfo(ModifierType type)
 		mti->requiredDataMask = explodeModifier_requiredDataMask;
 		mti->applyModifier = explodeModifier_applyModifier;
 
+		mti = INIT_TYPE(Shrinkwrap);
+		mti->type = eModifierTypeType_OnlyDeform;
+		mti->flags = eModifierTypeFlag_AcceptsMesh
+				| eModifierTypeFlag_AcceptsCVs
+				| eModifierTypeFlag_SupportsEditmode
+				| eModifierTypeFlag_EnableInEditmode;
+		mti->initData = shrinkwrapModifier_initData;
+		mti->copyData = shrinkwrapModifier_copyData;
+		mti->requiredDataMask = shrinkwrapModifier_requiredDataMask;
+		mti->isDisabled = shrinkwrapModifier_isDisabled;
+		mti->foreachObjectLink = shrinkwrapModifier_foreachObjectLink;
+		mti->deformVerts = shrinkwrapModifier_deformVerts;
+		mti->deformVertsEM = shrinkwrapModifier_deformVertsEM;
+		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
index 00000000000..292a800e9cd
--- /dev/null
+++ b/source/blender/blenkernel/intern/shrinkwrap.c
@@ -0,0 +1,588 @@
+/**
+ * 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 <memory.h>
+#include <stdio.h>
+#include <time.h>
+#include <assert.h>
+
+#include "DNA_object_types.h"
+#include "DNA_modifier_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_mesh_types.h"
+
+#include "BKE_shrinkwrap.h"
+#include "BKE_DerivedMesh.h"
+#include "BKE_utildefines.h"
+#include "BKE_deform.h"
+#include "BKE_cdderivedmesh.h"
+#include "BKE_displist.h"
+#include "BKE_global.h"
+#include "BKE_subsurf.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 !defined(_WIN32) && 0
+
+#include <sys/time.h>
+
+#define BENCH(a)	\
+	do {			\
+		double _t1, _t2;				\
+		struct timeval _tstart, _tend;	\
+		clock_t _clock_init = clock();	\
+		gettimeofday ( &_tstart, NULL);	\
+		(a);							\
+		gettimeofday ( &_tend, NULL);	\
+		_t1 = ( double ) _tstart.tv_sec + ( double ) _tstart.tv_usec/ ( 1000*1000 );	\
+		_t2 = ( double )   _tend.tv_sec + ( double )   _tend.tv_usec/ ( 1000*1000 );	\
+		printf("%s: %fs (real) %fs (cpu)\n", #a, _t2-_t1, (float)(clock()-_clock_init)/CLOCKS_PER_SEC);\
+	} while(0)
+
+#else
+
+#define BENCH(a)	(a)
+
+#endif
+
+typedef void ( *Shrinkwrap_ForeachVertexCallback) (DerivedMesh *target, float *co, float *normal);
+
+/* get derived mesh */
+//TODO is anyfunction that does this? returning the derivedFinal witouth we caring if its in edit mode or not?
+DerivedMesh *object_get_derived_final(Object *ob, CustomDataMask dataMask)
+{
+	if (ob==G.obedit)
+	{
+		DerivedMesh *final = NULL;
+		editmesh_get_derived_cage_and_final(&final, dataMask);
+		return final;
+	}
+	else
+		return mesh_get_derived_final(ob, dataMask);
+}
+
+/* Space transform */
+void space_transform_from_matrixs(SpaceTransform *data, float local[4][4], float target[4][4])
+{
+	float itarget[4][4];
+	Mat4Invert(itarget, target);
+	Mat4MulSerie(data->local2target, itarget, local, 0, 0, 0, 0, 0, 0);
+	Mat4Invert(data->target2local, data->local2target);
+}
+
+void space_transform_apply(const SpaceTransform *data, float *co)
+{
+	VecMat4MulVecfl(co, data->local2target, co);
+}
+
+void space_transform_invert(const SpaceTransform *data, float *co)
+{
+	VecMat4MulVecfl(co, data->target2local, co);
+}
+
+void space_transform_apply_normal(const SpaceTransform *data, float *no)
+{
+	Mat4Mul3Vecfl(data->local2target, no);
+	Normalize(no); // TODO: could we just determine de scale value from the matrix?
+}
+
+void space_transform_invert_normal(const SpaceTransform *data, float *no)
+{
+	Mat4Mul3Vecfl(data->target2local, no);
+	Normalize(no); // TODO: could we just determine de scale value from the matrix?
+}
+
+/*
+ * 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);
+}
+
+/* Main shrinkwrap function */
+void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, Object *ob, DerivedMesh *dm, float (*vertexCos)[3], int numVerts)
+{
+
+	ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;
+
+	//remove loop dependencies on derived meshs (TODO should this be done elsewhere?)
+	if(smd->target == ob) smd->target = NULL;
+	if(smd->auxTarget == ob) smd->auxTarget = NULL;
+
+
+	//Configure Shrinkwrap calc data
+	calc.smd = smd;
+	calc.ob = ob;
+	calc.original = dm;
+	calc.numVerts = numVerts;
+	calc.vertexCos = vertexCos;
+
+	if(smd->target)
+	{
+		//TODO currently we need a copy in case object_get_derived_final returns an emDM that does not defines getVertArray or getFace array
+		calc.target = CDDM_copy( object_get_derived_final(smd->target, CD_MASK_BAREMESH) );
+
+		//TODO there might be several "bugs" on non-uniform scales matrixs.. because it will no longer be nearest surface, not sphere projection
+		//because space has been deformed
+		space_transform_setup(&calc.local2target, ob, smd->target);
+
+		calc.keepDist = smd->keepDist;	//TODO: smd->keepDist is in global units.. must change to local
+	}
+
+
+	//Projecting target defined - lets work!
+	if(calc.target)
+	{
+		switch(smd->shrinkType)
+		{
+			case MOD_SHRINKWRAP_NEAREST_SURFACE:
+				BENCH(shrinkwrap_calc_nearest_surface_point(&calc));
+			break;
+
+			case MOD_SHRINKWRAP_PROJECT:
+				BENCH(shrinkwrap_calc_normal_projection(&calc));
+			break;
+
+			case MOD_SHRINKWRAP_NEAREST_VERTEX:
+				BENCH(shrinkwrap_calc_nearest_vertex(&calc));
+			break;
+		}
+	}
+
+	//free memory
+	if(calc.target)
+		calc.target->release( calc.target );
+}
+
+/*
+ * Shrinkwrap to the nearest vertex
+ *
+ * it builds a kdtree of vertexs we can attach to and then
+ * for each vertex performs a nearest vertex search on the tree
+ */
+void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
+{
+	int i;
+	const int vgroup		 = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
+	MDeformVert *const dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
+
+	BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
+	BVHTreeNearest  nearest  = NULL_BVHTreeNearest;
+
+
+	BENCH(bvhtree_from_mesh_verts(&treeData, calc->target, 0.0, 2, 6));
+	if(treeData.tree == NULL) return OUT_OF_MEMORY();
+
+	//Setup nearest
+	nearest.index = -1;
+	nearest.dist = FLT_MAX;
+
+#pragma omp parallel for default(none) private(i) firstprivate(nearest) shared(treeData,calc) schedule(static)
+	for(i = 0; i<calc->numVerts; ++i)
+	{
+		float *co = calc->vertexCos[i];
+		float tmp_co[3];
+		float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
+		if(weight == 0.0f) continue;
+
+		VECCOPY(tmp_co, co);
+		space_transform_apply(&calc->local2target, tmp_co); //Convert the coordinates to the tree coordinates
+
+		//Use local proximity heuristics (to reduce the nearest search)
+		//
+		//If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex
+		//so we can initiate the "nearest.dist" with the expected value to that last hit.
+		//This will lead in prunning of the search tree.
+		if(nearest.index != -1)
+			nearest.dist = squared_dist(tmp_co, nearest.co);
+		else
+			nearest.dist = FLT_MAX;
+
+		BLI_bvhtree_find_nearest(treeData.tree, tmp_co, &nearest, treeData.nearest_callback, &treeData);
+
+
+		//Found the nearest vertex
+		if(nearest.index != -1)
+		{
+			//Adjusting the vertex weight, so that after interpolating it keeps a certain distance from the nearest position
+			float dist = sasqrt(nearest.dist);
+			if(dist > FLT_EPSILON) weight *= (dist - calc->keepDist)/dist;
+
+			//Convert the coordinates back to mesh coordinates
+			VECCOPY(tmp_co, nearest.co);
+			space_transform_invert(&calc->local2target, tmp_co);
+
+			VecLerpf(co, co, tmp_co, weight);	//linear interpolation
+		}
+	}
+
+	free_bvhtree_from_mesh(&treeData);
+}
+
+/*
+ * This function raycast a single vertex and updates the hit if the "hit" is considered valid.
+ * Returns TRUE if "hit" was updated.
+ * Opts control whether an hit is valid or not
+ * Supported options are:
+ *	MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE (front faces hits are ignored)
+ *	MOD_SHRINKWRAP_CULL_TARGET_BACKFACE (back faces hits are ignored)
+ */
+int normal_projection_project_vertex(char options, const float *vert, const float *dir, const SpaceTransform *transf, BVHTree *tree, BVHTreeRayHit *hit, BVHTree_RayCastCallback callback, void *userdata)
+{
+	float tmp_co[3], tmp_no[3];
+	const float *co, *no;
+	BVHTreeRayHit hit_tmp;
+
+	//Copy from hit (we need to convert hit rays from one space coordinates to the other
+	memcpy( &hit_tmp, hit, sizeof(hit_tmp) );
+
+	//Apply space transform (TODO readjust dist)
+	if(transf)
+	{
+		VECCOPY( tmp_co, vert );
+		space_transform_apply( transf, tmp_co );
+		co = tmp_co;
+
+		VECCOPY( tmp_no, dir );
+		space_transform_apply_normal( transf, tmp_no );
+		no = tmp_no;
+
+		hit_tmp.dist *= Mat4ToScalef( transf->local2target );
+	}
+	else
+	{
+		co = vert;
+		no = dir;
+	}
+
+	hit_tmp.index = -1;
+
+	BLI_bvhtree_ray_cast(tree, co, no, &hit_tmp, callback, userdata);
+
+	if(hit_tmp.index != -1)
+	{
+		float dot = INPR( dir, hit_tmp.no);
+
+		if(((options & MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE) && dot <= 0.0f)
+		|| ((options & MOD_SHRINKWRAP_CULL_TARGET_BACKFACE) && dot >= 0.0f))
+			return FALSE; //Ignore hit
+
+
+		//Inverting space transform (TODO make coeherent with the initial dist readjust)
+		if(transf)
+		{
+			space_transform_invert( transf, hit_tmp.co );
+			space_transform_invert_normal( transf, hit_tmp.no );
+
+			hit_tmp.dist = VecLenf( vert, hit_tmp.co );
+		}
+
+		memcpy(hit, &hit_tmp, sizeof(hit_tmp) );
+		return TRUE;
+	}
+	return FALSE;
+}
+
+
+void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
+{
+	int i;
+
+	//Options about projection direction
+	const char use_normal    = calc->smd->shrinkOpts;
+	float proj_axis[3] = {0.0f, 0.0f, 0.0f};
+	MVert *vert  = NULL; //Needed in case of vertex normal
+	DerivedMesh* ss_mesh = NULL;
+
+	//Vertex group data
+	const int vgroup		   = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
+	const MDeformVert *dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
+
+
+	//Raycast and tree stuff
+	BVHTreeRayHit hit;
+	BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh; 	//target
+
+	//auxiliar target
+	DerivedMesh * aux_mesh = NULL;
+	BVHTreeFromMesh auxData= NULL_BVHTreeFromMesh;
+	SpaceTransform local2aux;
+
+do
+{
+
+	//Prepare data to retrieve the direction in which we should project each vertex
+	if(calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL)
+	{
+		//No Mvert information: jump to "free memory and return" part
+		if(calc->original == NULL) break;
+
+		if(calc->smd->subsurfLevels)
+		{
+			SubsurfModifierData smd;
+			memset(&smd, 0, sizeof(smd));
+			smd.subdivType = ME_CC_SUBSURF;			//catmull clark
+			smd.levels = calc->smd->subsurfLevels;	//levels
+
+			ss_mesh = subsurf_make_derived_from_derived(calc->original, &smd, FALSE, NULL, 0, 0);
+
+			if(ss_mesh)
+			{
+				vert = ss_mesh->getVertDataArray(ss_mesh, CD_MVERT);
+				if(vert)
+				{
+					//TRICKY: this code assumes subsurface will have the transformed original vertices
+					//in their original order at the end of the vert array.
+					vert = vert
+						 + ss_mesh->getNumVerts(ss_mesh)
+						 - calc->original->getNumVerts(calc->original);
+				}
+			}
+
+			//To make sure we are not letting any memory behind
+			assert(smd.emCache == NULL);
+			assert(smd.mCache == NULL);
+		}
+		else
+			vert = calc->original->getVertDataArray(calc->original, CD_MVERT);
+
+		//Not able to get vert information: jump to "free memory and return" part
+		if(vert == NULL) break;
+	}
+	else
+	{
+		//The code supports any axis that is a combination of X,Y,Z.. altought currently UI only allows to set the 3 diferent axis
+		if(calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS) proj_axis[0] = 1.0f;
+		if(calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS) proj_axis[1] = 1.0f;
+		if(calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS) proj_axis[2] = 1.0f;
+
+		Normalize(proj_axis);
+
+		//Invalid projection direction: jump to "free memory and return" part
+		if(INPR(proj_axis, proj_axis) < FLT_EPSILON) break; 
+	}
+
+	//If the user doesn't allows to project in any direction of projection axis... then theres nothing todo.
+	if((use_normal & (MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)) == 0)
+		break; //jump to "free memory and return" part
+
+
+	//Build target tree
+	BENCH(bvhtree_from_mesh_faces(&treeData, calc->target, calc->keepDist, 4, 6));
+	if(treeData.tree == NULL)
+		break; //jump to "free memory and return" part
+
+
+	//Build auxiliar target
+	if(calc->smd->auxTarget)
+	{
+		space_transform_setup( &local2aux, calc->ob, calc->smd->auxTarget);
+
+		aux_mesh = CDDM_copy( object_get_derived_final(calc->smd->auxTarget, CD_MASK_BAREMESH) ); 		//TODO currently we need a copy in case object_get_derived_final returns an emDM that does not defines getVertArray or getFace array
+		if(aux_mesh)
+			BENCH(bvhtree_from_mesh_faces(&auxData, aux_mesh, 0.0, 4, 6));
+		else
+			printf("Auxiliar target finalDerived mesh is null\n");
+	}
+
+
+	//Now, everything is ready to project the vertexs!
+#pragma omp parallel for private(i,hit) schedule(static)
+	for(i = 0; i<calc->numVerts; ++i)
+	{
+		float *co = calc->vertexCos[i];
+		float tmp_co[3], tmp_no[3];
+		float lim = 10000.0f; //TODO: we should use FLT_MAX here, but sweepsphere code isnt prepared for that
+		float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
+
+		if(weight == 0.0f) continue;
+
+		if(ss_mesh)
+		{
+			VECCOPY(tmp_co, vert[i].co);
+		}
+		else
+		{
+			VECCOPY(tmp_co, co);
+		}
+
+
+		if(vert)
+			NormalShortToFloat(tmp_no, vert[i].no);
+		else
+			VECCOPY( tmp_no, proj_axis );
+
+
+		hit.index = -1;
+		hit.dist = lim;
+
+
+		//Project over positive direction of axis
+		if(use_normal & MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR)
+		{
+
+			if(auxData.tree)
+				normal_projection_project_vertex(0, tmp_co, tmp_no, &local2aux, auxData.tree, &hit, auxData.raycast_callback, &auxData);
+
+			normal_projection_project_vertex(calc->smd->shrinkOpts, tmp_co, tmp_no, &calc->local2target, treeData.tree, &hit, treeData.raycast_callback, &treeData);
+		}
+
+		//Project over negative direction of axis
+		if(use_normal & MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)
+		{
+			float inv_no[3] = { -tmp_no[0], -tmp_no[1], -tmp_no[2] };
+
+
+			if(auxData.tree)
+				normal_projection_project_vertex(0, tmp_co, inv_no, &local2aux, auxData.tree, &hit, auxData.raycast_callback, &auxData);
+
+			normal_projection_project_vertex(calc->smd->shrinkOpts, tmp_co, inv_no, &calc->local2target, treeData.tree, &hit, treeData.raycast_callback, &treeData);
+		}
+
+
+		if(hit.index != -1)
+		{
+			VecLerpf(co, co, hit.co, weight);
+		}
+	}
+
+
+//Simple do{} while(0) structure to allow to easily jump to the "free memory and return" part
+} while(0);
+
+	//free data structures
+
+	free_bvhtree_from_mesh(&treeData);
+	free_bvhtree_from_mesh(&auxData);
+
+	if(aux_mesh)
+		aux_mesh->release(aux_mesh);
+
+	if(ss_mesh)
+		ss_mesh->release(ss_mesh);
+}
+
+/*
+ * 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;
+
+	const int vgroup = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
+	const MDeformVert *const dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
+
+	BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
+	BVHTreeNearest  nearest  = NULL_BVHTreeNearest;
+
+
+
+	//Create a bvh-tree of the given target
+	BENCH(bvhtree_from_mesh_faces( &treeData, calc->target, 0.0, 2, 6));
+	if(treeData.tree == NULL) return OUT_OF_MEMORY();
+
+	//Setup nearest
+	nearest.index = -1;
+	nearest.dist = FLT_MAX;
+
+
+	//Find the nearest vertex
+#pragma omp parallel for default(none) private(i) firstprivate(nearest) shared(calc,treeData) schedule(static)
+	for(i = 0; i<calc->numVerts; ++i)
+	{
+		float *co = calc->vertexCos[i];
+		float tmp_co[3];
+		float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
+		if(weight == 0.0f) continue;
+
+		//Convert the vertex to tree coordinates
+		VECCOPY(tmp_co, co);
+		space_transform_apply(&calc->local2target, tmp_co);
+
+		//Use local proximity heuristics (to reduce the nearest search)
+		//
+		//If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex
+		//so we can initiate the "nearest.dist" with the expected value to that last hit.
+		//This will lead in prunning of the search tree.
+		if(nearest.index != -1)
+			nearest.dist = squared_dist(tmp_co, nearest.co);
+		else
+			nearest.dist = FLT_MAX;
+
+		BLI_bvhtree_find_nearest(treeData.tree, tmp_co, &nearest, treeData.nearest_callback, &treeData);
+
+		//Found the nearest vertex
+		if(nearest.index != -1)
+		{
+			if(calc->smd->shrinkOpts & MOD_SHRINKWRAP_KEEP_ABOVE_SURFACE)
+			{
+				//Make the vertex stay on the front side of the face
+				VECADDFAC(tmp_co, nearest.co, nearest.no, calc->keepDist);
+			}
+			else
+			{
+				//Adjusting the vertex weight, so that after interpolating it keeps a certain distance from the nearest position
+				float dist = sasqrt( nearest.dist );
+				if(dist > FLT_EPSILON)
+					VecLerpf(tmp_co, tmp_co, nearest.co, (dist - calc->keepDist)/dist);	//linear interpolation
+				else
+					VECCOPY( tmp_co, nearest.co );
+			}
+
+			//Convert the coordinates back to mesh coordinates
+			space_transform_invert(&calc->local2target, tmp_co);
+			VecLerpf(co, co, tmp_co, weight);	//linear interpolation
+		}
+	}
+
+
+	free_bvhtree_from_mesh(&treeData);
+}
+
diff --git a/source/blender/blenlib/intern/BLI_kdopbvh.c b/source/blender/blenlib/intern/BLI_kdopbvh.c
index 49f3c3cc9e6..c41c3d8c3ab 100644
--- a/source/blender/blenlib/intern/BLI_kdopbvh.c
+++ b/source/blender/blenlib/intern/BLI_kdopbvh.c
@@ -46,6 +46,7 @@
 
 
 #define MAX_TREETYPE 32
+#define DEFAULT_FIND_NEAREST_HEAP_SIZE 1024
 
 typedef struct BVHNode
 {
@@ -119,6 +120,72 @@ static float KDOP_AXES[13][3] =
 {0, 1.0, -1.0}
 };
 
+/*
+ * Generic push and pop heap
+ */
+#define PUSH_HEAP_BODY(HEAP_TYPE,PRIORITY,heap,heap_size)	\
+{													\
+	HEAP_TYPE element = heap[heap_size-1];			\
+	int child = heap_size-1;						\
+	while(child != 0)								\
+	{												\
+		int parent = (child-1) / 2;					\
+		if(PRIORITY(element, heap[parent]))			\
+		{											\
+			heap[child] = heap[parent];				\
+			child = parent;							\
+		}											\
+		else break;									\
+	}												\
+	heap[child] = element;							\
+}
+
+#define POP_HEAP_BODY(HEAP_TYPE, PRIORITY,heap,heap_size)	\
+{													\
+	HEAP_TYPE element = heap[heap_size-1];			\
+	int parent = 0;									\
+	while(parent < (heap_size-1)/2 )				\
+	{												\
+		int child2 = (parent+1)*2;					\
+		if(PRIORITY(heap[child2-1], heap[child2]))	\
+			--child2;								\
+													\
+		if(PRIORITY(element, heap[child2]))			\
+			break;									\
+													\
+		heap[parent] = heap[child2];				\
+		parent = child2;							\
+	}												\
+	heap[parent] = element;							\
+}
+
+int ADJUST_MEMORY(void *local_memblock, void **memblock, int new_size, int *max_size, int size_per_item)
+{
+	int   new_max_size = *max_size * 2;
+	void *new_memblock = NULL;
+
+	if(new_size <= *max_size)
+		return TRUE;
+
+	if(*memblock == local_memblock)
+	{
+		new_memblock = malloc( size_per_item * new_max_size );
+		memcpy( new_memblock, *memblock, size_per_item * *max_size );
+	}
+	else
+		new_memblock = realloc(*memblock, size_per_item * new_max_size );
+
+	if(new_memblock)
+	{
+		*memblock = new_memblock;
+		*max_size = new_max_size;
+		return TRUE;
+	}
+	else
+		return FALSE;
+}
+
+
 //////////////////////////////////////////////////////////////////////////////////////////////////////
 // Introsort 
 // with permission deriven from the following Java code:
@@ -734,6 +801,11 @@ BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis)
 		return NULL;
 
 	tree = (BVHTree *)MEM_callocN(sizeof(BVHTree), "BVHTree");
+
+	//tree epsilon must be >= FLT_EPSILON
+	//so that tangent rays can still hit a bounding volume..
+	//this bug would show up when casting a ray aligned with a kdop-axis and with an edge of 2 faces
+	epsilon = MAX2(FLT_EPSILON, epsilon);
 	
 	if(tree)
 	{
@@ -1126,15 +1198,18 @@ static float calc_nearest_point(BVHNearestData *data, BVHNode *node, float *near
 }
 
 
-// TODO: use a priority queue to reduce the number of nodes looked on
-static void dfs_find_nearest(BVHNearestData *data, BVHNode *node)
+typedef struct NodeDistance
 {
-	int i;
-	float nearest[3], sdist;
+	BVHNode *node;
+	float dist;
 
-	sdist = calc_nearest_point(data, node, nearest);
-	if(sdist >= data->nearest.dist) return;
+} NodeDistance;
+
+#define NodeDistance_priority(a,b)	( (a).dist < (b).dist )
 
+// TODO: use a priority queue to reduce the number of nodes looked on
+static void dfs_find_nearest_dfs(BVHNearestData *data, BVHNode *node)
+{
 	if(node->totnode == 0)
 	{
 		if(data->callback)
@@ -1142,17 +1217,130 @@ static void dfs_find_nearest(BVHNearestData *data, BVHNode *node)
 		else
 		{
 			data->nearest.index	= node->index;
-			VECCOPY(data->nearest.co, nearest);
-			data->nearest.dist	= sdist;
+			data->nearest.dist	= calc_nearest_point(data, node, data->nearest.co);
 		}
 	}
 	else
 	{
-		for(i=0; i != node->totnode; i++)
-			dfs_find_nearest(data, node->children[i]);
+		//Better heuristic to pick the closest node to dive on
+		int i;
+		float nearest[3];
+
+		if(data->proj[ node->main_axis ] <= node->children[0]->bv[node->main_axis*2+1])
+		{
+
+			for(i=0; i != node->totnode; i++)
+			{
+				if( calc_nearest_point(data, node->children[i], nearest) >= data->nearest.dist) continue;
+				dfs_find_nearest_dfs(data, node->children[i]);
+			}
+		}
+		else
+		{
+			for(i=node->totnode-1; i >= 0 ; i--)
+			{
+				if( calc_nearest_point(data, node->children[i], nearest) >= data->nearest.dist) continue;
+				dfs_find_nearest_dfs(data, node->children[i]);
+			}
+		}
+	}
+}
+
+static void dfs_find_nearest_begin(BVHNearestData *data, BVHNode *node)
+{
+	int i;
+	float nearest[3], sdist;
+	sdist = calc_nearest_point(data, node, nearest);
+	if(sdist >= data->nearest.dist) return;
+	dfs_find_nearest_dfs(data, node);
+}
+
+
+static void NodeDistance_push_heap(NodeDistance *heap, int heap_size)
+PUSH_HEAP_BODY(NodeDistance, NodeDistance_priority, heap, heap_size)
+
+static void NodeDistance_pop_heap(NodeDistance *heap, int heap_size)
+POP_HEAP_BODY(NodeDistance, NodeDistance_priority, heap, heap_size)
+
+//NN function that uses an heap.. this functions leads to an optimal number of min-distance
+//but for normal tri-faces and BV 6-dop.. a simple dfs with local heuristics (as implemented
+//in source/blender/blenkernel/intern/shrinkwrap.c) works faster.
+//
+//It may make sense to use this function if the callback queries are very slow.. or if its impossible
+//to get a nice heuristic
+//
+//this function uses "malloc/free" instead of the MEM_* because it intends to be openmp safe
+static void bfs_find_nearest(BVHNearestData *data, BVHNode *node)
+{
+	int i;
+	NodeDistance default_heap[DEFAULT_FIND_NEAREST_HEAP_SIZE];
+	NodeDistance *heap=default_heap, current;
+	int heap_size = 0, max_heap_size = sizeof(default_heap)/sizeof(default_heap[0]);
+	float nearest[3];
+
+	int callbacks = 0, push_heaps = 0;
+
+	if(node->totnode == 0)
+	{
+		dfs_find_nearest_dfs(data, node);
+		return;
+	}
+
+	current.node = node;
+	current.dist = calc_nearest_point(data, node, nearest);
+
+	while(current.dist < data->nearest.dist)
+	{
+//		printf("%f : %f\n", current.dist, data->nearest.dist);
+		for(i=0; i< current.node->totnode; i++)
+		{
+			BVHNode *child = current.node->children[i];
+			if(child->totnode == 0)
+			{
+				callbacks++;
+				dfs_find_nearest_dfs(data, child);
+			}
+			else
+			{
+				//adjust heap size
+				if(heap_size >= max_heap_size
+				&& ADJUST_MEMORY(default_heap, (void**)&heap, heap_size+1, &max_heap_size, sizeof(heap[0])) == FALSE)
+				{
+					printf("WARNING: bvh_find_nearest got out of memory\n");
+
+					if(heap != default_heap)
+						free(heap);
+
+					return;
+				}
+
+				heap[heap_size].node = current.node->children[i];
+				heap[heap_size].dist = calc_nearest_point(data, current.node->children[i], nearest);
+
+				if(heap[heap_size].dist >= data->nearest.dist) continue;
+				heap_size++;
+
+				NodeDistance_push_heap(heap, heap_size);
+	//			PUSH_HEAP_BODY(NodeDistance, NodeDistance_priority, heap, heap_size);
+				push_heaps++;
+			}
+		}
+		
+		if(heap_size == 0) break;
+
+		current = heap[0];
+		NodeDistance_pop_heap(heap, heap_size);
+//		POP_HEAP_BODY(NodeDistance, NodeDistance_priority, heap, heap_size);
+		heap_size--;
 	}
+
+//	printf("hsize=%d, callbacks=%d, pushs=%d\n", heap_size, callbacks, push_heaps);
+
+	if(heap != default_heap)
+		free(heap);
 }
 
+
 int BLI_bvhtree_find_nearest(BVHTree *tree, const float *co, BVHTreeNearest *nearest, BVHTree_NearestPointCallback callback, void *userdata)
 {
 	int i;
@@ -1184,7 +1372,7 @@ int BLI_bvhtree_find_nearest(BVHTree *tree, const float *co, BVHTreeNearest *nea
 
 	//dfs search
 	if(root)
-		dfs_find_nearest(&data, root);
+		dfs_find_nearest_begin(&data, root);
 
 	//copy back results
 	if(nearest)
@@ -1224,7 +1412,7 @@ static float ray_nearest_hit(BVHRayCastData *data, BVHNode *node)
 			float ll = (bv[0] - data->ray.origin[i]) / data->ray_dot_axis[i];
 			float lu = (bv[1] - data->ray.origin[i]) / data->ray_dot_axis[i];
 
-			if(data->ray_dot_axis[i] > 0)
+			if(data->ray_dot_axis[i] > 0.0f)
 			{
 				if(ll > low)   low = ll;
 				if(lu < upper) upper = lu;
@@ -1264,7 +1452,7 @@ static void dfs_raycast(BVHRayCastData *data, BVHNode *node)
 	else
 	{
 		//pick loop direction to dive into the tree (based on ray direction and split axis)
-		if(data->ray_dot_axis[ node->main_axis ] > 0)
+		if(data->ray_dot_axis[ node->main_axis ] > 0.0f)
 		{
 			for(i=0; i != node->totnode; i++)
 			{
@@ -1301,7 +1489,7 @@ int BLI_bvhtree_ray_cast(BVHTree *tree, const float *co, const float *dir, BVHTr
 	{
 		data.ray_dot_axis[i] = INPR( data.ray.direction, KDOP_AXES[i]);
 
-		if(fabs(data.ray_dot_axis[i]) < 1e-7)
+		if(fabs(data.ray_dot_axis[i]) < FLT_EPSILON)
 			data.ray_dot_axis[i] = 0.0;
 	}
 
diff --git a/source/blender/makesdna/DNA_modifier_types.h b/source/blender/makesdna/DNA_modifier_types.h
index d8cc7633fb3..9599cc1d247 100644
--- a/source/blender/makesdna/DNA_modifier_types.h
+++ b/source/blender/makesdna/DNA_modifier_types.h
@@ -35,6 +35,7 @@ typedef enum ModifierType {
 	eModifierType_Cloth,
 	eModifierType_Collision,
 	eModifierType_Bevel,
+	eModifierType_Shrinkwrap,
 	NUM_MODIFIER_TYPES
 } ModifierType;
 
@@ -491,4 +492,45 @@ typedef struct ExplodeModifierData {
 	float protect;
 } ExplodeModifierData;
 
+typedef struct ShrinkwrapModifierData {
+	ModifierData modifier;
+
+	struct Object *target;	/* shrink target */
+	struct Object *auxTarget; /* additional shrink target */
+	char vgroup_name[32];	/* optional vertexgroup name */
+	float keepDist;			/* distance offset to keep from mesh/projection point */
+	short shrinkType;		/* shrink type projection */
+	short shrinkOpts;		/* shrink options */
+	char projAxis;			/* axis to project over */
+
+	/*
+	 * if using projection over vertex normal this controls the
+	 * the level of subsurface that must be done before getting the
+	 * vertex coordinates and normal
+	 */
+	char subsurfLevels;
+
+	char pad[6];
+
+} ShrinkwrapModifierData;
+
+/* Shrinkwrap->shrinkType */
+#define MOD_SHRINKWRAP_NEAREST_SURFACE	0
+#define MOD_SHRINKWRAP_PROJECT			1
+#define MOD_SHRINKWRAP_NEAREST_VERTEX	2
+
+/* Shrinkwrap->shrinkOpts */
+#define MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR	(1<<0)	/* allow shrinkwrap to move the vertex in the positive direction of axis */
+#define MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR	(1<<1)	/* allow shrinkwrap to move the vertex in the negative direction of axis */
+
+#define MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE	(1<<3)	/* ignore vertex moves if a vertex ends projected on a front face of the target */
+#define MOD_SHRINKWRAP_CULL_TARGET_BACKFACE		(1<<4)	/* ignore vertex moves if a vertex ends projected on a back face of the target */
+
+#define MOD_SHRINKWRAP_KEEP_ABOVE_SURFACE		(1<<5)	/* distance is measure to the front face of the target */
+
+#define MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS		(1<<0)
+#define MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS		(1<<1)
+#define MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS		(1<<2)
+#define MOD_SHRINKWRAP_PROJECT_OVER_NORMAL			0	/* projection over normal is used if no axis is selected */
+
 #endif
diff --git a/source/blender/src/buttons_editing.c b/source/blender/src/buttons_editing.c
index 27ddca1aaa5..76f244db6fc 100644
--- a/source/blender/src/buttons_editing.c
+++ b/source/blender/src/buttons_editing.c
@@ -1826,6 +1826,16 @@ 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 + 3;
+			if (smd->shrinkType == MOD_SHRINKWRAP_PROJECT)
+			{
+				height += 19*5;
+				if(smd->projAxis == 0) height += 19;
+			}
+			else if (smd->shrinkType == MOD_SHRINKWRAP_NEAREST_SURFACE)
+				height += 19;
 		}
 							/* 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 +2456,51 @@ 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|projection %x1|nearest vertex %x2";
+
+			uiDefIDPoinBut(block, modifier_testMeshObj, ID_OB, B_CHANGEDEP, "Ob: ",	lx, (cy-=19), buttonWidth,19, &smd->target, "Target to shrink to");
+
+			but=uiDefBut(block, TEX, B_MODIFIER_RECALC, "VGroup: ",		lx, (cy-=19), buttonWidth,19, &smd->vgroup_name, 0, 31, 0, 0, "Vertex Group name");
+			uiButSetCompleteFunc(but, autocomplete_vgroup, (void *)ob);
+
+			uiDefButF(block, NUM, B_MODIFIER_RECALC, "Offset:",	lx,(cy-=19),buttonWidth,19, &smd->keepDist, 0.0f, 100.0f, 1.0f, 0, "Specify distance to keep from the target");
+
+			cy -= 3;
+			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_PROJECT){
+
+
+				/* UI for projection axis */
+				uiBlockBeginAlign(block);
+				uiDefButC(block, ROW, B_MODIFIER_RECALC, "Normal"    , lx,(cy-=19),buttonWidth,19, &smd->projAxis, 18.0, MOD_SHRINKWRAP_PROJECT_OVER_NORMAL, 0, 0, "Projection over X axis");
+				if(smd->projAxis == 0)
+				{
+					uiDefButC(block, NUM, B_MODIFIER_RECALC, "SS Levels:",		lx, (cy-=19), buttonWidth,19, &smd->subsurfLevels, 0, 6, 0, 0, "This indicates the number of CCSubdivisions that must be performed before extracting vertexs positions and normals");
+				}
+
+				uiDefButBitC(block, TOG, MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS, B_MODIFIER_RECALC, "X",	lx+buttonWidth/3*0,(cy-=19),buttonWidth/3,19, &smd->projAxis, 0, 0, 0, 0, "Projection over X axis");
+				uiDefButBitC(block, TOG, MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS, B_MODIFIER_RECALC, "Y",	lx+buttonWidth/3*1,cy,buttonWidth/3,19, &smd->projAxis, 0, 0, 0, 0, "Projection over Y axis");
+				uiDefButBitC(block, TOG, MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS, B_MODIFIER_RECALC, "Z",	lx+buttonWidth/3*2,cy,buttonWidth/3,19, &smd->projAxis, 0, 0, 0, 0, "Projection over Z axis");
+
+
+				/* allowed directions of projection axis */
+				uiDefButBitS(block, TOG, MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR, B_MODIFIER_RECALC, "Negative",	lx,(cy-=19),buttonWidth/2,19, &smd->shrinkOpts, 0, 0, 0, 0, "Allows to move the vertex in the negative direction of axis");
+				uiDefButBitS(block, TOG, MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR, B_MODIFIER_RECALC, "Positive",	lx + buttonWidth/2,cy,buttonWidth/2,19, &smd->shrinkOpts, 0, 0, 0, 0, "Allows to move the vertex in the positive direction of axis");
+
+				uiDefButBitS(block, TOG, MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE, B_MODIFIER_RECALC, "Cull frontfaces",lx,(cy-=19),buttonWidth/2,19, &smd->shrinkOpts, 0, 0, 0, 0, "Controls whether a vertex can be projected to a front face on target");
+				uiDefButBitS(block, TOG, MOD_SHRINKWRAP_CULL_TARGET_BACKFACE,  B_MODIFIER_RECALC, "Cull backfaces",	lx+buttonWidth/2,cy,buttonWidth/2,19, &smd->shrinkOpts, 0, 0, 0, 0, "Controls whether a vertex can be projected to a back face on target");
+				uiDefIDPoinBut(block, modifier_testMeshObj, ID_OB, B_CHANGEDEP, "Ob2: ",	lx, (cy-=19), buttonWidth,19, &smd->auxTarget, "Aditional mesh to project over");
+			}
+			else if (smd->shrinkType == MOD_SHRINKWRAP_NEAREST_SURFACE){
+				uiDefButBitS(block, TOG, MOD_SHRINKWRAP_KEEP_ABOVE_SURFACE, B_MODIFIER_RECALC, "Above surface",	lx,(cy-=19),buttonWidth,19, &smd->shrinkOpts, 0, 0, 0, 0, "Vertices are kept on the front side of faces");
+			}
+
+			uiBlockEndAlign(block);
+
 		}
 
 		uiBlockEndAlign(block);