3dphotoshoot
view src/sanegl.c @ 11:ad49e1f9b627
foo
| author | John Tsiombikas <nuclear@member.fsf.org> |
|---|---|
| date | Sun, 31 May 2015 06:02:08 +0300 |
| parents | |
| children | a460b1e5af4a |
line source
1 /*
2 SaneGL - a small library to bring back sanity to OpenGL ES 2.x
3 Copyright (C) 2011-2013 John Tsiombikas <nuclear@member.fsf.org>
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <string.h>
22 #include <math.h>
23 #include <assert.h>
24 #include "sanegl.h"
26 #define MMODE_IDX(x) ((x) - GL_MODELVIEW)
27 #define MAT_STACK_SIZE 32
28 #define MAT_IDENT {1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1}
30 #define MAX_VERTS 512
32 static void gl_draw_immediate(void);
33 static void m4_transpose(double *res, double *m);
34 static double m4_determinant(double *m);
35 static void m4_adjoint(double *res, double *m);
36 static void m4_inverse(double *res, double *m);
39 typedef struct { float x, y; } vec2_t;
40 typedef struct { float x, y, z; } vec3_t;
41 typedef struct { float x, y, z, w; } vec4_t;
43 static int mm_idx = 0;
44 static float mat_stack[3][MAT_STACK_SIZE][16] = {{MAT_IDENT}, {MAT_IDENT}, {MAT_IDENT}};
45 static int stack_top[3];
46 static float mat_mvp[16];
47 static int mvp_valid;
48 static int prim = -1;
50 static vec3_t cur_normal;
51 static vec4_t cur_color, cur_attrib;
52 static vec2_t cur_texcoord;
54 static vec4_t *vert_arr, *col_arr, *attr_arr;
55 static vec3_t *norm_arr;
56 static vec2_t *texc_arr;
57 /*static unsigned int vbuf, cbuf, nbuf, tbuf, abuf;*/
58 static int vloc, nloc, cloc, tloc, aloc = -1;
60 static int num_verts, vert_calls;
61 static int cur_prog;
63 #ifdef GLDEF
64 #undef glEnable
65 #undef glDisable
66 #endif
68 void gl_enable(int state)
69 {
70 switch(state) {
71 case GL_TEXTURE_2D:
72 break;
74 default:
75 glEnable(state);
76 }
77 }
79 void gl_disable(int state)
80 {
81 switch(state) {
82 case GL_TEXTURE_2D:
83 glBindTexture(state, 0);
84 break;
86 default:
87 glDisable(state);
88 }
89 }
91 void gl_matrix_mode(int mm)
92 {
93 mm_idx = MMODE_IDX(mm);
94 }
96 void gl_push_matrix(void)
97 {
98 int top = stack_top[mm_idx];
100 memcpy(mat_stack[mm_idx][top + 1], mat_stack[mm_idx][top], 16 * sizeof(float));
101 stack_top[mm_idx]++;
102 mvp_valid = 0;
103 }
105 void gl_pop_matrix(void)
106 {
107 stack_top[mm_idx]--;
108 mvp_valid = 0;
109 }
111 void gl_load_identity(void)
112 {
113 static const float idmat[] = MAT_IDENT;
114 int top = stack_top[mm_idx];
115 float *mat = mat_stack[mm_idx][top];
117 memcpy(mat, idmat, sizeof idmat);
118 mvp_valid = 0;
119 }
121 void gl_load_matrixf(const float *m)
122 {
123 int top = stack_top[mm_idx];
124 float *mat = mat_stack[mm_idx][top];
126 memcpy(mat, m, 16 * sizeof *mat);
127 mvp_valid = 0;
128 }
130 void gl_load_matrixd(const double *m)
131 {
132 int i;
133 float mf[16];
135 for(i=0; i<16; i++) {
136 mf[i] = (float)m[i];
137 }
138 gl_load_matrixf(mf);
139 }
141 #define M(i, j) ((i << 2) + j)
143 void gl_mult_matrixf(const float *m2)
144 {
145 int i, j;
146 int top = stack_top[mm_idx];
147 float *m1 = mat_stack[mm_idx][top];
148 float res[16];
150 for(i=0; i<4; i++) {
151 for(j=0; j<4; j++) {
152 res[M(i,j)] = m2[M(i,0)] * m1[M(0,j)] +
153 m2[M(i,1)] * m1[M(1,j)] +
154 m2[M(i,2)] * m1[M(2,j)] +
155 m2[M(i,3)] * m1[M(3,j)];
156 }
157 }
159 memcpy(m1, res, sizeof res);
160 mvp_valid = 0;
161 }
163 void gl_mult_matrixd(const double *m)
164 {
165 int i;
166 float mf[16];
168 for(i=0; i<16; i++) {
169 mf[i] = (float)m[i];
170 }
171 gl_mult_matrixf(mf);
172 }
174 void gl_translatef(float x, float y, float z)
175 {
176 float mat[] = MAT_IDENT;
178 mat[12] = x;
179 mat[13] = y;
180 mat[14] = z;
182 gl_mult_matrixf(mat);
183 }
185 void gl_rotatef(float angle, float x, float y, float z)
186 {
187 float mat[] = MAT_IDENT;
189 float angle_rad = M_PI * angle / 180.0;
190 float sina = sin(angle_rad);
191 float cosa = cos(angle_rad);
192 float one_minus_cosa = 1.0 - cosa;
193 float nxsq = x * x;
194 float nysq = y * y;
195 float nzsq = z * z;
197 mat[0] = nxsq + (1.0 - nxsq) * cosa;
198 mat[4] = x * y * one_minus_cosa - z * sina;
199 mat[8] = x * z * one_minus_cosa + y * sina;
200 mat[1] = x * y * one_minus_cosa + z * sina;
201 mat[5] = nysq + (1.0 - nysq) * cosa;
202 mat[9] = y * z * one_minus_cosa - x * sina;
203 mat[2] = x * z * one_minus_cosa - y * sina;
204 mat[6] = y * z * one_minus_cosa + x * sina;
205 mat[10] = nzsq + (1.0 - nzsq) * cosa;
207 gl_mult_matrixf(mat);
208 }
210 void gl_scalef(float x, float y, float z)
211 {
212 float mat[] = MAT_IDENT;
214 mat[0] = x;
215 mat[5] = y;
216 mat[10] = z;
218 gl_mult_matrixf(mat);
219 }
221 void gl_ortho(float left, float right, float bottom, float top, float near, float far)
222 {
223 float mat[] = MAT_IDENT;
225 float dx = right - left;
226 float dy = top - bottom;
227 float dz = far - near;
229 float tx = -(right + left) / dx;
230 float ty = -(top + bottom) / dy;
231 float tz = -(far + near) / dz;
233 float sx = 2.0 / dx;
234 float sy = 2.0 / dy;
235 float sz = -2.0 / dz;
237 mat[0] = sx;
238 mat[5] = sy;
239 mat[10] = sz;
240 mat[12] = tx;
241 mat[13] = ty;
242 mat[14] = tz;
244 gl_mult_matrixf(mat);
245 }
247 void gl_frustum(float left, float right, float bottom, float top, float near, float far)
248 {
249 float mat[] = MAT_IDENT;
251 float dx = right - left;
252 float dy = top - bottom;
253 float dz = far - near;
255 float a = (right + left) / dx;
256 float b = (top + bottom) / dy;
257 float c = -(far + near) / dz;
258 float d = -2.0 * far * near / dz;
260 mat[0] = 2.0 * near / dx;
261 mat[5] = 2.0 * near / dy;
262 mat[8] = a;
263 mat[9] = b;
264 mat[10] = c;
265 mat[11] = -1.0;
266 mat[14] = d;
267 mat[15] = 0.0;
269 gl_mult_matrixf(mat);
270 }
272 void glu_perspective(float vfov, float aspect, float near, float far)
273 {
274 float vfov_rad = M_PI * vfov / 180.0;
275 float x = near * tan(vfov_rad / 2.0);
276 gl_frustum(-aspect * x, aspect * x, -x, x, near, far);
277 }
279 int glu_un_project(double winx, double winy, double winz,
280 const double *model, const double *proj, const int *viewp,
281 double *objx, double *objy, double *objz)
282 {
283 double mvp[16], inv_mvp[16];
285 double ndcx = 2.0 * (winx - viewp[0]) / viewp[2] - 1.0;
286 double ndcy = 2.0 * (winy - viewp[1]) / viewp[3] - 1.0;
287 double ndcz = 2.0 * winz - 1.0;
289 // calculate modelviewprojection
290 gl_matrix_mode(GL_MODELVIEW);
291 gl_push_matrix();
292 gl_load_matrixd(proj);
293 gl_mult_matrixd(model);
294 gl_get_doublev(GL_MODELVIEW_MATRIX, mvp);
295 gl_pop_matrix();
297 // invert modelviewprojection
298 m4_inverse(inv_mvp, mvp);
300 // transform ndc by modelview -> obj
301 /**objx = inv_mvp[0] * ndcx + inv_mvp[4] * ndcy + inv_mvp[8] * ndcz + inv_mvp[12];
302 *objy = inv_mvp[1] * ndcx + inv_mvp[5] * ndcy + inv_mvp[9] * ndcz + inv_mvp[13];
303 *objz = inv_mvp[2] * ndcx + inv_mvp[6] * ndcy + inv_mvp[10] * ndcz + inv_mvp[14];*/
304 *objx = inv_mvp[0] * ndcx + inv_mvp[1] * ndcy + inv_mvp[2] * ndcz + inv_mvp[3];
305 *objy = inv_mvp[4] * ndcx + inv_mvp[5] * ndcy + inv_mvp[6] * ndcz + inv_mvp[7];
306 *objz = inv_mvp[8] * ndcx + inv_mvp[9] * ndcy + inv_mvp[10] * ndcz + inv_mvp[11];
307 return 0;
308 }
310 void gl_apply_xform(unsigned int prog)
311 {
312 int loc, mvidx, pidx, tidx, mvtop, ptop, ttop;
314 mvidx = MMODE_IDX(GL_MODELVIEW);
315 pidx = MMODE_IDX(GL_PROJECTION);
316 tidx = MMODE_IDX(GL_TEXTURE);
318 mvtop = stack_top[mvidx];
319 ptop = stack_top[pidx];
320 ttop = stack_top[tidx];
322 assert(prog);
324 /*printf("APPLY XFORM\n");*/
326 CHECK_GLERROR;
327 if((loc = glGetUniformLocation(prog, "matrix_modelview")) != -1) {
328 CHECK_GLERROR;
329 /*printf(" MODELVIEW:\n");
330 for(i=0; i<16; i+=4) {
331 printf("%.2f %.2f %.2f %.2f\n", mat_stack[mvidx][mvtop][i], mat_stack[mvidx][mvtop][i + 1], mat_stack[mvidx][mvtop][i + 2], mat_stack[mvidx][mvtop][i + 3]);
332 }*/
333 glUniformMatrix4fv(loc, 1, 0, mat_stack[mvidx][mvtop]);
334 CHECK_GLERROR;
335 }
336 CHECK_GLERROR;
338 if((loc = glGetUniformLocation(prog, "matrix_projection")) != -1) {
339 CHECK_GLERROR;
340 /*printf(" PROJECTION:\n");
341 for(i=0; i<16; i+=4) {
342 printf("%.2f %.2f %.2f %.2f\n", mat_stack[pidx][ptop][i], mat_stack[pidx][ptop][i + 1], mat_stack[pidx][ptop][i + 2], mat_stack[pidx][ptop][i + 3]);
343 }*/
344 glUniformMatrix4fv(loc, 1, 0, mat_stack[pidx][ptop]);
345 CHECK_GLERROR;
346 }
347 CHECK_GLERROR;
349 if((loc = glGetUniformLocation(prog, "matrix_texture")) != -1) {
350 CHECK_GLERROR;
351 glUniformMatrix4fv(loc, 1, 0, mat_stack[tidx][ttop]);
352 CHECK_GLERROR;
353 }
354 CHECK_GLERROR;
356 if((loc = glGetUniformLocation(prog, "matrix_normal")) != -1) {
357 float nmat[9];
359 CHECK_GLERROR;
361 nmat[0] = mat_stack[mvidx][mvtop][0];
362 nmat[1] = mat_stack[mvidx][mvtop][1];
363 nmat[2] = mat_stack[mvidx][mvtop][2];
364 nmat[3] = mat_stack[mvidx][mvtop][4];
365 nmat[4] = mat_stack[mvidx][mvtop][5];
366 nmat[5] = mat_stack[mvidx][mvtop][6];
367 nmat[6] = mat_stack[mvidx][mvtop][8];
368 nmat[7] = mat_stack[mvidx][mvtop][9];
369 nmat[8] = mat_stack[mvidx][mvtop][10];
370 glUniformMatrix3fv(loc, 1, 0, nmat);
371 CHECK_GLERROR;
372 }
373 CHECK_GLERROR;
375 if((loc = glGetUniformLocation(prog, "matrix_modelview_projection")) != -1) {
376 CHECK_GLERROR;
377 if(!mvp_valid) {
378 /* TODO calc mvp */
379 }
380 glUniformMatrix4fv(loc, 1, 0, mat_mvp);
381 CHECK_GLERROR;
382 }
383 CHECK_GLERROR;
384 }
387 /* immediate mode rendering */
388 void gl_begin(int p)
389 {
390 if(!vert_arr) {
391 vert_arr = malloc(MAX_VERTS * sizeof *vert_arr);
392 norm_arr = malloc(MAX_VERTS * sizeof *norm_arr);
393 texc_arr = malloc(MAX_VERTS * sizeof *texc_arr);
394 col_arr = malloc(MAX_VERTS * sizeof *col_arr);
395 attr_arr = malloc(MAX_VERTS * sizeof *attr_arr);
396 assert(vert_arr && norm_arr && texc_arr && col_arr && attr_arr);
397 }
399 prim = p;
400 num_verts = vert_calls = 0;
402 glGetIntegerv(GL_CURRENT_PROGRAM, &cur_prog);
403 CHECK_GLERROR;
404 assert(cur_prog);
406 gl_apply_xform(cur_prog);
407 CHECK_GLERROR;
409 vloc = glGetAttribLocation(cur_prog, "attr_vertex");
410 CHECK_GLERROR;
411 nloc = glGetAttribLocation(cur_prog, "attr_normal");
412 CHECK_GLERROR;
413 cloc = glGetAttribLocation(cur_prog, "attr_color");
414 CHECK_GLERROR;
415 tloc = glGetAttribLocation(cur_prog, "attr_texcoord");
416 CHECK_GLERROR;
417 }
419 void gl_end(void)
420 {
421 if(num_verts > 0) {
422 gl_draw_immediate();
423 }
424 aloc = -1;
425 }
427 static void gl_draw_immediate(void)
428 {
429 int glprim;
431 if(vloc == -1) {
432 fprintf(stderr, "gl_draw_immediate call with vloc == -1\n");
433 return;
434 }
436 glprim = prim == GL_QUADS ? GL_TRIANGLES : prim;
438 CHECK_GLERROR;
439 glVertexAttribPointer(vloc, 4, GL_FLOAT, 0, 0, vert_arr);
440 CHECK_GLERROR;
441 glEnableVertexAttribArray(vloc);
442 CHECK_GLERROR;
444 if(nloc != -1) {
445 glVertexAttribPointer(nloc, 3, GL_FLOAT, 0, 0, norm_arr);
446 CHECK_GLERROR;
447 glEnableVertexAttribArray(nloc);
448 CHECK_GLERROR;
449 }
451 if(cloc != -1) {
452 glVertexAttribPointer(cloc, 4, GL_FLOAT, 1, 0, col_arr);
453 CHECK_GLERROR;
454 glEnableVertexAttribArray(cloc);
455 CHECK_GLERROR;
456 }
458 if(tloc != -1) {
459 glVertexAttribPointer(tloc, 2, GL_FLOAT, 0, 0, texc_arr);
460 CHECK_GLERROR;
461 glEnableVertexAttribArray(tloc);
462 CHECK_GLERROR;
463 }
465 if(aloc != -1) {
466 glVertexAttribPointer(aloc, 4, GL_FLOAT, 0, 0, attr_arr);
467 CHECK_GLERROR;
468 glEnableVertexAttribArray(aloc);
469 CHECK_GLERROR;
470 }
472 glDrawArrays(glprim, 0, num_verts);
473 CHECK_GLERROR;
475 glDisableVertexAttribArray(vloc);
476 CHECK_GLERROR;
477 if(nloc != -1) {
478 glDisableVertexAttribArray(nloc);
479 CHECK_GLERROR;
480 }
481 if(cloc != -1) {
482 glDisableVertexAttribArray(cloc);
483 CHECK_GLERROR;
484 }
485 if(tloc != -1) {
486 glDisableVertexAttribArray(tloc);
487 CHECK_GLERROR;
488 }
489 if(aloc != -1) {
490 glDisableVertexAttribArray(aloc);
491 CHECK_GLERROR;
492 }
493 }
496 void gl_vertex2f(float x, float y)
497 {
498 gl_vertex4f(x, y, 0.0f, 1.0f);
499 }
501 void gl_vertex3f(float x, float y, float z)
502 {
503 gl_vertex4f(x, y, z, 1.0f);
504 }
506 void gl_vertex4f(float x, float y, float z, float w)
507 {
508 int i, buffer_full;
510 if(prim == GL_QUADS && vert_calls % 4 == 3) {
511 for(i=0; i<2; i++) {
512 if(aloc != -1) {
513 attr_arr[num_verts] = attr_arr[num_verts - 3 + i];
514 }
515 if(cloc != -1) {
516 col_arr[num_verts] = col_arr[num_verts - 3 + i];
517 }
518 if(tloc != -1) {
519 texc_arr[num_verts] = texc_arr[num_verts - 3 + i];
520 }
521 if(nloc != -1) {
522 norm_arr[num_verts] = norm_arr[num_verts - 3 + i];
523 }
524 vert_arr[num_verts] = vert_arr[num_verts - 3 + i];
525 num_verts++;
526 }
527 }
529 vert_arr[num_verts].x = x;
530 vert_arr[num_verts].y = y;
531 vert_arr[num_verts].z = z;
532 vert_arr[num_verts].w = w;
534 if(cloc != -1) {
535 col_arr[num_verts] = cur_color;
536 }
537 if(nloc != -1) {
538 norm_arr[num_verts] = cur_normal;
539 }
540 if(tloc != -1) {
541 texc_arr[num_verts] = cur_texcoord;
542 }
543 if(aloc != -1) {
544 attr_arr[num_verts] = cur_attrib;
545 }
547 vert_calls++;
548 num_verts++;
550 if(prim == GL_QUADS) {
551 /* leave space for 6 more worst-case and don't allow flushes mid-quad */
552 buffer_full = num_verts >= MAX_VERTS - 6 && vert_calls % 4 == 0;
553 } else {
554 buffer_full = num_verts >= MAX_VERTS - prim;
555 }
557 if(buffer_full) {
558 gl_draw_immediate();
559 gl_begin(prim); /* reset everything */
560 }
561 }
564 void gl_normal3f(float x, float y, float z)
565 {
566 cur_normal.x = x;
567 cur_normal.y = y;
568 cur_normal.z = z;
569 }
572 void gl_color3f(float r, float g, float b)
573 {
574 cur_color.x = r;
575 cur_color.y = g;
576 cur_color.z = b;
577 cur_color.w = 1.0f;
578 }
580 void gl_color4f(float r, float g, float b, float a)
581 {
582 cur_color.x = r;
583 cur_color.y = g;
584 cur_color.z = b;
585 cur_color.w = a;
586 }
589 void gl_texcoord1f(float s)
590 {
591 cur_texcoord.x = s;
592 cur_texcoord.y = 0.0f;
593 }
595 void gl_texcoord2f(float s, float t)
596 {
597 cur_texcoord.x = s;
598 cur_texcoord.y = t;
599 }
601 void gl_vertex_attrib2f(int loc, float x, float y)
602 {
603 aloc = loc;
604 cur_attrib.x = x;
605 cur_attrib.y = y;
606 cur_attrib.z = 0.0f;
607 cur_attrib.w = 1.0f;
608 }
610 void gl_vertex_attrib3f(int loc, float x, float y, float z)
611 {
612 aloc = loc;
613 cur_attrib.x = x;
614 cur_attrib.y = y;
615 cur_attrib.z = z;
616 cur_attrib.w = 1.0f;
617 }
619 void gl_vertex_attrib4f(int loc, float x, float y, float z, float w)
620 {
621 aloc = loc;
622 cur_attrib.x = x;
623 cur_attrib.y = y;
624 cur_attrib.z = z;
625 cur_attrib.w = w;
626 }
628 #ifdef GLDEF
629 #undef glGetFloatv
630 #endif
632 void gl_get_floatv(int what, float *res)
633 {
634 int idx;
636 switch(what) {
637 case GL_MODELVIEW_MATRIX:
638 idx = MMODE_IDX(GL_MODELVIEW);
639 memcpy(res, mat_stack[idx][stack_top[idx]], 16 * sizeof *res);
640 break;
642 case GL_PROJECTION_MATRIX:
643 idx = MMODE_IDX(GL_PROJECTION);
644 memcpy(res, mat_stack[idx][stack_top[idx]], 16 * sizeof *res);
645 break;
647 default:
648 glGetFloatv(what, res);
649 }
650 }
652 void gl_get_doublev(int what, double *res)
653 {
654 int i, idx;
655 float tmp[16];
657 switch(what) {
658 case GL_MODELVIEW_MATRIX:
659 if(1) {
660 idx = MMODE_IDX(GL_MODELVIEW);
661 } else {
662 case GL_PROJECTION_MATRIX:
663 idx = MMODE_IDX(GL_PROJECTION);
664 }
665 for(i=0; i<16; i++) {
666 res[i] = mat_stack[idx][stack_top[idx]][i];
667 }
668 break;
670 default:
671 glGetFloatv(what, tmp);
672 for(i=0; i<16; i++) {
673 res[i] = tmp[i];
674 }
675 }
676 }
679 /* ---- matrix inversion stuff ---- */
680 static void m4_transpose(double *res, double *m)
681 {
682 int i, j;
683 double tmp[16];
685 if(res == m) {
686 memcpy(tmp, m, 16 * sizeof *m);
687 m = tmp;
688 }
690 for(i=0; i<4; i++) {
691 for(j=0; j<4; j++) {
692 res[M(i, j)] = m[M(j, i)];
693 }
694 }
695 }
697 static double m4_determinant(double *m)
698 {
699 double det11, det12, det13, det14;
701 det11 = (m[M(1, 1)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(3, 2)] * m[M(2, 3)])) -
702 (m[M(1, 2)] * (m[M(2, 1)] * m[M(3, 3)] - m[M(3, 1)] * m[M(2, 3)])) +
703 (m[M(1, 3)] * (m[M(2, 1)] * m[M(3, 2)] - m[M(3, 1)] * m[M(2, 2)]));
705 det12 = (m[M(1, 0)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(3, 2)] * m[M(2, 3)])) -
706 (m[M(1, 2)] * (m[M(2, 0)] * m[M(3, 3)] - m[M(3, 0)] * m[M(2, 3)])) +
707 (m[M(1, 3)] * (m[M(2, 0)] * m[M(3, 2)] - m[M(3, 0)] * m[M(2, 2)]));
709 det13 = (m[M(1, 0)] * (m[M(2, 1)] * m[M(3, 3)] - m[M(3, 1)] * m[M(2, 3)])) -
710 (m[M(1, 1)] * (m[M(2, 0)] * m[M(3, 3)] - m[M(3, 0)] * m[M(2, 3)])) +
711 (m[M(1, 3)] * (m[M(2, 0)] * m[M(3, 1)] - m[M(3, 0)] * m[M(2, 1)]));
713 det14 = (m[M(1, 0)] * (m[M(2, 1)] * m[M(3, 2)] - m[M(3, 1)] * m[M(2, 2)])) -
714 (m[M(1, 1)] * (m[M(2, 0)] * m[M(3, 2)] - m[M(3, 0)] * m[M(2, 2)])) +
715 (m[M(1, 2)] * (m[M(2, 0)] * m[M(3, 1)] - m[M(3, 0)] * m[M(2, 1)]));
717 return m[M(0, 0)] * det11 - m[M(0, 1)] * det12 + m[M(0, 2)] * det13 - m[M(0, 3)] * det14;
718 }
720 static void m4_adjoint(double *res, double *m)
721 {
722 int i, j;
723 double coef[16];
725 coef[M(0, 0)] = (m[M(1, 1)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(3, 2)] * m[M(2, 3)])) -
726 (m[M(1, 2)] * (m[M(2, 1)] * m[M(3, 3)] - m[M(3, 1)] * m[M(2, 3)])) +
727 (m[M(1, 3)] * (m[M(2, 1)] * m[M(3, 2)] - m[M(3, 1)] * m[M(2, 2)]));
728 coef[M(0, 1)] = (m[M(1, 0)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(3, 2)] * m[M(2, 3)])) -
729 (m[M(1, 2)] * (m[M(2, 0)] * m[M(3, 3)] - m[M(3, 0)] * m[M(2, 3)])) +
730 (m[M(1, 3)] * (m[M(2, 0)] * m[M(3, 2)] - m[M(3, 0)] * m[M(2, 2)]));
731 coef[M(0, 2)] = (m[M(1, 0)] * (m[M(2, 1)] * m[M(3, 3)] - m[M(3, 1)] * m[M(2, 3)])) -
732 (m[M(1, 1)] * (m[M(2, 0)] * m[M(3, 3)] - m[M(3, 0)] * m[M(2, 3)])) +
733 (m[M(1, 3)] * (m[M(2, 0)] * m[M(3, 1)] - m[M(3, 0)] * m[M(2, 1)]));
734 coef[M(0, 3)] = (m[M(1, 0)] * (m[M(2, 1)] * m[M(3, 2)] - m[M(3, 1)] * m[M(2, 2)])) -
735 (m[M(1, 1)] * (m[M(2, 0)] * m[M(3, 2)] - m[M(3, 0)] * m[M(2, 2)])) +
736 (m[M(1, 2)] * (m[M(2, 0)] * m[M(3, 1)] - m[M(3, 0)] * m[M(2, 1)]));
738 coef[M(1, 0)] = (m[M(0, 1)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(3, 2)] * m[M(2, 3)])) -
739 (m[M(0, 2)] * (m[M(2, 1)] * m[M(3, 3)] - m[M(3, 1)] * m[M(2, 3)])) +
740 (m[M(0, 3)] * (m[M(2, 1)] * m[M(3, 2)] - m[M(3, 1)] * m[M(2, 2)]));
741 coef[M(1, 1)] = (m[M(0, 0)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(3, 2)] * m[M(2, 3)])) -
742 (m[M(0, 2)] * (m[M(2, 0)] * m[M(3, 3)] - m[M(3, 0)] * m[M(2, 3)])) +
743 (m[M(0, 3)] * (m[M(2, 0)] * m[M(3, 2)] - m[M(3, 0)] * m[M(2, 2)]));
744 coef[M(1, 2)] = (m[M(0, 0)] * (m[M(2, 1)] * m[M(3, 3)] - m[M(3, 1)] * m[M(2, 3)])) -
745 (m[M(0, 1)] * (m[M(2, 0)] * m[M(3, 3)] - m[M(3, 0)] * m[M(2, 3)])) +
746 (m[M(0, 3)] * (m[M(2, 0)] * m[M(3, 1)] - m[M(3, 0)] * m[M(2, 1)]));
747 coef[M(1, 3)] = (m[M(0, 0)] * (m[M(2, 1)] * m[M(3, 2)] - m[M(3, 1)] * m[M(2, 2)])) -
748 (m[M(0, 1)] * (m[M(2, 0)] * m[M(3, 2)] - m[M(3, 0)] * m[M(2, 2)])) +
749 (m[M(0, 2)] * (m[M(2, 0)] * m[M(3, 1)] - m[M(3, 0)] * m[M(2, 1)]));
751 coef[M(2, 0)] = (m[M(0, 1)] * (m[M(1, 2)] * m[M(3, 3)] - m[M(3, 2)] * m[M(1, 3)])) -
752 (m[M(0, 2)] * (m[M(1, 1)] * m[M(3, 3)] - m[M(3, 1)] * m[M(1, 3)])) +
753 (m[M(0, 3)] * (m[M(1, 1)] * m[M(3, 2)] - m[M(3, 1)] * m[M(1, 2)]));
754 coef[M(2, 1)] = (m[M(0, 0)] * (m[M(1, 2)] * m[M(3, 3)] - m[M(3, 2)] * m[M(1, 3)])) -
755 (m[M(0, 2)] * (m[M(1, 0)] * m[M(3, 3)] - m[M(3, 0)] * m[M(1, 3)])) +
756 (m[M(0, 3)] * (m[M(1, 0)] * m[M(3, 2)] - m[M(3, 0)] * m[M(1, 2)]));
757 coef[M(2, 2)] = (m[M(0, 0)] * (m[M(1, 1)] * m[M(3, 3)] - m[M(3, 1)] * m[M(1, 3)])) -
758 (m[M(0, 1)] * (m[M(1, 0)] * m[M(3, 3)] - m[M(3, 0)] * m[M(1, 3)])) +
759 (m[M(0, 3)] * (m[M(1, 0)] * m[M(3, 1)] - m[M(3, 0)] * m[M(1, 1)]));
760 coef[M(2, 3)] = (m[M(0, 0)] * (m[M(1, 1)] * m[M(3, 2)] - m[M(3, 1)] * m[M(1, 2)])) -
761 (m[M(0, 1)] * (m[M(1, 0)] * m[M(3, 2)] - m[M(3, 0)] * m[M(1, 2)])) +
762 (m[M(0, 2)] * (m[M(1, 0)] * m[M(3, 1)] - m[M(3, 0)] * m[M(1, 1)]));
764 coef[M(3, 0)] = (m[M(0, 1)] * (m[M(1, 2)] * m[M(2, 3)] - m[M(2, 2)] * m[M(1, 3)])) -
765 (m[M(0, 2)] * (m[M(1, 1)] * m[M(2, 3)] - m[M(2, 1)] * m[M(1, 3)])) +
766 (m[M(0, 3)] * (m[M(1, 1)] * m[M(2, 2)] - m[M(2, 1)] * m[M(1, 2)]));
767 coef[M(3, 1)] = (m[M(0, 0)] * (m[M(1, 2)] * m[M(2, 3)] - m[M(2, 2)] * m[M(1, 3)])) -
768 (m[M(0, 2)] * (m[M(1, 0)] * m[M(2, 3)] - m[M(2, 0)] * m[M(1, 3)])) +
769 (m[M(0, 3)] * (m[M(1, 0)] * m[M(2, 2)] - m[M(2, 0)] * m[M(1, 2)]));
770 coef[M(3, 2)] = (m[M(0, 0)] * (m[M(1, 1)] * m[M(2, 3)] - m[M(2, 1)] * m[M(1, 3)])) -
771 (m[M(0, 1)] * (m[M(1, 0)] * m[M(2, 3)] - m[M(2, 0)] * m[M(1, 3)])) +
772 (m[M(0, 3)] * (m[M(1, 0)] * m[M(2, 1)] - m[M(2, 0)] * m[M(1, 1)]));
773 coef[M(3, 3)] = (m[M(0, 0)] * (m[M(1, 1)] * m[M(2, 2)] - m[M(2, 1)] * m[M(1, 2)])) -
774 (m[M(0, 1)] * (m[M(1, 0)] * m[M(2, 2)] - m[M(2, 0)] * m[M(1, 2)])) +
775 (m[M(0, 2)] * (m[M(1, 0)] * m[M(2, 1)] - m[M(2, 0)] * m[M(1, 1)]));
777 m4_transpose(res, coef);
779 for(i=0; i<4; i++) {
780 for(j=0; j<4; j++) {
781 res[M(i, j)] = j % 2 ? -res[M(i, j)] : res[M(i, j)];
782 if(i % 2) res[M(i, j)] = -res[M(i, j)];
783 }
784 }
785 }
787 static void m4_inverse(double *res, double *m)
788 {
789 int i, j;
790 double adj[16];
791 double det;
793 m4_adjoint(adj, m);
794 det = m4_determinant(m);
796 for(i=0; i<4; i++) {
797 for(j=0; j<4; j++) {
798 res[M(i, j)] = adj[M(i, j)] / det;
799 }
800 }
801 }