code cleanup: favor braces when blocks have mixed brace use.
[blender.git] / source / blender / blenlib / intern / math_geom.c
index 10a9a2df9f5092cd5e73c438482163c2b9baee2b..49be60dda36177bdaf68ea60cd9a3c8711da4f15 100644 (file)
@@ -39,9 +39,9 @@
 
 void cent_tri_v3(float cent[3], const float v1[3], const float v2[3], const float v3[3])
 {
 
 void cent_tri_v3(float cent[3], const float v1[3], const float v2[3], const float v3[3])
 {
-       cent[0] = 0.33333f * (v1[0] + v2[0] + v3[0]);
-       cent[1] = 0.33333f * (v1[1] + v2[1] + v3[1]);
-       cent[2] = 0.33333f * (v1[2] + v2[2] + v3[2]);
+       cent[0] = (v1[0] + v2[0] + v3[0]) / 3.0f;
+       cent[1] = (v1[1] + v2[1] + v3[1]) / 3.0f;
+       cent[2] = (v1[2] + v2[2] + v3[2]) / 3.0f;
 }
 
 void cent_quad_v3(float cent[3], const float v1[3], const float v2[3], const float v3[3], const float v4[3])
 }
 
 void cent_quad_v3(float cent[3], const float v1[3], const float v2[3], const float v3[3], const float v4[3])
@@ -106,12 +106,12 @@ float area_quad_v3(const float v1[3], const float v2[3], const float v3[3], cons
        sub_v3_v3v3(vec1, v2, v1);
        sub_v3_v3v3(vec2, v4, v1);
        cross_v3_v3v3(n, vec1, vec2);
        sub_v3_v3v3(vec1, v2, v1);
        sub_v3_v3v3(vec2, v4, v1);
        cross_v3_v3v3(n, vec1, vec2);
-       len = normalize_v3(n);
+       len = len_v3(n);
 
        sub_v3_v3v3(vec1, v4, v3);
        sub_v3_v3v3(vec2, v2, v3);
        cross_v3_v3v3(n, vec1, vec2);
 
        sub_v3_v3v3(vec1, v4, v3);
        sub_v3_v3v3(vec2, v2, v3);
        cross_v3_v3v3(n, vec1, vec2);
-       len += normalize_v3(n);
+       len += len_v3(n);
 
        return (len / 2.0f);
 }
 
        return (len / 2.0f);
 }
@@ -119,45 +119,68 @@ float area_quad_v3(const float v1[3], const float v2[3], const float v3[3], cons
 /* Triangles */
 float area_tri_v3(const float v1[3], const float v2[3], const float v3[3])
 {
 /* Triangles */
 float area_tri_v3(const float v1[3], const float v2[3], const float v3[3])
 {
-       float len, vec1[3], vec2[3], n[3];
+       float vec1[3], vec2[3], n[3];
 
        sub_v3_v3v3(vec1, v3, v2);
        sub_v3_v3v3(vec2, v1, v2);
        cross_v3_v3v3(n, vec1, vec2);
 
        sub_v3_v3v3(vec1, v3, v2);
        sub_v3_v3v3(vec2, v1, v2);
        cross_v3_v3v3(n, vec1, vec2);
-       len = normalize_v3(n);
 
 
-       return (len / 2.0f);
+       return len_v3(n) / 2.0f;
+}
+
+float area_tri_signed_v3(const float v1[3], const float v2[3], const float v3[3], const float normal[3])
+{
+       float area, vec1[3], vec2[3], n[3];
+
+       sub_v3_v3v3(vec1, v3, v2);
+       sub_v3_v3v3(vec2, v1, v2);
+       cross_v3_v3v3(n, vec1, vec2);
+       area = len_v3(n) / 2.0f;
+
+       /* negate area for flipped triangles */
+       if (dot_v3v3(n, normal) < 0.0f)
+               area = -area;
+
+       return area;
 }
 
 float area_poly_v3(int nr, float verts[][3], const float normal[3])
 {
 }
 
 float area_poly_v3(int nr, float verts[][3], const float normal[3])
 {
-       float x, y, z, area, max;
-       float *cur, *prev;
-       int a, px = 0, py = 1;
+       int a, px, py;
+       const float max = axis_dominant_v3_max(&px, &py, normal);
+       float area;
+       float *co_curr, *co_prev;
 
 
-       /* first: find dominant axis: 0==X, 1==Y, 2==Z
-        * don't use 'axis_dominant_v3()' because we need max axis too */
-       x = fabsf(normal[0]);
-       y = fabsf(normal[1]);
-       z = fabsf(normal[2]);
-       max = max_fff(x, y, z);
-       if (max == y) py = 2;
-       else if (max == x) {
-               px = 1;
-               py = 2;
+       /* The Trapezium Area Rule */
+       co_prev = verts[nr - 1];
+       co_curr = verts[0];
+       area = 0.0f;
+       for (a = 0; a < nr; a++) {
+               area += (co_curr[px] - co_prev[px]) * (co_curr[py] + co_prev[py]);
+               co_prev = verts[a];
+               co_curr = verts[a + 1];
        }
 
        }
 
+       return fabsf(0.5f * area / max);
+}
+
+float area_poly_v2(int nr, float verts[][2])
+{
+       int a;
+       float area;
+       float *co_curr, *co_prev;
+
        /* The Trapezium Area Rule */
        /* The Trapezium Area Rule */
-       prev = verts[nr - 1];
-       cur = verts[0];
-       area = 0;
+       co_prev = verts[nr - 1];
+       co_curr = verts[0];
+       area = 0.0f;
        for (a = 0; a < nr; a++) {
        for (a = 0; a < nr; a++) {
-               area += (cur[px] - prev[px]) * (cur[py] + prev[py]);
-               prev = verts[a];
-               cur = verts[a + 1];
+               area += (co_curr[0] - co_prev[0]) * (co_curr[1] + co_prev[1]);
+               co_prev = verts[a];
+               co_curr = verts[a + 1];
        }
 
        }
 
-       return fabsf(0.5f * area / max);
+       return fabsf(0.5f * area);
 }
 
 /********************************* Distance **********************************/
 }
 
 /********************************* Distance **********************************/
@@ -296,6 +319,97 @@ float dist_to_line_segment_v3(const float v1[3], const float v2[3], const float
        return len_v3v3(closest, v1);
 }
 
        return len_v3v3(closest, v1);
 }
 
