Bone Heat Weighting

===================

This is a new automatic vertex weighting method, next to the existing
envelope based method. The details are here:

http://www.blender.org/development/current-projects/changes-since-244/skinning/

This is based on section 4 of the paper:

"Automatic Rigging and Animation of 3D Characters"
Ilya Baran and Jovan Popovic, SIGGRAPH 2007

Implementation Notes:

- Generic code for making mesh laplacian matrices has been added, which
  is only used by bone heat weighting at the moment.
- Bone to vertex visibility checking is done with the raytracing code.
- Fixed an issue in the subsurf limit calculation function, where the
  position of vertices on boundary edges was wrong. It is still not the
  correct position, but at least it's in the neighbourhood now.

10 files changed, 1247 insertions, 250 deletions

diff --git a/source/blender/blenkernel/intern/subsurf_ccg.c b/source/blender/blenkernel/intern/subsurf_ccg.c
index 7f3c5ddf8a7..8fc92829258 100644
--- a/source/blender/blenkernel/intern/subsurf_ccg.c
+++ b/source/blender/blenkernel/intern/subsurf_ccg.c
@@ -2441,11 +2441,11 @@ struct DerivedMesh *subsurf_make_derived_from_derived(
 
 void subsurf_calculate_limit_positions(Mesh *me, float (*positions_r)[3]) 
 {
-		/* Finds the subsurf limit positions for the verts in a mesh 
-		 * and puts them in an array of floats. Please note that the 
-		 * calculated vert positions is incorrect for the verts 
-		 * on the boundary of the mesh.
-		 */
+	/* Finds the subsurf limit positions for the verts in a mesh 
+	 * and puts them in an array of floats. Please note that the 
+	 * calculated vert positions is incorrect for the verts 
+	 * on the boundary of the mesh.
+	 */
 	CCGSubSurf *ss = _getSubSurf(NULL, 1, 0, 1, 0);
 	float edge_sum[3], face_sum[3];
 	CCGVertIterator *vi;
@@ -2474,6 +2474,11 @@ void subsurf_calculate_limit_positions(Mesh *me, float (*positions_r)[3])
 			VecAddf(face_sum, face_sum, ccgSubSurf_getFaceCenterData(ss, f));
 		}
 
+		/* ad-hoc correction for boundary vertices, to at least avoid them
+		   moving completely out of place (brecht) */
+		if(numFaces && numFaces != N)
+			VecMulf(face_sum, (float)N/(float)numFaces);
+
 		co = ccgSubSurf_getVertData(ss, v);
 		positions_r[idx][0] = (co[0]*N*N + edge_sum[0]*4 + face_sum[0])/(N*(N+5));
 		positions_r[idx][1] = (co[1]*N*N + edge_sum[1]*4 + face_sum[1])/(N*(N+5));
diff --git a/source/blender/blenlib/BLI_arithb.h b/source/blender/blenlib/BLI_arithb.h
index 4ee01f7fe33..a732f09ca58 100644
--- a/source/blender/blenlib/BLI_arithb.h
+++ b/source/blender/blenlib/BLI_arithb.h
@@ -258,6 +258,8 @@ void euler_rot(float *beul, float ang, char axis);
 
 float DistVL2Dfl(float *v1, float *v2, float *v3);
 float PdistVL2Dfl(float *v1, float *v2, float *v3);
+float PdistVL3Dfl(float *v1, float *v2, float *v3);
+void PclosestVL3Dfl(float *closest, float *v1, float *v2, float *v3);
 float AreaF2Dfl(float *v1, float *v2, float *v3);
 float AreaQ3Dfl(float *v1, float *v2, float *v3, float *v4);
 float AreaT3Dfl(float *v1, float *v2, float *v3);
diff --git a/source/blender/blenlib/intern/arithb.c b/source/blender/blenlib/intern/arithb.c
index 125c24a7ada..e98ce92787f 100644
--- a/source/blender/blenlib/intern/arithb.c
+++ b/source/blender/blenlib/intern/arithb.c
@@ -3284,6 +3284,31 @@ int point_in_tri_prism(float p[3], float v1[3], float v2[3], float v3[3])
 	return 1;
 }
 
+/* point closest to v1 on line v2-v3 in 3D */
+void PclosestVL3Dfl(float *closest, float *v1, float *v2, float *v3)
+{
+	float lambda, cp[3], len;
+
+	lambda= lambda_cp_line_ex(v1, v2, v3, cp);
+
+	if(lambda <= 0.0f)
+		VecCopyf(closest, v2);
+	else if(lambda >= 1.0f)
+		VecCopyf(closest, v3);
+	else
+		VecCopyf(closest, cp);
+}
+
+/* distance v1 to line-piece v2-v3 in 3D */
+float PdistVL3Dfl(float *v1, float *v2, float *v3) 
+{
+	float closest[3];
+
+	PclosestVL3Dfl(closest, v1, v2, v3);
+
+	return VecLenf(closest, v1);
+}
+
 /********************************************************/
 
 /* make a 4x4 matrix out of 3 transform components */
@@ -3331,3 +3356,4 @@ void LocQuatSizeToMat4(float mat[][4], float loc[3], float quat[4], float size[3
 	mat[3][1] = loc[1];
 	mat[3][2] = loc[2];
 }
+
diff --git a/source/blender/include/BIF_meshlaplacian.h b/source/blender/include/BIF_meshlaplacian.h
new file mode 100644
index 00000000000..1e10336723f
--- /dev/null
+++ b/source/blender/include/BIF_meshlaplacian.h
@@ -0,0 +1,81 @@
+/**
+ * $Id$
+ *
+ * ***** 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 use 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) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): none yet.
+ *
+ * ***** END GPL/BL DUAL LICENSE BLOCK *****
+ * BIF_meshlaplacian.h: Algorithms using the mesh laplacian.
+ */
+
+#ifndef BIF_MESHLAPLACIAN_H
+#define BIF_MESHLAPLACIAN_H
+
+//#define RIGID_DEFORM
+
+struct Object;
+struct Mesh;
+struct bDeformGroup;
+
+#ifdef RIGID_DEFORM
+struct EditMesh;
+#endif
+
+/* Laplacian System */
+
+struct LaplacianSystem;
+typedef struct LaplacianSystem LaplacianSystem;
+
+LaplacianSystem *laplacian_construct_begin(int totvert, int totface);
+
+void laplacian_add_vertex(LaplacianSystem *sys, float *co, int pinned);
+void laplacian_add_triangle(LaplacianSystem *sys, int v1, int v2, int v3);
+
+void laplacian_construct_end(LaplacianSystem *sys);
+void laplacian_delete(LaplacianSystem *sys);
+
+void laplacian_begin_solve(LaplacianSystem *sys, int index);
+void laplacian_add_right_hand_side(LaplacianSystem *sys, int v, float value);
+int laplacian_system_solve(LaplacianSystem *sys);
+float laplacian_system_get_solution(int v);
+
+/* Heat Weighting */
+
+void heat_bone_weighting(struct Object *ob, struct Mesh *me, float (*verts)[3],
+	int numbones, struct bDeformGroup **dgrouplist,
+	struct bDeformGroup **dgroupflip, float (*root)[3], float (*tip)[3],
+	int *selected);
+
+#ifdef RIGID_DEFORM
+/* As-Rigid-As-Possible Deformation */
+
+void rigid_deform_begin(struct EditMesh *em);
+void rigid_deform_iteration(void);
+void rigid_deform_end(int cancel);
+#endif
+
+#endif
+
diff --git a/source/blender/include/BIF_poseobject.h b/source/blender/include/BIF_poseobject.h
index 835b808c9ba..a640d3abe84 100644
--- a/source/blender/include/BIF_poseobject.h
+++ b/source/blender/include/BIF_poseobject.h
@@ -58,7 +58,7 @@ void free_posebuf(void);
 void copy_posebuf (void);
 void paste_posebuf (int flip);
 
-void pose_adds_vgroups(struct Object *meshobj);
+void pose_adds_vgroups(struct Object *meshobj, int heatweights);
 
 void pose_calculate_path(struct Object *ob);
 void pose_clear_paths(struct Object *ob);
diff --git a/source/blender/src/editarmature.c b/source/blender/src/editarmature.c
index 949decb1171..01488836570 100644
--- a/source/blender/src/editarmature.c
+++ b/source/blender/src/editarmature.c
@@ -67,6 +67,7 @@
 #include "BKE_constraint.h"
 #include "BKE_deform.h"
 #include "BKE_depsgraph.h"
+#include "BKE_derivedmesh.h"
 #include "BKE_global.h"
 #include "BKE_main.h"
 #include "BKE_object.h"
@@ -82,6 +83,8 @@
 #include "BIF_gl.h"
 #include "BIF_graphics.h"
 #include "BIF_interface.h"
+#include "BIF_meshlaplacian.h"
+#include "BIF_meshtools.h"
 #include "BIF_poseobject.h"
 #include "BIF_mywindow.h"
 #include "BIF_resources.h"
@@ -2334,14 +2337,16 @@ static int bone_skinnable(Object *ob, Bone *bone, void *data)
      */
     Bone ***hbone;
 
-    if (!(bone->flag & BONE_NO_DEFORM)) {
-		if (data != NULL) {
-			hbone = (Bone ***) data;
-            **hbone = bone;
-            ++*hbone;
-        }
-        return 1;
-    }
+	if(!(G.f & G_WEIGHTPAINT) || !(bone->flag & BONE_HIDDEN_P)) {
+		if (!(bone->flag & BONE_NO_DEFORM)) {
+			if (data != NULL) {
+				hbone = (Bone ***) data;
+				**hbone = bone;
+				++*hbone;
+			}
+			return 1;
+		}
+	}
     return 0;
 }
 
