Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
path: root/source/blender/blenlib/intern/math_geom.c
diff options
context:
space:
mode:
Diffstat (limited to 'source/blender/blenlib/intern/math_geom.c')
-rw-r--r--source/blender/blenlib/intern/math_geom.c533
1 files changed, 428 insertions, 105 deletions
diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c
index 931025606db..579f397f833 100644
--- a/source/blender/blenlib/intern/math_geom.c
+++ b/source/blender/blenlib/intern/math_geom.c
@@ -39,9 +39,9 @@
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])
@@ -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);
- 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);
- len += normalize_v3(n);
+ len += len_v3(n);
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])
{
- 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);
- 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 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 = MAX3(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 */
- 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++) {
- 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 **********************************/
@@ -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);
}
+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 */
@@ -614,6 +728,86 @@ static short IsectLLPt2Df(const float x0, const float y0, const float x1, const
return 1;
}
+/* point in polygon (keep float and int versions in sync) */
+bool isect_point_poly_v2(const float pt[2], const float verts[][2], const int nr)
+{
+ /* we do the angle rule, define that all added angles should be about zero or (2 * PI) */
+ float angletot = 0.0;
+ float fp1[2], fp2[2];
+ int i;
+ const float *p1, *p2;
+
+ p1 = verts[nr - 1];
+ p2 = verts[0];
+
+ /* first vector */
+ fp1[0] = (float)(p1[0] - pt[0]);
+ fp1[1] = (float)(p1[1] - pt[1]);
+ normalize_v2(fp1);
+
+ for (i = 0; i < nr; i++) {
+ float dot, ang, cross;
+ /* second vector */
+ fp2[0] = (float)(p2[0] - pt[0]);
+ fp2[1] = (float)(p2[1] - pt[1]);
+ normalize_v2(fp2);
+
+ /* dot and angle and cross */
+ dot = dot_v2v2(fp1, fp2);
+ ang = fabsf(saacos(dot));
+ cross = (float)((p1[1] - p2[1]) * (p1[0] - pt[0]) + (p2[0] - p1[0]) * (p1[1] - pt[1]));
+
+ if (cross < 0.0f) angletot -= ang;
+ else angletot += ang;
+
+ /* circulate */
+ copy_v2_v2(fp1, fp2);
+ p1 = p2;
+ p2 = verts[i + 1];
+ }
+
+ return (fabsf(angletot) > 4.0f);
+}
+bool isect_point_poly_v2_int(const int pt[2], const int verts[][2], const int nr)
+{
+ /* we do the angle rule, define that all added angles should be about zero or (2 * PI) */
+ float angletot = 0.0;
+ float fp1[2], fp2[2];
+ int i;
+ const int *p1, *p2;
+
+ p1 = verts[nr - 1];
+ p2 = verts[0];
+
+ /* first vector */
+ fp1[0] = (float)(p1[0] - pt[0]);
+ fp1[1] = (float)(p1[1] - pt[1]);
+ normalize_v2(fp1);
+
+ for (i = 0; i < nr; i++) {
+ float dot, ang, cross;
+ /* second vector */
+ fp2[0] = (float)(p2[0] - pt[0]);
+ fp2[1] = (float)(p2[1] - pt[1]);
+ normalize_v2(fp2);
+
+ /* dot and angle and cross */
+ dot = dot_v2v2(fp1, fp2);
+ ang = fabsf(saacos(dot));
+ cross = (float)((p1[1] - p2[1]) * (p1[0] - pt[0]) + (p2[0] - p1[0]) * (p1[1] - pt[1]));
+
+ if (cross < 0.0f) angletot -= ang;
+ else angletot += ang;
+
+ /* circulate */
+ copy_v2_v2(fp1, fp2);
+ p1 = p2;
+ p2 = verts[i + 1];
+ }
+
+ return (fabsf(angletot) > 4.0f);
+}
+
/* point in tri */
/* only single direction */
@@ -893,7 +1087,7 @@ int isect_ray_tri_threshold_v3(const float p1[3], const float d[3],
*/
int isect_line_plane_v3(float out[3],
const float l1[3], const float l2[3],
- const float plane_co[3], const float plane_no[3], const short no_flip)
+ const float plane_co[3], const float plane_no[3], const bool no_flip)
{
float l_vec[3]; /* l1 -> l2 normalized vector */
float p_no[3]; /* 'plane_no' normalized */
@@ -1037,7 +1231,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;
- /* clamp to [0,1] */
+ /* clamp to [0, 1] */
CLAMP(t0, 0.0f, 1.0f);
CLAMP(t1, 0.0f, 1.0f);
@@ -1157,10 +1351,10 @@ int isect_sweeping_sphere_tri_v3(const float p1[3], const float p2[3], const flo
}
/*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);
@@ -1195,13 +1389,13 @@ 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
- 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 (MIN3(v0[a1],v1[a1],v2[a1]) > p1[a1]) return 0;
- if (MIN3(v0[a2],v1[a2],v2[a2]) > p1[a2]) return 0;
- if (MAX3(v0[a1],v1[a1],v2[a1]) < p1[a1]) return 0;
- if (MAX3(v0[a2],v1[a2],v2[a2]) < p1[a2]) return 0;
+ if (min_fff(v0[a1], v1[a1], v2[a1]) > p1[a1]) return 0;
+ if (min_fff(v0[a2], v1[a2], v2[a2]) > p1[a2]) return 0;
+ if (max_fff(v0[a1], v1[a1], v2[a1]) < p1[a1]) return 0;
+ if (max_fff(v0[a2], v1[a2], v2[a2]) < p1[a2]) return 0;
/* then a full intersection test */
#endif
@@ -1415,8 +1609,8 @@ int isect_ray_aabb(const IsectRayAABBData *data, const float bb_min[3],
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])
{
@@ -1702,9 +1896,9 @@ static int point_in_slice(const float p[3], const float v1[3], const float l1[3]
/*
* 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
- * 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'
@@ -1735,7 +1929,7 @@ static int point_in_slice_as(float p[3], float origin[3], float normal[3])
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];
@@ -1803,7 +1997,7 @@ int clip_line_plane(float p1[3], float p2[3], const float plane[4])
}
}
-void plot_line_v2v2i(const int p1[2], const int p2[2], int (*callback)(int, int, void *), void *userData)
+void plot_line_v2v2i(const int p1[2], const int p2[2], bool (*callback)(int, int, void *), void *userData)
{
int x1 = p1[0];
int y1 = p1[1];
@@ -1867,20 +2061,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 */
-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) { *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; }
+ 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]);
+
+ 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]));
@@ -2051,6 +2298,13 @@ void barycentric_weights_v2_quad(const float v1[2], const float v2[2], const flo
len_v2(dirs[3]),
};
+ /* variable 'area' is just for storage,
+ * the order its initialized doesn't matter */
+#ifdef __clang__
+# pragma clang diagnostic push
+# pragma clang diagnostic ignored "-Wunsequenced"
+#endif
+
/* inline mean_value_half_tan four times here */
float t[4] = {
MEAN_VALUE_HALF_TAN_V2(area, 0, 1),
@@ -2059,6 +2313,10 @@ void barycentric_weights_v2_quad(const float v1[2], const float v2[2], const flo
MEAN_VALUE_HALF_TAN_V2(area, 3, 0),
};
+#ifdef __clang__
+# pragma clang diagnostic pop
+#endif
+
#undef MEAN_VALUE_HALF_TAN_V2
w[0] = (t[3] + t[0]) / lens[0];
@@ -2252,59 +2510,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;
- int i;
+ int i, i_next, i_curr;
+ bool edge_interp = false;
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];
- 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);
- w[i] = (len != 0.0f)? (t1 + t2) / len: 0.0f;
+ w[i] = (len != 0.0f) ? (t1 + t2) / len: 0.0f;
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;
- int i;
+ int i, i_next, i_curr;
+ bool edge_interp = false;
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];
- 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);
- w[i] = (len != 0.0f)? (t1 + t2) / len: 0.0f;
+ w[i] = (len != 0.0f) ? (t1 + t2) / len: 0.0f;
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];
@@ -2349,7 +2654,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);
}
- else zero_v2(r_uv);
+ else {
+ zero_v2(r_uv);
+ }
}
/* bilinear reverse */
@@ -2791,8 +3098,8 @@ void tangent_from_uv(float uv1[2], float uv2[2], float uv3[3], float co1[3], flo
/****************************** 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
* (
@@ -2848,7 +3155,9 @@ void vcloud_estimate_transform(int list_size, float (*pos)[3], float *weight, fl
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];
@@ -2857,8 +3166,9 @@ void vcloud_estimate_transform(int list_size, float (*pos)[3], float *weight, fl
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;
@@ -2881,9 +3191,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);
- /* 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);
- /* 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];
@@ -2956,9 +3266,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]);
}
-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];
@@ -3317,8 +3627,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;
@@ -3364,11 +3674,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;
- 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;
}
@@ -3376,44 +3686,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])
{
- 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! */
- 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 */
- 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 */
- return (isect_line_line_v2(v1, v3, v2, v4) > 0) ? TRUE : FALSE;
+ return (isect_line_line_v2(v1, v3, v2, v4) > 0);
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-08-12 07:37:52 +0300