diff options
Diffstat (limited to 'source/blender/blenlib/intern/math_matrix.c')
| -rw-r--r-- | source/blender/blenlib/intern/math_matrix.c | 174 |
1 files changed, 155 insertions, 19 deletions
diff --git a/source/blender/blenlib/intern/math_matrix.c b/source/blender/blenlib/intern/math_matrix.c index 33d0fb87aca..4bf6d162970 100644 --- a/source/blender/blenlib/intern/math_matrix.c +++ b/source/blender/blenlib/intern/math_matrix.c @@ -516,7 +516,8 @@ void mul_v3_mat3_m4v3(float r[3], float mat[4][4], const float vec[3]) void mul_project_m4_v3(float mat[4][4], float vec[3]) { - const float w = mul_project_m4_v3_zfac(mat, vec); + /* absolute value to not flip the frustum upside down behind the camera */ + const float w = fabsf(mul_project_m4_v3_zfac(mat, vec)); mul_m4_v3(mat, vec); vec[0] /= w; @@ -526,7 +527,7 @@ void mul_project_m4_v3(float mat[4][4], float vec[3]) void mul_v3_project_m4_v3(float r[3], float mat[4][4], const float vec[3]) { - const float w = mul_project_m4_v3_zfac(mat, vec); + const float w = fabsf(mul_project_m4_v3_zfac(mat, vec)); mul_v3_m4v3(r, mat, vec); r[0] /= w; @@ -536,7 +537,7 @@ void mul_v3_project_m4_v3(float r[3], float mat[4][4], const float vec[3]) void mul_v2_project_m4_v3(float r[2], float mat[4][4], const float vec[3]) { - const float w = mul_project_m4_v3_zfac(mat, vec); + const float w = fabsf(mul_project_m4_v3_zfac(mat, vec)); mul_v2_m4v3(r, mat, vec); r[0] /= w; @@ -1247,36 +1248,74 @@ bool is_uniform_scaled_m4(float m[4][4]) return is_uniform_scaled_m3(t); } +void normalize_m3_ex(float mat[3][3], float r_scale[3]) +{ + int i; + for (i = 0; i < 3; i++) { + r_scale[i] = normalize_v3(mat[i]); + } +} void normalize_m3(float mat[3][3]) { - normalize_v3(mat[0]); - normalize_v3(mat[1]); - normalize_v3(mat[2]); + int i; + for (i = 0; i < 3; i++) { + normalize_v3(mat[i]); + } } +void normalize_m3_m3_ex(float rmat[3][3], float mat[3][3], float r_scale[3]) +{ + int i; + for (i = 0; i < 3; i++) { + r_scale[i] = normalize_v3_v3(rmat[i], mat[i]); + } +} void normalize_m3_m3(float rmat[3][3], float mat[3][3]) { - normalize_v3_v3(rmat[0], mat[0]); - normalize_v3_v3(rmat[1], mat[1]); - normalize_v3_v3(rmat[2], mat[2]); + int i; + for (i = 0; i < 3; i++) { + normalize_v3_v3(rmat[i], mat[i]); + } } +void normalize_m4_ex(float mat[4][4], float r_scale[3]) +{ + int i; + for (i = 0; i < 3; i++) { + r_scale[i] = normalize_v3(mat[i]); + if (r_scale[i] != 0.0f) { + mat[i][3] /= r_scale[i]; + } + } +} void normalize_m4(float mat[4][4]) { - float len; - - len = normalize_v3(mat[0]); - if (len != 0.0f) mat[0][3] /= len; - len = normalize_v3(mat[1]); - if (len != 0.0f) mat[1][3] /= len; - len = normalize_v3(mat[2]); - if (len != 0.0f) mat[2][3] /= len; + int i; + for (i = 0; i < 3; i++) { + float len = normalize_v3(mat[i]); + if (len != 0.0f) { + mat[i][3] /= len; + } + } } +void normalize_m4_m4_ex(float rmat[4][4], float mat[4][4], float r_scale[3]) +{ + int i; + for (i = 0; i < 3; i++) { + r_scale[i] = normalize_v3_v3(rmat[i], mat[i]); + rmat[i][3] = (r_scale[i] != 0.0f) ? (mat[i][3] / r_scale[i]) : mat[i][3]; + } + copy_v4_v4(rmat[3], mat[3]); +} void normalize_m4_m4(float rmat[4][4], float mat[4][4]) { - copy_m4_m4(rmat, mat); - normalize_m4(rmat); + int i; + for (i = 0; i < 3; i++) { + float len = normalize_v3_v3(rmat[i], mat[i]); + rmat[i][3] = (len != 0.0f) ? (mat[i][3] / len) : mat[i][3]; + } + copy_v4_v4(rmat[3], mat[3]); } void adjoint_m2_m2(float m1[2][2], float m[2][2]) @@ -1521,6 +1560,34 @@ void mat4_decompose(float loc[3], float quat[4], float size[3], float wmat[4][4] mat3_to_quat(quat, rot); } +/** + * Right polar decomposition: + * M = UP + * + * U is the 'rotation'-like component, the closest orthogonal matrix to M. + * P is the 'scaling'-like component, defined in U space. + * + * See https://en.wikipedia.org/wiki/Polar_decomposition for more. + */ +void mat3_polar_decompose(float mat3[3][3], float r_U[3][3], float r_P[3][3]) +{ + /* From svd decomposition (M = WSV*), we have: + * U = WV* + * P = VSV* + */ + float W[3][3], S[3][3], V[3][3], Vt[3][3]; + float sval[3]; + + BLI_svd_m3(mat3, W, sval, V); + + size_to_mat3(S, sval); + + transpose_m3_m3(Vt, V); + mul_m3_m3m3(r_U, W, Vt); + mul_m3_series(r_P, V, S, Vt); +} + + void scale_m3_fl(float m[3][3], float scale) { m[0][0] = m[1][1] = m[2][2] = scale; @@ -1660,6 +1727,75 @@ void blend_m4_m4m4(float out[4][4], float dst[4][4], float src[4][4], const floa loc_quat_size_to_mat4(out, floc, fquat, fsize); } +/** + * A polar-decomposition-based interpolation between matrix A and matrix B. + * + * \note This code is about five times slower as the 'naive' interpolation done by \a blend_m3_m3m3 + * (it typically remains below 2 usec on an average i74700, while \a blend_m3_m3m3 remains below 0.4 usec). + * However, it gives expected results even with non-uniformaly scaled matrices, see T46418 for an example. + * + * Based on "Matrix Animation and Polar Decomposition", by Ken Shoemake & Tom Duff + * + * @return R the interpolated matrix. + * @param A the intput matrix which is totally effective with \a t = 0.0. + * @param B the intput matrix which is totally effective with \a t = 1.0. + * @param t the interpolation factor. + */ +void interp_m3_m3m3(float R[3][3], float A[3][3], float B[3][3], const float t) +{ + /* 'Rotation' component ('U' part of polar decomposition, the closest orthogonal matrix to M3 rot/scale + * transformation matrix), spherically interpolated. */ + float U_A[3][3], U_B[3][3], U[3][3]; + float quat_A[4], quat_B[4], quat[4]; + /* 'Scaling' component ('P' part of polar decomposition, i.e. scaling in U-defined space), linearly interpolated. */ + float P_A[3][3], P_B[3][3], P[3][3]; + + int i; + + mat3_polar_decompose(A, U_A, P_A); + mat3_polar_decompose(B, U_B, P_B); + + mat3_to_quat(quat_A, U_A); + mat3_to_quat(quat_B, U_B); + interp_qt_qtqt(quat, quat_A, quat_B, t); + quat_to_mat3(U, quat); + + for (i = 0; i < 3; i++) { + interp_v3_v3v3(P[i], P_A[i], P_B[i], t); + } + + /* And we reconstruct rot/scale matrix from interpolated polar components */ + mul_m3_m3m3(R, U, P); +} + +/** + * Complete transform matrix interpolation, based on polar-decomposition-based interpolation from interp_m3_m3m3. + * + * @return R the interpolated matrix. + * @param A the intput matrix which is totally effective with \a t = 0.0. + * @param B the intput matrix which is totally effective with \a t = 1.0. + * @param t the interpolation factor. + */ +void interp_m4_m4m4(float R[4][4], float A[4][4], float B[4][4], const float t) +{ + float A3[3][3], B3[3][3], R3[3][3]; + + /* Location component, linearly interpolated. */ + float loc_A[3], loc_B[3], loc[3]; + + copy_v3_v3(loc_A, A[3]); + copy_v3_v3(loc_B, B[3]); + interp_v3_v3v3(loc, loc_A, loc_B, t); + + copy_m3_m4(A3, A); + copy_m3_m4(B3, B); + + interp_m3_m3m3(R3, A3, B3, t); + + copy_m4_m3(R, R3); + copy_v3_v3(R[3], loc); +} + bool is_negative_m3(float mat[3][3]) { float vec[3]; |