curve twist

* added new twist method - "Tangent", suggested by Martin.
  the nice thing about this is its stable no matter how you rotate the data, rotation is local to each segment.
* added smooth option that smooths the twisting (before applying user twist), to workaround Z-Up and Tangent's ugly curve twisting. Id prefer not to have this however it makes tangent much nicer. Possibly tangent can be improved some other way and this can be removed.
  A smooth value of 1.0 will iterate over and smooth the twisting by the resolution value of the spline.
* Minimum-Twist method now corrects for cyclic twist by taking the roll difference between first and last, then increasingly counter rotate each segment over the entire curve. Previously it calculated from both directions and blended them.

details
* BevPoints use quats rather then 3x3 matrix.
* added BevPoint direction "dir" and tangent "tan" used only for 3D curves.
* don't calculate BevPoint->cosa, BevPoint->sina for 3D curves.
* split bevel tilt calculation into functions.
* nurbs curves currently don't generate tangents and wont work with tangent twist method.
* some of the use of quats should be optimized.
* smoothing is not animation safe, the higher the smoothing the higher the likelyhood of flipping.

6 files changed, 474 insertions, 198 deletions

diff --git a/source/blender/blenkernel/BKE_utildefines.h b/source/blender/blenkernel/BKE_utildefines.h
index 76e0da98f69..4d43518901e 100644
--- a/source/blender/blenkernel/BKE_utildefines.h
+++ b/source/blender/blenkernel/BKE_utildefines.h
@@ -107,6 +107,7 @@
 #define VECSUB(v1,v2,v3) 	{*(v1)= *(v2) - *(v3); *(v1+1)= *(v2+1) - *(v3+1); *(v1+2)= *(v2+2) - *(v3+2);}
 #define VECSUB2D(v1,v2,v3) 	{*(v1)= *(v2) - *(v3); *(v1+1)= *(v2+1) - *(v3+1);}
 #define VECADDFAC(v1,v2,v3,fac) {*(v1)= *(v2) + *(v3)*(fac); *(v1+1)= *(v2+1) + *(v3+1)*(fac); *(v1+2)= *(v2+2) + *(v3+2)*(fac);}
+#define VECSUBFAC(v1,v2,v3,fac) {*(v1)= *(v2) - *(v3)*(fac); *(v1+1)= *(v2+1) - *(v3+1)*(fac); *(v1+2)= *(v2+2) - *(v3+2)*(fac);}
 #define QUATADDFAC(v1,v2,v3,fac) {*(v1)= *(v2) + *(v3)*(fac); *(v1+1)= *(v2+1) + *(v3+1)*(fac); *(v1+2)= *(v2+2) + *(v3+2)*(fac); *(v1+3)= *(v2+3) + *(v3+3)*(fac);}
 
 #define INPR(v1, v2)		( (v1)[0]*(v2)[0] + (v1)[1]*(v2)[1] + (v1)[2]*(v2)[2] )
diff --git a/source/blender/blenkernel/intern/curve.c b/source/blender/blenkernel/intern/curve.c
index 802809b11e5..2b8c91fcbed 100644
--- a/source/blender/blenkernel/intern/curve.c
+++ b/source/blender/blenkernel/intern/curve.c
@@ -975,7 +975,26 @@ void forward_diff_bezier(float q0, float q1, float q2, float q3, float *p, int i
  		q1+= q2;
  		q2+= q3;
  	}
-}	
+}
+
+void forward_diff_bezier_cotangent(float *p0, float *p1, float *p2, float *p3, float *p, int it, int stride)
+{
+	/* note that these are not purpendicular to the curve
+	 * they need to be rotated for this,
+	 *
+	 * This could also be optimized like forward_diff_bezier */
+	int a;
+  	for(a=0; a<=it; a++) {
+		float t = (float)a / (float)it;
+
+		int i;
+		for(i=0; i<3; i++) {
+			p[i]= (-6*t + 6)*p0[i] + (18*t - 12)*p1[i] + (-18*t + 6)*p2[i] + (6*t)*p3[i];
+		}
+		Normalize(p);
+		p = (float *)(((char *)p)+stride);
+ 	}
+}
 
 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
 
@@ -1545,6 +1564,348 @@ static void alfa_bezpart(BezTriple *prevbezt, BezTriple *bezt, Nurb *nu, float *
 	}
 }
 
