Quaternion Deform Interpolation

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

This is a new armature deform interpolation method using Dual Quaternions,
which reduces the artifacts of linear blend skinning:

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

Based on the paper and provided code:

Skinning with Dual Quaternions
Ladislav Kavan, Steven Collins, Jiri Zara, Carol O'Sullivan.
Symposium on Interactive 3D Graphics and Games, 2007.

6 files changed, 371 insertions, 70 deletions

diff --git a/source/blender/blenkernel/intern/armature.c b/source/blender/blenkernel/intern/armature.c
index db384181982..f75174e47a2 100644
--- a/source/blender/blenkernel/intern/armature.c
+++ b/source/blender/blenkernel/intern/armature.c
@@ -499,11 +499,12 @@ Mat4 *b_bone_spline_setup(bPoseChannel *pchan)
 
 /* ************ Armature Deform ******************* */
 
-static void pchan_b_bone_defmats(bPoseChannel *pchan)
+static void pchan_b_bone_defmats(bPoseChannel *pchan, int use_quaternion)
 {
 	Bone *bone= pchan->bone;
 	Mat4 *b_bone= b_bone_spline_setup(pchan);
 	Mat4 *b_bone_mats;
+	DualQuat *b_bone_dual_quats= NULL;
 	float tmat[4][4];
 	int a;
 	
@@ -511,6 +512,11 @@ static void pchan_b_bone_defmats(bPoseChannel *pchan)
 	b_bone_mats= MEM_mallocN((1+bone->segments)*sizeof(Mat4), "BBone defmats");
 	pchan->b_bone_mats= b_bone_mats;
 
+	if(use_quaternion) {
+		b_bone_dual_quats= MEM_mallocN((bone->segments)*sizeof(DualQuat), "BBone dqs");
+		pchan->b_bone_dual_quats= b_bone_dual_quats;
+	}
+	
 	/* first matrix is the inverse arm_mat, to bring points in local bone space
 	   for finding out which segment it belongs to */
 	Mat4Invert(b_bone_mats[0].mat, bone->arm_mat);
@@ -527,10 +533,13 @@ static void pchan_b_bone_defmats(bPoseChannel *pchan)
 
 		Mat4MulSerie(b_bone_mats[a+1].mat, pchan->chan_mat, bone->arm_mat,
 			b_bone[a].mat, tmat, b_bone_mats[0].mat, NULL, NULL, NULL);
+
+		if(use_quaternion)
+			Mat4ToDQuat(bone->arm_mat, b_bone_mats[a+1].mat, &b_bone_dual_quats[a]);
 	}
 }
 