@@ -2387,57 +2392,97 @@ static int dgroup_skinnable(Object *ob, Bone *bone, void *data)
      */
     bDeformGroup ***hgroup, *defgroup;
 
-   if (!(bone->flag & BONE_NO_DEFORM)) {
-        if ( !(defgroup = get_named_vertexgroup(ob, bone->name)) ) {
-            defgroup = add_defgroup_name(ob, bone->name);
-        }
-
-        if (data != NULL) {
-            hgroup = (bDeformGroup ***) data;
-            **hgroup = defgroup;
-            ++*hgroup;
-        }
-        return 1;
-    }
+	if(!(G.f & G_WEIGHTPAINT) || !(bone->flag & BONE_HIDDEN_P)) {
+	   if (!(bone->flag & BONE_NO_DEFORM)) {
+			if ( !(defgroup = get_named_vertexgroup(ob, bone->name)) ) {
+				defgroup = add_defgroup_name(ob, bone->name);
+			}
+
+			if (data != NULL) {
+				hgroup = (bDeformGroup ***) data;
+				**hgroup = defgroup;
+				++*hgroup;
+			}
+			return 1;
+		}
+	}
     return 0;
 }
 
-static void add_verts_to_closest_dgroup(Object *ob, Object *par)
-{
-    /* This function implements a crude form of 
-     * auto-skinning: vertices are assigned to the
-     * deformation groups associated with bones based
-     * on thier proximity to a bone. Every vert is
-     * given a weight of 1.0 to the weight group
-     * cooresponding to the bone that it is
-     * closest to. The vertex may also be assigned to
-     * a deformation group associated to a bone
-     * that is within 10% of the mninimum distance
-     * between the bone and the nearest vert -- the
-     * cooresponding weight will fall-off to zero
-     * as the distance approaches the 10% tolerance mark.
-	 * If the mesh has subsurf enabled then the verts
-	 * on the subsurf limit surface is used to generate 
-	 * the weights rather than the verts on the cage
-	 * mesh.
+static void add_vgroups__mapFunc(void *userData, int index, float *co, float *no_f, short *no_s)
+{
+	/* DerivedMesh mapFunc for getting final coords in weight paint mode */
+
+	float (*verts)[3] = userData;
+	VECCOPY(verts[index], co);
+}
+
+static void envelope_bone_weighting(Object *ob, Mesh *mesh, float (*verts)[3], int numbones, Bone **bonelist, bDeformGroup **dgrouplist, bDeformGroup **dgroupflip, float (*root)[3], float (*tip)[3], int *selected)
+{
+	/* Create vertex group weights from envelopes */
+
+	Bone *bone;
+	bDeformGroup *dgroup;
+	float distance;
+	int i, iflip, j;
+
+	/* for each vertex in the mesh */
+	for (i=0; i < mesh->totvert; i++) {
+		iflip = (dgroupflip)? mesh_get_x_mirror_vert(ob, i): 0;
+
+		/* for each skinnable bone */
+		for (j=0; j < numbones; ++j) {
+			if(!selected[j])
+				continue;
+
+			bone = bonelist[j];
+			dgroup = dgrouplist[j];
+
+			/* store the distance-factor from the vertex to the bone */
+			distance = distfactor_to_bone (verts[i], root[j], tip[j],
+				bone->rad_head, bone->rad_tail, bone->dist);
+
+			/* add the vert to the deform group if weight!=0.0 */
+			if (distance!=0.0)
+				add_vert_to_defgroup (ob, dgroup, i, distance, WEIGHT_REPLACE);
+			else
+				remove_vert_defgroup (ob, dgroup, i);
+
+			/* do same for mirror */
+			if (dgroupflip && dgroupflip[j] && iflip >= 0) {
+				if (distance!=0.0)
+					add_vert_to_defgroup (ob, dgroupflip[j], iflip, distance,
+						WEIGHT_REPLACE);
+				else
+					remove_vert_defgroup (ob, dgroupflip[j], iflip);
+			}
+		}
+	}
+}
+
+void add_verts_to_dgroups(Object *ob, Object *par, int heat, int mirror)
+{
+	/* This functions implements the automatic computation of vertex group
+	 * weights, either through envelopes or using a heat equilibrium.
 	 *
-	 * ("Limit surface" = same amount of vertices as mesh, but vertices 
-     *   moved to the subsurfed position, like for 'optimal').
+	 * This function can be called both when parenting a mesh to an armature,
+	 * or in weightpaint + posemode. In the latter case selection is taken
+	 * into account and vertex weights can be mirrored.
+	 *
+	 * The mesh vertex positions used are either the final deformed coords
+	 * from the derivedmesh in weightpaint mode, the final subsurf coords
+	 * when parenting, or simply the original mesh coords.
      */
 
     bArmature *arm;
     Bone **bonelist, **bonehandle, *bone;
-    bDeformGroup **dgrouplist, **dgrouphandle, *defgroup;
-    float *distance;
-    float   root[3];
-    float   tip[3];
-    float real_co[3];
-	float *subverts = NULL;
-    float *subvert;
-    Mesh  *mesh;
-    MVert *vert;
-
-    int numbones, i, j;
+    bDeformGroup **dgrouplist, **dgroupflip, **dgrouphandle;
+	bDeformGroup *dgroup, *curdg;
+    Mesh *mesh;
+    float (*root)[3], (*tip)[3], (*verts)[3];
+	int *selected;
+    int numbones, vertsfilled = 0, i, j;
+	int wpmode = (G.f & G_WEIGHTPAINT);
 
     /* If the parent object is not an armature exit */
     arm = get_armature(par);
@@ -2445,97 +2490,113 @@ static void add_verts_to_closest_dgroup(Object *ob, Object *par)
         return;
 
     /* count the number of skinnable bones */
-    numbones = bone_looper(ob, arm->bonebase.first, NULL,
-                                  bone_skinnable);
+    numbones = bone_looper(ob, arm->bonebase.first, NULL, bone_skinnable);
+	
+	if (numbones == 0)
+		return;
 	
     /* create an array of pointer to bones that are skinnable
-     * and fill it with all of the skinnable bones
-     */
-    bonelist = MEM_mallocN(numbones*sizeof(Bone *), "bonelist");
+     * and fill it with all of the skinnable bones */
+    bonelist = MEM_callocN(numbones*sizeof(Bone *), "bonelist");
     bonehandle = bonelist;
-    bone_looper(ob, arm->bonebase.first, &bonehandle,
-                       bone_skinnable);
+    bone_looper(ob, arm->bonebase.first, &bonehandle, bone_skinnable);
 
     /* create an array of pointers to the deform groups that
      * coorespond to the skinnable bones (creating them
-     * as necessary.
-     */
-    dgrouplist = MEM_mallocN(numbones*sizeof(bDeformGroup *), "dgrouplist");
+     * as necessary. */
+    dgrouplist = MEM_callocN(numbones*sizeof(bDeformGroup *), "dgrouplist");
+    dgroupflip = MEM_callocN(numbones*sizeof(bDeformGroup *), "dgroupflip");
+
     dgrouphandle = dgrouplist;
-    bone_looper(ob, arm->bonebase.first, &dgrouphandle,
-                       dgroup_skinnable);
+    bone_looper(ob, arm->bonebase.first, &dgrouphandle, dgroup_skinnable);
 
-    /* create an array of floats that will be used for each vert
-     * to hold the distance-factor to each bone.
-     */
-    distance = MEM_mallocN(numbones*sizeof(float), "distance");
+    /* create an array of root and tip positions transformed into
+	 * global coords */
+    root = MEM_callocN(numbones*sizeof(float)*3, "root");
+    tip = MEM_callocN(numbones*sizeof(float)*3, "tip");
+	selected = MEM_callocN(numbones*sizeof(int), "selected");
 
-    mesh = (Mesh*)ob->data;
+	for (j=0; j < numbones; ++j) {
+   		bone = bonelist[j];
+		dgroup = dgrouplist[j];
 
-	/* Is subsurf on? Lets use the verts on the limit surface then */
-	if (modifiers_findByType(ob, eModifierType_Subsurf)) {
-		subverts = MEM_mallocN(3*mesh->totvert*sizeof(float), "subverts");
-		subsurf_calculate_limit_positions(mesh, (void *)subverts);	/* (ton) made void*, dunno how to cast */
-	}
+		/* compute root and tip */
+		VECCOPY(root[j], bone->arm_head);
+		Mat4MulVecfl(par->obmat, root[j]);
 
-    /* for each vertex in the mesh ...
-     */
-    for ( i=0 ; i < mesh->totvert ; ++i ) {
-        /* get the vert in global coords
-         */
-		
-		if (subverts) {
-			subvert = subverts + i*3;
-			VECCOPY (real_co, subvert);
+		VECCOPY(tip[j], bone->arm_tail);
+		Mat4MulVecfl(par->obmat, tip[j]);
+
+		/* set selected */
+		if(wpmode) {
+			if ((arm->layer & bone->layer) && (bone->flag & BONE_SELECTED))
+				selected[j] = 1;
 		}
-		else {
-			vert = mesh->mvert + i;
-			VECCOPY (real_co, vert->co);
+		else
+			selected[j] = 1;
+
+		/* find flipped group */
+		if(mirror) {
+			char name[32];
+			
+			BLI_strncpy(name, dgroup->name, 32);
+			// 0 = don't strip off number extensions
+			bone_flip_name(name, 0);
+
+			for (curdg = ob->defbase.first; curdg; curdg=curdg->next)
+				if (!strcmp(curdg->name, name))
+					break;
+
+			dgroupflip[j] = curdg;
 		}
-        Mat4MulVecfl(ob->obmat, real_co);
+	}
 
+	/* create verts */
+    mesh = (Mesh*)ob->data;
+	verts = MEM_callocN(mesh->totvert*sizeof(*verts), "closestboneverts");
 
-        /* for each skinnable bone ...
-         */
-        for (j=0; j < numbones; ++j) {
-            bone = bonelist[j];
-
-            /* get the root of the bone in global coords
-             */
-			VECCOPY(root, bone->arm_head);
-			Mat4MulVecfl(par->obmat, root);
-
-            /* get the tip of the bone in global coords
-             */
-			VECCOPY(tip, bone->arm_tail);
-            Mat4MulVecfl(par->obmat, tip);
-
-            /* store the distance-factor from the vertex to
-             * the bone
-             */
-			distance[j]= distfactor_to_bone (real_co, root, tip, bone->rad_head, bone->rad_tail, bone->dist);
-        }
-
-        /* for each deform group ...
-         */
-        for (j=0; j < numbones; ++j) {
-            defgroup = dgrouplist[j];
-
-            /* add the vert to the deform group if weight!=0.0
-             */
-            if (distance[j]!=0.0)
-                add_vert_to_defgroup (ob, defgroup, i, distance[j], WEIGHT_REPLACE);
-            else
-                remove_vert_defgroup (ob, defgroup, i);
-        }
-    }
+	if (wpmode) {
+		/* if in weight paint mode, use final verts from derivedmesh */
+		DerivedMesh *dm = mesh_get_derived_final(ob, CD_MASK_BAREMESH);
 
-    /* free the memory allocated
-     */
+		if(dm->foreachMappedVert) {
+			dm->foreachMappedVert(dm, add_vgroups__mapFunc, (void*)verts);
+			vertsfilled = 1;
+		}
+
+		dm->release(dm);
+	}
+	else if (modifiers_findByType(ob, eModifierType_Subsurf)) {
+		/* is subsurf on? Lets use the verts on the limit surface then.
+	 	 * = same amount of vertices as mesh, but vertices  moved to the
+		 * subsurfed position, like for 'optimal'. */
+		subsurf_calculate_limit_positions(mesh, verts);
+		vertsfilled = 1;
+	}
+
+	/* transform verts to global space */
+	for (i=0; i < mesh->totvert; i++) {
+		if (!vertsfilled)
+			VECCOPY(verts[i], mesh->mvert[i].co)
+		Mat4MulVecfl(ob->obmat, verts[i]);
+	}
+
+	/* compute the weights based on gathered vertices and bones */
+	if (heat)
+		heat_bone_weighting(ob, mesh, verts, numbones, dgrouplist, dgroupflip,
+			root, tip, selected);
+	else
+		envelope_bone_weighting(ob, mesh, verts, numbones, bonelist, dgrouplist,
+			dgroupflip, root, tip, selected);
+
+    /* free the memory allocated */
     MEM_freeN(bonelist);
     MEM_freeN(dgrouplist);
-    MEM_freeN(distance);
-	if (subverts) MEM_freeN(subverts);
+	MEM_freeN(dgroupflip);
+	MEM_freeN(root);
+	MEM_freeN(tip);
+	MEM_freeN(selected);
+	MEM_freeN(verts);
 }
 
 void create_vgroups_from_armature(Object *ob, Object *par)
@@ -2557,7 +2618,8 @@ void create_vgroups_from_armature(Object *ob, Object *par)
     mode= pupmenu("Create Vertex Groups? %t|"
 				  "Don't Create Groups %x1|"
 				  "Name Groups %x2|"
-                  "Create From Closest Bones %x3");
+                  "Create From Envelopes %x3|"
+				  "Create From Bone Heat %x4|");
 	switch (mode){
 	case 2:
 		/* Traverse the bone list, trying to create empty vertex 
@@ -2571,11 +2633,12 @@ void create_vgroups_from_armature(Object *ob, Object *par)
 		break;
 
 	case 3:
+	case 4:
 		/* Traverse the bone list, trying to create vertex groups 
 		 * that are populated with the vertices for which the
 		 * bone is closest.
 		 */
-		add_verts_to_closest_dgroup(ob, par);
+		add_verts_to_dgroups(ob, par, (mode == 4), 0);
 		break;
 
 	}
diff --git a/source/blender/src/editobject.c b/source/blender/src/editobject.c
index 22ca7a9494d..c4736a84837 100644
--- a/source/blender/src/editobject.c
+++ b/source/blender/src/editobject.c
@@ -2263,12 +2263,12 @@ void special_editmenu(void)
 		}
 		else if(G.f & G_WEIGHTPAINT) {
 			Object *par= modifiers_isDeformedByArmature(ob);
+
 			if(par && (par->flag & OB_POSEMODE)) {
-				nr= pupmenu("Specials%t|Apply Bone Envelopes to VertexGroups %x1");
-				if(nr==1) {
-					pose_adds_vgroups(ob);
-					BIF_undo_push("Apply Bone Envelopes to VertexGroups");
-				}
+				nr= pupmenu("Specials%t|Apply Bone Envelopes to Vertex Groups %x1|Apply Bone Heat Weights to Vertex Groups %x2");
+
+				if(nr==1 || nr==2)
+					pose_adds_vgroups(ob, (nr == 2));
 			}
 		}
 		else {
diff --git a/source/blender/src/header_view3d.c b/source/blender/src/header_view3d.c
index 4dc8759f859..05b4df137ed 100644
--- a/source/blender/src/header_view3d.c
+++ b/source/blender/src/header_view3d.c
@@ -4129,10 +4129,9 @@ static uiBlock *view3d_tpaintmenu(void *arg_unused)
 static void do_view3d_wpaintmenu(void *arg, int event)
 {
 	Object *ob= OBACT;
-	Object *par= modifiers_isDeformedByArmature(ob);
 	
-	/* events >= 2 are registered bpython scripts */
-	if (event >= 3) BPY_menu_do_python(PYMENU_WEIGHTPAINT, event - 3);
+	/* events >= 3 are registered bpython scripts */
+	if (event >= 4) BPY_menu_do_python(PYMENU_WEIGHTPAINT, event - 4);
 	
 	switch(event) {
 	case 0: /* undo weight painting */
@@ -4142,13 +4141,11 @@ static void do_view3d_wpaintmenu(void *arg, int event)
 		clear_wpaint_selectedfaces();
 		break;
 	case 2: /* vgroups from envelopes */
-		if(par && (par->flag & OB_POSEMODE))  {
-			pose_adds_vgroups(ob);
-			BIF_undo_push("Apply Bone Envelopes to VertexGroups");
-		} else {
-			error("The active object must have a deforming armature in pose mode");
-		}
-	break;
+		pose_adds_vgroups(ob, 0);
+		break;
+	case 3: /* vgroups from bone heat */
+		pose_adds_vgroups(ob, 1);
+		break;
 	}
 	allqueue(REDRAWVIEW3D, 0);
 }
@@ -4164,6 +4161,10 @@ static uiBlock *view3d_wpaintmenu(void *arg_unused)
 	uiBlockSetButmFunc(block, do_view3d_wpaintmenu, NULL);
 	
 	uiDefIconTextBut(block, BUTM, 1, ICON_BLANK1, "Undo Weight Painting|U",		0, yco-=20, menuwidth, 19, NULL, 0.0, 0.0, 1, 0, "");
+
+	uiDefBut(block, SEPR, 0, "",				0, yco-=6, menuwidth, 6, NULL, 0.0, 0.0, 0, 0, "");
+
+	uiDefIconTextBut(block, BUTM, 1, ICON_BLANK1, "Apply Bone Heat Weights to Vertex Groups|W, 2",		0, yco-=20, menuwidth, 19, NULL, 0.0, 0.0, 1, 3, "");
 	uiDefIconTextBut(block, BUTM, 1, ICON_BLANK1, "Apply Bone Envelopes to Vertex Groups|W, 1",		0, yco-=20, menuwidth, 19, NULL, 0.0, 0.0, 1, 2, "");
 	
 	uiDefBut(block, SEPR, 0, "",				0, yco-=6, menuwidth, 6, NULL, 0.0, 0.0, 0, 0, "");
@@ -4174,11 +4175,11 @@ static uiBlock *view3d_wpaintmenu(void *arg_unused)
 		menunr++;
 	}
 	
