barycentric transform utility geometry function.

From 2 triangles and 1 point, the relative position between the point and the first triangle is applied to the second triangle to find the target point.
the barycentric weights are calculated in 2D space with a signed area so values outside the triangle bounds are supported.

wrapped by python:
 pt_to = Geometry.BarycentricTransform(pt_from, t1a, t1b, t1c, t2a, t1b, t1c)

NOTE: 
- moved some barycentric weight functions out of projection painting into the math lib.
- ended up making some of the math functions use const args.
TODO:
- support exceptional cases. zero area tries and similar.

9 files changed, 176 insertions, 82 deletions

diff --git a/source/blender/blenlib/BLI_math_geom.h b/source/blender/blenlib/BLI_math_geom.h
index 49d335f6c5c..e676001fba5 100644
--- a/source/blender/blenlib/BLI_math_geom.h
+++ b/source/blender/blenlib/BLI_math_geom.h
@@ -37,12 +37,13 @@ extern "C" {
 void cent_tri_v3(float r[3], float a[3], float b[3], float c[3]);
 void cent_quad_v3(float r[3], float a[3], float b[3], float c[3], float d[3]);
 
-float normal_tri_v3(float r[3], float a[3], float b[3], float c[3]);
-float normal_quad_v3(float r[3], float a[3], float b[3], float c[3], float d[3]);
+float normal_tri_v3(float r[3], const float a[3], const float b[3], const float c[3]);
+float normal_quad_v3(float r[3], const float a[3], const float b[3], const float c[3], const float d[3]);
 
-float area_tri_v2(float a[2], float b[2], float c[2]);
-float area_tri_v3(float a[3], float b[3], float c[3]);
-float area_quad_v3(float a[3], float b[3], float c[3], float d[3]);
+float area_tri_v2(const float a[2], const float b[2], const float c[2]);
+float area_tri_signed_v2(const float v1[2], const float v2[2], const float v3[2]);
+float area_tri_v3(const float a[3], const float b[3], const float c[3]);
+float area_quad_v3(const float a[3], const float b[3], const float c[3], const float d[3]);
 float area_poly_v3(int nr, float verts[][3], float normal[3]);
 
 /********************************* Distance **********************************/
@@ -120,6 +121,13 @@ void interp_weights_poly_v3(float w[], float v[][3], int n, float p[3]);
 void interp_cubic_v3(float x[3], float v[3],
 	float x1[3], float v1[3], float x2[3], float v2[3], float t);
 
+void barycentric_transform(float pt_tar[3], float const pt_src[3],
+	const float tri_tar_p1[3], const float tri_tar_p2[3], const float tri_tar_p3[3],
+	const float tri_src_p1[3], const float tri_src_p2[3], const float tri_src_p3[3]);
+
+void barycentric_weights_v2(const float v1[2], const float v2[2], const float v3[2],
+	const float co[2], float w[3]);
+
 /***************************** View & Projection *****************************/
 
 void lookat_m4(float mat[4][4], float vx, float vy, 
diff --git a/source/blender/blenlib/BLI_math_rotation.h b/source/blender/blenlib/BLI_math_rotation.h
index b221d89487f..b640cf6dc3c 100644
--- a/source/blender/blenlib/BLI_math_rotation.h
+++ b/source/blender/blenlib/BLI_math_rotation.h
@@ -70,7 +70,8 @@ void mat3_to_quat(float q[4], float mat[3][3]);
 void mat4_to_quat(float q[4], float mat[4][4]);
 void tri_to_quat(float q[4], float a[3], float b[3], float c[3]);
 void vec_to_quat(float q[4], float vec[3], short axis, short upflag);
-void rotation_between_vecs_to_quat(float q[4], float v1[3], float v2[3]);
+/* note: v1 and v2 must be normalized */
+void rotation_between_vecs_to_quat(float q[4], const float v1[3], const float v2[3]);
 
 /* TODO: don't what this is, but it's not the same as mat3_to_quat */
 void mat3_to_quat_is_ok(float q[4], float mat[3][3]);
diff --git a/source/blender/blenlib/BLI_math_vector.h b/source/blender/blenlib/BLI_math_vector.h
index e915a9a85f3..26e7ff5abe9 100644
--- a/source/blender/blenlib/BLI_math_vector.h
+++ b/source/blender/blenlib/BLI_math_vector.h
@@ -51,8 +51,8 @@ extern "C" {
 MINLINE void zero_v2(float r[2]);
 MINLINE void zero_v3(float r[3]);
 
-MINLINE void copy_v2_v2(float r[2], float a[2]);
-MINLINE void copy_v3_v3(float r[3], float a[3]);
+MINLINE void copy_v2_v2(float r[2], const float a[2]);
+MINLINE void copy_v3_v3(float r[3], const float a[3]);
 
 MINLINE void swap_v2_v2(float a[2], float b[2]);
 MINLINE void swap_v3_v3(float a[3], float b[3]);
@@ -67,7 +67,7 @@ MINLINE void add_v3_v3v3(float r[3], float a[3], float b[3]);
 MINLINE void sub_v2_v2(float r[2], float a[2]);
 MINLINE void sub_v2_v2v2(float r[2], float a[2], float b[2]);
 MINLINE void sub_v3_v3(float r[3], float a[3]);
-MINLINE void sub_v3_v3v3(float r[3], float a[3], float b[3]);
+MINLINE void sub_v3_v3v3(float r[3], const float a[3], const float b[3]);
 
 MINLINE void mul_v2_fl(float r[2], float f);
 MINLINE void mul_v3_fl(float r[3], float f);
@@ -87,7 +87,7 @@ MINLINE float dot_v2v2(float a[2], float b[2]);
 MINLINE float dot_v3v3(float a[3], float b[3]);
 
 MINLINE float cross_v2v2(float a[2], float b[2]);
-MINLINE void cross_v3_v3v3(float r[3], float a[3], float b[3]);
+MINLINE void cross_v3_v3v3(float r[3], const float a[3], const float b[3]);
 
 MINLINE void star_m3_v3(float R[3][3],float a[3]);
 
@@ -103,11 +103,11 @@ MINLINE float normalize_v3(float r[3]);
 
 /******************************* Interpolation *******************************/
 
-void interp_v2_v2v2(float r[2], float a[2], float b[2], float t);
-void interp_v2_v2v2v2(float r[2], float a[2], float b[2], float c[3], float t[3]);
-void interp_v3_v3v3(float r[3], float a[3], float b[3], float t);
-void interp_v3_v3v3v3(float p[3], float v1[3], float v2[3], float v3[3], float w[3]);
-void interp_v3_v3v3v3v3(float p[3], float v1[3], float v2[3], float v3[3], float v4[3], float w[4]);
+void interp_v2_v2v2(float r[2], const float a[2], const float b[2], const float t);
+void interp_v2_v2v2v2(float r[2], const float a[2], const float b[2], const float c[3], const float t[3]);
+void interp_v3_v3v3(float r[3], const float a[3], const float b[3], const float t);
+void interp_v3_v3v3v3(float p[3], const float v1[3], const float v2[3], const float v3[3], const float w[3]);
+void interp_v3_v3v3v3v3(float p[3], const float v1[3], const float v2[3], const float v3[3], const float v4[3], const float w[4]);
 
 void mid_v3_v3v3(float r[3], float a[3], float b[3]);
 
@@ -130,7 +130,7 @@ float angle_v2v2v2(float a[2], float b[2], float c[2]);
 float angle_normalized_v2v2(float a[2], float b[2]);
 float angle_v3v3(float a[2], float b[2]);
 float angle_v3v3v3(float a[2], float b[2], float c[2]);
-float angle_normalized_v3v3(float a[3], float b[3]);
+float angle_normalized_v3v3(const float v1[3], const float v2[3]);
 void angle_tri_v3(float angles[3], const float v1[3], const float v2[3], const float v3[3]);
 void angle_quad_v3(float angles[4], const float v1[3], const float v2[3], const float v3[3], const float v4[3]);
 
diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c
index 3e714b384a7..dd773cb00c1 100644
--- a/source/blender/blenlib/intern/math_geom.c
+++ b/source/blender/blenlib/intern/math_geom.c
@@ -49,7 +49,7 @@ void cent_quad_v3(float *cent, float *v1, float *v2, float *v3, float *v4)
 	cent[2]= 0.25f*(v1[2]+v2[2]+v3[2]+v4[2]);
 }
 
-float normal_tri_v3(float *n, float *v1, float *v2, float *v3)
+float normal_tri_v3(float n[3], const float v1[3], const float v2[3], const float v3[3])
 {
 	float n1[3],n2[3];
 
@@ -65,7 +65,7 @@ float normal_tri_v3(float *n, float *v1, float *v2, float *v3)
 	return normalize_v3(n);
 }
 
-float normal_quad_v3(float *n, float *v1, float *v2, float *v3, float *v4)
+float normal_quad_v3(float n[3], const float v1[3], const float v2[3], const float v3[3], const float v4[3])
 {
 	/* real cross! */
 	float n1[3],n2[3];
@@ -85,13 +85,17 @@ float normal_quad_v3(float *n, float *v1, float *v2, float *v3, float *v4)
 	return normalize_v3(n);
 }
 
-float area_tri_v2(float *v1, float *v2, float *v3)
+float area_tri_v2(const float v1[2], const float v2[2], const float v3[2])
 {
-	return (float)(0.5*fabs((v1[0]-v2[0])*(v2[1]-v3[1]) + (v1[1]-v2[1])*(v3[0]-v2[0])));
+	return (float)(0.5f*fabs((v1[0]-v2[0])*(v2[1]-v3[1]) + (v1[1]-v2[1])*(v3[0]-v2[0])));
 }
 
+float area_tri_signed_v2(const float v1[2], const float v2[2], const float v3[2])
+{
+   return (float)(0.5f*((v1[0]-v2[0])*(v2[1]-v3[1]) + (v1[1]-v2[1])*(v3[0]-v2[0])));
+}
 
-float area_quad_v3(float *v1, float *v2, float *v3,  float *v4)  /* only convex Quadrilaterals */
+float area_quad_v3(const float v1[3], const float v2[3], const float v3[3], const float v4[3])  /* only convex Quadrilaterals */
 {
 	float len, vec1[3], vec2[3], n[3];
 
@@ -108,7 +112,7 @@ float area_quad_v3(float *v1, float *v2, float *v3,  float *v4)  /* only convex
 	return (len/2.0f);
 }
 
-float area_tri_v3(float *v1, float *v2, float *v3)  /* Triangles */
+float area_tri_v3(const float v1[3], const float v2[3], const float v3[3])  /* Triangles */
 {
 	float len, vec1[3], vec2[3], n[3];
 
@@ -1334,6 +1338,77 @@ void interp_weights_face_v3(float *w,float *v1, float *v2, float *v3, float *v4,
 	}
 }
 
+/* used by projection painting
+ * note: using area_tri_signed_v2 means locations outside the triangle are correctly weighted */
+void barycentric_weights_v2(const float v1[2], const float v2[2], const float v3[2], const float co[2], float w[3])
+{
+   float wtot_inv, wtot;
+
+   w[0] = area_tri_signed_v2(v2, v3, co);
+   w[1] = area_tri_signed_v2(v3, v1, co);
+   w[2] = area_tri_signed_v2(v1, v2, co);
+   wtot = w[0]+w[1]+w[2];
+
+   if (wtot != 0.0f) {
+       wtot_inv = 1.0f/wtot;
+
+       w[0] = w[0]*wtot_inv;
+       w[1] = w[1]*wtot_inv;
+       w[2] = w[2]*wtot_inv;
+   }
+   else /* dummy values for zero area face */
+       w[0] = w[1] = w[2] = 1.0f/3.0f;
+}
+
+/* given 2 triangles in 3D space, and a point in relation to the first triangle.
+ * calculate the location of a point in relation to the second triangle.
+ * Useful for finding relative positions with geometry */
+void barycentric_transform(float pt_tar[3], float const pt_src[3],
+		const float tri_tar_p1[3], const float tri_tar_p2[3], const float tri_tar_p3[3],
+		const float tri_src_p1[3], const float tri_src_p2[3], const float tri_src_p3[3])
+{
+	/* this works by moving the source triangle so its normal is pointing on the Z
+	 * axis where its barycentric wights can be calculated in 2D and its Z offset can
+	 *  be re-applied. The weights are applied directly to the targets 3D points and the
+	 *  z-depth is used to scale the targets normal as an offset.
+	 * This saves transforming the target into its Z-Up orientation and back (which could also work) */
+	const float z_up[3] = {0, 0, 1};
+	float no_tar[3], no_src[3];
+	float quat_src[4];
+	float pt_src_xy[3];
+	float  tri_xy_src[3][3];
+	float w_src[3];
+	float area_tar, area_src;
+	float z_ofs_src;
+
+	normal_tri_v3(no_tar, tri_tar_p1, tri_tar_p2, tri_tar_p3);
+	normal_tri_v3(no_src, tri_src_p1, tri_src_p2, tri_src_p3);
+
+	rotation_between_vecs_to_quat(quat_src, no_src, z_up);
+	normalize_qt(quat_src);
+
+	copy_v3_v3(pt_src_xy, pt_src);
+	copy_v3_v3(tri_xy_src[0], tri_src_p1);
+	copy_v3_v3(tri_xy_src[1], tri_src_p2);
+	copy_v3_v3(tri_xy_src[2], tri_src_p3);
+
+	/* make the source tri xy space */
+	mul_qt_v3(quat_src, pt_src_xy);
+	mul_qt_v3(quat_src, tri_xy_src[0]);
+	mul_qt_v3(quat_src, tri_xy_src[1]);
+	mul_qt_v3(quat_src, tri_xy_src[2]);
+
+	barycentric_weights_v2(tri_xy_src[0], tri_xy_src[1], tri_xy_src[2], pt_src_xy, w_src);
+	interp_v3_v3v3v3(pt_tar, tri_tar_p1, tri_tar_p2, tri_tar_p3, w_src);
+
+	area_tar= sqrtf(area_tri_v3(tri_tar_p1, tri_tar_p2, tri_tar_p3));
+	area_src= sqrtf(area_tri_v2(tri_xy_src[0], tri_xy_src[1], tri_xy_src[2]));
+
+	z_ofs_src= tri_xy_src[0][2] - pt_src_xy[2];
+	madd_v3_v3v3fl(pt_tar, pt_tar, no_tar, (z_ofs_src / area_src) * area_tar);
+}
+
+
 /* Mean value weights - smooth interpolation weights for polygons with
  * more than 3 vertices */
 static float mean_value_half_tan(float *v1, float *v2, float *v3)
@@ -1789,3 +1864,4 @@ void vcloud_estimate_transform(int list_size, float (*pos)[3], float *weight,flo
 		}
 	}
 }
+
diff --git a/source/blender/blenlib/intern/math_rotation.c b/source/blender/blenlib/intern/math_rotation.c
index 084db725409..61d5eeaaa2c 100644
--- a/source/blender/blenlib/intern/math_rotation.c
+++ b/source/blender/blenlib/intern/math_rotation.c
@@ -327,7 +327,7 @@ void normalize_qt(float *q)
 	}
 }
 