-static void b_bone_deform(bPoseChannel *pchan, Bone *bone, float *co, float defmat[][3])
+static void b_bone_deform(bPoseChannel *pchan, Bone *bone, float *co, DualQuat *dq, float defmat[][3])
 {
 	Mat4 *b_bone= pchan->b_bone_mats;
 	float (*mat)[4]= b_bone[0].mat;
@@ -548,10 +557,15 @@ static void b_bone_deform(bPoseChannel *pchan, Bone *bone, float *co, float defm
 	   straight joints in restpos. */
 	CLAMP(a, 0, bone->segments-1);
 
-	Mat4MulVecfl(b_bone[a+1].mat, co);
+	if(dq) {
+		DQuatCpyDQuat(dq, &((DualQuat*)pchan->b_bone_dual_quats)[a]);
+	}
+	else {
+		Mat4MulVecfl(b_bone[a+1].mat, co);
 
-	if(defmat)
-		Mat3CpyMat4(defmat, b_bone[a+1].mat);
+		if(defmat)
+			Mat3CpyMat4(defmat, b_bone[a+1].mat);
+	}
 }
 
 /* using vec with dist to bone b1 - b2 */
@@ -618,11 +632,12 @@ static void pchan_deform_mat_add(bPoseChannel *pchan, float weight, float bbonem
 	Mat3AddMat3(mat, mat, wmat);
 }
 
-static float dist_bone_deform(bPoseChannel *pchan, float *vec, float mat[][3], float *co)
+static float dist_bone_deform(bPoseChannel *pchan, float *vec, DualQuat *dq, float mat[][3], float *co)
 {
 	Bone *bone= pchan->bone;
 	float fac, contrib=0.0;
 	float cop[3], bbonemat[3][3];
+	DualQuat bbonedq;
 
 	if(bone==NULL) return 0.0f;
 	
@@ -635,46 +650,67 @@ static float dist_bone_deform(bPoseChannel *pchan, float *vec, float mat[][3], f
 		fac*=bone->weight;
 		contrib= fac;
 		if(contrib>0.0) {
-			if(bone->segments>1)
-				// applies on cop and bbonemat
-				b_bone_deform(pchan, bone, cop, (mat)?bbonemat:NULL);
-			else
-				Mat4MulVecfl(pchan->chan_mat, cop);
+			if(vec) {
+				if(bone->segments>1)
+					// applies on cop and bbonemat
+					b_bone_deform(pchan, bone, cop, NULL, (mat)?bbonemat:NULL);
+				else
+					Mat4MulVecfl(pchan->chan_mat, cop);
 
-			//	Make this a delta from the base position
-			VecSubf (cop, cop, co);
-			cop[0]*=fac; cop[1]*=fac; cop[2]*=fac;
-			VecAddf (vec, vec, cop);
+				//	Make this a delta from the base position
+				VecSubf (cop, cop, co);
+				cop[0]*=fac; cop[1]*=fac; cop[2]*=fac;
+				VecAddf (vec, vec, cop);
 
-			if(mat)
-				pchan_deform_mat_add(pchan, fac, bbonemat, mat);
+				if(mat)
+					pchan_deform_mat_add(pchan, fac, bbonemat, mat);
+			}
+			else {
+				if(bone->segments>1) {
+					b_bone_deform(pchan, bone, cop, &bbonedq, NULL);
+					DQuatAddWeighted(dq, &bbonedq, fac);
+				}
+				else
+					DQuatAddWeighted(dq, pchan->dual_quat, fac);
+			}
 		}
 	}
 	
 	return contrib;
 }
 
-static void pchan_bone_deform(bPoseChannel *pchan, float weight, float *vec, float mat[][3], float *co, float *contrib)
+static void pchan_bone_deform(bPoseChannel *pchan, float weight, float *vec, DualQuat *dq, float mat[][3], float *co, float *contrib)
 {
 	float cop[3], bbonemat[3][3];
+	DualQuat bbonedq;
 
 	if (!weight)
 		return;
 
 	VECCOPY(cop, co);
 
-	if(pchan->bone->segments>1)
-		// applies on cop and bbonemat
-		b_bone_deform(pchan, pchan->bone, cop, (mat)?bbonemat:NULL);
-	else
-		Mat4MulVecfl(pchan->chan_mat, cop);
-	
-	vec[0]+=(cop[0]-co[0])*weight;
-	vec[1]+=(cop[1]-co[1])*weight;
-	vec[2]+=(cop[2]-co[2])*weight;
+	if(vec) {
+		if(pchan->bone->segments>1)
+			// applies on cop and bbonemat
+			b_bone_deform(pchan, pchan->bone, cop, NULL, (mat)?bbonemat:NULL);
+		else
+			Mat4MulVecfl(pchan->chan_mat, cop);
+		
+		vec[0]+=(cop[0]-co[0])*weight;
+		vec[1]+=(cop[1]-co[1])*weight;
+		vec[2]+=(cop[2]-co[2])*weight;
 
-	if(mat)
-		pchan_deform_mat_add(pchan, weight, bbonemat, mat);
+		if(mat)
+			pchan_deform_mat_add(pchan, weight, bbonemat, mat);
+	}
+	else {
+		if(pchan->bone->segments>1) {
+			b_bone_deform(pchan, pchan->bone, cop, &bbonedq, NULL);
+			DQuatAddWeighted(dq, &bbonedq, weight);
+		}
+		else
+			DQuatAddWeighted(dq, pchan->dual_quat, weight);
+	}
 
 	(*contrib)+=weight;
 }
@@ -686,12 +722,15 @@ void armature_deform_verts(Object *armOb, Object *target, DerivedMesh *dm,
 	bPoseChannel *pchan, **defnrToPC = NULL;
 	MDeformVert *dverts = NULL;
 	bDeformGroup *dg;
+	DualQuat *dualquats= NULL;
 	float obinv[4][4], premat[4][4], postmat[4][4];
 	int use_envelope = deformflag & ARM_DEF_ENVELOPE;
+	int use_quaternion = deformflag & ARM_DEF_QUATERNION;
 	int numGroups = 0;		/* safety for vertexgroup index overflow */
 	int i, target_totvert = 0;	/* safety for vertexgroup overflow */
 	int use_dverts = 0;
 	int armature_def_nr = -1;
+	int totchan;
 
 	if(armOb == G.obedit) return;
 	
@@ -703,10 +742,23 @@ void armature_deform_verts(Object *armOb, Object *target, DerivedMesh *dm,
 	/* bone defmats are already in the channels, chan_mat */
 	
 	/* initialize B_bone matrices and dual quaternions */
-	for(pchan = armOb->pose->chanbase.first; pchan; pchan = pchan->next)
-		if(!(pchan->bone->flag & BONE_NO_DEFORM))
+	if(use_quaternion) {
+		totchan= BLI_countlist(&armOb->pose->chanbase);
+		dualquats= MEM_callocN(sizeof(DualQuat)*totchan, "dualquats");
+	}
+
+	totchan= 0;
+	for(pchan = armOb->pose->chanbase.first; pchan; pchan = pchan->next) {
+		if(!(pchan->bone->flag & BONE_NO_DEFORM)) {
 			if(pchan->bone->segments > 1)
-				pchan_b_bone_defmats(pchan);
+				pchan_b_bone_defmats(pchan, use_quaternion);
+
+			if(use_quaternion) {
+				pchan->dual_quat= &dualquats[totchan++];
+				Mat4ToDQuat(pchan->bone->arm_mat, pchan->chan_mat, pchan->dual_quat);
+			}
+		}
+	}
 
 	/* get the def_nr for the overall armature vertex group if present */
 	for(i = 0, dg = target->defbase.first; dg; i++, dg = dg->next)
@@ -754,18 +806,26 @@ void armature_deform_verts(Object *armOb, Object *target, DerivedMesh *dm,
 
 	for(i = 0; i < numVerts; i++) {
 		MDeformVert *dvert;
+		DualQuat sumdq, *dq = NULL;
 		float *co = vertexCos[i];
-		float summat[3][3], (*smat)[3] = NULL;
-		float vec[3];
+		float sumvec[3], summat[3][3];
+		float *vec = NULL, (*smat)[3] = NULL;
 		float contrib = 0.0f;
 		float armature_weight = 1.0f; /* default to 1 if no overall def group */
 		int	  j;
 
-		vec[0] = vec[1] = vec[2] = 0.0f;
+		if(use_quaternion) {
+			memset(&sumdq, 0, sizeof(DualQuat));
+			dq= &sumdq;
+		}
+		else {
+			sumvec[0] = sumvec[1] = sumvec[2] = 0.0f;
+			vec= sumvec;
 
-		if(defMats) {
-			Mat3Clr((float*)summat);
-			smat = summat;
+			if(defMats) {
+				Mat3Clr((float*)summat);
+				smat = summat;
+			}
 		}
 
 		if(use_dverts || armature_def_nr >= 0) {
@@ -810,7 +870,7 @@ void armature_deform_verts(Object *armOb, Object *target, DerivedMesh *dm,
 						                             bone->rad_tail,
 						                             bone->dist);
 					}
-					pchan_bone_deform(pchan, weight, vec, smat, co, &contrib);
+					pchan_bone_deform(pchan, weight, vec, dq, smat, co, &contrib);
 				}
 			}
 			/* if there are vertexgroups but not groups with bones
@@ -820,7 +880,7 @@ void armature_deform_verts(Object *armOb, Object *target, DerivedMesh *dm,
 				for(pchan = armOb->pose->chanbase.first; pchan;
 				    pchan = pchan->next) {
 					if(!(pchan->bone->flag & BONE_NO_DEFORM))
-						contrib += dist_bone_deform(pchan, vec, smat, co);
+						contrib += dist_bone_deform(pchan, vec, dq, smat, co);
 				}
 			}
 		}
@@ -828,16 +888,21 @@ void armature_deform_verts(Object *armOb, Object *target, DerivedMesh *dm,
 			for(pchan = armOb->pose->chanbase.first; pchan;
 			    pchan = pchan->next) {
 				if(!(pchan->bone->flag & BONE_NO_DEFORM))
-					contrib += dist_bone_deform(pchan, vec, smat, co);
+					contrib += dist_bone_deform(pchan, vec, dq, smat, co);
 			}
 		}
 
 		/* actually should be EPSILON? weight values and contrib can be like 10e-39 small */
 		if(contrib > 0.0001f) {
-			float scale = armature_weight/contrib;
-
-			VecMulf(vec, scale);
-			VecAddf(co, vec, co);
+			if(use_quaternion) {
+				DQuatNormalize(dq, contrib, armature_weight);
+				DQuatMulVecfl(dq, co, (defMats)? summat: NULL);
+				smat = summat;
+			}
+			else {
+				VecMulf(vec, armature_weight/contrib);
+				VecAddf(co, vec, co);
+			}
 
 			if(defMats) {
 				float pre[3][3], post[3][3], tmpmat[3][3];
@@ -846,17 +911,19 @@ void armature_deform_verts(Object *armOb, Object *target, DerivedMesh *dm,
 				Mat3CpyMat4(post, postmat);
 				Mat3CpyMat3(tmpmat, defMats[i]);
 
-				Mat3MulFloat((float*)smat, scale);
+				if(!use_quaternion) /* quaternion already is scale corrected */
+					Mat3MulFloat((float*)smat, armature_weight/contrib);
+
 				Mat3MulSerie(defMats[i], tmpmat, pre, smat, post,
 					NULL, NULL, NULL, NULL);
 			}
-
 		}
 		
 		/* always, check above code */
 		Mat4MulVecfl(postmat, co);
 	}
 
+	if(dualquats) MEM_freeN(dualquats);
 	if(defnrToPC) MEM_freeN(defnrToPC);
 	
 	/* free B_bone matrices */
@@ -865,6 +932,12 @@ void armature_deform_verts(Object *armOb, Object *target, DerivedMesh *dm,
 			MEM_freeN(pchan->b_bone_mats);
 			pchan->b_bone_mats = NULL;
 		}
+		if(pchan->b_bone_dual_quats) {
+			MEM_freeN(pchan->b_bone_dual_quats);
+			pchan->b_bone_dual_quats = NULL;
+		}
+
+		pchan->dual_quat = NULL;
 	}
 }
 
diff --git a/source/blender/blenlib/BLI_arithb.h b/source/blender/blenlib/BLI_arithb.h
index 7a1bedf5c14..1adc840dfe4 100644
--- a/source/blender/blenlib/BLI_arithb.h
+++ b/source/blender/blenlib/BLI_arithb.h
@@ -128,6 +128,7 @@ void QuatConj(float *q);
 void QuatInv(float *q);
 void QuatMulf(float *q, float f);
 float QuatDot(float *q1, float *q2);
+void QuatCopy(float *q1, float *q2);
 
 void printquat(char *str, float q[4]);
 
@@ -353,6 +354,21 @@ void spheremap(float x, float y, float z, float *u, float *v);
 int LineIntersectsTriangle(float p1[3], float p2[3], float v0[3], float v1[3], float v2[3], float *lambda);
 int point_in_tri_prism(float p[3], float v1[3], float v2[3], float v3[3]);
 
+typedef struct DualQuat {
+	float quat[4];
+	float trans[4];
+
+	float scale[4][4];
+	float scale_weight;
+} DualQuat;
+
+void Mat4ToDQuat(float basemat[][4], float mat[][4], DualQuat *dq);
+void DQuatToMat4(DualQuat *dq, float mat[][4]);
+void DQuatAddWeighted(DualQuat *dqsum, DualQuat *dq, float weight);
+void DQuatNormalize(DualQuat *dq, float totweight, float factor);
+void DQuatMulVecfl(DualQuat *dq, float *co, float mat[][3]);
+void DQuatCpyDQuat(DualQuat *dq1, DualQuat *dq2);
+			  
 #ifdef __cplusplus
 }
 #endif
diff --git a/source/blender/blenlib/intern/arithb.c b/source/blender/blenlib/intern/arithb.c
index 387e2227797..237b446bf2c 100644
--- a/source/blender/blenlib/intern/arithb.c
+++ b/source/blender/blenlib/intern/arithb.c
@@ -936,7 +936,7 @@ void Mat4MulFloat(float *m, float f)
 {
 	int i;
 
-	for(i=0;i<12;i++) m[i]*=f;	/* count to 12: without vector component */
+	for(i=0;i<16;i++) m[i]*=f;	/* count to 12: without vector component */
 }
 
 
@@ -1597,6 +1597,221 @@ void QuatAdd(float *result, float *quat1, float *quat2, float t)
 	result[3]= quat1[3] + t*quat2[3];
 }
 
