matrix4.c 7.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299
  1. /******************************************************************************
  2. Copyright (C) 2013 by Hugh Bailey <[email protected]>
  3. This program is free software: you can redistribute it and/or modify
  4. it under the terms of the GNU General Public License as published by
  5. the Free Software Foundation, either version 2 of the License, or
  6. (at your option) any later version.
  7. This program is distributed in the hope that it will be useful,
  8. but WITHOUT ANY WARRANTY; without even the implied warranty of
  9. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  10. GNU General Public License for more details.
  11. You should have received a copy of the GNU General Public License
  12. along with this program. If not, see <http://www.gnu.org/licenses/>.
  13. ******************************************************************************/
  14. #include "math-defs.h"
  15. #include "matrix4.h"
  16. #include "matrix3.h"
  17. #include "quat.h"
  18. void matrix4_from_matrix3(struct matrix4 *dst, const struct matrix3 *m)
  19. {
  20. dst->x.m = m->x.m;
  21. dst->y.m = m->y.m;
  22. dst->z.m = m->z.m;
  23. dst->t.m = m->t.m;
  24. dst->t.w = 1.0f;
  25. }
  26. void matrix4_from_quat(struct matrix4 *dst, const struct quat *q)
  27. {
  28. float norm = quat_dot(q, q);
  29. float s = (norm > 0.0f) ? (2.0f / norm) : 0.0f;
  30. float xx = q->x * q->x * s;
  31. float yy = q->y * q->y * s;
  32. float zz = q->z * q->z * s;
  33. float xy = q->x * q->y * s;
  34. float xz = q->x * q->z * s;
  35. float yz = q->y * q->z * s;
  36. float wx = q->w * q->x * s;
  37. float wy = q->w * q->y * s;
  38. float wz = q->w * q->z * s;
  39. vec4_set(&dst->x, 1.0f - (yy + zz), xy + wz, xz - wy, 0.0f);
  40. vec4_set(&dst->y, xy - wz, 1.0f - (xx + zz), yz + wx, 0.0f);
  41. vec4_set(&dst->z, xz + wy, yz - wx, 1.0f - (xx + yy), 0.0f);
  42. vec4_set(&dst->t, 0.0f, 0.0f, 0.0f, 1.0f);
  43. }
  44. void matrix4_from_axisang(struct matrix4 *dst, const struct axisang *aa)
  45. {
  46. struct quat q;
  47. quat_from_axisang(&q, aa);
  48. matrix4_from_quat(dst, &q);
  49. }
  50. void matrix4_mul(struct matrix4 *dst, const struct matrix4 *m1,
  51. const struct matrix4 *m2)
  52. {
  53. const struct vec4 *m1v = (const struct vec4 *)m1;
  54. const float *m2f = (const float *)m2;
  55. struct vec4 out[4];
  56. int i, j;
  57. for (i = 0; i < 4; i++) {
  58. for (j = 0; j < 4; j++) {
  59. struct vec4 temp;
  60. vec4_set(&temp, m2f[j], m2f[j + 4], m2f[j + 8],
  61. m2f[j + 12]);
  62. out[i].ptr[j] = vec4_dot(&m1v[i], &temp);
  63. }
  64. }
  65. matrix4_copy(dst, (struct matrix4 *)out);
  66. }
  67. static inline void get_3x3_submatrix(float *dst, const struct matrix4 *m, int i,
  68. int j)
  69. {
  70. const float *mf = (const float *)m;
  71. int ti, tj, idst, jdst;
  72. for (ti = 0; ti < 4; ti++) {
  73. if (ti < i)
  74. idst = ti;
  75. else if (ti > i)
  76. idst = ti - 1;
  77. else
  78. continue;
  79. for (tj = 0; tj < 4; tj++) {
  80. if (tj < j)
  81. jdst = tj;
  82. else if (tj > j)
  83. jdst = tj - 1;
  84. else
  85. continue;
  86. dst[(idst * 3) + jdst] = mf[(ti * 4) + tj];
  87. }
  88. }
  89. }
  90. static inline float get_3x3_determinant(const float *m)
  91. {
  92. return (m[0] * ((m[4] * m[8]) - (m[7] * m[5]))) -
  93. (m[1] * ((m[3] * m[8]) - (m[6] * m[5]))) +
  94. (m[2] * ((m[3] * m[7]) - (m[6] * m[4])));
  95. }
  96. float matrix4_determinant(const struct matrix4 *m)
  97. {
  98. const float *mf = (const float *)m;
  99. float det, result = 0.0f, i = 1.0f;
  100. float m3x3[9];
  101. int n;
  102. for (n = 0; n < 4; n++, i = -i) { // NOLINT(clang-tidy-cert-flp30-c)
  103. get_3x3_submatrix(m3x3, m, 0, n);
  104. det = get_3x3_determinant(m3x3);
  105. result += mf[n] * det * i;
  106. }
  107. return result;
  108. }
  109. void matrix4_translate3v(struct matrix4 *dst, const struct matrix4 *m,
  110. const struct vec3 *v)
  111. {
  112. struct matrix4 temp;
  113. vec4_set(&temp.x, 1.0f, 0.0f, 0.0f, 0.0f);
  114. vec4_set(&temp.y, 0.0f, 1.0f, 0.0f, 0.0f);
  115. vec4_set(&temp.z, 0.0f, 0.0f, 1.0f, 0.0f);
  116. vec4_from_vec3(&temp.t, v);
  117. matrix4_mul(dst, m, &temp);
  118. }
  119. void matrix4_translate4v(struct matrix4 *dst, const struct matrix4 *m,
  120. const struct vec4 *v)
  121. {
  122. struct matrix4 temp;
  123. vec4_set(&temp.x, 1.0f, 0.0f, 0.0f, 0.0f);
  124. vec4_set(&temp.y, 0.0f, 1.0f, 0.0f, 0.0f);
  125. vec4_set(&temp.z, 0.0f, 0.0f, 1.0f, 0.0f);
  126. vec4_copy(&temp.t, v);
  127. matrix4_mul(dst, m, &temp);
  128. }
  129. void matrix4_rotate(struct matrix4 *dst, const struct matrix4 *m,
  130. const struct quat *q)
  131. {
  132. struct matrix4 temp;
  133. matrix4_from_quat(&temp, q);
  134. matrix4_mul(dst, m, &temp);
  135. }
  136. void matrix4_rotate_aa(struct matrix4 *dst, const struct matrix4 *m,
  137. const struct axisang *aa)
  138. {
  139. struct matrix4 temp;
  140. matrix4_from_axisang(&temp, aa);
  141. matrix4_mul(dst, m, &temp);
  142. }
  143. void matrix4_scale(struct matrix4 *dst, const struct matrix4 *m,
  144. const struct vec3 *v)
  145. {
  146. struct matrix4 temp;
  147. vec4_set(&temp.x, v->x, 0.0f, 0.0f, 0.0f);
  148. vec4_set(&temp.y, 0.0f, v->y, 0.0f, 0.0f);
  149. vec4_set(&temp.z, 0.0f, 0.0f, v->z, 0.0f);
  150. vec4_set(&temp.t, 0.0f, 0.0f, 0.0f, 1.0f);
  151. matrix4_mul(dst, m, &temp);
  152. }
  153. void matrix4_translate3v_i(struct matrix4 *dst, const struct vec3 *v,
  154. const struct matrix4 *m)
  155. {
  156. struct matrix4 temp;
  157. vec4_set(&temp.x, 1.0f, 0.0f, 0.0f, 0.0f);
  158. vec4_set(&temp.y, 0.0f, 1.0f, 0.0f, 0.0f);
  159. vec4_set(&temp.z, 0.0f, 0.0f, 1.0f, 0.0f);
  160. vec4_from_vec3(&temp.t, v);
  161. matrix4_mul(dst, &temp, m);
  162. }
  163. void matrix4_translate4v_i(struct matrix4 *dst, const struct vec4 *v,
  164. const struct matrix4 *m)
  165. {
  166. struct matrix4 temp;
  167. vec4_set(&temp.x, 1.0f, 0.0f, 0.0f, 0.0f);
  168. vec4_set(&temp.y, 0.0f, 1.0f, 0.0f, 0.0f);
  169. vec4_set(&temp.z, 0.0f, 0.0f, 1.0f, 0.0f);
  170. vec4_copy(&temp.t, v);
  171. matrix4_mul(dst, &temp, m);
  172. }
  173. void matrix4_rotate_i(struct matrix4 *dst, const struct quat *q,
  174. const struct matrix4 *m)
  175. {
  176. struct matrix4 temp;
  177. matrix4_from_quat(&temp, q);
  178. matrix4_mul(dst, &temp, m);
  179. }
  180. void matrix4_rotate_aa_i(struct matrix4 *dst, const struct axisang *aa,
  181. const struct matrix4 *m)
  182. {
  183. struct matrix4 temp;
  184. matrix4_from_axisang(&temp, aa);
  185. matrix4_mul(dst, &temp, m);
  186. }
  187. void matrix4_scale_i(struct matrix4 *dst, const struct vec3 *v,
  188. const struct matrix4 *m)
  189. {
  190. struct matrix4 temp;
  191. vec4_set(&temp.x, v->x, 0.0f, 0.0f, 0.0f);
  192. vec4_set(&temp.y, 0.0f, v->y, 0.0f, 0.0f);
  193. vec4_set(&temp.z, 0.0f, 0.0f, v->z, 0.0f);
  194. vec4_set(&temp.t, 0.0f, 0.0f, 0.0f, 1.0f);
  195. matrix4_mul(dst, &temp, m);
  196. }
  197. bool matrix4_inv(struct matrix4 *dst, const struct matrix4 *m)
  198. {
  199. struct vec4 *dstv;
  200. float det;
  201. float m3x3[9];
  202. int i, j, sign;
  203. if (dst == m) {
  204. struct matrix4 temp = *m;
  205. return matrix4_inv(dst, &temp);
  206. }
  207. dstv = (struct vec4 *)dst;
  208. det = matrix4_determinant(m);
  209. if (fabs(det) < 0.0005f)
  210. return false;
  211. for (i = 0; i < 4; i++) {
  212. for (j = 0; j < 4; j++) {
  213. sign = 1 - ((i + j) % 2) * 2;
  214. get_3x3_submatrix(m3x3, m, i, j);
  215. dstv[j].ptr[i] =
  216. get_3x3_determinant(m3x3) * (float)sign / det;
  217. }
  218. }
  219. return true;
  220. }
  221. void matrix4_transpose(struct matrix4 *dst, const struct matrix4 *m)
  222. {
  223. if (dst == m) {
  224. struct matrix4 temp = *m;
  225. matrix4_transpose(dst, &temp);
  226. return;
  227. }
  228. #ifdef NO_INTRINSICS
  229. dst->x.x = m->x.x;
  230. dst->x.y = m->y.x;
  231. dst->x.z = m->z.x;
  232. dst->x.w = m->t.x;
  233. dst->y.x = m->x.y;
  234. dst->y.y = m->y.y;
  235. dst->y.z = m->z.y;
  236. dst->y.w = m->t.y;
  237. dst->z.x = m->x.z;
  238. dst->z.y = m->y.z;
  239. dst->z.z = m->z.z;
  240. dst->z.w = m->t.z;
  241. dst->t.x = m->x.w;
  242. dst->t.y = m->y.w;
  243. dst->t.z = m->z.w;
  244. dst->t.w = m->t.w;
  245. #else
  246. __m128 a0 = _mm_unpacklo_ps(m->x.m, m->z.m);
  247. __m128 a1 = _mm_unpacklo_ps(m->y.m, m->t.m);
  248. __m128 a2 = _mm_unpackhi_ps(m->x.m, m->z.m);
  249. __m128 a3 = _mm_unpackhi_ps(m->y.m, m->t.m);
  250. dst->x.m = _mm_unpacklo_ps(a0, a1);
  251. dst->y.m = _mm_unpackhi_ps(a0, a1);
  252. dst->z.m = _mm_unpacklo_ps(a2, a3);
  253. dst->t.m = _mm_unpackhi_ps(a2, a3);
  254. #endif
  255. }