-void rotation_between_vecs_to_quat(float *q, float v1[3], float v2[3])
+void rotation_between_vecs_to_quat(float *q, const float v1[3], const float v2[3])
 {
 	float axis[3];
 	float angle;
diff --git a/source/blender/blenlib/intern/math_vector.c b/source/blender/blenlib/intern/math_vector.c
index 3188587a96c..d99d96d28ff 100644
--- a/source/blender/blenlib/intern/math_vector.c
+++ b/source/blender/blenlib/intern/math_vector.c
@@ -34,7 +34,7 @@
 
 //******************************* Interpolation *******************************/
 
-void interp_v2_v2v2(float *target, float *a, float *b, float t)
+void interp_v2_v2v2(float *target, const float *a, const float *b, const float t)
 {
 	float s = 1.0f-t;
 
@@ -44,13 +44,13 @@ void interp_v2_v2v2(float *target, float *a, float *b, float t)
 
 /* weight 3 2D vectors,
  * 'w' must be unit length but is not a vector, just 3 weights */
-void interp_v2_v2v2v2(float p[2], float v1[2], float v2[2], float v3[2], float w[3])
+void interp_v2_v2v2v2(float p[2], const float v1[2], const float v2[2], const float v3[2], const float w[3])
 {
 	p[0] = v1[0]*w[0] + v2[0]*w[1] + v3[0]*w[2];
 	p[1] = v1[1]*w[0] + v2[1]*w[1] + v3[1]*w[2];
 }
 
-void interp_v3_v3v3(float *target, float *a, float *b, float t)
+void interp_v3_v3v3(float target[3], const float a[3], const float b[3], const float t)
 {
 	float s = 1.0f-t;
 
@@ -61,7 +61,7 @@ void interp_v3_v3v3(float *target, float *a, float *b, float t)
 
 /* weight 3 vectors,
  * 'w' must be unit length but is not a vector, just 3 weights */
-void interp_v3_v3v3v3(float p[3], float v1[3], float v2[3], float v3[3], float w[3])
+void interp_v3_v3v3v3(float p[3], const float v1[3], const float v2[3], const float v3[3], const float w[3])
 {
 	p[0] = v1[0]*w[0] + v2[0]*w[1] + v3[0]*w[2];
 	p[1] = v1[1]*w[0] + v2[1]*w[1] + v3[1]*w[2];
@@ -70,7 +70,7 @@ void interp_v3_v3v3v3(float p[3], float v1[3], float v2[3], float v3[3], float w
 
 /* weight 3 vectors,
  * 'w' must be unit length but is not a vector, just 3 weights */
-void interp_v3_v3v3v3v3(float p[3], float v1[3], float v2[3], float v3[3], float v4[3], float w[4])
+void interp_v3_v3v3v3v3(float p[3], const float v1[3], const float v2[3], const float v3[3], const float v4[3], const float w[4])
 {
 	p[0] = v1[0]*w[0] + v2[0]*w[1] + v3[0]*w[2] + v4[0]*w[3];
 	p[1] = v1[1]*w[0] + v2[1]*w[1] + v3[1]*w[2] + v4[1]*w[3];
@@ -191,7 +191,7 @@ float angle_v2v2(float *v1, float *v2)
 	return angle_normalized_v2v2(vec1, vec2);
 }
 
-float angle_normalized_v3v3(float *v1, float *v2)
+float angle_normalized_v3v3(const float v1[3], const float v2[3])
 {
 	/* this is the same as acos(dot_v3v3(v1, v2)), but more accurate */
 	if (dot_v3v3(v1, v2) < 0.0f) {
diff --git a/source/blender/blenlib/intern/math_vector_inline.c b/source/blender/blenlib/intern/math_vector_inline.c
index cb629312712..16f35dbc5fa 100644
--- a/source/blender/blenlib/intern/math_vector_inline.c
+++ b/source/blender/blenlib/intern/math_vector_inline.c
@@ -45,13 +45,13 @@ MINLINE void zero_v3(float r[3])
 	r[2]= 0.0f;
 }
 
-MINLINE void copy_v2_v2(float r[2], float a[2])
+MINLINE void copy_v2_v2(float r[2], const float a[2])
 {
 	r[0]= a[0];
 	r[1]= a[1];
 }
 
-MINLINE void copy_v3_v3(float r[3], float a[3])
+MINLINE void copy_v3_v3(float r[3], const float a[3])
 {
 	r[0]= a[0];
 	r[1]= a[1];
@@ -118,7 +118,7 @@ MINLINE void sub_v3_v3(float *r, float *a)
 	r[2] -= a[2];
 }
 
-MINLINE void sub_v3_v3v3(float *r, float *a, float *b)
+MINLINE void sub_v3_v3v3(float r[3], const float a[3], const float b[3])
 {
 	r[0]= a[0] - b[0];
 	r[1]= a[1] - b[1];
@@ -216,7 +216,7 @@ MINLINE float cross_v2v2(float a[2], float b[2])
 	 return a[0]*b[1] - a[1]*b[0];
 }
 
-MINLINE void cross_v3_v3v3(float r[3], float a[3], float b[3])
+MINLINE void cross_v3_v3v3(float r[3], const float a[3], const float b[3])
 {
 	r[0]= a[1]*b[2] - a[2]*b[1];
 	r[1]= a[2]*b[0] - a[0]*b[2];
diff --git a/source/blender/editors/sculpt_paint/paint_image.c b/source/blender/editors/sculpt_paint/paint_image.c
index 058dab4ba15..526eb2d6661 100644
--- a/source/blender/editors/sculpt_paint/paint_image.c
+++ b/source/blender/editors/sculpt_paint/paint_image.c
@@ -464,40 +464,14 @@ static int project_bucket_offset_safe(const ProjPaintState *ps, const float proj
 
 #define SIDE_OF_LINE(pa, pb, pp)	((pa[0]-pp[0])*(pb[1]-pp[1]))-((pb[0]-pp[0])*(pa[1]-pp[1]))
 
-static float AreaSignedF2Dfl(float *v1, float *v2, float *v3)
-{
-   return (float)(0.5f*((v1[0]-v2[0])*(v2[1]-v3[1]) +
-(v1[1]-v2[1])*(v3[0]-v2[0])));
-}
-
-static void BarycentricWeights2f(float pt[2], float v1[2], float v2[2], float v3[2], float w[3])
-{
-   float wtot_inv, wtot;
-
-   w[0] = AreaSignedF2Dfl(v2, v3, pt);
-   w[1] = AreaSignedF2Dfl(v3, v1, pt);
-   w[2] = AreaSignedF2Dfl(v1, v2, pt);
-   wtot = w[0]+w[1]+w[2];
-
-   if (wtot != 0.0f) {
-       wtot_inv = 1.0f/wtot;
-
-       w[0] = w[0]*wtot_inv;
-       w[1] = w[1]*wtot_inv;
-       w[2] = w[2]*wtot_inv;
-   }
-   else /* dummy values for zero area face */
-       w[0] = w[1] = w[2] = 1.0f/3.0f;
-}
-
 /* still use 2D X,Y space but this works for verts transformed by a perspective matrix, using their 4th component as a weight */
-static void BarycentricWeightsPersp2f(float pt[2], float v1[4], float v2[4], float v3[4], float w[3])
+static void barycentric_weights_v2_persp(float v1[4], float v2[4], float v3[4], float co[2], float w[3])
 {
    float wtot_inv, wtot;
 
-   w[0] = AreaSignedF2Dfl(v2, v3, pt) / v1[3];
-   w[1] = AreaSignedF2Dfl(v3, v1, pt) / v2[3];
-   w[2] = AreaSignedF2Dfl(v1, v2, pt) / v3[3];
+   w[0] = area_tri_signed_v2(v2, v3, co) / v1[3];
+   w[1] = area_tri_signed_v2(v3, v1, co) / v2[3];
+   w[2] = area_tri_signed_v2(v1, v2, co) / v3[3];
    wtot = w[0]+w[1]+w[2];
 
    if (wtot != 0.0f) {
@@ -513,13 +487,13 @@ static void BarycentricWeightsPersp2f(float pt[2], float v1[4], float v2[4], flo
 
 static float VecZDepthOrtho(float pt[2], float v1[3], float v2[3], float v3[3], float w[3])
 {
-	BarycentricWeights2f(pt, v1, v2, v3, w);
+	barycentric_weights_v2(v1, v2, v3, pt, w);
 	return (v1[2]*w[0]) + (v2[2]*w[1]) + (v3[2]*w[2]);
 }
 
 static float VecZDepthPersp(float pt[2], float v1[3], float v2[3], float v3[3], float w[3])
 {
-	BarycentricWeightsPersp2f(pt, v1, v2, v3, w);
+	barycentric_weights_v2_persp(v1, v2, v3, pt, w);
 	return (v1[2]*w[0]) + (v2[2]*w[1]) + (v3[2]*w[2]);
 }
 
@@ -738,8 +712,8 @@ static int project_paint_occlude_ptv_clip(
 		return ret;
 
 	if (ret==1) { /* weights not calculated */
-		if (ps->is_ortho)	BarycentricWeights2f(pt, v1, v2, v3, w);
-		else				BarycentricWeightsPersp2f(pt, v1, v2, v3, w);
+		if (ps->is_ortho)	barycentric_weights_v2(v1, v2, v3, pt, w);
+		else				barycentric_weights_v2_persp(v1, v2, v3, pt, w);
 	}
 
 	/* Test if we're in the clipped area, */
@@ -1187,7 +1161,7 @@ static void screen_px_from_ortho(
 		float pixelScreenCo[4],
 		float w[3])
 {
-	BarycentricWeights2f(uv, uv1co, uv2co, uv3co, w);
+	barycentric_weights_v2(uv1co, uv2co, uv3co, uv, w);
 	interp_v3_v3v3v3(pixelScreenCo, v1co, v2co, v3co, w);
 }
 
@@ -1202,7 +1176,7 @@ static void screen_px_from_persp(
 {
 
 	float wtot_inv, wtot;
-	BarycentricWeights2f(uv, uv1co, uv2co, uv3co, w);
+	barycentric_weights_v2(uv1co, uv2co, uv3co, uv, w);
 	
 	/* re-weight from the 4th coord of each screen vert */
 	w[0] *= v1co[3];
@@ -1774,26 +1748,26 @@ static void rect_to_uvspace_ortho(
 	/* get the UV space bounding box */
 	uv[0] = bucket_bounds->xmax;
 	uv[1] = bucket_bounds->ymin;
-	BarycentricWeights2f(uv, v1coSS, v2coSS, v3coSS, w);
+	barycentric_weights_v2(v1coSS, v2coSS, v3coSS, uv, w);
 	interp_v2_v2v2v2(bucket_bounds_uv[flip?3:0], uv1co, uv2co, uv3co, w);
 
 	//uv[0] = bucket_bounds->xmax; // set above
 	uv[1] = bucket_bounds->ymax;
-	BarycentricWeights2f(uv, v1coSS, v2coSS, v3coSS, w);
+	barycentric_weights_v2(v1coSS, v2coSS, v3coSS, uv, w);
 	interp_v2_v2v2v2(bucket_bounds_uv[flip?2:1], uv1co, uv2co, uv3co, w);
 
 	uv[0] = bucket_bounds->xmin;
 	//uv[1] = bucket_bounds->ymax; // set above
-	BarycentricWeights2f(uv, v1coSS, v2coSS, v3coSS, w);
+	barycentric_weights_v2(v1coSS, v2coSS, v3coSS, uv, w);
 	interp_v2_v2v2v2(bucket_bounds_uv[flip?1:2], uv1co, uv2co, uv3co, w);
 
 	//uv[0] = bucket_bounds->xmin; // set above
 	uv[1] = bucket_bounds->ymin;
-	BarycentricWeights2f(uv, v1coSS, v2coSS, v3coSS, w);
+	barycentric_weights_v2(v1coSS, v2coSS, v3coSS, uv, w);
 	interp_v2_v2v2v2(bucket_bounds_uv[flip?0:3], uv1co, uv2co, uv3co, w);
 }
 
-/* same as above but use BarycentricWeightsPersp2f */
+/* same as above but use barycentric_weights_v2_persp */
 static void rect_to_uvspace_persp(
 		rctf *bucket_bounds,
 		float *v1coSS, float *v2coSS, float *v3coSS,
@@ -1808,22 +1782,22 @@ static void rect_to_uvspace_persp(
 	/* get the UV space bounding box */
 	uv[0] = bucket_bounds->xmax;
 	uv[1] = bucket_bounds->ymin;
-	BarycentricWeightsPersp2f(uv, v1coSS, v2coSS, v3coSS, w);
+	barycentric_weights_v2_persp(v1coSS, v2coSS, v3coSS, uv, w);
 	interp_v2_v2v2v2(bucket_bounds_uv[flip?3:0], uv1co, uv2co, uv3co, w);
 
 	//uv[0] = bucket_bounds->xmax; // set above
 	uv[1] = bucket_bounds->ymax;
-	BarycentricWeightsPersp2f(uv, v1coSS, v2coSS, v3coSS, w);
+	barycentric_weights_v2_persp(v1coSS, v2coSS, v3coSS, uv, w);
 	interp_v2_v2v2v2(bucket_bounds_uv[flip?2:1], uv1co, uv2co, uv3co, w);
 
 	uv[0] = bucket_bounds->xmin;
 	//uv[1] = bucket_bounds->ymax; // set above
-	BarycentricWeightsPersp2f(uv, v1coSS, v2coSS, v3coSS, w);
+	barycentric_weights_v2_persp(v1coSS, v2coSS, v3coSS, uv, w);
 	interp_v2_v2v2v2(bucket_bounds_uv[flip?1:2], uv1co, uv2co, uv3co, w);
 
 	//uv[0] = bucket_bounds->xmin; // set above
 	uv[1] = bucket_bounds->ymin;
-	BarycentricWeightsPersp2f(uv, v1coSS, v2coSS, v3coSS, w);
+	barycentric_weights_v2_persp(v1coSS, v2coSS, v3coSS, uv, w);
 	interp_v2_v2v2v2(bucket_bounds_uv[flip?0:3], uv1co, uv2co, uv3co, w);
 }
 
@@ -2067,13 +2041,13 @@ static void project_bucket_clip_face(
 		
 		if (is_ortho) {
 			for(i=0; i<(*tot); i++) {
-				BarycentricWeights2f(isectVCosSS[i], v1coSS, v2coSS, v3coSS, w);
+				barycentric_weights_v2(v1coSS, v2coSS, v3coSS, isectVCosSS[i], w);
 				interp_v2_v2v2v2(bucket_bounds_uv[i], uv1co, uv2co, uv3co, w);
 			}
 		}
 		else {
 			for(i=0; i<(*tot); i++) {
-				BarycentricWeightsPersp2f(isectVCosSS[i], v1coSS, v2coSS, v3coSS, w);
+				barycentric_weights_v2_persp(v1coSS, v2coSS, v3coSS, isectVCosSS[i], w);
 				interp_v2_v2v2v2(bucket_bounds_uv[i], uv1co, uv2co, uv3co, w);
 			}
 		}
@@ -2521,10 +2495,10 @@ static void project_paint_face_init(const ProjPaintState *ps, const int thread_i
 #if 0
 												/* This is not QUITE correct since UV is not inside the UV's but good enough for seams */
 												if (side) {
-													BarycentricWeights2f(uv, tf_uv_pxoffset[0], tf_uv_pxoffset[2], tf_uv_pxoffset[3], w);
+													barycentric_weights_v2(tf_uv_pxoffset[0], tf_uv_pxoffset[2], tf_uv_pxoffset[3], uv, w);
 												}
 												else {
-													BarycentricWeights2f(uv, tf_uv_pxoffset[0], tf_uv_pxoffset[1], tf_uv_pxoffset[2], w);
+													barycentric_weights_v2(tf_uv_pxoffset[0], tf_uv_pxoffset[1], tf_uv_pxoffset[2], uv, w);
 												}
 #endif
 #if 1
diff --git a/source/blender/python/generic/Geometry.c b/source/blender/python/generic/Geometry.c
index 0d59df6ceca..579c9d987cf 100644
--- a/source/blender/python/generic/Geometry.c
+++ b/source/blender/python/generic/Geometry.c
@@ -55,6 +55,7 @@ static PyObject *M_Geometry_PointInTriangle2D( PyObject * self, PyObject * args
 static PyObject *M_Geometry_PointInQuad2D( PyObject * self, PyObject * args );
 static PyObject *M_Geometry_BoxPack2D( PyObject * self, PyObject * args );
 static PyObject *M_Geometry_BezierInterp( PyObject * self, PyObject * args );
+static PyObject *M_Geometry_BarycentricTransform( PyObject * self, PyObject * args );
 
 
 /*-------------------------DOC STRINGS ---------------------------*/
@@ -75,6 +76,7 @@ struct PyMethodDef M_Geometry_methods[] = {
 	{"PointInQuad2D", ( PyCFunction ) M_Geometry_PointInQuad2D, METH_VARARGS, M_Geometry_PointInQuad2D_doc},
 	{"BoxPack2D", ( PyCFunction ) M_Geometry_BoxPack2D, METH_O, M_Geometry_BoxPack2D_doc},
 	{"BezierInterp", ( PyCFunction ) M_Geometry_BezierInterp, METH_VARARGS, M_Geometry_BezierInterp_doc},
+	{"BarycentricTransform", ( PyCFunction ) M_Geometry_BarycentricTransform, METH_VARARGS, NULL},
 	{NULL, NULL, 0, NULL}
 };
 
@@ -534,3 +536,36 @@ static PyObject *M_Geometry_BezierInterp( PyObject * self, PyObject * args )
 	MEM_freeN(coord_array);
 	return list;
 }
+
+static PyObject *M_Geometry_BarycentricTransform(PyObject * self, PyObject * args)
+{
+	VectorObject *vec_pt;
+	VectorObject *vec_t1_tar, *vec_t2_tar, *vec_t3_tar;
+	VectorObject *vec_t1_src, *vec_t2_src, *vec_t3_src;
+	float vec[3];
+
+	if( !PyArg_ParseTuple ( args, "O!O!O!O!O!O!O!",
+	  &vector_Type, &vec_pt,
+	  &vector_Type, &vec_t1_src,
+	  &vector_Type, &vec_t2_src,
+	  &vector_Type, &vec_t3_src,
+	  &vector_Type, &vec_t1_tar,
+	  &vector_Type, &vec_t2_tar,
+	  &vector_Type, &vec_t3_tar) || (	vec_pt->size != 3 ||
+										vec_t1_src->size != 3 ||
+										vec_t2_src->size != 3 ||
+										vec_t3_src->size != 3 ||
+										vec_t1_tar->size != 3 ||
+										vec_t2_tar->size != 3 ||
+										vec_t3_tar->size != 3)
+	) {
+		PyErr_SetString( PyExc_TypeError, "expected 7, 3D vector types\n" );
+		return NULL;
+	}
+
+	barycentric_transform(vec, vec_pt->vec,
+			vec_t1_tar->vec, vec_t2_tar->vec, vec_t3_tar->vec,
+			vec_t1_src->vec, vec_t2_src->vec, vec_t3_src->vec);
+
+	return newVectorObject(vec, 3, Py_NEW, NULL);
+}