+void QuatCopy(float *q1, float *q2)
+{
+	q1[0]= q2[0];
+	q1[1]= q2[1];
+	q1[2]= q2[2];
+	q1[3]= q2[3];
+}
+
+/* **************** DUAL QUATERNIONS ************** */
+
+/*
+   Conversion routines between (regular quaternion, translation) and
+   dual quaternion.
+
+   Version 1.0.0, February 7th, 2007
+
+   Copyright (C) 2006-2007 University of Dublin, Trinity College, All Rights 
+   Reserved
+
+   This software is provided 'as-is', without any express or implied
+   warranty.  In no event will the author(s) be held liable for any damages
+   arising from the use of this software.
+
+   Permission is granted to anyone to use this software for any purpose,
+   including commercial applications, and to alter it and redistribute it
+   freely, subject to the following restrictions:
+
+   1. The origin of this software must not be misrepresented; you must not
+      claim that you wrote the original software. If you use this software
+      in a product, an acknowledgment in the product documentation would be
+      appreciated but is not required.
+   2. Altered source versions must be plainly marked as such, and must not be
+      misrepresented as being the original software.
+   3. This notice may not be removed or altered from any source distribution.
+
+   Author: Ladislav Kavan, kavanl@cs.tcd.ie
+
+   Changes for Blender:
+   - renaming, style changes and optimizations
+   - added support for scaling
+*/
+
+void Mat4ToDQuat(float basemat[][4], float mat[][4], DualQuat *dq)
+{
+	float *t, *q, dscale[3], scale[3], basequat[4];
+	float baseRS[4][4], baseinv[4][4], baseR[4][4], baseRinv[4][4];
+	float R[4][4], S[4][4];
+
+	/* split scaling and rotation, there is probably a faster way to do
+	   this, it's done like this now to correctly get negative scaling */
+	Mat4MulMat4(baseRS, basemat, mat);
+	Mat4ToSize(baseRS, scale);
+
+	VecCopyf(dscale, scale);
+	dscale[0] -= 1.0f; dscale[1] -= 1.0f; dscale[2] -= 1.0f;
+
+	if((Det4x4(mat) < 0.0f) || VecLength(dscale) > 1e-4) {
+		/* extract R and S  */
+		Mat4ToQuat(baseRS, basequat);
+		QuatToMat4(basequat, baseR);
+		VecCopyf(baseR[3], baseRS[3]);
+
+		Mat4Invert(baseinv, basemat);
+		Mat4MulMat4(R, baseinv, baseR);
+
+		Mat4Invert(baseRinv, baseR);
+		Mat4MulMat4(S, baseRS, baseRinv);
+
+		/* set scaling part */
+		Mat4MulSerie(dq->scale, basemat, S, baseinv, 0, 0, 0, 0, 0);
+		dq->scale_weight= 1.0f;
+	}
+	else {
+		/* matrix does not contain scaling */
+		Mat4CpyMat4(R, mat);
+		dq->scale_weight= 0.0f;
+	}
+
+	/* non-dual part */
+	Mat4ToQuat(R, dq->quat);
+
+	/* dual part */
+	t= R[3];
+	q= dq->quat;
+	dq->trans[0]= -0.5f*( t[0]*q[1] + t[1]*q[2] + t[2]*q[3]);
+	dq->trans[1]=  0.5f*( t[0]*q[0] + t[1]*q[3] - t[2]*q[2]);
+	dq->trans[2]=  0.5f*(-t[0]*q[3] + t[1]*q[0] + t[2]*q[1]);
+	dq->trans[3]=  0.5f*( t[0]*q[2] - t[1]*q[1] + t[2]*q[0]);
+}
+
+void DQuatToMat4(DualQuat *dq, float mat[][4])
+{
+	float len, *t, q0[4];
+	
+	/* regular quaternion */
+	QuatCopy(q0, dq->quat);
+
+	/* normalize */
+	len= sqrt(QuatDot(q0, q0)); 
+	if(len != 0.0f)
+		QuatMulf(q0, 1.0f/len);
+	
+	/* rotation */
+	QuatToMat4(q0, mat);
+
+	/* translation */
+	t= dq->trans;
+	mat[3][0]= 2.0*(-t[0]*q0[1] + t[1]*q0[0] - t[2]*q0[3] + t[3]*q0[2]);
+	mat[3][1]= 2.0*(-t[0]*q0[2] + t[1]*q0[3] + t[2]*q0[0] - t[3]*q0[1]);
+	mat[3][2]= 2.0*(-t[0]*q0[3] - t[1]*q0[2] + t[2]*q0[1] + t[3]*q0[0]);
+
+	/* note: this does not handle scaling */
+}	
+
+void DQuatAddWeighted(DualQuat *dqsum, DualQuat *dq, float weight)
+{
+	/* make sure we interpolate quats in the right direction */
+	if (QuatDot(dq->quat, dqsum->quat) < 0)
+		weight = -weight;
+
+	/* interpolate rotation and translation */
+	dqsum->quat[0] += weight*dq->quat[0];
+	dqsum->quat[1] += weight*dq->quat[1];
+	dqsum->quat[2] += weight*dq->quat[2];
+	dqsum->quat[3] += weight*dq->quat[3];
+
+	dqsum->trans[0] += weight*dq->trans[0];
+	dqsum->trans[1] += weight*dq->trans[1];
+	dqsum->trans[2] += weight*dq->trans[2];
+	dqsum->trans[3] += weight*dq->trans[3];
+
+	/* interpolate scale - but only if needed */
+	if (dq->scale_weight) {
+		float wmat[4][4];
+
+		Mat4CpyMat4(wmat, dq->scale);
+		Mat4MulFloat((float*)wmat, weight);
+		Mat4AddMat4(dqsum->scale, dqsum->scale, wmat);
+		dqsum->scale_weight += weight;
+	}
+}
+
+void DQuatNormalize(DualQuat *dq, float totweight, float factor)
+{
+	float scale= factor/totweight;
+
+	QuatMulf(dq->quat, scale);
+	QuatMulf(dq->trans, scale);
+	
+	if(dq->scale_weight) {
+		float addweight= totweight - dq->scale_weight;
+
+		if(addweight) {
+			dq->scale[0][0] += addweight;
+			dq->scale[1][1] += addweight;
+			dq->scale[2][2] += addweight;
+			dq->scale[3][3] += addweight;
+		}
+
+		Mat4MulFloat((float*)dq->scale, scale);
+		dq->scale_weight= 1.0f;
+	}
+}
+
+void DQuatMulVecfl(DualQuat *dq, float *co, float mat[][3])
+{	
+	float M[3][3], t[3], scalemat[3][3], len2;
+	float w= dq->quat[0], x= dq->quat[1], y= dq->quat[2], z= dq->quat[3];
+	float t0= dq->trans[0], t1= dq->trans[1], t2= dq->trans[2], t3= dq->trans[3];
+	
+	/* rotation matrix */
+	M[0][0]= w*w + x*x - y*y - z*z;
+	M[1][0]= 2*(x*y - w*z);
+	M[2][0]= 2*(x*z + w*y);
+
+	M[0][1]= 2*(x*y + w*z);
+	M[1][1]= w*w + y*y - x*x - z*z;
+	M[2][1]= 2*(y*z - w*x); 
+
+	M[0][2]= 2*(x*z - w*y);
+	M[1][2]= 2*(y*z + w*x);
+	M[2][2]= w*w + z*z - x*x - y*y;
+	
+	len2= QuatDot(dq->quat, dq->quat);
+	if(len2 > 0.0f)
+		len2= 1.0f/len2;
+	
+	/* translation */
+	t[0]= 2*(-t0*x + w*t1 - t2*z + y*t3);
+	t[1]= 2*(-t0*y + t1*z - x*t3 + w*t2);
+	t[2]= 2*(-t0*z + x*t2 + w*t3 - t1*y);
+
+	/* apply scaling */
+	if(dq->scale_weight)
+		Mat4MulVecfl(dq->scale, co);
+	
+	/* apply rotation and translation */
+	Mat3MulVecfl(M, co);
+	co[0]= (co[0] + t[0])*len2;
+	co[1]= (co[1] + t[1])*len2;
+	co[2]= (co[2] + t[2])*len2;
+
+	/* compute crazyspace correction mat */
+	if(mat) {
+		Mat3CpyMat4(scalemat, dq->scale);
+		Mat3MulMat3(mat, M, scalemat);
+		Mat3MulFloat((float*)mat, len2);
+	}
+}
+
+void DQuatCpyDQuat(DualQuat *dq1, DualQuat *dq2)
+{
+	memcpy(dq1, dq2, sizeof(DualQuat));
+}
+
 /* **************** VIEW / PROJECTION ********************************  */
 
 