+/* make_bevel_list_3D_* funcs, at a minimum these must
+ * fill in the bezp->quat and bezp->dir values */
+
+/* correct non-cyclic cases by copying direction and rotation
+ * values onto the first & last end-points */
+static void bevel_list_cyclic_fix(BevList *bl)
+{
+	BevPoint *bevp, *bevp1;
+
+	bevp= (BevPoint *)(bl+1);
+	bevp1= bevp+1;
+	QUATCOPY(bevp->quat, bevp1->quat);
+	VECCOPY(bevp->dir, bevp1->dir);
+	VECCOPY(bevp->tan, bevp1->tan);
+	bevp= (BevPoint *)(bl+1);
+	bevp+= (bl->nr-1);
+	bevp1= bevp-1;
+	QUATCOPY(bevp->quat, bevp1->quat);
+	VECCOPY(bevp->dir, bevp1->dir);
+	VECCOPY(bevp->tan, bevp1->tan);
+}
+/* utility for make_bevel_list_3D_* funcs */
+static void bevel_list_calc_bisect(BevList *bl)
+{
+	BevPoint *bevp2, *bevp1, *bevp0;
+	int nr;
+
+	bevp2= (BevPoint *)(bl+1);
+	bevp1= bevp2+(bl->nr-1);
+	bevp0= bevp1-1;
+
+	nr= bl->nr;
+	while(nr--) {
+		/* totally simple */
+		VecBisect3(bevp1->dir, bevp0->vec, bevp1->vec, bevp2->vec);
+
+		bevp0= bevp1;
+		bevp1= bevp2;
+		bevp2++;
+	}
+}
+static void bevel_list_flip_tangents(BevList *bl)
+{
+	BevPoint *bevp2, *bevp1, *bevp0;
+	int nr;
+
+	bevp2= (BevPoint *)(bl+1);
+	bevp1= bevp2+(bl->nr-1);
+	bevp0= bevp1-1;
+
+	nr= bl->nr;
+	while(nr--) {
+		if(VecAngle2(bevp0->tan, bevp1->tan) > 90)
+			VecNegf(bevp1->tan);
+
+		bevp0= bevp1;
+		bevp1= bevp2;
+		bevp2++;
+	}
+}
+/* apply user tilt */
+static void bevel_list_apply_tilt(BevList *bl)
+{
+	BevPoint *bevp2, *bevp1, *bevp0;
+	int nr;
+	float q[4];
+
+	bevp2= (BevPoint *)(bl+1);
+	bevp1= bevp2+(bl->nr-1);
+	bevp0= bevp1-1;
+
+	nr= bl->nr;
+	while(nr--) {
+		AxisAngleToQuat(q, bevp1->dir, bevp1->alfa);
+		QuatMul(bevp1->quat, q, bevp1->quat);
+		NormalQuat(bevp1->quat);
+
+		bevp0= bevp1;
+		bevp1= bevp2;
+		bevp2++;
+	}
+}
+/* smooth quats, this function should be optimized, it can get slow with many iterations. */
+static void bevel_list_smooth(BevList *bl, int smooth_iter)
+{
+	BevPoint *bevp2, *bevp1, *bevp0;
+	int nr;
+
+	float q[4];
+	float bevp0_quat[4];
+	int a;
+
+	for(a=0; a < smooth_iter; a++) {
+
+		bevp2= (BevPoint *)(bl+1);
+		bevp1= bevp2+(bl->nr-1);
+		bevp0= bevp1-1;
+
+		nr= bl->nr;
+
+		if(bl->poly== -1) { /* check its not cyclic */
+			/* skip the first point */
+			bevp0= bevp1;
+			bevp1= bevp2;
+			bevp2++;
+			nr--;
+
+			bevp0= bevp1;
+			bevp1= bevp2;
+			bevp2++;
+			nr--;
+
+		}
+
+		QUATCOPY(bevp0_quat, bevp0->quat);
+
+		while(nr--) {
+			/* interpolate quats */
+			float zaxis[3] = {0,0,1}, cross[3], q2[4];
+			QuatInterpol(q, bevp0_quat, bevp2->quat, 0.5);
+			NormalQuat(q);
+
+			QuatMulVecf(q, zaxis);
+			Crossf(cross, zaxis, bevp1->dir);
+			AxisAngleToQuat(q2, cross, NormalizedVecAngle2(zaxis, bevp1->dir));
+			NormalQuat(q2);
+
+			QUATCOPY(bevp0_quat, bevp1->quat);
+			QuatMul(q, q2, q);
+			QuatInterpol(bevp1->quat, bevp1->quat, q, 0.5);
+			NormalQuat(bevp1->quat);
+
+
+			bevp0= bevp1;
+			bevp1= bevp2;
+			bevp2++;
+		}
+	}
+}
+
+static void make_bevel_list_3D_zup(BevList *bl)
+{
+	BevPoint *bevp2, *bevp1, *bevp0; /* standard for all make_bevel_list_3D_* funcs */
+	int nr;
+
+	bevp2= (BevPoint *)(bl+1);
+	bevp1= bevp2+(bl->nr-1);
+	bevp0= bevp1-1;
+
+	nr= bl->nr;
+	while(nr--) {
+		/* totally simple */
+		VecBisect3(bevp1->dir, bevp0->vec, bevp1->vec, bevp2->vec);
+		vectoquat(bevp1->dir, 5, 1, bevp1->quat);
+
+		bevp0= bevp1;
+		bevp1= bevp2;
+		bevp2++;
+	}
+}
+
+static void make_bevel_list_3D_minimum_twist(BevList *bl)
+{
+	BevPoint *bevp2, *bevp1, *bevp0; /* standard for all make_bevel_list_3D_* funcs */
+	int nr;
+	float q[4];
+
+	float cross_tmp[3];
+
+	bevel_list_calc_bisect(bl);
+
+	bevp2= (BevPoint *)(bl+1);
+	bevp1= bevp2+(bl->nr-1);
+	bevp0= bevp1-1;
+
+	nr= bl->nr;
+	while(nr--) {
+
+		if(nr+4 > bl->nr) { /* first time and second time, otherwise first point adjusts last */
+			vectoquat(bevp1->dir, 5, 1, bevp1->quat);
+		}
+		else {
+			float angle= NormalizedVecAngle2(bevp0->dir, bevp1->dir);
+
+			if(angle > 0.0f) { /* otherwise we can keep as is */
+				Crossf(cross_tmp, bevp0->dir, bevp1->dir);
+				AxisAngleToQuat(q, cross_tmp, angle);
+				QuatMul(bevp1->quat, q, bevp0->quat);
+			}
+			else {
+				QUATCOPY(bevp1->quat, bevp0->quat);
+			}
+		}
+
+		bevp0= bevp1;
+		bevp1= bevp2;
+		bevp2++;
+	}
+
+	if(bl->poly != -1) { /* check for cyclic */
+
+		/* Need to correct for the start/end points not matching
+		 * do this by calculating the tilt angle difference, then apply
+		 * the rotation gradually over the entire curve
+		 *
+		 * note that the split is between last and second last, rather then first/last as youd expect.
+		 *
+		 * real order is like this
+		 * 0,1,2,3,4 --> 1,2,3,4,0
+		 *
+		 * this is why we compare last with second last
+		 * */
+		float vec_1[3]= {0,1,0}, vec_2[3]= {0,1,0}, angle, ang_fac, cross_tmp[3];
+
+		BevPoint *bevp_first;
+		BevPoint *bevp_last;
+
+
+		bevp_first= (BevPoint *)(bl+1);
+		bevp_first+= bl->nr-1;
+		bevp_last = bevp_first;
+		bevp_last--;
+
+		/* quats and vec's are normalized, should not need to re-normalize */
+		QuatMulVecf(bevp_first->quat, vec_1);
+		QuatMulVecf(bevp_last->quat, vec_2);
+		Normalize(vec_1);
+		Normalize(vec_2);
+
+		/* align the vector, can avoid this and it looks 98% OK but
+		 * better to align the angle quat roll's before comparing */
+		{
+			Crossf(cross_tmp, bevp_last->dir, bevp_first->dir);
+			angle = NormalizedVecAngle2(bevp_first->dir, bevp_last->dir);
+			AxisAngleToQuat(q, cross_tmp, angle);
+			QuatMulVecf(q, vec_2);
+		}
+
+		angle= NormalizedVecAngle2(vec_1, vec_2);
+
+		/* flip rotation if needs be */
+		Crossf(cross_tmp, vec_1, vec_2);
+		Normalize(cross_tmp);
+		if(NormalizedVecAngle2(bevp_first->dir, cross_tmp) < 90/(180.0/M_PI))
+			angle = -angle;
+
+		bevp2= (BevPoint *)(bl+1);
+		bevp1= bevp2+(bl->nr-1);
+		bevp0= bevp1-1;
+
+		nr= bl->nr;
+		while(nr--) {
+			ang_fac= angle * (1.0f-((float)nr/bl->nr)); /* also works */
+
+			AxisAngleToQuat(q, bevp1->dir, ang_fac);
+			QuatMul(bevp1->quat, q, bevp1->quat);
+
+			bevp0= bevp1;
+			bevp1= bevp2;
+			bevp2++;
+		}
+	}
+}
+
+static void make_bevel_list_3D_tangent(BevList *bl)
+{
+	BevPoint *bevp2, *bevp1, *bevp0; /* standard for all make_bevel_list_3D_* funcs */
+	int nr;
+
+	float bevp0_tan[3], cross_tmp[3];
+
+	bevel_list_calc_bisect(bl);
+	if(bl->poly== -1) /* check its not cyclic */
+		bevel_list_cyclic_fix(bl); // XXX - run this now so tangents will be right before doing the flipping
+	bevel_list_flip_tangents(bl);
+
+	/* correct the tangents */
+	bevp2= (BevPoint *)(bl+1);
+	bevp1= bevp2+(bl->nr-1);
+	bevp0= bevp1-1;
+
+	nr= bl->nr;
+	while(nr--) {
+
+		Crossf(cross_tmp, bevp1->tan, bevp1->dir);
+		Crossf(bevp1->tan, cross_tmp, bevp1->dir);
+		Normalize(bevp1->tan);
+
+		bevp0= bevp1;
+		bevp1= bevp2;
+		bevp2++;
+	}
+
+
+	/* now for the real twist calc */
+	bevp2= (BevPoint *)(bl+1);
+	bevp1= bevp2+(bl->nr-1);
+	bevp0= bevp1-1;
+
+	VECCOPY(bevp0_tan, bevp0->tan);
+
+	nr= bl->nr;
+	while(nr--) {
+
+		/* make perpendicular, modify tan in place, is ok */
+		float cross_tmp[3];
+		float zero[3] = {0,0,0};
+
+		Crossf(cross_tmp, bevp1->tan, bevp1->dir);
+		Normalize(cross_tmp);
+		triatoquat(zero, cross_tmp, bevp1->tan, bevp1->quat); /* XXX - could be faster */
+
+		bevp0= bevp1;
+		bevp1= bevp2;
+		bevp2++;
+	}
+}
+
+void make_bevel_list_3D(BevList *bl, int smooth_iter, int twist_mode)
+{
+	switch(twist_mode) {
+	case CU_TWIST_TANGENT:
+		make_bevel_list_3D_tangent(bl);
+		break;
+	case CU_TWIST_MINIMUM:
+		make_bevel_list_3D_minimum_twist(bl);
+		break;
+	default: /* CU_TWIST_Z_UP default, pre 2.49c */
+		make_bevel_list_3D_zup(bl);
+	}
+
+	if(bl->poly== -1) /* check its not cyclic */
+		bevel_list_cyclic_fix(bl);
+
+	if(smooth_iter)
+		bevel_list_smooth(bl, smooth_iter);
+
+	bevel_list_apply_tilt(bl);
+}
+
+
+
 void makeBevelList(Object *ob)
 {
 	/*
@@ -1559,7 +1920,7 @@ void makeBevelList(Object *ob)
 	BPoint *bp;
 	BevList *bl, *blnew, *blnext;
 	BevPoint *bevp, *bevp2, *bevp1 = NULL, *bevp0;
-	float min, inp, x1, x2, y1, y2, vec[3], vec_prev[3], q[4], quat[4], quat_prev[4], cross[3];
+	float min, inp, x1, x2, y1, y2;
 	struct bevelsort *sortdata, *sd, *sd1;
 	int a, b, nr, poly, resolu, len=0;
 	int do_tilt, do_radius;
@@ -1664,7 +2025,15 @@ void makeBevelList(Object *ob)
 										 do_tilt	? &bevp->alfa : NULL,
 										 do_radius	? &bevp->radius : NULL,
 										 resolu, sizeof(BevPoint));
+
 						
+						if(cu->twist_mode==CU_TWIST_TANGENT) {
+							forward_diff_bezier_cotangent(
+													prevbezt->vec[1],	prevbezt->vec[2],
+													bezt->vec[0],		bezt->vec[1],
+													bevp->tan, resolu, sizeof(BevPoint));
+						}
+
 						/* indicate with handlecodes double points */
 						if(prevbezt->h1==prevbezt->h2) {
 							if(prevbezt->h1==0 || prevbezt->h1==HD_VECT) bevp->split_tag= TRUE;
@@ -1843,208 +2212,90 @@ void makeBevelList(Object *ob)
 		MEM_freeN(sortdata);
 	}
 
-	/* STEP 4: COSINES */
-	bl= cu->bev.first;
-	while(bl) {
-	
-		if(bl->nr < 2) {
-			/* do nothing */
-		}
-		else if(bl->nr==2) {	/* 2 pnt, treat separate */
-			bevp2= (BevPoint *)(bl+1);
-			bevp1= bevp2+1;
-
-			x1= bevp1->vec[0]- bevp2->vec[0];
-			y1= bevp1->vec[1]- bevp2->vec[1];
-
-			calc_bevel_sin_cos(x1, y1, -x1, -y1, &(bevp1->sina), &(bevp1->cosa));
-			bevp2->sina= bevp1->sina;
-			bevp2->cosa= bevp1->cosa;
+	/* STEP 4: 2D-COSINES or 3D ORIENTATION */
+	if((cu->flag & CU_3D)==0) { /* 3D */
+		bl= cu->bev.first;
+		while(bl) {
 
-			if(cu->flag & CU_3D) {	/* 3D */
-				float quat[4], q[4];
-			
-				VecSubf(vec, bevp1->vec, bevp2->vec);
-				
-				vectoquat(vec, 5, 1, quat);
-				
-				AxisAngleToQuat(q, vec, bevp1->alfa);
-				QuatMul(quat, q, quat);
-				
-				QuatToMat3(quat, bevp1->mat);
-				Mat3CpyMat3(bevp2->mat, bevp1->mat);
+			if(bl->nr < 2) {
+				/* do nothing */
 			}
+			else if(bl->nr==2) {	/* 2 pnt, treat separate */
+				bevp2= (BevPoint *)(bl+1);
+				bevp1= bevp2+1;
 
-		}	/* this has to be >2 points */
-		else if(cu->flag & CU_NO_TWIST && cu->flag & CU_3D && bl->poly != -1) {
-
-			/* Special case, cyclic curve with no twist. tricky... */
-
-			float quat[4], q[4], cross[3];
-
-			/* correcting a cyclic curve is more complicated, need to be corrected from both ends */
-			float *quat_tmp1, *quat_tmp2; /* store a quat in the matrix temporarily */
-			int iter_dir;
-			BevPoint *bevp_start= (BevPoint *)(bl+1);
-
-			/* loop over the points twice, once up, once back, accumulate the quat rotations
-			 * in both directions, then blend them in the 3rd loop and apply the tilt */
-			for(iter_dir = 0; iter_dir < 2; iter_dir++) {
+				x1= bevp1->vec[0]- bevp2->vec[0];
+				y1= bevp1->vec[1]- bevp2->vec[1];
 
+				calc_bevel_sin_cos(x1, y1, -x1, -y1, &(bevp1->sina), &(bevp1->cosa));
+				bevp2->sina= bevp1->sina;
+				bevp2->cosa= bevp1->cosa;
+			}
+			else {
 				bevp2= (BevPoint *)(bl+1);
 				bevp1= bevp2+(bl->nr-1);
 				bevp0= bevp1-1;
 
 				nr= bl->nr;
 				while(nr--) {
-	
-					/* Normalizes */
-					VecBisect3(vec, bevp0->vec, bevp1->vec, bevp2->vec);
+					x1= bevp1->vec[0]- bevp0->vec[0];
+					x2= bevp1->vec[0]- bevp2->vec[0];
+					y1= bevp1->vec[1]- bevp0->vec[1];
+					y2= bevp1->vec[1]- bevp2->vec[1];
 
-					if(bl->nr==nr+1) { /* first time */
-						vectoquat(vec, 5, 1, quat);
-					}
-					else {
-						float angle = NormalizedVecAngle2(vec_prev, vec);
-					
-						if(angle > 0.0f) { /* otherwise we can keep as is */
-							Crossf(cross, vec_prev, vec);
-							AxisAngleToQuat(q, cross, angle);
-							QuatMul(quat, q, quat_prev);
-						}
-						else {
-							QUATCOPY(quat, quat_prev);
-						}
-					}
-					QUATCOPY(quat_prev, quat); /* quat_prev can't have the tilt applied */
-					VECCOPY(vec_prev, vec);
-
-					if(iter_dir==0) { /* up, first time */
-						quat_tmp1= (float *)bevp1->mat;
+					calc_bevel_sin_cos(x1, y1, x2, y2, &(bevp1->sina), &(bevp1->cosa));
 
-						bevp0= bevp1;
-						bevp1= bevp2;
-						bevp2++;
-					}
-					else { /* down second time */
-						quat_tmp1= ((float *)bevp1->mat)+4;
-
-						bevp2= bevp1;
-						bevp1= bevp0;
-						bevp0--;
-
-						/* wrap around */
-						if (bevp0 < bevp_start)
-							bevp0= bevp_start+(bl->nr-1);
-					}
-
-					QUATCOPY(quat_tmp1, quat);
+					bevp0= bevp1;
+					bevp1= bevp2;
+					bevp2++;
 				}
-			}
-
-			/* Now interpolate the 2 quats and apply tilt */
 
-			bevp2= (BevPoint *)(bl+1);
-			bevp1= bevp2+(bl->nr-1);
-			bevp0= bevp1-1;
-
-			nr= bl->nr;
-			while(nr--) {
-
-				VecBisect3(vec, bevp0->vec, bevp1->vec, bevp2->vec);
-
-				quat_tmp1= (float *)bevp1->mat;
-				quat_tmp2= quat_tmp1+4;
-
-				/* blend the 2 rotations gathered from both directions */
-				QuatInterpol(quat, quat_tmp1, quat_tmp2, 1.0 - (((float)nr)/bl->nr));
-
-				AxisAngleToQuat(q, vec, bevp1->alfa);
-				QuatMul(quat, q, quat);
-				QuatToMat3(quat, bevp1->mat);
-				
-				/* generic */
-				x1= bevp1->vec[0]- bevp0->vec[0];
-				x2= bevp1->vec[0]- bevp2->vec[0];
-				y1= bevp1->vec[1]- bevp0->vec[1];
-				y2= bevp1->vec[1]- bevp2->vec[1];
-			
-				calc_bevel_sin_cos(x1, y1, x2, y2, &(bevp1->sina), &(bevp1->cosa));
-				
-				bevp0= bevp1;
-				bevp1= bevp2;
-				bevp2++;
-			}
-		}
-		else {
-			/* Any curve with 3 or more points */
-
-			bevp2= (BevPoint *)(bl+1);
-			bevp1= bevp2+(bl->nr-1);
-			bevp0= bevp1-1;
-
-			nr= bl->nr;
-			while(nr--) {
-
-				if(cu->flag & CU_3D) {	/* 3D */
-
-					/* Normalizes */
-					VecBisect3(vec, bevp0->vec, bevp1->vec, bevp2->vec);
-
-					if(bl->nr==nr+1 || !(cu->flag & CU_NO_TWIST)) { /* first time */
-						vectoquat(vec, 5, 1, quat);
-					}
-					else {
-						float angle = NormalizedVecAngle2(vec_prev, vec);
-
-						if(angle > 0.0f) { /* otherwise we can keep as is */
-							Crossf(cross, vec_prev, vec);
-							AxisAngleToQuat(q, cross, angle);
-							QuatMul(quat, q, quat_prev);
-						}
-						else {
-							QUATCOPY(quat, quat_prev);
-						}
-					}
-					QUATCOPY(quat_prev, quat); /* quat_prev can't have the tilt applied */
-					VECCOPY(vec_prev, vec);
-					
-					AxisAngleToQuat(q, vec, bevp1->alfa);
-					QuatMul(quat, q, quat);
-					QuatToMat3(quat, bevp1->mat);
-				}
-
-				x1= bevp1->vec[0]- bevp0->vec[0];
-				x2= bevp1->vec[0]- bevp2->vec[0];
-				y1= bevp1->vec[1]- bevp0->vec[1];
-				y2= bevp1->vec[1]- bevp2->vec[1];
-
-				calc_bevel_sin_cos(x1, y1, x2, y2, &(bevp1->sina), &(bevp1->cosa));
-
-
-				bevp0= bevp1;
-				bevp1= bevp2;
-				bevp2++;
-			}
-
-			/* correct non-cyclic cases */
-			if(bl->poly== -1) {
-				if(bl->nr>2) {
+				/* correct non-cyclic cases */
+				if(bl->poly== -1) {
 					bevp= (BevPoint *)(bl+1);
 					bevp1= bevp+1;
 					bevp->sina= bevp1->sina;
 					bevp->cosa= bevp1->cosa;
-					Mat3CpyMat3(bevp->mat, bevp1->mat);
 					bevp= (BevPoint *)(bl+1);
 					bevp+= (bl->nr-1);
 					bevp1= bevp-1;
 					bevp->sina= bevp1->sina;
 					bevp->cosa= bevp1->cosa;
-					Mat3CpyMat3(bevp->mat, bevp1->mat);
 				}
 			}
+			bl= bl->next;
+		}
+	}
+	else { /* 3D Curves */
+		bl= cu->bev.first;
+		while(bl) {
+
+			if(bl->nr < 2) {
+				/* do nothing */
+			}
+			else if(bl->nr==2) {	/* 2 pnt, treat separate */
+				float q[4];
+
+				bevp2= (BevPoint *)(bl+1);
+				bevp1= bevp2+1;
+
+				/* simple quat/dir */
+				VecSubf(bevp1->dir, bevp1->vec, bevp2->vec);
+				Normalize(bevp1->dir);
+				
+				vectoquat(bevp1->dir, 5, 1, bevp1->quat);
+				
+				AxisAngleToQuat(q, bevp1->dir, bevp1->alfa);
+				QuatMul(bevp1->quat, q, bevp1->quat);
+				NormalQuat(bevp1->quat);
+				VECCOPY(bevp2->dir, bevp1->dir);
+				QUATCOPY(bevp2->quat, bevp1->quat);
+			}
+			else {
+				make_bevel_list_3D(bl, (int)(resolu*cu->twist_smooth), cu->twist_mode);
+			}
+			bl= bl->next;
 		}
-		bl= bl->next;
 	}
 }
 
diff --git a/source/blender/blenkernel/intern/displist.c b/source/blender/blenkernel/intern/displist.c
index 99baa17b582..d87dbc833c5 100644
--- a/source/blender/blenkernel/intern/displist.c
+++ b/source/blender/blenkernel/intern/displist.c
@@ -1599,7 +1599,7 @@ void makeDispListCurveTypes(Scene *scene, Object *ob, int forOrco)
 										vec[1]= fp1[2];
 										vec[2]= 0.0;
 										
-										Mat3MulVecfl(bevp->mat, vec);
+										QuatMulVecf(bevp->quat, vec);
 										
 										data[0]= bevp->vec[0] + fac*vec[0];
 										data[1]= bevp->vec[1] + fac*vec[1];
diff --git a/source/blender/editors/space_view3d/drawobject.c b/source/blender/editors/space_view3d/drawobject.c
index 53fd6972641..69bc8e1f67e 100644
--- a/source/blender/editors/space_view3d/drawobject.c
+++ b/source/blender/editors/space_view3d/drawobject.c
@@ -4129,21 +4129,21 @@ static void drawnurb(Scene *scene, View3D *v3d, RegionView3D *rv3d, Base *base,
 			
 			while (nr-->0) { /* accounts for empty bevel lists */
 				float fac= bevp->radius * ts->normalsize;
-				float ox,oy,oz; // Offset perpendicular to the curve
-				float dx,dy,dz; // Delta along the curve
+				float vec_a[3] = { fac,0, 0}; // Offset perpendicular to the curve
+				float vec_b[3] = {-fac,0, 0}; // Delta along the curve
+
+				QuatMulVecf(bevp->quat, vec_a);
+				QuatMulVecf(bevp->quat, vec_b);
+				VecAddf(vec_a, vec_a, bevp->vec);
+				VecAddf(vec_b, vec_b, bevp->vec);
 				
-				ox = fac*bevp->mat[0][0];
-				oy = fac*bevp->mat[0][1];
-				oz = fac*bevp->mat[0][2];
-			
-				dx = fac*bevp->mat[2][0];
-				dy = fac*bevp->mat[2][1];
-				dz = fac*bevp->mat[2][2];
+				VECSUBFAC(vec_a, vec_a, bevp->dir, fac);
+				VECSUBFAC(vec_b, vec_b, bevp->dir, fac);
 
 				glBegin(GL_LINE_STRIP);
-				glVertex3f(bevp->vec[0] - ox - dx, bevp->vec[1] - oy - dy, bevp->vec[2] - oz - dz);
+				glVertex3fv(vec_a);
 				glVertex3fv(bevp->vec);
-				glVertex3f(bevp->vec[0] + ox - dx, bevp->vec[1] + oy - dy, bevp->vec[2] + oz - dz);
+				glVertex3fv(vec_b);
 				glEnd();
 				
 				bevp += skip+1;
diff --git a/source/blender/makesdna/DNA_curve_types.h b/source/blender/makesdna/DNA_curve_types.h
index 3a7cf874dbd..013491187ae 100644
--- a/source/blender/makesdna/DNA_curve_types.h
+++ b/source/blender/makesdna/DNA_curve_types.h
@@ -71,7 +71,9 @@ typedef struct BevList {
 #
 #
 typedef struct BevPoint {
-	float vec[3], mat[3][3], alfa, radius, sina, cosa;
+	float vec[3], alfa, radius;
+	float sina, cosa;				/* 2D Only */
+	float dir[3], tan[3], quat[4];	/* 3D Only */
 	short split_tag, dupe_tag;
 } BevPoint;
 
@@ -165,7 +167,8 @@ typedef struct Curve {
 
 	int texflag; /* keep an int because of give_obdata_texspace() */
 
-	short drawflag, pad[3];
+	short drawflag, twist_mode,  pad[2];
+	float twist_smooth, pad2;
 
 	short pathlen, totcol;
 	short flag, bevresol;
@@ -174,7 +177,7 @@ typedef struct Curve {
 	/* default */
 	short resolu, resolv;
 	short resolu_ren, resolv_ren;
-	
+
 	/* edit, index in nurb list */
 	int actnu;
 	/* edit, last selected bpoint */
@@ -231,7 +234,12 @@ typedef struct Curve {
 #define CU_RETOPO               1024
 #define CU_DS_EXPAND	2048
 
-#define CU_NO_TWIST		4096
+/* twist mode */
+#define CU_TWIST_Z_UP			0
+// #define CU_TWIST_Y_UP			1 // not used yet
+// #define CU_TWIST_X_UP			2
+#define CU_TWIST_MINIMUM		3
+#define CU_TWIST_TANGENT		4
 
 /* spacemode */
 #define CU_LEFT			0
diff --git a/source/blender/makesrna/intern/rna_curve.c b/source/blender/makesrna/intern/rna_curve.c
index 91628bab376..d09546231d5 100644
--- a/source/blender/makesrna/intern/rna_curve.c
+++ b/source/blender/makesrna/intern/rna_curve.c
@@ -664,6 +664,12 @@ static void rna_def_curve(BlenderRNA *brna)
 	StructRNA *srna;
 	PropertyRNA *prop;
 	
+	static EnumPropertyItem curve_twist_mode_items[] = {
+			{CU_TWIST_Z_UP, "Z_UP", 0, "Z-Up", "Use Z-Up axis to calculate the curve twist at each point"},
+			{CU_TWIST_MINIMUM, "MINIMUM", 0, "Minimum", "Use the least twist over the entire curve"},
+			{CU_TWIST_TANGENT, "TANGENT", 0, "Tangent", "Use the tangent to calculate twist"},
+			{0, NULL, 0, NULL, NULL}};
+
 	srna= RNA_def_struct(brna, "Curve", "ID");
 	RNA_def_struct_ui_text(srna, "Curve", "Curve datablock storing curves, splines and NURBS.");
 	RNA_def_struct_ui_icon(srna, ICON_CURVE_DATA);
@@ -782,12 +788,22 @@ static void rna_def_curve(BlenderRNA *brna)
 	RNA_def_property_boolean_sdna(prop, NULL, "flag", CU_BACK);
 	RNA_def_property_ui_text(prop, "Back", "Draw filled back for extruded/beveled curves.");
 	RNA_def_property_update(prop, 0, "rna_Curve_update_data");
-	
-	prop= RNA_def_property(srna, "use_twist_correction", PROP_BOOLEAN, PROP_NONE);
-	RNA_def_property_boolean_sdna(prop, NULL, "flag", CU_NO_TWIST);
-	RNA_def_property_ui_text(prop, "Minimal Twist", "Correct for twisting.");
+
+	prop= RNA_def_property(srna, "twist_mode", PROP_ENUM, PROP_NONE);
+	RNA_def_property_enum_sdna(prop, NULL, "twist_mode");
+	RNA_def_property_enum_items(prop, curve_twist_mode_items);
+	RNA_def_property_ui_text(prop, "Twist Method", "The type of tilt calculation for 3D Curves.");
+	RNA_def_property_update(prop, 0, "rna_Curve_update_data");
+
+	// XXX - would be nice to have a better way to do this, only add for testing.
+	prop= RNA_def_property(srna, "twist_smooth", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "twist_smooth");
+	RNA_def_property_ui_range(prop, 0, 100.0, 0.1, 0);
+	RNA_def_property_ui_text(prop, "Twist Smooth", "Smoothing iteration for tangents");
 	RNA_def_property_update(prop, 0, "rna_Curve_update_data");
 
+
+
 	prop= RNA_def_property(srna, "retopo", PROP_BOOLEAN, PROP_NONE);
 	RNA_def_property_boolean_sdna(prop, NULL, "flag", CU_RETOPO);
 	RNA_def_property_ui_text(prop, "Retopo", "Turn on the re-topology tool.");