+float dist_to_line_v3(const float v1[3], const float v2[3], const float v3[3])
+{
+       float closest[3];
+
+       closest_to_line_v3(closest, v1, v2, v3);
+
+       return len_v3v3(closest, v1);
+}
+
+/* Adapted from "Real-Time Collision Detection" by Christer Ericson,
+ * published by Morgan Kaufmann Publishers, copyright 2005 Elsevier Inc.
+ * 
+ * Set 'r' to the point in triangle (a, b, c) closest to point 'p' */
+void closest_on_tri_to_point_v3(float r[3], const float p[3],
+                                const float a[3], const float b[3], const float c[3])
+{
+       float ab[3], ac[3], ap[3], d1, d2;
+       float bp[3], d3, d4, vc, cp[3], d5, d6, vb, va;
+       float denom, v, w;
+
+       /* Check if P in vertex region outside A */
+       sub_v3_v3v3(ab, b, a);
+       sub_v3_v3v3(ac, c, a);
+       sub_v3_v3v3(ap, p, a);
+       d1 = dot_v3v3(ab, ap);
+       d2 = dot_v3v3(ac, ap);
+       if (d1 <= 0.0f && d2 <= 0.0f) {
+               /* barycentric coordinates (1,0,0) */
+               copy_v3_v3(r, a);
+               return;
+       }
+
+       /* Check if P in vertex region outside B */
+       sub_v3_v3v3(bp, p, b);
+       d3 = dot_v3v3(ab, bp);
+       d4 = dot_v3v3(ac, bp);
+       if (d3 >= 0.0f && d4 <= d3) {
+               /* barycentric coordinates (0,1,0) */
+               copy_v3_v3(r, b);
+               return;
+       }
+       /* Check if P in edge region of AB, if so return projection of P onto AB */
+       vc = d1 * d4 - d3 * d2;
+       if (vc <= 0.0f && d1 >= 0.0f && d3 <= 0.0f) {
+               float v = d1 / (d1 - d3);
+               /* barycentric coordinates (1-v,v,0) */
+               madd_v3_v3v3fl(r, a, ab, v);
+               return;
+       }
+       /* Check if P in vertex region outside C */
+       sub_v3_v3v3(cp, p, c);
+       d5 = dot_v3v3(ab, cp);
+       d6 = dot_v3v3(ac, cp);
+       if (d6 >= 0.0f && d5 <= d6) {
+               /* barycentric coordinates (0,0,1) */
+               copy_v3_v3(r, c);
+               return;
+       }
+       /* Check if P in edge region of AC, if so return projection of P onto AC */
+       vb = d5 * d2 - d1 * d6;
+       if (vb <= 0.0f && d2 >= 0.0f && d6 <= 0.0f) {
+               float w = d2 / (d2 - d6);
+               /* barycentric coordinates (1-w,0,w) */
+               madd_v3_v3v3fl(r, a, ac, w);
+               return;
+       }
+       /* Check if P in edge region of BC, if so return projection of P onto BC */
+       va = d3 * d6 - d5 * d4;
+       if (va <= 0.0f && (d4 - d3) >= 0.0f && (d5 - d6) >= 0.0f) {
+               float w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
+               /* barycentric coordinates (0,1-w,w) */
+               sub_v3_v3v3(r, c, b);
+               mul_v3_fl(r, w);
+               add_v3_v3(r, b);
+               return;
+       }
+
+       /* P inside face region. Compute Q through its barycentric coordinates (u,v,w) */
+       denom = 1.0f / (va + vb + vc);
+       v = vb * denom;
+       w = vc * denom;
+
+       /* = u*a + v*b + w*c, u = va * denom = 1.0f - v - w */
+       /* ac * w */
+       mul_v3_fl(ac, w);
+       /* a + ab * v */
+       madd_v3_v3v3fl(r, a, ab, v);
+       /* a + ab * v + ac * w */
+       add_v3_v3(r, ac);
+}
+
 /******************************* Intersection ********************************/
 
 /* intersect Line-Line, shorts */
 /******************************* Intersection ********************************/
 
 /* intersect Line-Line, shorts */