@@ -2635,28 +2850,16 @@ void SizeToMat4( float *size, float mat[][4])
 
 void Mat3ToSize( float mat[][3], float *size)
 {
-	float vec[3];
-
-	VecCopyf(vec, mat[0]);
-	size[0]= Normalize(vec);
-	VecCopyf(vec, mat[1]);
-	size[1]= Normalize(vec);
-	VecCopyf(vec, mat[2]);
-	size[2]= Normalize(vec);
-
+	size[0]= VecLength(mat[0]);
+	size[1]= VecLength(mat[1]);
+	size[2]= VecLength(mat[2]);
 }
 
 void Mat4ToSize( float mat[][4], float *size)
 {
-	float vec[3];
-	
-
-	VecCopyf(vec, mat[0]);
-	size[0]= Normalize(vec);
-	VecCopyf(vec, mat[1]);
-	size[1]= Normalize(vec);
-	VecCopyf(vec, mat[2]);
-	size[2]= Normalize(vec);
+	size[0]= VecLength(mat[0]);
+	size[1]= VecLength(mat[1]);
+	size[2]= VecLength(mat[2]);
 }
 
 /* ************* SPECIALS ******************* */
diff --git a/source/blender/makesdna/DNA_action_types.h b/source/blender/makesdna/DNA_action_types.h
index f0640dc5178..c12a4bf4022 100644
--- a/source/blender/makesdna/DNA_action_types.h
+++ b/source/blender/makesdna/DNA_action_types.h
@@ -61,7 +61,12 @@ typedef struct bPoseChannel {
 	struct bPoseChannel *parent;	/* set on read file or rebuild pose */
 	struct bPoseChannel *child;		/* set on read file or rebuild pose, the 'ik' child, for b-bones */
 	struct ListBase		 iktree;		/* only while evaluating pose */
-	void				*b_bone_mats;	/* only while deform, stores precalculated b_bone deform mats */
+	
+	/* only while deform, stores precalculated b_bone deform mats,
+	   dual quaternions */
+	void				*b_bone_mats;	
+	void				*dual_quat;
+	void				*b_bone_dual_quats;
 	
 	float		loc[3];				/* written in by actions or transform */
 	float		size[3];
diff --git a/source/blender/makesdna/DNA_armature_types.h b/source/blender/makesdna/DNA_armature_types.h
index b3235fc71f3..53a5f26033c 100644
--- a/source/blender/makesdna/DNA_armature_types.h
+++ b/source/blender/makesdna/DNA_armature_types.h
@@ -113,6 +113,7 @@ typedef struct bArmature {
 /* armature->deformflag */
 #define		ARM_DEF_VGROUP		1
 #define		ARM_DEF_ENVELOPE	2
+#define		ARM_DEF_QUATERNION	4
 
 /* armature->pathflag */
 #define		ARM_PATH_FNUMS	0x001
diff --git a/source/blender/src/buttons_editing.c b/source/blender/src/buttons_editing.c
index 05d08e4c9fc..65b299b9f96 100644
--- a/source/blender/src/buttons_editing.c
+++ b/source/blender/src/buttons_editing.c
@@ -1592,7 +1592,7 @@ static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco
 			if(wmd->flag & MOD_WAVE_NORM)
 				height += 19;
 		} else if (md->type==eModifierType_Armature) {
-			height = 67;
+			height = 86;
 		} else if (md->type==eModifierType_Hook) {
 			HookModifierData *hmd = (HookModifierData*) md;
 			height = 86;
@@ -1913,6 +1913,8 @@ static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco
 			uiButSetCompleteFunc(but, autocomplete_vgroup, (void *)ob);
 			uiDefButBitS(block, TOG, ARM_DEF_VGROUP, B_ARM_RECALCDATA, "Vert.Groups",	lx,cy-=19,buttonWidth/2,20, &amd->deformflag, 0, 0, 0, 0, "Enable VertexGroups defining deform");
 			uiDefButBitS(block, TOG, ARM_DEF_ENVELOPE, B_ARM_RECALCDATA, "Envelopes",	lx+buttonWidth/2,cy,(buttonWidth + 1)/2,20, &amd->deformflag, 0, 0, 0, 0, "Enable Bone Envelopes defining deform");
+			uiDefButBitS(block, TOG, ARM_DEF_QUATERNION, B_ARM_RECALCDATA, "Quaternion",	lx,(cy-=19),buttonWidth + 1,20, &amd->deformflag, 0, 0, 0, 0, "Enable deform rotation interpolation with Quaternions");
+
 			
 		} else if (md->type==eModifierType_Hook) {
 			HookModifierData *hmd = (HookModifierData*) md;
@@ -3624,8 +3626,9 @@ static void editing_panel_armature_type(Object *ob, bArmature *arm)
 	
 	uiDefBut(block, LABEL, 0, "Deform Options", 10,40,150,20, 0, 0, 0, 0, 0, "");
 	uiBlockBeginAlign(block);
-	uiDefButBitS(block, TOG, ARM_DEF_VGROUP, B_ARM_RECALCDATA, "Vertex Groups",	10, 20,150,20, &arm->deformflag, 0, 0, 0, 0, "Enable VertexGroups defining deform (not for Modifiers)");
-	uiDefButBitS(block, TOG, ARM_DEF_ENVELOPE, B_ARM_RECALCDATA, "Envelopes",	160,20,150,20, &arm->deformflag, 0, 0, 0, 0, "Enable Bone Envelopes defining deform (not for Modifiers)");
+	uiDefButBitS(block, TOG, ARM_DEF_VGROUP, B_ARM_RECALCDATA, "Vertex Groups",	10, 20,100,20, &arm->deformflag, 0, 0, 0, 0, "Enable VertexGroups defining deform (not for Modifiers)");
+	uiDefButBitS(block, TOG, ARM_DEF_ENVELOPE, B_ARM_RECALCDATA, "Envelopes",	110,20,100,20, &arm->deformflag, 0, 0, 0, 0, "Enable Bone Envelopes defining deform (not for Modifiers)");
+	uiDefButBitS(block, TOG, ARM_DEF_QUATERNION, B_ARM_RECALCDATA, "Quaternion", 210,20,100,20, &arm->deformflag, 0, 0, 0, 0, "Enable deform rotation interpolation with Quaternions (not for Modifiers)");
 	uiDefButBitI(block, TOG, ARM_RESTPOS, B_ARM_RECALCDATA,"Rest Position",		10,0,150,20, &arm->flag, 0, 0, 0, 0, "Show armature rest position, no posing possible");
 	uiDefButBitI(block, TOG, ARM_DELAYDEFORM, REDRAWVIEW3D, "Delay Deform",		160,0,150,20, &arm->flag, 0, 0, 0, 0, "Don't deform children when manipulating bones in pose mode");
 	uiBlockEndAlign(block);