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(!(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;
-    }
+			if (data != NULL) {
+				hgroup = (bDeformGroup ***) data;
+				**hgroup = defgroup;
+				++*hgroup;
+			}
+			return 1;
+		}
+	}
     return 0;
 }
 
-static void add_verts_to_closest_dgroup(Object *ob, Object *par)
+static void add_vgroups__mapFunc(void *userData, int index, float *co, float *no_f, short *no_s)
 {
-    /* 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.
+	/* 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]);
+
+		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
+			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;
+		}
 	}
 
-    /* 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);
+	/* create verts */
+    mesh = (Mesh*)ob->data;
+	verts = MEM_callocN(mesh->totvert*sizeof(*verts), "closestboneverts");
+
+	if (wpmode) {
+		/* if in weight paint mode, use final verts from derivedmesh */
+		DerivedMesh *dm = mesh_get_derived_final(ob, CD_MASK_BAREMESH);
+
+		if(dm->foreachMappedVert) {
+			dm->foreachMappedVert(dm, add_vgroups__mapFunc, (void*)verts);
+			vertsfilled = 1;
 		}
-		else {
-			vert = mesh->mvert + i;
-			VECCOPY (real_co, vert->co);
-		}
-        Mat4MulVecfl(ob->obmat, real_co);
 
+		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;
+	}
 
-        /* for each skinnable bone ...
-         */
-        for (j=0; j < numbones; ++j) {
-            bone = bonelist[j];
+	/* 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]);
+	}
 
-            /* get the root of the bone in global coords
-             */
-			VECCOPY(root, bone->arm_head);
-			Mat4MulVecfl(par->obmat, root);
+	/* 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);
 
-            /* 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);
-        }
-    }
-
-    /* free the memory allocated
-     */
+    /* 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);
 }
 
 /* ********************************************** */