-	/* note that we account for the 3 previous entries with i+3:
+	/* note that we account for the 4 previous entries with i+4:
 	even if the last item isnt displayed, it dosent matter */
 	for (pym = BPyMenuTable[PYMENU_WEIGHTPAINT]; pym; pym = pym->next, i++) {
 		uiDefIconTextBut(block, BUTM, 1, ICON_PYTHON, pym->name, 0, yco-=20,
-			menuwidth, 19, NULL, 0.0, 0.0, 1, i+3,
+			menuwidth, 19, NULL, 0.0, 0.0, 1, i+4,
 			pym->tooltip?pym->tooltip:pym->filename);
 	}
 	
diff --git a/source/blender/src/meshlaplacian.c b/source/blender/src/meshlaplacian.c
new file mode 100644
index 00000000000..a660b42996e
--- /dev/null
+++ b/source/blender/src/meshlaplacian.c
@@ -0,0 +1,906 @@
+/**
+ * $Id$
+ *
+ * ***** 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 use 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) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): none yet.
+ *
+ * ***** END GPL/BL DUAL LICENSE BLOCK *****
+ * meshlaplacian.c: Algorithms using the mesh laplacian.
+ */
+
+#include <string.h>
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_listBase.h"
+#include "DNA_object_types.h"
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+
+#include "BLI_arithb.h"
+#include "BLI_edgehash.h"
+
+#include "BKE_utildefines.h"
+
+#include "BIF_editdeform.h"
+#include "BIF_meshlaplacian.h"
+#include "BIF_meshtools.h"
+#include "BIF_toolbox.h"
+
+#ifdef RIGID_DEFORM
+#include "BLI_editVert.h"
+#include "BLI_polardecomp.h"
+#endif
+
+#include "RE_raytrace.h"
+
+#include "ONL_opennl.h"
+
+/************************** Laplacian System *****************************/
+
+struct LaplacianSystem {
+	NLContext context;	/* opennl context */
+
+	int totvert, totface;
+
+	float **verts;			/* vertex coordinates */
+	float *varea;			/* vertex weights for laplacian computation */
+	char *vpinned;			/* vertex pinning */
+	int (*faces)[3];		/* face vertex indices */
+	float (*fweights)[3];	/* cotangent weights per face */
+
+	int areaweights;		/* use area in cotangent weights? */
+	int storeweights;		/* store cotangent weights in fweights */
+	int nlbegun;			/* nlBegin(NL_SYSTEM/NL_MATRIX) done */
+
+	EdgeHash *edgehash;		/* edge hash for construction */
+
+	struct HeatWeighting {
+		Mesh *mesh;
+		float (*verts)[3];	/* vertex coordinates */
+		float (*vnors)[3];	/* vertex normals */
+
+		float (*root)[3];	/* bone root */
+		float (*tip)[3];	/* bone tip */
+		int numbones;
+
+		float *H;			/* diagonal H matrix */
+		float *p;			/* values from all p vectors */
+		float *mindist;		/* minimum distance to a bone for all vertices */
+		
+		RayTree *raytree;	/* ray tracing acceleration structure */
+		MFace **vface;		/* a face that the vertex belongs to */
+	} heat;
+
+#ifdef RIGID_DEFORM
+	struct RigidDeformation {
+		EditMesh *mesh;
+
+		float (*R)[3][3];
+		float (*rhs)[3];
+		float (*origco)[3];
+		int thrownerror;
+	} rigid;
+#endif
+};
+
+/* Laplacian matrix construction */
+
+/* Computation of these weights for the laplacian is based on:
+   "Discrete Differential-Geometry Operators for Triangulated 2-Manifolds",
+   Meyer et al, 2002. Section 3.5, formula (8).
+   
+   We do it a bit different by going over faces instead of going over each
+   vertex and adjacent faces, since we don't store this adjacency. Also, the
+   formulas are tweaked a bit to work for non-manifold meshes. */
+
+static void laplacian_increase_edge_count(EdgeHash *edgehash, int v1, int v2)
+{
+	void **p = BLI_edgehash_lookup_p(edgehash, v1, v2);
+
+	if(p)
+		*p = (void*)((long)*p + (long)1);
+	else
+		BLI_edgehash_insert(edgehash, v1, v2, (void*)(long)1);
+}
+
+static int laplacian_edge_count(EdgeHash *edgehash, int v1, int v2)
+{
+	return (int)(long)BLI_edgehash_lookup(edgehash, v1, v2);
+}
+
+static float cotan_weight(float *v1, float *v2, float *v3)
+{
+	float a[3], b[3], c[3], clen;
+
+	VecSubf(a, v2, v1);
+	VecSubf(b, v3, v1);
+	Crossf(c, a, b);
+
+	clen = VecLength(c);
+
+	if (clen == 0.0f)
+		return 0.0f;
+	
+	return Inpf(a, b)/clen;
+}
+
+static void laplacian_triangle_area(LaplacianSystem *sys, int i1, int i2, int i3)
+{
+	float t1, t2, t3, len1, len2, len3, area;
+	float *varea= sys->varea, *v1, *v2, *v3;
+	int obtuse = 0;
+
+	v1= sys->verts[i1];
+	v2= sys->verts[i2];
+	v3= sys->verts[i3];
+
+	t1= cotan_weight(v1, v2, v3);
+	t2= cotan_weight(v2, v3, v1);
+	t3= cotan_weight(v3, v1, v2);
+
+	if(VecAngle3(v2, v1, v3) > 90) obtuse= 1;
+	else if(VecAngle3(v1, v2, v3) > 90) obtuse= 2;
+	else if(VecAngle3(v1, v3, v2) > 90) obtuse= 3;
+
+	if (obtuse > 0) {
+		area= AreaT3Dfl(v1, v2, v3);
+
+		varea[i1] += (obtuse == 1)? area: area*0.5;
+		varea[i2] += (obtuse == 2)? area: area*0.5;
+		varea[i3] += (obtuse == 3)? area: area*0.5;
+	}
+	else {
+		len1= VecLenf(v2, v3);
+		len2= VecLenf(v1, v3);
+		len3= VecLenf(v1, v2);
+
+		t1 *= len1*len1;
+		t2 *= len2*len2;
+		t3 *= len3*len3;
+
+		varea[i1] += (t2 + t3)*0.25f;
+		varea[i2] += (t1 + t3)*0.25f;
+		varea[i3] += (t1 + t2)*0.25f;
+	}
+}
+
+static void laplacian_triangle_weights(LaplacianSystem *sys, int f, int i1, int i2, int i3)
+{
+	float t1, t2, t3;
+	float *varea= sys->varea, *v1, *v2, *v3;
+
+	v1= sys->verts[i1];
+	v2= sys->verts[i2];
+	v3= sys->verts[i3];
+
+	/* instead of *0.5 we divided by the number of faces of the edge, it still
+	   needs to be varified that this is indeed the correct thing to do! */
+	t1= cotan_weight(v1, v2, v3)/laplacian_edge_count(sys->edgehash, i2, i3);
+	t2= cotan_weight(v2, v3, v1)/laplacian_edge_count(sys->edgehash, i3, i1);
+	t3= cotan_weight(v3, v1, v2)/laplacian_edge_count(sys->edgehash, i1, i2);
+
+	nlMatrixAdd(i1, i1, (t2+t3)*varea[i1]);
+	nlMatrixAdd(i2, i2, (t1+t3)*varea[i2]);
+	nlMatrixAdd(i3, i3, (t1+t2)*varea[i3]);
+
+	nlMatrixAdd(i1, i2, -t3*varea[i1]);
+	nlMatrixAdd(i2, i1, -t3*varea[i2]);
+
+	nlMatrixAdd(i2, i3, -t1*varea[i2]);
+	nlMatrixAdd(i3, i2, -t1*varea[i3]);
+
+	nlMatrixAdd(i3, i1, -t2*varea[i3]);
+	nlMatrixAdd(i1, i3, -t2*varea[i1]);
+
+	if(sys->storeweights) {
+		sys->fweights[f][0]= t1*varea[i1];
+		sys->fweights[f][1]= t2*varea[i2];
+		sys->fweights[f][2]= t3*varea[i3];
+	}
+}
+
+LaplacianSystem *laplacian_system_construct_begin(int totvert, int totface)
+{
+	LaplacianSystem *sys;
+
+	sys= MEM_callocN(sizeof(LaplacianSystem), "LaplacianSystem");
+
+	sys->verts= MEM_callocN(sizeof(float*)*totvert, "LaplacianSystemVerts");
+	sys->vpinned= MEM_callocN(sizeof(char)*totvert, "LaplacianSystemVpinned");
+	sys->faces= MEM_callocN(sizeof(int)*3*totface, "LaplacianSystemFaces");
+
+	sys->totvert= 0;
+	sys->totface= 0;
+
+	sys->areaweights= 1;
+	sys->storeweights= 0;
+
+	/* create opennl context */
+	nlNewContext();
+	nlSolverParameteri(NL_NB_VARIABLES, totvert);
+
+	sys->context= nlGetCurrent();
+
+	return sys;
+}
+
+void laplacian_add_vertex(LaplacianSystem *sys, float *co, int pinned)
+{
+	sys->verts[sys->totvert]= co;
+	sys->vpinned[sys->totvert]= pinned;
+	sys->totvert++;
+}
+
+void laplacian_add_triangle(LaplacianSystem *sys, int v1, int v2, int v3)
+{
+	sys->faces[sys->totface][0]= v1;
+	sys->faces[sys->totface][1]= v2;
+	sys->faces[sys->totface][2]= v3;
+	sys->totface++;
+}
+
+void laplacian_system_construct_end(LaplacianSystem *sys)
+{
+	int (*face)[3];
+	int a, totvert=sys->totvert, totface=sys->totface;
+
+	laplacian_begin_solve(sys, 0);
+
+	sys->varea= MEM_callocN(sizeof(float)*totvert, "LaplacianSystemVarea");
+
+	sys->edgehash= BLI_edgehash_new();
+	for(a=0, face=sys->faces; a<sys->totface; a++, face++) {
+		laplacian_increase_edge_count(sys->edgehash, (*face)[0], (*face)[1]);
+		laplacian_increase_edge_count(sys->edgehash, (*face)[1], (*face)[2]);
+		laplacian_increase_edge_count(sys->edgehash, (*face)[2], (*face)[0]);
+	}
+
+	if(sys->areaweights)
+		for(a=0, face=sys->faces; a<sys->totface; a++, face++)
+			laplacian_triangle_area(sys, (*face)[0], (*face)[1], (*face)[2]);
+	
+	for(a=0; a<totvert; a++) {
+		if(sys->areaweights) {
+			if(sys->varea[a] != 0.0f)
+				sys->varea[a]= 0.5f/sys->varea[a];
+		}
+		else
+			sys->varea[a]= 1.0f;
+
+		/* for heat weighting */
+		if(sys->heat.H)
+			nlMatrixAdd(a, a, sys->heat.H[a]);
+	}
+
+	if(sys->storeweights)
+		sys->fweights= MEM_callocN(sizeof(float)*3*totface, "LaplacianFWeight");
+	
+	for(a=0, face=sys->faces; a<totface; a++, face++)
+		laplacian_triangle_weights(sys, a, (*face)[0], (*face)[1], (*face)[2]);
+
+	MEM_freeN(sys->faces);
+	sys->faces= NULL;
+
+	if(sys->varea) {
+		MEM_freeN(sys->varea);
+		sys->varea= NULL;
+	}
+
+	BLI_edgehash_free(sys->edgehash, NULL);
+	sys->edgehash= NULL;
+}
+
+void laplacian_system_delete(LaplacianSystem *sys)
+{
+	if(sys->verts) MEM_freeN(sys->verts);
+	if(sys->varea) MEM_freeN(sys->varea);
+	if(sys->vpinned) MEM_freeN(sys->vpinned);
+	if(sys->faces) MEM_freeN(sys->faces);
+	if(sys->fweights) MEM_freeN(sys->fweights);
+
+	nlDeleteContext(sys->context);
+	MEM_freeN(sys);
+}
+
+void laplacian_begin_solve(LaplacianSystem *sys, int index)
+{
+	int a;
+
+	if (!sys->nlbegun) {
+		nlBegin(NL_SYSTEM);
+
+		if(index >= 0) {
+			for(a=0; a<sys->totvert; a++) {
+				if(sys->vpinned[a]) {
+					nlSetVariable(a, sys->verts[a][index]);
+					nlLockVariable(a);
+				}
+			}
+		}
+
+		nlBegin(NL_MATRIX);
+		sys->nlbegun = 1;
+	}
+}
+
+void laplacian_add_right_hand_side(LaplacianSystem *sys, int v, float value)
+{
+	nlRightHandSideAdd(v, value);
+}
+
+int laplacian_system_solve(LaplacianSystem *sys)
+{
+	nlEnd(NL_MATRIX);
+	nlEnd(NL_SYSTEM);
+	sys->nlbegun = 0;
+
+	//nlPrintMatrix();
+
+	return nlSolveAdvanced(NULL, NL_TRUE);
+}
+
+float laplacian_system_get_solution(int v)
+{
+	return nlGetVariable(v);
+}
+
+/************************* Heat Bone Weighting ******************************/
+/* From "Automatic Rigging and Animation of 3D Characters"
+         Ilya Baran and Jovan Popovic, SIGGRAPH 2007 */
+
+#define C_WEIGHT			1.0f
+#define WEIGHT_LIMIT		0.05f
+#define DISTANCE_EPSILON	1e-4f
+
+/* Raytracing for vertex to bone visibility */
+
+static LaplacianSystem *HeatSys = NULL;
+
+static void heat_ray_coords_func(RayFace *face, float **v1, float **v2, float **v3, float **v4)
+{
+	MFace *mface= (MFace*)face;
+	float (*verts)[3]= HeatSys->heat.verts;
+
+	*v1= verts[mface->v1];
+	*v2= verts[mface->v2];
+	*v3= verts[mface->v3];
+	*v4= (mface->v4)? verts[mface->v4]: NULL;
+}
+
+static int heat_ray_check_func(Isect *is, RayFace *face)
+{
+	float *v1, *v2, *v3, *v4, nor[3];
+
+	/* don't intersect if the ray faces along the face normal */
+	heat_ray_coords_func(face, &v1, &v2, &v3, &v4);
+
+	if(v4) CalcNormFloat4(v1, v2, v3, v4, nor);
+	else CalcNormFloat(v1, v2, v3, nor);
+	
+	return (INPR(nor, is->vec) < 0);
+}
+
+static void heat_ray_tree_create(LaplacianSystem *sys)
+{
+	Mesh *me = sys->heat.mesh;
+	RayTree *tree;
+	MFace *mface;
+	float min[3], max[3];
+	int a;
+
+	/* create a raytrace tree from the mesh */
+	INIT_MINMAX(min, max);
+
+	for(a=0; a<me->totvert; a++)
+		DO_MINMAX(sys->heat.verts[a], min, max);
+
+	tree= RE_ray_tree_create(64, me->totface, min, max,
+		heat_ray_coords_func, heat_ray_check_func);
+	
+	sys->heat.vface= MEM_callocN(sizeof(MFace*)*me->totvert, "HeatVFaces");
+
+	HeatSys= sys;
+
+	for(a=0, mface=me->mface; a<me->totface; a++, mface++) {
+		RE_ray_tree_add_face(tree, mface);
+
+		sys->heat.vface[mface->v1]= mface;
+		sys->heat.vface[mface->v2]= mface;
+		sys->heat.vface[mface->v3]= mface;
+		if(mface->v4) sys->heat.vface[mface->v4]= mface;
+	}
+
+	HeatSys= NULL;
+	
+	RE_ray_tree_done(tree);
+
+	sys->heat.raytree= tree;
+}
+
+static int heat_ray_bone_visible(LaplacianSystem *sys, int vertex, int bone)
+{
+	Isect isec;
+	MFace *mface;
+	float dir[3];
+	int visible;
+
+	mface= sys->heat.vface[vertex];
+	if(!mface)
+		return 1;
+
+	/* setup isec */
+	isec.mode= RE_RAY_SHADOW;
+	isec.lay= -1;
+	isec.face_last= NULL;
+	isec.faceorig= mface;
+
+	VECCOPY(isec.start, sys->heat.verts[vertex]);
+	PclosestVL3Dfl(isec.end, isec.start,
+		sys->heat.root[bone], sys->heat.tip[bone]);
+
+	/* add an extra offset to the start position to avoid self intersection */
+	VECSUB(dir, isec.end, isec.start);
+	Normalize(dir);
+	VecMulf(dir, 1e-5);
+	VecAddf(isec.start, isec.start, dir);
+	
+	HeatSys= sys;
+	visible= !RE_ray_tree_intersect(sys->heat.raytree, &isec);
+	HeatSys= NULL;
+
+	return visible;
+}
+
+static float heat_bone_distance(LaplacianSystem *sys, int vertex, int bone)
+{
+	float closest[3], d[3], dist, cosine;
+	
+	/* compute euclidian distance */
+	PclosestVL3Dfl(closest, sys->heat.verts[vertex],
+		sys->heat.root[bone], sys->heat.tip[bone]);
+
+	VecSubf(d, sys->heat.verts[vertex], closest);
+	dist= Normalize(d);
+
+	/* if the vertex normal does not point along the bone, increase distance */
+	cosine= INPR(d, sys->heat.vnors[vertex]);
+
+	return dist/(0.5f*(cosine + 1.001f));
+}
+
+static int heat_bone_closest(LaplacianSystem *sys, int vertex, int bone)
+{
+	float dist;
+	
+	dist= heat_bone_distance(sys, vertex, bone);
+
+	if(dist <= sys->heat.mindist[vertex]*(1.0f + DISTANCE_EPSILON))
+		if(heat_ray_bone_visible(sys, vertex, bone))
+			return 1;
+	
+	return 0;
+}
+
+static void heat_set_H(LaplacianSystem *sys, int vertex)
+{
+	float dist, mindist, h;
+	int j, numclosest = 0;
+
+	mindist= 1e10;
+
+	/* compute minimum distance */
+	for(j=0; j<sys->heat.numbones; j++) {
+		dist= heat_bone_distance(sys, vertex, j);
+
+		if(dist < mindist)
+			mindist= dist;
+	}
+
+	sys->heat.mindist[vertex]= mindist;
+
+	/* count number of bones with approximately this minimum distance */
+	for(j=0; j<sys->heat.numbones; j++)
+		if(heat_bone_closest(sys, vertex, j))
+			numclosest++;
+
+	sys->heat.p[vertex]= (numclosest > 0)? 1.0f/numclosest: 0.0f;
+
+	/* compute H entry */
+	if(numclosest > 0) {
+		if(mindist > 1e-5)
+			h= numclosest*C_WEIGHT/(mindist*mindist);
+		else
+			h= 1e10f;
+	}
+	else
+		h= 0.0f;
+	
+	sys->heat.H[vertex]= h;
+}
+
+void heat_calc_vnormals(LaplacianSystem *sys)
+{
+	float fnor[3];
+	int a, v1, v2, v3, (*face)[3];
+
+	sys->heat.vnors= MEM_callocN(sizeof(float)*3*sys->totvert, "HeatVNors");
+
+	for(a=0, face=sys->faces; a<sys->totface; a++, face++) {
+		v1= (*face)[0];
+		v2= (*face)[1];
+		v3= (*face)[2];
+
+		CalcNormFloat(sys->verts[v1], sys->verts[v2], sys->verts[v3], fnor);
+		
+		VecAddf(sys->heat.vnors[v1], sys->heat.vnors[v1], fnor);
+		VecAddf(sys->heat.vnors[v2], sys->heat.vnors[v2], fnor);
+		VecAddf(sys->heat.vnors[v3], sys->heat.vnors[v3], fnor);
+	}
+
+	for(a=0; a<sys->totvert; a++)
+		Normalize(sys->heat.vnors[a]);
+}
+
+static void heat_laplacian_create(LaplacianSystem *sys)
+{
+	Mesh *me = sys->heat.mesh;
+	MFace *mface;
+	int a;
+
+	/* heat specific definitions */
+	sys->heat.mindist= MEM_callocN(sizeof(float)*me->totvert, "HeatMinDist");
+	sys->heat.H= MEM_callocN(sizeof(float)*me->totvert, "HeatH");
+	sys->heat.p= MEM_callocN(sizeof(float)*me->totvert, "HeatP");
+
+	/* add verts and faces to laplacian */
+	for(a=0; a<me->totvert; a++)
+		laplacian_add_vertex(sys, sys->heat.verts[a], 0);
+
+	for(a=0, mface=me->mface; a<me->totface; a++, mface++) {
+		laplacian_add_triangle(sys, mface->v1, mface->v2, mface->v3);
+		if(mface->v4)
+			laplacian_add_triangle(sys, mface->v1, mface->v3, mface->v4);
+	}
+
+	/* for distance computation in set_H */
+	heat_calc_vnormals(sys);
+
+	for(a=0; a<me->totvert; a++)
+		heat_set_H(sys, a);
+}
+
+void heat_bone_weighting(Object *ob, Mesh *me, float (*verts)[3], int numbones, bDeformGroup **dgrouplist, bDeformGroup **dgroupflip, float (*root)[3], float (*tip)[3], int *selected)
+{
+	LaplacianSystem *sys;
+	MFace *mface;
+	float solution;
+	int a, aflip, totface, j, thrownerror = 0;
+
+	/* count triangles */
+	for(totface=0, a=0, mface=me->mface; a<me->totface; a++, mface++) {
+		totface++;
+		if(mface->v4) totface++;
+	}
+
+	/* create laplacian */
+	sys = laplacian_system_construct_begin(me->totvert, totface);
+
+	sys->heat.mesh= me;
+	sys->heat.verts= verts;
+	sys->heat.root= root;
+	sys->heat.tip= tip;
+	sys->heat.numbones= numbones;
+
+	heat_ray_tree_create(sys);
+	heat_laplacian_create(sys);
+
+	laplacian_system_construct_end(sys);
+
+	/* compute weights per bone */
+	for(j=0; j<numbones; j++) {
+		if(!selected[j])
+			continue;
+
+		laplacian_begin_solve(sys, -1);
+
+		for(a=0; a<me->totvert; a++)
+			if(heat_bone_closest(sys, a, j))
+				laplacian_add_right_hand_side(sys, a,
+					sys->heat.H[a]*sys->heat.p[a]);
+
+		if(laplacian_system_solve(sys)) {
+			for(a=0; a<me->totvert; a++) {
+				solution= laplacian_system_get_solution(a);
+
+				if(solution > WEIGHT_LIMIT)
+					add_vert_to_defgroup(ob, dgrouplist[j], a, solution,
+						WEIGHT_REPLACE);
+				else
+					remove_vert_defgroup(ob, dgrouplist[j], a);
+
+				/* do same for mirror */
+				aflip = (dgroupflip)? mesh_get_x_mirror_vert(ob, a): 0;
+				if (dgroupflip && dgroupflip[j] && aflip >= 0) {
+					if(solution > WEIGHT_LIMIT)
+						add_vert_to_defgroup(ob, dgroupflip[j], aflip,
+							solution, WEIGHT_REPLACE);
+					else
+						remove_vert_defgroup(ob, dgroupflip[j], aflip);
+				}
+			}
+		}
+		else if(!thrownerror) {
+			error("Bone Heat Weighting:"
+				" failed to find solution for one or more bones");
+			thrownerror= 1;
+			break;
+		}
+	}
+
+	/* free */
+	RE_ray_tree_free(sys->heat.raytree);
+	MEM_freeN(sys->heat.vface);
+
+	MEM_freeN(sys->heat.mindist);
+	MEM_freeN(sys->heat.H);
+	MEM_freeN(sys->heat.p);
+	MEM_freeN(sys->heat.vnors);
+
+	laplacian_system_delete(sys);
+}
+
+#ifdef RIGID_DEFORM
+/********************** As-Rigid-As-Possible Deformation ******************/
+/* From "As-Rigid-As-Possible Surface Modeling",
+        Olga Sorkine and Marc Alexa, ESGP 2007. */
+
+/* investigate:
+   - transpose R in orthogonal
+   - flipped normals and per face adding
+   - move cancelling to transform, make origco pointer
+*/
+
+static LaplacianSystem *RigidDeformSystem = NULL;
+
+static void rigid_add_half_edge_to_R(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
+{
+	float e[3], e_[3];
+	int i;
+
+	VecSubf(e, sys->rigid.origco[v1->tmp.l], sys->rigid.origco[v2->tmp.l]);
+	VecSubf(e_, v1->co, v2->co);
+
+	/* formula (5) */
+	for (i=0; i<3; i++) {
+		sys->rigid.R[v1->tmp.l][i][0] += w*e[0]*e_[i];
+		sys->rigid.R[v1->tmp.l][i][1] += w*e[1]*e_[i];
+		sys->rigid.R[v1->tmp.l][i][2] += w*e[2]*e_[i];
+	}
+}
+
+static void rigid_add_edge_to_R(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
+{
+	rigid_add_half_edge_to_R(sys, v1, v2, w);
+	rigid_add_half_edge_to_R(sys, v2, v1, w);
+}
+
+static void rigid_orthogonalize_R(float R[][3])
+{
+	HMatrix M, Q, S;
+
+	Mat4CpyMat3(M, R);
+	polar_decomp(M, Q, S);
+	Mat3CpyMat4(R, Q);
+}
+
+static void rigid_add_half_edge_to_rhs(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
+{
+	/* formula (8) */
+	float Rsum[3][3], rhs[3];
+
+	if (sys->vpinned[v1->tmp.l])
+		return;
+
+	Mat3AddMat3(Rsum, sys->rigid.R[v1->tmp.l], sys->rigid.R[v2->tmp.l]);
+	Mat3Transp(Rsum);
+
+	VecSubf(rhs, sys->rigid.origco[v1->tmp.l], sys->rigid.origco[v2->tmp.l]);
+	Mat3MulVecfl(Rsum, rhs);
+	VecMulf(rhs, 0.5f);
+	VecMulf(rhs, w);
+
+	VecAddf(sys->rigid.rhs[v1->tmp.l], sys->rigid.rhs[v1->tmp.l], rhs);
+}
+
+static void rigid_add_edge_to_rhs(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
+{
+	rigid_add_half_edge_to_rhs(sys, v1, v2, w);
+	rigid_add_half_edge_to_rhs(sys, v2, v1, w);
+}
+
+void rigid_deform_iteration()
+{
+	LaplacianSystem *sys= RigidDeformSystem;
+	EditMesh *em;
+	EditVert *eve;
+	EditFace *efa;
+	int a, i;
+
+	if(!sys)
+		return;
+	
+	nlMakeCurrent(sys->context);
+	em= sys->rigid.mesh;
+
+	/* compute R */
+	memset(sys->rigid.R, 0, sizeof(float)*3*3*sys->totvert);
+	memset(sys->rigid.rhs, 0, sizeof(float)*3*sys->totvert);
+
+	for(a=0, efa=em->faces.first; efa; efa=efa->next, a++) {
+		rigid_add_edge_to_R(sys, efa->v1, efa->v2, sys->fweights[a][2]);
+		rigid_add_edge_to_R(sys, efa->v2, efa->v3, sys->fweights[a][0]);
+		rigid_add_edge_to_R(sys, efa->v3, efa->v1, sys->fweights[a][1]);
+
+		if(efa->v4) {
+			a++;
+			rigid_add_edge_to_R(sys, efa->v1, efa->v3, sys->fweights[a][2]);
+			rigid_add_edge_to_R(sys, efa->v3, efa->v4, sys->fweights[a][0]);
+			rigid_add_edge_to_R(sys, efa->v4, efa->v1, sys->fweights[a][1]);
+		}
+	}
+
+	for(a=0, eve=em->verts.first; eve; eve=eve->next, a++) {
+		rigid_orthogonalize_R(sys->rigid.R[a]);
+		eve->tmp.l= a;
+	}
+	
+	/* compute right hand sides for solving */
+	for(a=0, efa=em->faces.first; efa; efa=efa->next, a++) {
+		rigid_add_edge_to_rhs(sys, efa->v1, efa->v2, sys->fweights[a][2]);
+		rigid_add_edge_to_rhs(sys, efa->v2, efa->v3, sys->fweights[a][0]);
+		rigid_add_edge_to_rhs(sys, efa->v3, efa->v1, sys->fweights[a][1]);
+
+		if(efa->v4) {
+			a++;
+			rigid_add_edge_to_rhs(sys, efa->v1, efa->v3, sys->fweights[a][2]);
+			rigid_add_edge_to_rhs(sys, efa->v3, efa->v4, sys->fweights[a][0]);
+			rigid_add_edge_to_rhs(sys, efa->v4, efa->v1, sys->fweights[a][1]);
+		}
+	}
+
+	/* solve for positions, for X,Y and Z separately */
+	for(i=0; i<3; i++) {
+		laplacian_begin_solve(sys, i);
+
+		for(a=0; a<sys->totvert; a++)
+			if(!sys->vpinned[a]) {
+				/*if (i==0)
+					printf("rhs %f\n", sys->rigid.rhs[a][0]);*/
+				laplacian_add_right_hand_side(sys, a, sys->rigid.rhs[a][i]);
+			}
+
+		if(laplacian_system_solve(sys)) {
+			for(a=0, eve=em->verts.first; eve; eve=eve->next, a++)
+				eve->co[i]= laplacian_system_get_solution(a);
+		}
+		else {
+			if(!sys->rigid.thrownerror) {
+				error("RigidDeform: failed to find solution.");
+				sys->rigid.thrownerror= 1;
+			}
+			break;
+		}
+	}
+
+	/*printf("\n--------------------------------------------\n\n");*/
+}
+
+static void rigid_laplacian_create(LaplacianSystem *sys)
+{
+	EditMesh *em = sys->rigid.mesh;
+	EditVert *eve;
+	EditFace *efa;
+	int a;
+
+	/* add verts and faces to laplacian */
+	for(a=0, eve=em->verts.first; eve; eve=eve->next, a++) {
+		laplacian_add_vertex(sys, eve->co, eve->pinned);
+		eve->tmp.l= a;
+	}
+
+	for(efa=em->faces.first; efa; efa=efa->next) {
+		laplacian_add_triangle(sys,
+			efa->v1->tmp.l, efa->v2->tmp.l, efa->v3->tmp.l);
+		if(efa->v4)
+			laplacian_add_triangle(sys,
+				efa->v1->tmp.l, efa->v3->tmp.l, efa->v4->tmp.l);
+	}
+}
+
+void rigid_deform_begin(EditMesh *em)
+{
+	LaplacianSystem *sys;
+	EditVert *eve;
+	EditFace *efa;
+	int a, totvert, totface;
+
+	/* count vertices, triangles */
+	for(totvert=0, eve=em->verts.first; eve; eve=eve->next)
+		totvert++;
+
+	for(totface=0, efa=em->faces.first; efa; efa=efa->next) {
+		totface++;
+		if(efa->v4) totface++;
+	}
+
+	/* create laplacian */
+	sys = laplacian_system_construct_begin(totvert, totface);
+
+	sys->rigid.mesh= em;
+	sys->rigid.R = MEM_callocN(sizeof(float)*3*3*totvert, "RigidDeformR");
+	sys->rigid.rhs = MEM_callocN(sizeof(float)*3*totvert, "RigidDeformRHS");
+	sys->rigid.origco = MEM_callocN(sizeof(float)*3*totvert, "RigidDeformCo");
+
+	for(a=0, eve=em->verts.first; eve; eve=eve->next, a++)
+		VecCopyf(sys->rigid.origco[a], eve->co);
+
+	sys->areaweights= 0;
+	sys->storeweights= 1;
+
+	rigid_laplacian_create(sys);
+
+	laplacian_system_construct_end(sys);
+
+	RigidDeformSystem = sys;
+}
+
+void rigid_deform_end(int cancel)
+{
+	LaplacianSystem *sys = RigidDeformSystem;
+
+	if(sys) {
+		EditMesh *em = sys->rigid.mesh;
+		EditVert *eve;
+		int a;
+
+		if(cancel)
+			for(a=0, eve=em->verts.first; eve; eve=eve->next, a++)
+				if(!eve->pinned)
+					VecCopyf(eve->co, sys->rigid.origco[a]);
+
+		if(sys->rigid.R) MEM_freeN(sys->rigid.R);
+		if(sys->rigid.rhs) MEM_freeN(sys->rigid.rhs);
+		if(sys->rigid.origco) MEM_freeN(sys->rigid.origco);
+
+		/* free */
+		laplacian_system_delete(sys);
+	}
+
+	RigidDeformSystem = NULL;
+}
+#endif
+
diff --git a/source/blender/src/poseobject.c b/source/blender/src/poseobject.c
index 0b78601ab9b..dda173d0d62 100644
--- a/source/blender/src/poseobject.c
+++ b/source/blender/src/poseobject.c
@@ -51,7 +51,6 @@
 #include "BKE_constraint.h"
 #include "BKE_deform.h"
 #include "BKE_depsgraph.h"
-#include "BKE_DerivedMesh.h"
 #include "BKE_displist.h"
 #include "BKE_global.h"
 #include "BKE_modifier.h"
@@ -66,7 +65,6 @@
 #include "BIF_graphics.h"
 #include "BIF_interface.h"
 #include "BIF_poseobject.h"
-#include "BIF_meshtools.h"
 #include "BIF_space.h"
 #include "BIF_toolbox.h"
 #include "BIF_screen.h"
@@ -754,114 +752,29 @@ void paste_posebuf (int flip)
 
 /* ********************************************** */
 
-struct vgroup_map {
-	float head[3], tail[3];
-	Bone *bone;
-	bDeformGroup *dg, *dgflip;
-	Object *meshobj;
-};
-
-static void pose_adds_vgroups__mapFunc(void *userData, int index, float *co, float *no_f, short *no_s)
-{
-	struct vgroup_map *map= userData;
-	float vec[3], fac;
-	
-	VECCOPY(vec, co);
-	Mat4MulVecfl(map->meshobj->obmat, vec);
-		
-	/* get the distance-factor from the vertex to bone */
-	fac= distfactor_to_bone (vec, map->head, map->tail, map->bone->rad_head, map->bone->rad_tail, map->bone->dist);
-	
-	/* add to vgroup. this call also makes me->dverts */
-	if(fac!=0.0f) 
-		add_vert_to_defgroup (map->meshobj, map->dg, index, fac, WEIGHT_REPLACE);
-	else
-		remove_vert_defgroup (map->meshobj, map->dg, index);
-	
-	if(map->dgflip) {
-		int j= mesh_get_x_mirror_vert(map->meshobj, index);
-		if(j>=0) {
-			if(fac!=0.0f) 
-				add_vert_to_defgroup (map->meshobj, map->dgflip, j, fac, WEIGHT_REPLACE);
-			else
-				remove_vert_defgroup (map->meshobj, map->dgflip, j);
-		}
-	}
-}
-
 /* context weightpaint and deformer in posemode */
-void pose_adds_vgroups(Object *meshobj)
+void pose_adds_vgroups(Object *meshobj, int heatweights)
 {
 	extern VPaint Gwp;         /* from vpaint */
-	struct vgroup_map map;
-	DerivedMesh *dm;
 	Object *poseobj= modifiers_isDeformedByArmature(meshobj);
-	bArmature *arm= poseobj->data;
-	bPoseChannel *pchan;
-	Bone *bone;
-	bDeformGroup *dg, *curdef;
-	
-	if(poseobj==NULL || (poseobj->flag & OB_POSEMODE)==0) return;
-	
-	dm = mesh_get_derived_final(meshobj, CD_MASK_BAREMESH);
-	
-	map.meshobj= meshobj;
-	
-	for(pchan= poseobj->pose->chanbase.first; pchan; pchan= pchan->next) {
-		bone= pchan->bone;
-		if(arm->layer & pchan->bone->layer) {
-			if(bone->flag & (BONE_SELECTED)) {
-				
-				/* check if mesh has vgroups */
-				dg= get_named_vertexgroup(meshobj, bone->name);
-				if(dg==NULL)
-					dg= add_defgroup_name(meshobj, bone->name);
-				
-				/* flipped bone */
-				if(Gwp.flag & VP_MIRROR_X) {
-					char name[32];
-					
-					BLI_strncpy(name, dg->name, 32);
-					bone_flip_name(name, 0);		// 0 = don't strip off number extensions
-					
-					for (curdef = meshobj->defbase.first; curdef; curdef=curdef->next)
-						if (!strcmp(curdef->name, name))
-							break;
-					map.dgflip= curdef;
-				}
-				else map.dgflip= NULL;
-				
-				/* get the root of the bone in global coords */
-				VECCOPY(map.head, bone->arm_head);
-				Mat4MulVecfl(poseobj->obmat, map.head);
-				
-				/* get the tip of the bone in global coords */
-				VECCOPY(map.tail, bone->arm_tail);
-				Mat4MulVecfl(poseobj->obmat, map.tail);
-				
-				/* use the optimal vertices instead of mverts */
-				map.dg= dg;
-				map.bone= bone;
-				if(dm->foreachMappedVert) 
-					dm->foreachMappedVert(dm, pose_adds_vgroups__mapFunc, (void*) &map);
-				else {
-					Mesh *me= meshobj->data;
-					int i;
-					for(i=0; i<me->totvert; i++) 
-						pose_adds_vgroups__mapFunc(&map, i, (me->mvert+i)->co, NULL, NULL);
-				}
-				
-			}
-		}
+
+	if(poseobj==NULL || (poseobj->flag & OB_POSEMODE)==0) {
+		error("The active object must have a deforming armature in pose mode");
+		return;
 	}
-	
-	dm->release(dm);
+
+	add_verts_to_dgroups(meshobj, poseobj, heatweights, (Gwp.flag & VP_MIRROR_X));
+
+	if(heatweights)
+		BIF_undo_push("Apply Bone Heat Weights to Vertex Groups");
+	else
+		BIF_undo_push("Apply Bone Envelopes to Vertex Groups");
 
 	allqueue(REDRAWVIEW3D, 0);
 	allqueue(REDRAWBUTSEDIT, 0);
 	
-	DAG_object_flush_update(G.scene, meshobj, OB_RECALC_DATA);	// and all its relations
-
+	// and all its relations
+	DAG_object_flush_update(G.scene, meshobj, OB_RECALC_DATA);
 }
 
 /* ********************************************** */