@@ -1037,7 +1151,7 @@ int isect_sweeping_sphere_tri_v3(const float p1[3], const float p2[3], const flo
 
                if (t0 > 1.0f || t1 < 0.0f) return 0;
 
 
                if (t0 > 1.0f || t1 < 0.0f) return 0;
 
-               /* clamp to [0,1] */
+               /* clamp to [0, 1] */
                CLAMP(t0, 0.0f, 1.0f);
                CLAMP(t1, 0.0f, 1.0f);
 
                CLAMP(t0, 0.0f, 1.0f);
                CLAMP(t1, 0.0f, 1.0f);
 
@@ -1157,10 +1271,10 @@ int isect_sweeping_sphere_tri_v3(const float p1[3], const float p2[3], const flo
        }
 
        /*e3*/
        }
 
        /*e3*/
-       /* sub_v3_v3v3(bv,v0,p1); */ /* UNUSED */
-       /* elen2 = dot_v3v3(e1,e1); */ /* UNUSED */
-       /* edotv = dot_v3v3(e1,vel); */ /* UNUSED */
-       /* edotbv = dot_v3v3(e1,bv); */ /* UNUSED */
+       /* sub_v3_v3v3(bv, v0, p1); */ /* UNUSED */
+       /* elen2 = dot_v3v3(e1, e1); */ /* UNUSED */
+       /* edotv = dot_v3v3(e1, vel); */ /* UNUSED */
+       /* edotbv = dot_v3v3(e1, bv); */ /* UNUSED */
 
        sub_v3_v3v3(bv, v1, p1);
        elen2 = dot_v3v3(e3, e3);
 
        sub_v3_v3v3(bv, v1, p1);
        elen2 = dot_v3v3(e3, e3);
@@ -1195,7 +1309,7 @@ int isect_axial_line_tri_v3(const int axis, const float p1[3], const float p2[3]
        int a0 = axis, a1 = (axis + 1) % 3, a2 = (axis + 2) % 3;
 
 #if 0
        int a0 = axis, a1 = (axis + 1) % 3, a2 = (axis + 2) % 3;
 
 #if 0
-       return isect_line_tri_v3(p1,p2,v0,v1,v2,lambda);
+       return isect_line_tri_v3(p1, p2, v0, v1, v2, lambda);
 
        /* first a simple bounding box test */
        if (min_fff(v0[a1], v1[a1], v2[a1]) > p1[a1]) return 0;
 
        /* first a simple bounding box test */
        if (min_fff(v0[a1], v1[a1], v2[a1]) > p1[a1]) return 0;
@@ -1415,8 +1529,8 @@ int isect_ray_aabb(const IsectRayAABBData *data, const float bb_min[3],
        return TRUE;
 }
 
        return TRUE;
 }
 
