diff options
Diffstat (limited to 'source/blender/blenlib/intern/math_geom.c')
| -rw-r--r-- | source/blender/blenlib/intern/math_geom.c | 299 |
1 files changed, 255 insertions, 44 deletions
diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c index 3cf26ccf904..ee6a3dcc9b3 100644 --- a/source/blender/blenlib/intern/math_geom.c +++ b/source/blender/blenlib/intern/math_geom.c @@ -572,8 +572,8 @@ float dist_squared_to_ray_v3( } /** * Find the closest point in a seg to a ray and return the distance squared. - * \param r_point : Is the point on segment closest to ray (or to ray_origin if the ray and the segment are parallel). - * \param depth: the distance of r_point projection on ray to the ray_origin. + * \param r_point: Is the point on segment closest to ray (or to ray_origin if the ray and the segment are parallel). + * \param r_depth: the distance of r_point projection on ray to the ray_origin. */ float dist_squared_ray_to_seg_v3( const float ray_origin[3], const float ray_direction[3], @@ -619,6 +619,152 @@ float dist_squared_ray_to_seg_v3( return len_squared_v3(t) - SQUARE(*r_depth); } +/* -------------------------------------------------------------------- */ +/** \name dist_squared_to_ray_to_aabb and helpers + * \{ */ + +void dist_squared_ray_to_aabb_v3_precalc( + struct DistRayAABB_Precalc *neasrest_precalc, + const float ray_origin[3], const float ray_direction[3]) +{ + copy_v3_v3(neasrest_precalc->ray_origin, ray_origin); + copy_v3_v3(neasrest_precalc->ray_direction, ray_direction); + + for (int i = 0; i < 3; i++) { + neasrest_precalc->ray_inv_dir[i] = + (neasrest_precalc->ray_direction[i] != 0.0f) ? + (1.0f / neasrest_precalc->ray_direction[i]) : FLT_MAX; + neasrest_precalc->sign[i] = (neasrest_precalc->ray_inv_dir[i] < 0.0f); + } +} + +/** + * Returns the distance from a ray to a bound-box (projected on ray) + */ +float dist_squared_ray_to_aabb_v3( + const struct DistRayAABB_Precalc *data, + const float bb_min[3], const float bb_max[3], + float r_point[3], float *r_depth) +{ + // bool r_axis_closest[3]; + float local_bvmin[3], local_bvmax[3]; + if (data->sign[0]) { + local_bvmin[0] = bb_max[0]; + local_bvmax[0] = bb_min[0]; + } + else { + local_bvmin[0] = bb_min[0]; + local_bvmax[0] = bb_max[0]; + } + if (data->sign[1]) { + local_bvmin[1] = bb_max[1]; + local_bvmax[1] = bb_min[1]; + } + else { + local_bvmin[1] = bb_min[1]; + local_bvmax[1] = bb_max[1]; + } + if (data->sign[2]) { + local_bvmin[2] = bb_max[2]; + local_bvmax[2] = bb_min[2]; + } + else { + local_bvmin[2] = bb_min[2]; + local_bvmax[2] = bb_max[2]; + } + + const float tmin[3] = { + (local_bvmin[0] - data->ray_origin[0]) * data->ray_inv_dir[0], + (local_bvmin[1] - data->ray_origin[1]) * data->ray_inv_dir[1], + (local_bvmin[2] - data->ray_origin[2]) * data->ray_inv_dir[2], + }; + const float tmax[3] = { + (local_bvmax[0] - data->ray_origin[0]) * data->ray_inv_dir[0], + (local_bvmax[1] - data->ray_origin[1]) * data->ray_inv_dir[1], + (local_bvmax[2] - data->ray_origin[2]) * data->ray_inv_dir[2], + }; + /* `va` and `vb` are the coordinates of the AABB edge closest to the ray */ + float va[3], vb[3]; + /* `rtmin` and `rtmax` are the minimum and maximum distances of the ray hits on the AABB */ + float rtmin, rtmax; + int main_axis; + + if ((tmax[0] <= tmax[1]) && (tmax[0] <= tmax[2])) { + rtmax = tmax[0]; + va[0] = vb[0] = local_bvmax[0]; + main_axis = 3; + // r_axis_closest[0] = data->sign[0]; + } + else if ((tmax[1] <= tmax[0]) && (tmax[1] <= tmax[2])) { + rtmax = tmax[1]; + va[1] = vb[1] = local_bvmax[1]; + main_axis = 2; + // r_axis_closest[1] = data->sign[1]; + } + else { + rtmax = tmax[2]; + va[2] = vb[2] = local_bvmax[2]; + main_axis = 1; + // r_axis_closest[2] = data->sign[2]; + } + + if ((tmin[0] >= tmin[1]) && (tmin[0] >= tmin[2])) { + rtmin = tmin[0]; + va[0] = vb[0] = local_bvmin[0]; + main_axis -= 3; + // r_axis_closest[0] = !data->sign[0]; + } + else if ((tmin[1] >= tmin[0]) && (tmin[1] >= tmin[2])) { + rtmin = tmin[1]; + va[1] = vb[1] = local_bvmin[1]; + main_axis -= 1; + // r_axis_closest[1] = !data->sign[1]; + } + else { + rtmin = tmin[2]; + va[2] = vb[2] = local_bvmin[2]; + main_axis -= 2; + // r_axis_closest[2] = !data->sign[2]; + } + if (main_axis < 0) { + main_axis += 3; + } + + /* if rtmin <= rtmax, ray intersect `AABB` */ + if (rtmin <= rtmax) { + float dvec[3]; + copy_v3_v3(r_point, local_bvmax); + sub_v3_v3v3(dvec, local_bvmax, data->ray_origin); + *r_depth = dot_v3v3(dvec, data->ray_direction); + return 0.0f; + } + + if (data->sign[main_axis]) { + va[main_axis] = local_bvmax[main_axis]; + vb[main_axis] = local_bvmin[main_axis]; + } + else { + va[main_axis] = local_bvmin[main_axis]; + vb[main_axis] = local_bvmax[main_axis]; + } + + return dist_squared_ray_to_seg_v3( + data->ray_origin, data->ray_direction, va, vb, + r_point, r_depth); +} + +float dist_squared_ray_to_aabb_v3_simple( + const float ray_origin[3], const float ray_direction[3], + const float bbmin[3], const float bbmax[3], + float r_point[3], float *r_depth) +{ + struct DistRayAABB_Precalc data; + dist_squared_ray_to_aabb_v3_precalc(&data, ray_origin, ray_direction); + return dist_squared_ray_to_aabb_v3(&data, bbmin, bbmax, r_point, r_depth); +} +/** \} */ + + /* Adapted from "Real-Time Collision Detection" by Christer Ericson, * published by Morgan Kaufmann Publishers, copyright 2005 Elsevier Inc. * @@ -765,18 +911,29 @@ int isect_seg_seg_v2(const float v1[2], const float v2[2], const float v3[2], co return ISECT_LINE_LINE_NONE; } -/* get intersection point of two 2D segments and return intersection type: - * -1: collinear - * 1: intersection +/** + * Get intersection point of two 2D segments. + * + * \param endpoint_bias: Bias to use when testing for end-point overlap. + * A positive value considers intersections that extend past the endpoints, + * negative values contract the endpoints. + * Note the bias is applied to a 0-1 factor, not scaled to the length of segments. + * + * \returns intersection type: + * - -1: collinear. + * - 1: intersection. + * - 0: no intersection. */ -int isect_seg_seg_v2_point( +int isect_seg_seg_v2_point_ex( const float v0[2], const float v1[2], const float v2[2], const float v3[2], + const float endpoint_bias, float r_vi[2]) { float s10[2], s32[2], s30[2], d; const float eps = 1e-6f; - const float eps_sq = eps * eps; + const float endpoint_min = -endpoint_bias; + const float endpoint_max = endpoint_bias + 1.0f; sub_v2_v2v2(s10, v1, v0); sub_v2_v2v2(s32, v3, v2); @@ -790,8 +947,8 @@ int isect_seg_seg_v2_point( u = cross_v2v2(s30, s32) / d; v = cross_v2v2(s10, s30) / d; - if ((u >= -eps && u <= 1.0f + eps) && - (v >= -eps && v <= 1.0f + eps)) + if ((u >= endpoint_min && u <= endpoint_max) && + (v >= endpoint_min && v <= endpoint_max)) { /* intersection */ float vi_test[2]; @@ -810,7 +967,7 @@ int isect_seg_seg_v2_point( sub_v2_v2v2(s_vi_v2, vi_test, v2); v = (dot_v2v2(s32, s_vi_v2) / dot_v2v2(s32, s32)); #endif - if (v >= -eps && v <= 1.0f + eps) { + if (v >= endpoint_min && v <= endpoint_max) { copy_v2_v2(r_vi, vi_test); return 1; } @@ -828,7 +985,7 @@ int isect_seg_seg_v2_point( float u_a, u_b; if (equals_v2v2(v0, v1)) { - if (len_squared_v2v2(v2, v3) > eps_sq) { + if (len_squared_v2v2(v2, v3) > SQUARE(eps)) { /* use non-point segment as basis */ SWAP(const float *, v0, v2); SWAP(const float *, v1, v3); @@ -855,7 +1012,7 @@ int isect_seg_seg_v2_point( if (u_a > u_b) SWAP(float, u_a, u_b); - if (u_a > 1.0f + eps || u_b < -eps) { + if (u_a > endpoint_max || u_b < endpoint_min) { /* non-overlapping segments */ return -1; } @@ -871,6 +1028,15 @@ int isect_seg_seg_v2_point( } } +int isect_seg_seg_v2_point( + const float v0[2], const float v1[2], + const float v2[2], const float v3[2], + float r_vi[2]) +{ + const float endpoint_bias = 1e-6f; + return isect_seg_seg_v2_point_ex(v0, v1, v2, v3, endpoint_bias, r_vi); +} + bool isect_seg_seg_v2_simple(const float v1[2], const float v2[2], const float v3[2], const float v4[2]) { #define CCW(A, B, C) \ @@ -1828,7 +1994,7 @@ bool isect_tri_tri_epsilon_v3( (range[0].max < range[1].min)) == 0) { if (r_i1 && r_i2) { - project_plane_v3_v3v3(plane_co, plane_co, plane_no); + project_plane_normalized_v3_v3v3(plane_co, plane_co, plane_no); madd_v3_v3v3fl(r_i1, plane_co, plane_no, max_ff(range[0].min, range[1].min)); madd_v3_v3v3fl(r_i2, plane_co, plane_no, min_ff(range[0].max, range[1].max)); } @@ -2851,8 +3017,13 @@ bool barycentric_coords_v2(const float v1[2], const float v2[2], const float v3[ /** * \note: using #cross_tri_v2 means locations outside the triangle are correctly weighted + * + * \note This is *exactly* the same calculation as #resolve_tri_uv_v2, + * although it has double precision and is used for texture baking, so keep both. */ -void barycentric_weights_v2(const float v1[2], const float v2[2], const float v3[2], const float co[2], float w[3]) +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; @@ -2870,10 +3041,35 @@ void barycentric_weights_v2(const float v1[2], const float v2[2], const float v3 } /** + * A version of #barycentric_weights_v2 that doesn't allow negative weights. + * Useful when negative values cause problems and points are only ever slightly outside of the triangle. + */ +void barycentric_weights_v2_clamped( + const float v1[2], const float v2[2], const float v3[2], + const float co[2], float w[3]) +{ + float wtot; + + w[0] = max_ff(cross_tri_v2(v2, v3, co), 0.0f); + w[1] = max_ff(cross_tri_v2(v3, v1, co), 0.0f); + w[2] = max_ff(cross_tri_v2(v1, v2, co), 0.0f); + wtot = w[0] + w[1] + w[2]; + + if (wtot != 0.0f) { + mul_v3_fl(w, 1.0f / wtot); + } + else { /* dummy values for zero area face */ + copy_v3_fl(w, 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 */ -void barycentric_weights_v2_persp(const float v1[4], const float v2[4], const float v3[4], const float co[2], float w[3]) +void barycentric_weights_v2_persp( + const float v1[4], const float v2[4], const float v3[4], + const float co[2], float w[3]) { float wtot; @@ -2890,11 +3086,14 @@ void barycentric_weights_v2_persp(const float v1[4], const float v2[4], const fl } } -/* same as #barycentric_weights_v2 but works with a quad, +/** + * same as #barycentric_weights_v2 but works with a quad, * note: untested for values outside the quad's bounds - * this is #interp_weights_poly_v2 expanded for quads only */ -void barycentric_weights_v2_quad(const float v1[2], const float v2[2], const float v3[2], const float v4[2], - const float co[2], float w[4]) + * this is #interp_weights_poly_v2 expanded for quads only + */ +void barycentric_weights_v2_quad( + const float v1[2], const float v2[2], const float v3[2], const float v4[2], + const float co[2], float w[4]) { /* note: fabsf() here is not needed for convex quads (and not used in interp_weights_poly_v2). * but in the case of concave/bow-tie quads for the mask rasterizer it gives unreliable results @@ -3345,6 +3544,8 @@ void interp_cubic_v3(float x[3], float v[3], const float x1[3], const float v1[3 * Barycentric reverse * * Compute coordinates (u, v) for point \a st with respect to triangle (\a st0, \a st1, \a st2) + * + * \note same basic result as #barycentric_weights_v2, see it's comment for details. */ void resolve_tri_uv_v2(float r_uv[2], const float st[2], const float st0[2], const float st1[2], const float st2[2]) @@ -3575,7 +3776,7 @@ void orthographic_m4(float matrix[4][4], const float left, const float right, co matrix[3][0] = -(right + left) / Xdelta; matrix[1][1] = 2.0f / Ydelta; matrix[3][1] = -(top + bottom) / Ydelta; - matrix[2][2] = -2.0f / Zdelta; /* note: negate Z */ + matrix[2][2] = -2.0f / Zdelta; /* note: negate Z */ matrix[3][2] = -(farClip + nearClip) / Zdelta; } @@ -3594,7 +3795,7 @@ void perspective_m4(float mat[4][4], const float left, const float right, const } mat[0][0] = nearClip * 2.0f / Xdelta; mat[1][1] = nearClip * 2.0f / Ydelta; - mat[2][0] = (right + left) / Xdelta; /* note: negate Z */ + mat[2][0] = (right + left) / Xdelta; /* note: negate Z */ mat[2][1] = (top + bottom) / Ydelta; mat[2][2] = -(farClip + nearClip) / Zdelta; mat[2][3] = -1.0f; @@ -3880,7 +4081,7 @@ void map_to_plane_axis_angle_v2_v3v3fl(float r_co[2], const float co[3], const f /********************************* Normals **********************************/ -void accumulate_vertex_normals_tri( +void accumulate_vertex_normals_tri_v3( float n1[3], float n2[3], float n3[3], const float f_no[3], const float co1[3], const float co2[3], const float co3[3]) @@ -3914,7 +4115,7 @@ void accumulate_vertex_normals_tri( } } -void accumulate_vertex_normals( +void accumulate_vertex_normals_v3( float n1[3], float n2[3], float n3[3], float n4[3], const float f_no[3], const float co1[3], const float co2[3], const float co3[3], const float co4[3]) @@ -3958,7 +4159,7 @@ void accumulate_vertex_normals( /* Add weighted face normal component into normals of the face vertices. * Caller must pass pre-allocated vdiffs of nverts length. */ -void accumulate_vertex_normals_poly(float **vertnos, const float polyno[3], +void accumulate_vertex_normals_poly_v3(float **vertnos, const float polyno[3], const float **vertcos, float vdiffs[][3], const int nverts) { int i; @@ -3989,7 +4190,7 @@ void accumulate_vertex_normals_poly(float **vertnos, const float polyno[3], /********************************* Tangents **********************************/ -void tangent_from_uv( +void tangent_from_uv_v3( const float uv1[2], const float uv2[2], const float uv3[3], const float co1[3], const float co2[3], const float co3[3], const float n[3], @@ -4031,30 +4232,28 @@ void tangent_from_uv( /****************************** 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]) - * +/** * input - * ( - * int list_size - * 4 lists as pointer to array[list_size] - * 1. current pos array of 'new' positions - * 2. current weight array of 'new'weights (may be NULL pointer if you have no weights ) - * 3. reference rpos array of 'old' positions - * 4. reference rweight array of 'old'weights (may be NULL pointer if you have no weights ) - * ) + * + * \param list_size: 4 lists as pointer to array[list_size] + * \param pos: current pos array of 'new' positions + * \param weight: current weight array of 'new'weights (may be NULL pointer if you have no weights) + * \param rpos: Reference rpos array of 'old' positions + * \param rweight: Reference rweight array of 'old'weights (may be NULL pointer if you have no weights). + * * output - * ( - * float lloc[3] center of mass pos - * float rloc[3] center of mass rpos - * float lrot[3][3] rotation matrix - * float lscale[3][3] scale matrix + * + * \param lloc: Center of mass pos. + * \param rloc: Center of mass rpos. + * \param lrot: Rotation matrix. + * \param lscale: Scale matrix. + * * pointers may be NULL if not needed - * ) */ -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_v3( + const int list_size, const float (*pos)[3], const float *weight, const float (*rpos)[3], const float *rweight, + float lloc[3], float rloc[3], float lrot[3][3], float lscale[3][3]) { float accu_com[3] = {0.0f, 0.0f, 0.0f}, accu_rcom[3] = {0.0f, 0.0f, 0.0f}; float accu_weight = 0.0f, accu_rweight = 0.0f; @@ -4732,6 +4931,18 @@ int is_quad_flip_v3(const float v1[3], const float v2[3], const float v3[3], con return ret; } +bool is_quad_flip_v3_first_third_fast(const float v1[3], const float v2[3], const float v3[3], const float v4[3]) +{ + float d_12[3], d_13[3], d_14[3]; + float cross_a[3], cross_b[3]; + sub_v3_v3v3(d_12, v2, v1); + sub_v3_v3v3(d_13, v3, v1); + sub_v3_v3v3(d_14, v4, v1); + cross_v3_v3v3(cross_a, d_12, d_13); + cross_v3_v3v3(cross_b, d_14, d_13); + return dot_v3v3(cross_a, cross_b) > 0.0f; +} + /** * Return the value which the distance between points will need to be scaled by, * to define a handle, given both points are on a perfect circle. |