-/* find closest point to p on line through l1,l2 and return lambda,
- * where (0 <= lambda <= 1) when cp is in the line segment l1,l2
+/* find closest point to p on line through (l1, l2) and return lambda,
+ * where (0 <= lambda <= 1) when cp is in the line segment (l1, l2)
  */
 float closest_to_line_v3(float cp[3], const float p[3], const float l1[3], const float l2[3])
 {
  */
 float closest_to_line_v3(float cp[3], const float p[3], const float l1[3], const float l2[3])
 {
@@ -1702,9 +1816,9 @@ static int point_in_slice(const float p[3], const float v1[3], const float l1[3]
        /*
         * what is a slice ?
         * some maths:
        /*
         * what is a slice ?
         * some maths:
-        * a line including l1,l2 and a point not on the line
+        * a line including (l1, l2) and a point not on the line
         * define a subset of R3 delimited by planes parallel to the line and orthogonal
         * define a subset of R3 delimited by planes parallel to the line and orthogonal
-        * to the (point --> line) distance vector,one plane on the line one on the point,
+        * to the (point --> line) distance vector, one plane on the line one on the point,
         * the room inside usually is rather small compared to R3 though still infinite
         * useful for restricting (speeding up) searches
         * e.g. all points of triangular prism are within the intersection of 3 'slices'
         * the room inside usually is rather small compared to R3 though still infinite
         * useful for restricting (speeding up) searches
         * e.g. all points of triangular prism are within the intersection of 3 'slices'
@@ -1735,7 +1849,7 @@ static int point_in_slice_as(float p[3], float origin[3], float normal[3])
        return 1;
 }
 
        return 1;
 }
 
-/*mama (knowing the squared length of the normal)*/
+/*mama (knowing the squared length of the normal) */
 static int point_in_slice_m(float p[3], float origin[3], float normal[3], float lns)
 {
        float h, rp[3];
 static int point_in_slice_m(float p[3], float origin[3], float normal[3], float lns)
 {
        float h, rp[3];
@@ -1867,20 +1981,73 @@ void plot_line_v2v2i(const int p1[2], const int p2[2], int (*callback)(int, int,
        }
 }
 
        }
 }
 
-/****************************** Interpolation ********************************/
+/****************************** Axis Utils ********************************/
+
+/**
+ * \brief Normal to x,y matrix
+ *
+ * Creates a 3x3 matrix from a normal.
+ * This matrix can be applied to vectors so their 'z' axis runs along \a normal.
+ * In practice it means you can use x,y as 2d coords. \see
+ *
+ * \param r_mat The matrix to return.
+ * \param normal A unit length vector.
+ */
+bool axis_dominant_v3_to_m3(float r_mat[3][3], const float normal[3])
+{
+       float up[3] = {0.0f, 0.0f, 1.0f};
+       float axis[3];
+       float angle;
+
+       /* double check they are normalized */
+       BLI_ASSERT_UNIT_V3(normal);
+
+       cross_v3_v3v3(axis, normal, up);
+       angle = saacos(dot_v3v3(normal, up));
+
+       if (angle >= FLT_EPSILON) {
+               if (len_squared_v3(axis) < FLT_EPSILON) {
+                       axis[0] = 0.0f;
+                       axis[1] = 1.0f;
+                       axis[2] = 0.0f;
+               }
+
+               axis_angle_to_mat3(r_mat, axis, angle);
+               return true;
+       }
+       else {
+               unit_m3(r_mat);
+               return false;
+       }
+}
 
 /* get the 2 dominant axis values, 0==X, 1==Y, 2==Z */
 
 /* get the 2 dominant axis values, 0==X, 1==Y, 2==Z */
-void axis_dominant_v3(int *axis_a, int *axis_b, const float axis[3])
+void axis_dominant_v3(int *r_axis_a, int *r_axis_b, const float axis[3])
+{
+       const float xn = fabsf(axis[0]);
+       const float yn = fabsf(axis[1]);
+       const float zn = fabsf(axis[2]);
+
+       if      (zn >= xn && zn >= yn) { *r_axis_a = 0; *r_axis_b = 1; }
+       else if (yn >= xn && yn >= zn) { *r_axis_a = 0; *r_axis_b = 2; }
+       else                           { *r_axis_a = 1; *r_axis_b = 2; }
+}
+
+/* same as axis_dominant_v3 but return the max value */
+float axis_dominant_v3_max(int *r_axis_a, int *r_axis_b, const float axis[3])
 {
        const float xn = fabsf(axis[0]);
        const float yn = fabsf(axis[1]);
        const float zn = fabsf(axis[2]);
 
 {
        const float xn = fabsf(axis[0]);
        const float yn = fabsf(axis[1]);
        const float zn = fabsf(axis[2]);
 
-       if      (zn >= xn && zn >= yn) { *axis_a = 0; *axis_b = 1; }
-       else if (yn >= xn && yn >= zn) { *axis_a = 0; *axis_b = 2; }
-       else                           { *axis_a = 1; *axis_b = 2; }
+       if      (zn >= xn && zn >= yn) { *r_axis_a = 0; *r_axis_b = 1; return zn; }
+       else if (yn >= xn && yn >= zn) { *r_axis_a = 0; *r_axis_b = 2; return yn; }
+       else                           { *r_axis_a = 1; *r_axis_b = 2; return xn; }
 }
 
 }
 
+
+/****************************** Interpolation ********************************/
+
 static float tri_signed_area(const float v1[3], const float v2[3], const float v3[3], const int i, const int j)
 {
        return 0.5f * ((v1[i] - v2[i]) * (v2[j] - v3[j]) + (v1[j] - v2[j]) * (v3[i] - v2[i]));
 static float tri_signed_area(const float v1[3], const float v2[3], const float v3[3], const int i, const int j)
 {
        return 0.5f * ((v1[i] - v2[i]) * (v2[j] - v3[j]) + (v1[j] - v2[j]) * (v3[i] - v2[i]));
@@ -2252,59 +2419,106 @@ void interp_weights_poly_v3(float *w, float v[][3], const int n, const float co[
 {
        /* TODO: t1 and t2 overlap each iter, we could call this only once per iter and reuse previous value */
        float totweight, t1, t2, len, *vmid, *vprev, *vnext;
 {
        /* TODO: t1 and t2 overlap each iter, we could call this only once per iter and reuse previous value */
        float totweight, t1, t2, len, *vmid, *vprev, *vnext;
-       int i;
+       int i, i_next, i_curr;
+       bool edge_interp = false;
 
        totweight = 0.0f;
 
        for (i = 0; i < n; i++) {
 
        totweight = 0.0f;
 
        for (i = 0; i < n; i++) {
+               i_curr = i;
+               i_next = (i == n - 1) ? 0 : i + 1;
+
                vmid = v[i];
                vprev = (i == 0) ? v[n - 1] : v[i - 1];
                vmid = v[i];
                vprev = (i == 0) ? v[n - 1] : v[i - 1];
-               vnext = (i == n - 1) ? v[0] : v[i + 1];
+               vnext = v[i_next];
+
+               /* Mark Mayer et al algorithm that is used here does not operate well if vertex is close
+                * to borders of face. In that case, do simple linear interpolation between the two edge vertices */
+               if (dist_to_line_segment_v3(co, vmid, vnext) < 10 * FLT_EPSILON) {
+                       edge_interp = true;
+                       break;
+               }
 
                t1 = mean_value_half_tan_v3(co, vprev, vmid);
                t2 = mean_value_half_tan_v3(co, vmid, vnext);
 
                len = len_v3v3(co, vmid);
 
                t1 = mean_value_half_tan_v3(co, vprev, vmid);
                t2 = mean_value_half_tan_v3(co, vmid, vnext);
 
                len = len_v3v3(co, vmid);
-               w[i] = (len != 0.0f)? (t1 + t2) / len: 0.0f;
+               w[i] = (len != 0.0f) ? (t1 + t2) / len: 0.0f;
                totweight += w[i];
        }
 
                totweight += w[i];
        }
 
-       if (totweight != 0.0f) {
-               for (i = 0; i < n; i++) {
-                       w[i] /= totweight;
+       if (edge_interp) {
+               float len_curr = len_v3v3(co, vmid);
+               float len_next = len_v3v3(co, vnext);
+               float edge_len = len_curr + len_next;
+               for (i = 0; i < n; i++)
+                       w[i] = 0.0;
+
+               w[i_curr] = len_next / edge_len;
+               w[i_next] = len_curr / edge_len;
+       }
+       else {
+               if (totweight != 0.0f) {
+                       for (i = 0; i < n; i++) {
+                               w[i] /= totweight;
+                       }
                }
        }
 }
 
                }
        }
 }
 
+
 void interp_weights_poly_v2(float *w, float v[][2], const int n, const float co[2])
 {
        /* TODO: t1 and t2 overlap each iter, we could call this only once per iter and reuse previous value */
        float totweight, t1, t2, len, *vmid, *vprev, *vnext;
 void interp_weights_poly_v2(float *w, float v[][2], const int n, const float co[2])
 {
        /* TODO: t1 and t2 overlap each iter, we could call this only once per iter and reuse previous value */
        float totweight, t1, t2, len, *vmid, *vprev, *vnext;
-       int i;
+       int i, i_next, i_curr;
+       bool edge_interp = false;
 
        totweight = 0.0f;
 
        for (i = 0; i < n; i++) {
 
        totweight = 0.0f;
 
        for (i = 0; i < n; i++) {
+               i_curr = i;
+               i_next = (i == n - 1) ? 0 : i + 1;
+
                vmid = v[i];
                vprev = (i == 0) ? v[n - 1] : v[i - 1];
                vmid = v[i];
                vprev = (i == 0) ? v[n - 1] : v[i - 1];
-               vnext = (i == n - 1) ? v[0] : v[i + 1];
+               vnext = v[i_next];
+
+               /* Mark Mayer et al algorithm that is used here does not operate well if vertex is close
+                * to borders of face. In that case, do simple linear interpolation between the two edge vertices */
+               if (dist_to_line_segment_v2(co, vmid, vnext) < 10 * FLT_EPSILON) {
+                       edge_interp = true;
+                       break;
+               }
 
                t1 = mean_value_half_tan_v2(co, vprev, vmid);
                t2 = mean_value_half_tan_v2(co, vmid, vnext);
 
                len = len_v2v2(co, vmid);
 
                t1 = mean_value_half_tan_v2(co, vprev, vmid);
                t2 = mean_value_half_tan_v2(co, vmid, vnext);
 
                len = len_v2v2(co, vmid);
-               w[i] = (len != 0.0f)? (t1 + t2) / len: 0.0f;
+               w[i] = (len != 0.0f) ? (t1 + t2) / len: 0.0f;
                totweight += w[i];
        }
 
                totweight += w[i];
        }
 
-       if (totweight != 0.0f) {
-               for (i = 0; i < n; i++) {
-                       w[i] /= totweight;
+       if (edge_interp) {
+               float len_curr = len_v2v2(co, vmid);
+               float len_next = len_v2v2(co, vnext);
+               float edge_len = len_curr + len_next;
+               for (i = 0; i < n; i++)
+                       w[i] = 0.0;
+
+               w[i_curr] = len_next / edge_len;
+               w[i_next] = len_curr / edge_len;
+       }
+       else {
+               if (totweight != 0.0f) {
+                       for (i = 0; i < n; i++) {
+                               w[i] /= totweight;
+                       }
                }
        }
 }
 
                }
        }
 }
 
-/* (x1,v1)(t1=0)------(x2,v2)(t2=1), 0<t<1 --> (x,v)(t) */
+/* (x1, v1)(t1=0)------(x2, v2)(t2=1), 0<t<1 --> (x, v)(t) */
 void interp_cubic_v3(float x[3], float v[3], const float x1[3], const float v1[3], const float x2[3], const float v2[3], const float t)
 {
        float a[3], b[3];
 void interp_cubic_v3(float x[3], float v[3], const float x1[3], const float v1[3], const float x2[3], const float v2[3], const float t)
 {
        float a[3], b[3];
@@ -2349,7 +2563,9 @@ void resolve_tri_uv(float r_uv[2], const float st[2], const float st0[2], const
                r_uv[0] = (float)((d * x[0] - b * x[1]) / det);
                r_uv[1] = (float)(((-c) * x[0] + a * x[1]) / det);
        }
                r_uv[0] = (float)((d * x[0] - b * x[1]) / det);
                r_uv[1] = (float)(((-c) * x[0] + a * x[1]) / det);
        }
-       else zero_v2(r_uv);
+       else {
+               zero_v2(r_uv);
+       }
 }
 
 /* bilinear reverse */
 }
 
 /* bilinear reverse */
@@ -2791,8 +3007,8 @@ void tangent_from_uv(float uv1[2], float uv2[2], float uv3[3], float co1[3], flo
 /****************************** Vector Clouds ********************************/
 
 /* vector clouds */
 /****************************** Vector Clouds ********************************/
 
 /* vector clouds */
-/* void vcloud_estimate_transform(int list_size, float (*pos)[3], float *weight,float (*rpos)[3], float *rweight,
- *                                float lloc[3],float rloc[3],float lrot[3][3],float lscale[3][3])
+/* void vcloud_estimate_transform(int list_size, float (*pos)[3], float *weight, float (*rpos)[3], float *rweight,
+ *                                float lloc[3], float rloc[3], float lrot[3][3], float lscale[3][3])
  *
  * input
  * (
  *
  * input
  * (
@@ -2848,7 +3064,9 @@ void vcloud_estimate_transform(int list_size, float (*pos)[3], float *weight, fl
                                add_v3_v3(accu_com, v);
                                accu_weight += weight[a];
                        }
                                add_v3_v3(accu_com, v);
                                accu_weight += weight[a];
                        }
-                       else add_v3_v3(accu_com, pos[a]);
+                       else {
+                               add_v3_v3(accu_com, pos[a]);
+                       }
 
                        if (rweight) {
                                float v[3];
 
                        if (rweight) {
                                float v[3];
@@ -2857,8 +3075,9 @@ void vcloud_estimate_transform(int list_size, float (*pos)[3], float *weight, fl
                                add_v3_v3(accu_rcom, v);
                                accu_rweight += rweight[a];
                        }
                                add_v3_v3(accu_rcom, v);
                                accu_rweight += rweight[a];
                        }
-                       else add_v3_v3(accu_rcom, rpos[a]);
-
+                       else {
+                               add_v3_v3(accu_rcom, rpos[a]);
+                       }
                }
                if (!weight || !rweight) {
                        accu_weight = accu_rweight = list_size;
                }
                if (!weight || !rweight) {
                        accu_weight = accu_rweight = list_size;
@@ -2881,9 +3100,9 @@ void vcloud_estimate_transform(int list_size, float (*pos)[3], float *weight, fl
                        /* build 'projection' matrix */
                        for (a = 0; a < list_size; a++) {
                                sub_v3_v3v3(va, rpos[a], accu_rcom);
                        /* build 'projection' matrix */
                        for (a = 0; a < list_size; a++) {
                                sub_v3_v3v3(va, rpos[a], accu_rcom);
-                               /* mul_v3_fl(va,bp->mass);  mass needs renormalzation here ?? */
+                               /* mul_v3_fl(va, bp->mass);  mass needs renormalzation here ?? */
                                sub_v3_v3v3(vb, pos[a], accu_com);
                                sub_v3_v3v3(vb, pos[a], accu_com);
-                               /* mul_v3_fl(va,rp->mass); */
+                               /* mul_v3_fl(va, rp->mass); */
                                m[0][0] += va[0] * vb[0];
                                m[0][1] += va[0] * vb[1];
                                m[0][2] += va[0] * vb[2];
                                m[0][0] += va[0] * vb[0];
                                m[0][1] += va[0] * vb[1];
                                m[0][2] += va[0] * vb[2];
@@ -2956,9 +3175,9 @@ static void vec_add_dir(float r[3], const float v1[3], const float v2[3], const
        r[2] = v1[2] + fac * (v2[2] - v1[2]);
 }
 
        r[2] = v1[2] + fac * (v2[2] - v1[2]);
 }
 
-static int ff_visible_quad(const float p[3], const float n[3],
-                           const float v0[3], const float v1[3], const float v2[3],
-                           float q0[3], float q1[3], float q2[3], float q3[3])
+int form_factor_visible_quad(const float p[3], const float n[3],
+                             const float v0[3], const float v1[3], const float v2[3],
+                             float q0[3], float q1[3], float q2[3], float q3[3])
 {
        static const float epsilon = 1e-6f;
        float c, sd[3];
 {
        static const float epsilon = 1e-6f;
        float c, sd[3];
@@ -3317,8 +3536,8 @@ static void ff_normalize(float n[3])
        }
 }
 
        }
 }
 
-static float ff_quad_form_factor(const float p[3], const float n[3],
-                                 const float q0[3], const float q1[3], const float q2[3], const float q3[3])
+float form_factor_quad(const float p[3], const float n[3],
+                       const float q0[3], const float q1[3], const float q2[3], const float q3[3])
 {
        float r0[3], r1[3], r2[3], r3[3], g0[3], g1[3], g2[3], g3[3];
        float a1, a2, a3, a4, dot1, dot2, dot3, dot4, result;
 {
        float r0[3], r1[3], r2[3], r3[3], g0[3], g1[3], g2[3], g3[3];
        float a1, a2, a3, a4, dot1, dot2, dot3, dot4, result;
@@ -3364,11 +3583,11 @@ float form_factor_hemi_poly(float p[3], float n[3], float v1[3], float v2[3], fl
         * covered by a quad or triangle, cosine weighted */
        float q0[3], q1[3], q2[3], q3[3], contrib = 0.0f;
 
         * covered by a quad or triangle, cosine weighted */
        float q0[3], q1[3], q2[3], q3[3], contrib = 0.0f;
 
-       if (ff_visible_quad(p, n, v1, v2, v3, q0, q1, q2, q3))
-               contrib += ff_quad_form_factor(p, n, q0, q1, q2, q3);
+       if (form_factor_visible_quad(p, n, v1, v2, v3, q0, q1, q2, q3))
+               contrib += form_factor_quad(p, n, q0, q1, q2, q3);
 
 
-       if (v4 && ff_visible_quad(p, n, v1, v3, v4, q0, q1, q2, q3))
-               contrib += ff_quad_form_factor(p, n, q0, q1, q2, q3);
+       if (v4 && form_factor_visible_quad(p, n, v1, v3, v4, q0, q1, q2, q3))
+               contrib += form_factor_quad(p, n, q0, q1, q2, q3);
 
        return contrib;
 }
 
        return contrib;
 }
@@ -3376,44 +3595,57 @@ float form_factor_hemi_poly(float p[3], float n[3], float v1[3], float v2[3], fl
 /* evaluate if entire quad is a proper convex quad */
 int is_quad_convex_v3(const float v1[3], const float v2[3], const float v3[3], const float v4[3])
 {
 /* evaluate if entire quad is a proper convex quad */
 int is_quad_convex_v3(const float v1[3], const float v2[3], const float v3[3], const float v4[3])
 {
-       float nor[3], nor1[3], nor2[3], vec[4][2];
-       int axis_a, axis_b;
+       float nor[3], nor_a[3], nor_b[3], vec[4][2];
+       float mat[3][3];
+       const bool is_ok_a = (normal_tri_v3(nor_a, v1, v2, v3) > FLT_EPSILON);
+       const bool is_ok_b = (normal_tri_v3(nor_b, v1, v3, v4) > FLT_EPSILON);
 
        /* define projection, do both trias apart, quad is undefined! */
 
 
        /* define projection, do both trias apart, quad is undefined! */
 
-       normal_tri_v3(nor1, v1, v2, v3);
-       normal_tri_v3(nor2, v1, v3, v4);
+       /* check normal length incase one size is zero area */
+       if (is_ok_a) {
+               if (is_ok_b) {
+                       /* use both, most common outcome */
+
+                       /* when the face is folded over as 2 tris we probably don't want to create
+                        * a quad from it, but go ahead with the intersection test since this
+                        * isn't a function for degenerate faces */
+                       if (UNLIKELY(dot_v3v3(nor_a, nor_b) < 0.0f)) {
+                               /* flip so adding normals in the opposite direction
+                                * doesn't give a zero length vector */
+                               negate_v3(nor_b);
+                       }
 
 
-       /* when the face is folded over as 2 tris we probably don't want to create
-        * a quad from it, but go ahead with the intersection test since this
-        * isn't a function for degenerate faces */
-       if (UNLIKELY(dot_v3v3(nor1, nor2) < 0.0f)) {
-               /* flip so adding normals in the opposite direction
-                * doesnt give a zero length vector */
-               negate_v3(nor2);
+                       add_v3_v3v3(nor, nor_a, nor_b);
+                       normalize_v3(nor);
+               }
+               else {
+                       copy_v3_v3(nor, nor_a);  /* only 'a' */
+               }
+       }
+       else {
+               if (is_ok_b) {
+                       copy_v3_v3(nor, nor_b);  /* only 'b' */
+               }
+               else {
+                       return false;  /* both zero, we can't do anything useful here */
+               }
        }
 
        }
 
-       add_v3_v3v3(nor, nor1, nor2);
-
-       axis_dominant_v3(&axis_a, &axis_b, nor);
-
-       vec[0][0] = v1[axis_a];
-       vec[0][1] = v1[axis_b];
-       vec[1][0] = v2[axis_a];
-       vec[1][1] = v2[axis_b];
+       axis_dominant_v3_to_m3(mat, nor);
 
 
-       vec[2][0] = v3[axis_a];
-       vec[2][1] = v3[axis_b];
-       vec[3][0] = v4[axis_a];
-       vec[3][1] = v4[axis_b];
+       mul_v2_m3v3(vec[0], mat, v1);
+       mul_v2_m3v3(vec[1], mat, v2);
+       mul_v2_m3v3(vec[2], mat, v3);
+       mul_v2_m3v3(vec[3], mat, v4);
 
        /* linetests, the 2 diagonals have to instersect to be convex */
 
        /* linetests, the 2 diagonals have to instersect to be convex */
-       return (isect_line_line_v2(vec[0], vec[2], vec[1], vec[3]) > 0) ? TRUE : FALSE;
+       return (isect_line_line_v2(vec[0], vec[2], vec[1], vec[3]) > 0);
 }
 
 int is_quad_convex_v2(const float v1[2], const float v2[2], const float v3[2], const float v4[2])
 {
        /* linetests, the 2 diagonals have to instersect to be convex */
 }
 
 int is_quad_convex_v2(const float v1[2], const float v2[2], const float v3[2], const float v4[2])
 {
        /* linetests, the 2 diagonals have to instersect to be convex */
-       return (isect_line_line_v2(v1, v3, v2, v4) > 0) ? TRUE : FALSE;
+       return (isect_line_line_v2(v1, v3, v2, v4) > 0);
 }
 
 }