BRL-CAD
mat.c
Go to the documentation of this file.
1 /* M A T . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1996-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 
21 /** @addtogroup mat */
22 /** @{ */
23 /** @file libbn/mat.c
24  *
25  * TODO: need a better way to control tolerancing, either via
26  * additional tolerance parameters or perhaps providing a global
27  * tolerance state interface.
28  */
29 
30 
31 #include "common.h"
32 
33 #include <math.h>
34 #include <string.h>
35 #include "bio.h"
36 
37 #include "bu/debug.h"
38 #include "bu/log.h"
39 #include "bu/malloc.h"
40 #include "bu/str.h"
41 #include "vmath.h"
42 #include "bn/mat.h"
43 #include "bn/plane_calc.h"
44 
45 
46 const mat_t bn_mat_identity = MAT_INIT_IDN;
47 
48 
49 void
51  const char *title,
52  const mat_t m,
53  char *obuf,
54  int len)
55 {
56  register int i;
57  register char *cp;
58 
59  snprintf(obuf, len, "MATRIX %s:\n ", title);
60  cp = obuf + strlen(obuf);
61  if (!m) {
62  bu_strlcat(obuf, "(Identity)", len);
63  } else {
64  for (i = 0; i < 16; i++) {
65  snprintf(cp, len - (cp - obuf), " %8.3f", m[i]);
66  cp += strlen(cp);
67  if (i == 15) {
68  break;
69  } else if ((i & 3) == 3) {
70  *cp++ = '\n';
71  *cp++ = ' ';
72  *cp++ = ' ';
73  }
74  }
75  *cp++ = '\0';
76  }
77 }
78 
79 
80 void
81 bn_mat_print(const char *title, const mat_t m)
82 {
83  char obuf[1024]; /* snprintf may be non-PARALLEL */
84 
85  bn_mat_print_guts(title, m, obuf, 1024);
86  bu_log("%s\n", obuf);
87 }
88 
89 
90 void
92  const char *title,
93  const mat_t m,
94  struct bu_vls *vls)
95 {
96  char obuf[1024];
97 
98  bn_mat_print_guts(title, m, obuf, 1024);
99  bu_vls_printf(vls, "%s\n", obuf);
100 }
101 
102 
103 double
104 bn_atan2(double y, double x)
105 {
106  if (x > -1.0e-20 && x < 1.0e-20) {
107  /* X is equal to zero, check Y */
108  if (y < -1.0e-20) {
109  return -M_PI_2;
110  }
111  if (y > 1.0e-20) {
112  return M_PI_2;
113  }
114  return 0.0;
115  }
116  return atan2(y, x);
117 }
118 
119 
120 void
121 bn_mat_mul(register mat_t o, register const mat_t a, register const mat_t b)
122 {
123  o[ 0] = a[ 0] * b[ 0] + a[ 1] * b[ 4] + a[ 2] * b[ 8] + a[ 3] * b[12];
124  o[ 1] = a[ 0] * b[ 1] + a[ 1] * b[ 5] + a[ 2] * b[ 9] + a[ 3] * b[13];
125  o[ 2] = a[ 0] * b[ 2] + a[ 1] * b[ 6] + a[ 2] * b[10] + a[ 3] * b[14];
126  o[ 3] = a[ 0] * b[ 3] + a[ 1] * b[ 7] + a[ 2] * b[11] + a[ 3] * b[15];
127 
128  o[ 4] = a[ 4] * b[ 0] + a[ 5] * b[ 4] + a[ 6] * b[ 8] + a[ 7] * b[12];
129  o[ 5] = a[ 4] * b[ 1] + a[ 5] * b[ 5] + a[ 6] * b[ 9] + a[ 7] * b[13];
130  o[ 6] = a[ 4] * b[ 2] + a[ 5] * b[ 6] + a[ 6] * b[10] + a[ 7] * b[14];
131  o[ 7] = a[ 4] * b[ 3] + a[ 5] * b[ 7] + a[ 6] * b[11] + a[ 7] * b[15];
132 
133  o[ 8] = a[ 8] * b[ 0] + a[ 9] * b[ 4] + a[10] * b[ 8] + a[11] * b[12];
134  o[ 9] = a[ 8] * b[ 1] + a[ 9] * b[ 5] + a[10] * b[ 9] + a[11] * b[13];
135  o[10] = a[ 8] * b[ 2] + a[ 9] * b[ 6] + a[10] * b[10] + a[11] * b[14];
136  o[11] = a[ 8] * b[ 3] + a[ 9] * b[ 7] + a[10] * b[11] + a[11] * b[15];
137 
138  o[12] = a[12] * b[ 0] + a[13] * b[ 4] + a[14] * b[ 8] + a[15] * b[12];
139  o[13] = a[12] * b[ 1] + a[13] * b[ 5] + a[14] * b[ 9] + a[15] * b[13];
140  o[14] = a[12] * b[ 2] + a[13] * b[ 6] + a[14] * b[10] + a[15] * b[14];
141  o[15] = a[12] * b[ 3] + a[13] * b[ 7] + a[14] * b[11] + a[15] * b[15];
142 }
143 
144 
145 void
146 bn_mat_mul2(register const mat_t i, register mat_t o)
147 {
148  mat_t temp;
149 
150  bn_mat_mul(temp, i, o);
151  MAT_COPY(o, temp);
152 }
153 
154 
155 void
156 bn_mat_mul3(mat_t o, const mat_t a, const mat_t b, const mat_t c)
157 {
158  mat_t t;
159 
160  bn_mat_mul(t, b, c);
161  bn_mat_mul(o, a, t);
162 }
163 
164 
165 void
167  mat_t ao,
168  const mat_t a,
169  const mat_t b,
170  const mat_t c,
171  const mat_t d)
172 {
173  mat_t t, u;
174 
175  bn_mat_mul(u, c, d);
176  bn_mat_mul(t, a, b);
177  bn_mat_mul(ao, t, u);
178 }
179 
180 
181 void
183  register vect_t ov,
184  register const mat_t im,
185  register const vect_t iv)
186 {
187  register int eo = 0; /* Position in output vector */
188  register int em = 0; /* Position in input matrix */
189  register int ei; /* Position in input vector */
190 
191  /* For each element in the output array... */
192  for (; eo < 4; eo++) {
193 
194  ov[eo] = 0; /* Start with zero in output */
195 
196  for (ei = 0; ei < 4; ei++) {
197  ov[eo] += im[em++] * iv[ei];
198  }
199  }
200 }
201 
202 
203 void
204 bn_mat_inv(register mat_t output, const mat_t input)
205 {
206  if (bn_mat_inverse(output, input) == 0) {
207 
208  if (bu_debug & BU_DEBUG_MATH) {
209  bu_log("bn_mat_inv: error!");
210  bn_mat_print("singular matrix", input);
211  }
212  bu_bomb("bn_mat_inv: singular matrix\n");
213  /* NOTREACHED */
214  }
215 }
216 
217 
218 int
219 bn_mat_inverse(register mat_t output, const mat_t input)
220 {
221  register int i, j; /* Indices */
222  int k; /* Indices */
223  int z[4]; /* Temporary */
224  fastf_t b[4]; /* Temporary */
225  fastf_t c[4]; /* Temporary */
226 
227  MAT_COPY(output, input); /* Duplicate */
228 
229  /* Initialization */
230  for (j = 0; j < 4; j++) {
231  z[j] = j;
232  }
233 
234  /* Main Loop */
235  for (i = 0; i < 4; i++) {
236  register fastf_t y; /* local temporary */
237 
238  k = i;
239  y = output[i * 4 + i];
240  for (j = i + 1; j < 4; j++) {
241  register fastf_t w; /* local temporary */
242 
243  w = output[i * 4 + j];
244  if (fabs(w) > fabs(y)) {
245  k = j;
246  y = w;
247  }
248  }
249 
250  if (ZERO(y)) {
251  /* SINGULAR */
252  return 0;
253  }
254  y = 1.0 / y;
255 
256  for (j = 0; j < 4; j++) {
257  register fastf_t temp; /* Local */
258 
259  c[j] = output[j * 4 + k];
260  output[j * 4 + k] = output[j * 4 + i];
261  output[j * 4 + i] = - c[j] * y;
262  temp = output[i * 4 + j] * y;
263  b[j] = temp;
264  output[i * 4 + j] = temp;
265  }
266 
267  output[i * 4 + i] = y;
268  j = z[i];
269  z[i] = z[k];
270  z[k] = j;
271  for (k = 0; k < 4; k++) {
272  if (k == i) {
273  continue;
274  }
275  for (j = 0; j < 4; j++) {
276  if (j == i) {
277  continue;
278  }
279  output[k * 4 + j] = output[k * 4 + j] - b[j] * c[k];
280  }
281  }
282  }
283 
284  /* Second Loop */
285  for (i = 0; i < 4; i++) {
286  while ((k = z[i]) != i) {
287  int p; /* Local temp */
288 
289  for (j = 0; j < 4; j++) {
290  register fastf_t w; /* Local temp */
291 
292  w = output[i * 4 + j];
293  output[i * 4 + j] = output[k * 4 + j];
294  output[k * 4 + j] = w;
295  }
296  p = z[i];
297  z[i] = z[k];
298  z[k] = p;
299  }
300  }
301 
302  return 1;
303 }
304 
305 
306 void
307 bn_vtoh_move(register vect_t h, register const vect_t v)
308 {
309  VMOVE(h, v);
310  h[W] = 1.0;
311 }
312 
313 
314 void
315 bn_htov_move(register vect_t v, register const vect_t h)
316 {
317  register fastf_t inv;
318 
319  if (ZERO(h[3] - 1.0)) {
320  VMOVE(v, h);
321  } else {
322  if (ZERO(h[W])) {
323  bu_log("bn_htov_move: divide by %f!\n", h[W]);
324  return;
325  }
326  inv = 1.0 / h[W];
327  VSCALE(v, h, inv);
328  }
329 }
330 
331 
332 void
333 bn_mat_trn(mat_t om, register const mat_t im)
334 {
335  MAT_TRANSPOSE(om, im);
336 }
337 
338 
339 void
340 bn_mat_ae(register fastf_t *m, double azimuth, double elev)
341 {
342  double sin_az, sin_el;
343  double cos_az, cos_el;
344 
345  azimuth *= DEG2RAD;
346  elev *= DEG2RAD;
347 
348  sin_az = sin(azimuth);
349  cos_az = cos(azimuth);
350  sin_el = sin(elev);
351  cos_el = cos(elev);
352 
353  m[0] = cos_el * cos_az;
354  m[1] = -sin_az;
355  m[2] = -sin_el * cos_az;
356  m[3] = 0;
357 
358  m[4] = cos_el * sin_az;
359  m[5] = cos_az;
360  m[6] = -sin_el * sin_az;
361  m[7] = 0;
362 
363  m[8] = sin_el;
364  m[9] = 0;
365  m[10] = cos_el;
366  m[11] = 0;
367 
368  m[12] = m[13] = m[14] = 0;
369  m[15] = 1.0;
370 }
371 
372 
373 void
374 bn_ae_vec(fastf_t *azp, fastf_t *elp, const vect_t v)
375 {
376  register fastf_t az;
377 
378  if ((az = bn_atan2(v[Y], v[X]) * RAD2DEG) < 0) {
379  *azp = 360 + az;
380  } else if (az >= 360) {
381  *azp = az - 360;
382  } else {
383  *azp = az;
384  }
385  *elp = bn_atan2(v[Z], hypot(v[X], v[Y])) * RAD2DEG;
386 }
387 
388 
389 void
391  fastf_t *az,
392  fastf_t *el,
393  fastf_t *twist,
394  fastf_t *vec_ae,
395  fastf_t *vec_twist,
396  fastf_t accuracy)
397 {
398  vect_t zero_twist, ninety_twist;
399  vect_t z_dir;
400 
401  /* Get az and el as usual */
402  bn_ae_vec(az, el, vec_ae);
403 
404  /* stabilize fluctuation between 0 and 360
405  * change azimuth near 360 to 0 */
406  if (NEAR_EQUAL(*az, 360.0, accuracy)) {
407  *az = 0.0;
408  }
409 
410  /* if elevation is +/-90 set twist to zero and calculate azimuth */
411  if (NEAR_EQUAL(*el, 90.0, accuracy) || NEAR_ZERO(*el + 90.0, accuracy)) {
412  *twist = 0.0;
413  *az = bn_atan2(-vec_twist[X], vec_twist[Y]) * RAD2DEG;
414  } else {
415  /* Calculate twist from vec_twist */
416  VSET(z_dir, 0, 0, 1);
417  VCROSS(zero_twist, z_dir, vec_ae);
418  VUNITIZE(zero_twist);
419  VCROSS(ninety_twist, vec_ae, zero_twist);
420  VUNITIZE(ninety_twist);
421 
422  *twist = bn_atan2(VDOT(vec_twist, ninety_twist), VDOT(vec_twist, zero_twist)) * RAD2DEG;
423 
424  /* stabilize flutter between +/- 180 */
425  if (NEAR_EQUAL(*twist, -180.0, accuracy)) {
426  *twist = 180.0;
427  }
428  }
429 }
430 
431 
432 void
433 bn_vec_ae(vect_t vect, fastf_t az, fastf_t el)
434 {
435  fastf_t vx, vy, vz, rtemp;
436  vz = sin(el);
437  rtemp = cos(el);
438  vy = rtemp * sin(az);
439  vx = rtemp * cos(az);
440  VSET(vect, vx, vy , vz);
441  VUNITIZE(vect);
442 }
443 
444 
445 void
446 bn_vec_aed(vect_t vect, fastf_t az, fastf_t el, fastf_t distance)
447 {
448  fastf_t vx, vy, vz, rtemp;
449  vz = distance * sin(el);
450  rtemp = sqrt((distance * distance) - (vz * vz));
451  vy = rtemp * sin(az);
452  vx = rtemp * cos(az);
453  VSET(vect, vx, vy , vz);
454 }
455 
456 
457 void
459  register fastf_t *mat,
460  double alpha_in,
461  double beta_in,
462  double ggamma_in)
463 {
464  double alpha, beta, ggamma;
465  double calpha, cbeta, cgamma;
466  double salpha, sbeta, sgamma;
467 
468  if (NEAR_ZERO(alpha_in, 0.0) && NEAR_ZERO(beta_in, 0.0) && NEAR_ZERO(ggamma_in, 0.0)) {
469  MAT_IDN(mat);
470  return;
471  }
472 
473  alpha = alpha_in * DEG2RAD;
474  beta = beta_in * DEG2RAD;
475  ggamma = ggamma_in * DEG2RAD;
476 
477  calpha = cos(alpha);
478  cbeta = cos(beta);
479  cgamma = cos(ggamma);
480 
481  /* sine of "180*DEG2RAD" will not be exactly zero and will
482  * result in errors when some codes try to convert this back to
483  * azimuth and elevation. do_frame() uses this technique!!!
484  */
485  if (ZERO(alpha_in - 180.0)) {
486  salpha = 0.0;
487  } else {
488  salpha = sin(alpha);
489  }
490 
491  if (ZERO(beta_in - 180.0)) {
492  sbeta = 0.0;
493  } else {
494  sbeta = sin(beta);
495  }
496 
497  if (ZERO(ggamma_in - 180.0)) {
498  sgamma = 0.0;
499  } else {
500  sgamma = sin(ggamma);
501  }
502 
503  mat[0] = cbeta * cgamma;
504  mat[1] = -cbeta * sgamma;
505  mat[2] = sbeta;
506  mat[3] = 0.0;
507 
508  mat[4] = salpha * sbeta * cgamma + calpha * sgamma;
509  mat[5] = -salpha * sbeta * sgamma + calpha * cgamma;
510  mat[6] = -salpha * cbeta;
511  mat[7] = 0.0;
512 
513  mat[8] = salpha * sgamma - calpha * sbeta * cgamma;
514  mat[9] = salpha * cgamma + calpha * sbeta * sgamma;
515  mat[10] = calpha * cbeta;
516  mat[11] = 0.0;
517  mat[12] = mat[13] = mat[14] = 0.0;
518  mat[15] = 1.0;
519 }
520 
521 
522 void
524  register mat_t mat,
525  double alpha,
526  double beta,
527  double ggamma)
528 {
529  double calpha, cbeta, cgamma;
530  double salpha, sbeta, sgamma;
531 
532  if (ZERO(alpha) && ZERO(beta) && ZERO(ggamma)) {
533  MAT_IDN(mat);
534  return;
535  }
536 
537  calpha = cos(alpha);
538  cbeta = cos(beta);
539  cgamma = cos(ggamma);
540 
541  salpha = sin(alpha);
542  sbeta = sin(beta);
543  sgamma = sin(ggamma);
544 
545  mat[0] = cbeta * cgamma;
546  mat[1] = -cbeta * sgamma;
547  mat[2] = sbeta;
548  mat[3] = 0.0;
549 
550  mat[4] = salpha * sbeta * cgamma + calpha * sgamma;
551  mat[5] = -salpha * sbeta * sgamma + calpha * cgamma;
552  mat[6] = -salpha * cbeta;
553  mat[7] = 0.0;
554 
555  mat[8] = salpha * sgamma - calpha * sbeta * cgamma;
556  mat[9] = salpha * cgamma + calpha * sbeta * sgamma;
557  mat[10] = calpha * cbeta;
558  mat[11] = 0.0;
559  mat[12] = mat[13] = mat[14] = 0.0;
560  mat[15] = 1.0;
561 }
562 
563 
564 void
566  fastf_t *val1,
567  fastf_t *val2,
568  fastf_t *vec1,
569  fastf_t *vec2,
570  fastf_t a,
571  fastf_t b,
572  fastf_t c)
573 {
574  fastf_t d, root;
575  fastf_t v1, v2;
576 
577  d = 0.5 * (c - a);
578 
579  /* Check for diagonal matrix */
580  if (NEAR_ZERO(b, 1.0e-10)) {
581  /* smaller mag first */
582  if (fabs(c) < fabs(a)) {
583  *val1 = c;
584  VSET(vec1, 0.0, 1.0, 0.0);
585  *val2 = a;
586  VSET(vec2, -1.0, 0.0, 0.0);
587  } else {
588  *val1 = a;
589  VSET(vec1, 1.0, 0.0, 0.0);
590  *val2 = c;
591  VSET(vec2, 0.0, 1.0, 0.0);
592  }
593  return;
594  }
595 
596  root = sqrt(d * d + b * b);
597  v1 = 0.5 * (c + a) - root;
598  v2 = 0.5 * (c + a) + root;
599 
600  /* smaller mag first */
601  if (fabs(v1) < fabs(v2)) {
602  *val1 = v1;
603  *val2 = v2;
604  VSET(vec1, b, d - root, 0.0);
605  } else {
606  *val1 = v2;
607  *val2 = v1;
608  VSET(vec1, root - d, b, 0.0);
609  }
610  VUNITIZE(vec1);
611  VSET(vec2, -vec1[Y], vec1[X], 0.0); /* vec1 X vec2 = +Z */
612 }
613 
614 
615 void
616 bn_vec_perp(vect_t new_vec, const vect_t old_vec)
617 {
618  register int i;
619  vect_t another_vec; /* Another vector, different */
620 
621  i = X;
622  if (fabs(old_vec[Y]) < fabs(old_vec[i])) {
623  i = Y;
624  }
625  if (fabs(old_vec[Z]) < fabs(old_vec[i])) {
626  i = Z;
627  }
628  VSETALL(another_vec, 0);
629  another_vec[i] = 1.0;
630  if (ZERO(old_vec[X]) && ZERO(old_vec[Y]) && ZERO(old_vec[Z])) {
631  VMOVE(new_vec, another_vec);
632  } else {
633  VCROSS(new_vec, another_vec, old_vec);
634  }
635 }
636 
637 
638 void
640  mat_t m,
641  const fastf_t *from,
642  const fastf_t *to,
643  const struct bn_tol *tol)
644 {
645  vect_t test_to;
646  vect_t unit_from, unit_to;
647  fastf_t dot;
648  point_t origin = VINIT_ZERO;
649 
650  /**
651  * The method used here is from Graphics Gems, A. Glassner, ed.
652  * page 531, "The Use of Coordinate Frames in Computer Graphics",
653  * by Ken Turkowski, Example 6.
654  */
655  mat_t Q, Qt;
656  mat_t R;
657  mat_t A;
658  mat_t temp;
659  vect_t N, M;
660  vect_t w_prime; /* image of "to" ("w") in Qt */
661 
662  VMOVE(unit_from, from);
663  VUNITIZE(unit_from); /* aka "v" */
664  VMOVE(unit_to, to);
665  VUNITIZE(unit_to); /* aka "w" */
666 
667  /* If from and to are the same or opposite, special handling is
668  * needed, because the cross product isn't defined. asin(0.00001)
669  * = 0.0005729 degrees (1/2000 degree)
670  */
671  if (bn_lseg3_lseg3_parallel(origin, from, origin, to, tol)) {
672  dot = VDOT(from, to);
673  if (dot > SMALL_FASTF) {
674  MAT_IDN(m);
675  return;
676  } else {
677  bn_vec_perp(N, unit_from);
678  }
679  } else {
680  VCROSS(N, unit_from, unit_to);
681  VUNITIZE(N); /* should be unnecessary */
682  }
683  VCROSS(M, N, unit_from);
684  VUNITIZE(M); /* should be unnecessary */
685 
686  /* Almost everything here is done with pre-multiplys: vector * mat */
687  MAT_IDN(Q);
688  VMOVE(&Q[0], unit_from);
689  VMOVE(&Q[4], M);
690  VMOVE(&Q[8], N);
691  bn_mat_trn(Qt, Q);
692 
693  /* w_prime = w * Qt */
694  MAT4X3VEC(w_prime, Q, unit_to); /* post-multiply by transpose */
695 
696  MAT_IDN(R);
697  VMOVE(&R[0], w_prime);
698  VSET(&R[4], -w_prime[Y], w_prime[X], w_prime[Z]);
699  VSET(&R[8], 0, 0, 1); /* is unnecessary */
700 
701  bn_mat_mul(temp, R, Q);
702  bn_mat_mul(A, Qt, temp);
703  bn_mat_trn(m, A); /* back to post-multiply style */
704 
705  /* Verify that it worked */
706  MAT4X3VEC(test_to, m, unit_from);
707  if (UNLIKELY(!bn_lseg3_lseg3_parallel(origin, unit_to, origin, test_to, tol))) {
708  dot = VDOT(unit_to, test_to);
709  bu_log("bn_mat_fromto() ERROR! from (%g, %g, %g) to (%g, %g, %g) went to (%g, %g, %g), dot=%g?\n",
710  V3ARGS(from), V3ARGS(to), V3ARGS(test_to), dot);
711  }
712 }
713 
714 
715 void
716 bn_mat_xrot(fastf_t *m, double sinx, double cosx)
717 {
718  m[0] = 1.0;
719  m[1] = 0.0;
720  m[2] = 0.0;
721  m[3] = 0.0;
722 
723  m[4] = 0.0;
724  m[5] = cosx;
725  m[6] = -sinx;
726  m[7] = 0.0;
727 
728  m[8] = 0.0;
729  m[9] = sinx;
730  m[10] = cosx;
731  m[11] = 0.0;
732 
733  m[12] = m[13] = m[14] = 0.0;
734  m[15] = 1.0;
735 }
736 
737 
738 void
739 bn_mat_yrot(fastf_t *m, double siny, double cosy)
740 {
741  m[0] = cosy;
742  m[1] = 0.0;
743  m[2] = -siny;
744  m[3] = 0.0;
745 
746  m[4] = 0.0;
747  m[5] = 1.0;
748  m[6] = 0.0;
749  m[7] = 0.0;
750 
751  m[8] = siny;
752  m[9] = 0.0;
753  m[10] = cosy;
754 
755  m[11] = 0.0;
756  m[12] = m[13] = m[14] = 0.0;
757  m[15] = 1.0;
758 }
759 
760 
761 void
762 bn_mat_zrot(fastf_t *m, double sinz, double cosz)
763 {
764  m[0] = cosz;
765  m[1] = -sinz;
766  m[2] = 0.0;
767  m[3] = 0.0;
768 
769  m[4] = sinz;
770  m[5] = cosz;
771  m[6] = 0.0;
772  m[7] = 0.0;
773 
774  m[8] = 0.0;
775  m[9] = 0.0;
776  m[10] = 1.0;
777  m[11] = 0.0;
778 
779  m[12] = m[13] = m[14] = 0.0;
780  m[15] = 1.0;
781 }
782 
783 
784 void
785 bn_mat_lookat(mat_t rot, const vect_t dir, int yflip)
786 {
787  mat_t first;
788  mat_t second;
789  mat_t prod12;
790  mat_t third;
791  vect_t x;
792  vect_t z;
793  vect_t t1;
794  fastf_t hypot_xy;
795  vect_t xproj;
796  vect_t zproj;
797 
798  /* First, rotate D around Z axis to match +X axis (azimuth) */
799  hypot_xy = hypot(dir[X], dir[Y]);
800  bn_mat_zrot(first, -dir[Y] / hypot_xy, dir[X] / hypot_xy);
801 
802  /* Next, rotate D around Y axis to match -Z axis (elevation) */
803  bn_mat_yrot(second, -hypot_xy, -dir[Z]);
804  bn_mat_mul(prod12, second, first);
805 
806  /* Produce twist correction, by re-orienting projection of X axis */
807  VSET(x, 1, 0, 0);
808  MAT4X3VEC(xproj, prod12, x);
809  hypot_xy = hypot(xproj[X], xproj[Y]);
810  if (hypot_xy < 1.0e-10) {
811  bu_log("Warning: bn_mat_lookat: unable to twist correct, hypot=%g\n", hypot_xy);
812  VPRINT("xproj", xproj);
813  MAT_COPY(rot, prod12);
814  return;
815  }
816  bn_mat_zrot(third, -xproj[Y] / hypot_xy, xproj[X] / hypot_xy);
817  bn_mat_mul(rot, third, prod12);
818 
819  if (yflip) {
820  VSET(z, 0, 0, 1);
821  MAT4X3VEC(zproj, rot, z);
822  /* If original Z inverts sign, flip sign on resulting Y */
823  if (zproj[Y] < 0.0) {
824  MAT_COPY(prod12, rot);
825  MAT_IDN(third);
826  third[5] = -1;
827  bn_mat_mul(rot, third, prod12);
828  }
829  }
830 
831  /* Check the final results */
832  MAT4X3VEC(t1, rot, dir);
833  if (t1[Z] > -0.98) {
834  bu_log("Error: bn_mat_lookat final= (%g, %g, %g)\n", V3ARGS(t1));
835  }
836 }
837 
838 
839 void
840 bn_vec_ortho(register vect_t out, register const vect_t in)
841 {
842  register int j, k;
843  register fastf_t f;
844  register int i;
845 
846  if (UNLIKELY(NEAR_ZERO(MAGSQ(in), SQRT_SMALL_FASTF))) {
847  bu_log("bn_vec_ortho(): zero magnitude input vector %g %g %g\n", V3ARGS(in));
848  VSETALL(out, 0);
849  return;
850  }
851 
852  /* Find component closest to zero */
853  f = fabs(in[X]);
854  i = X;
855  j = Y;
856  k = Z;
857  if (fabs(in[Y]) < f) {
858  f = fabs(in[Y]);
859  i = Y;
860  j = Z;
861  k = X;
862  }
863  if (fabs(in[Z]) < f) {
864  i = Z;
865  j = X;
866  k = Y;
867  }
868  f = hypot(in[j], in[k]);
869  if (UNLIKELY(ZERO(f))) {
870  bu_log("bn_vec_ortho(): zero hypot on %g %g %g\n", V3ARGS(in));
871  VSETALL(out, 0);
872  return;
873  }
874  f = 1.0 / f;
875  out[i] = 0.0;
876  out[j] = -in[k] * f;
877  out[k] = in[j] * f;
878 
879  return;
880 }
881 
882 
883 int
884 bn_mat_scale_about_pt(mat_t mat, const point_t pt, const double scale)
885 {
886  mat_t xlate;
887  mat_t s;
888  mat_t tmp;
889 
890  MAT_IDN(xlate);
891  MAT_DELTAS_VEC_NEG(xlate, pt);
892 
893  MAT_IDN(s);
894  if (ZERO(scale)) {
895  MAT_ZERO(mat);
896  return -1; /* ERROR */
897  }
898  s[15] = 1 / scale;
899 
900  bn_mat_mul(tmp, s, xlate);
901 
902  MAT_DELTAS_VEC(xlate, pt);
903  bn_mat_mul(mat, xlate, tmp);
904  return 0; /* OK */
905 }
906 
907 
908 void
909 bn_mat_xform_about_pt(mat_t mat, const mat_t xform, const point_t pt)
910 {
911  mat_t xlate;
912  mat_t tmp;
913 
914  MAT_IDN(xlate);
915  MAT_DELTAS_VEC_NEG(xlate, pt);
916 
917  bn_mat_mul(tmp, xform, xlate);
918 
919  MAT_DELTAS_VEC(xlate, pt);
920  bn_mat_mul(mat, xlate, tmp);
921 }
922 
923 
924 int
925 bn_mat_is_equal(const mat_t a, const mat_t b, const struct bn_tol *tol)
926 {
927  register int i;
928  register double f;
929  register double tdist, tperp;
930 
931  BN_CK_TOL(tol);
932 
933  tdist = tol->dist;
934  tperp = tol->perp;
935 
936  /* First, check that the translation part of the matrix (dist) are
937  * within the distance tolerance. Because most instancing
938  * involves translation and no rotation, doing this first should
939  * detect most non-equal cases rapidly.
940  */
941  for (i = 3; i < 12; i += 4) {
942  f = a[i] - b[i];
943  if (!NEAR_ZERO(f, tdist)) {
944  return 0;
945  }
946  }
947 
948  /* Check that the rotation part of the matrix (cosines) are within
949  * the perpendicular tolerance.
950  */
951  for (i = 0; i < 16; i += 4) {
952  f = a[i] - b[i];
953  if (!NEAR_ZERO(f, tperp)) {
954  return 0;
955  }
956  f = a[i + 1] - b[i + 1];
957  if (!NEAR_ZERO(f, tperp)) {
958  return 0;
959  }
960  f = a[i + 2] - b[i + 2];
961  if (!NEAR_ZERO(f, tperp)) {
962  return 0;
963  }
964  }
965 
966  /* Check that the scale part of the matrix (ratio) is within the
967  * perpendicular tolerance. There is no ratio tolerance so we use
968  * the tighter of dist or perp.
969  */
970  f = a[15] - b[15];
971  if (!NEAR_ZERO(f, tperp)) {
972  return 0;
973  }
974 
975  return 1;
976 }
977 
978 
979 int
980 bn_mat_is_identity(const mat_t m)
981 {
982  return ! memcmp(m, bn_mat_identity, sizeof(mat_t));
983 }
984 
985 
986 void
988  mat_t m,
989  const point_t pt,
990  const vect_t dir,
991  const fastf_t ang)
992 {
993  mat_t tran1, tran2, rot;
994  double cos_ang, sin_ang, one_m_cosang;
995  double n1_sq, n2_sq, n3_sq;
996  double n1_n2, n1_n3, n2_n3;
997 
998  if (ZERO(ang)) {
999  MAT_IDN(m);
1000  return;
1001  }
1002 
1003  MAT_IDN(tran1);
1004  MAT_IDN(tran2);
1005 
1006  /* construct translation matrix to pt */
1007  tran1[MDX] = (-pt[X]);
1008  tran1[MDY] = (-pt[Y]);
1009  tran1[MDZ] = (-pt[Z]);
1010 
1011  /* construct translation back from pt */
1012  tran2[MDX] = pt[X];
1013  tran2[MDY] = pt[Y];
1014  tran2[MDZ] = pt[Z];
1015 
1016  /* construct rotation matrix */
1017  cos_ang = cos(ang);
1018  sin_ang = sin(ang);
1019  one_m_cosang = 1.0 - cos_ang;
1020  n1_sq = dir[X] * dir[X];
1021  n2_sq = dir[Y] * dir[Y];
1022  n3_sq = dir[Z] * dir[Z];
1023  n1_n2 = dir[X] * dir[Y];
1024  n1_n3 = dir[X] * dir[Z];
1025  n2_n3 = dir[Y] * dir[Z];
1026 
1027  MAT_IDN(rot);
1028  rot[0] = n1_sq + (1.0 - n1_sq) * cos_ang;
1029  rot[1] = n1_n2 * one_m_cosang - dir[Z] * sin_ang;
1030  rot[2] = n1_n3 * one_m_cosang + dir[Y] * sin_ang;
1031 
1032  rot[4] = n1_n2 * one_m_cosang + dir[Z] * sin_ang;
1033  rot[5] = n2_sq + (1.0 - n2_sq) * cos_ang;
1034  rot[6] = n2_n3 * one_m_cosang - dir[X] * sin_ang;
1035 
1036  rot[8] = n1_n3 * one_m_cosang - dir[Y] * sin_ang;
1037  rot[9] = n2_n3 * one_m_cosang + dir[X] * sin_ang;
1038  rot[10] = n3_sq + (1.0 - n3_sq) * cos_ang;
1039 
1040  bn_mat_mul(m, rot, tran1);
1041  bn_mat_mul2(tran2, m);
1042 }
1043 
1044 
1045 matp_t
1046 bn_mat_dup(const mat_t in)
1047 {
1048  mat_t *out;
1049 
1050  BU_ALLOC(out, mat_t);
1051  MAT_COPY(*out, in);
1052 
1053  return (matp_t)out;
1054 }
1055 
1056 
1057 int
1058 bn_mat_ck(const char *title, const mat_t m)
1059 {
1060  vect_t A, B, C;
1061  fastf_t fx, fy, fz;
1062 
1063  if (!m) {
1064  return 0; /* implies identity matrix */
1065  }
1066 
1067  /* Validate that matrix preserves perpendicularity of axis by
1068  * checking that A.B == 0, B.C == 0, A.C == 0 XXX these vectors
1069  * should just be grabbed out of the matrix
1070  */
1071  VMOVE(A, &m[0]);
1072  VMOVE(B, &m[4]);
1073  VMOVE(C, &m[8]);
1074 
1075  fx = VDOT(A, B);
1076  fy = VDOT(B, C);
1077  fz = VDOT(A, C);
1078 
1079  /* NOTE: this tolerance cannot be any more tight than 0.00001 due
1080  * to default calculation tolerancing used by models. Matrices
1081  * exported to disk outside of tolerance will fail import if
1082  * set too restrictive.
1083  */
1084  if (!NEAR_ZERO(fx, 0.00001)
1085  || !NEAR_ZERO(fy, 0.00001)
1086  || !NEAR_ZERO(fz, 0.00001)
1087  || NEAR_ZERO(m[15], VDIVIDE_TOL)) {
1088  if (bu_debug & BU_DEBUG_MATH) {
1089  bu_log("bn_mat_ck(%s): bad matrix, does not preserve axis perpendicularity.\n X.Y=%g, Y.Z=%g, X.Z=%g, s=%g\n", title, fx, fy, fz, m[15]);
1090  bn_mat_print("bn_mat_ck() bad matrix", m);
1091  }
1092 
1093  if (bu_debug & (BU_DEBUG_MATH | BU_DEBUG_COREDUMP)) {
1095  bu_bomb("bn_mat_ck() bad matrix\n");
1096  }
1097  return -1; /* FAIL */
1098  }
1099  return 0; /* OK */
1100 }
1101 
1102 
1103 fastf_t
1104 bn_mat_det3(const mat_t m)
1105 {
1106  register fastf_t sum;
1107 
1108  sum = m[0] * (m[5] * m[10] - m[6] * m[9])
1109  - m[1] * (m[4] * m[10] - m[6] * m[8])
1110  + m[2] * (m[4] * m[9] - m[5] * m[8]);
1111 
1112  return sum;
1113 }
1114 
1115 
1116 fastf_t
1117 bn_mat_determinant(const mat_t m)
1118 {
1119  fastf_t det[4];
1120  fastf_t sum;
1121 
1122  det[0] = m[5] * (m[10] * m[15] - m[11] * m[14])
1123  - m[6] * (m[ 9] * m[15] - m[11] * m[13])
1124  + m[7] * (m[ 9] * m[14] - m[10] * m[13]);
1125 
1126  det[1] = m[4] * (m[10] * m[15] - m[11] * m[14])
1127  - m[6] * (m[ 8] * m[15] - m[11] * m[12])
1128  + m[7] * (m[ 8] * m[14] - m[10] * m[12]);
1129 
1130  det[2] = m[4] * (m[ 9] * m[15] - m[11] * m[13])
1131  - m[5] * (m[ 8] * m[15] - m[11] * m[12])
1132  + m[7] * (m[ 8] * m[13] - m[ 9] * m[12]);
1133 
1134  det[3] = m[4] * (m[ 9] * m[14] - m[10] * m[13])
1135  - m[5] * (m[ 8] * m[14] - m[10] * m[12])
1136  + m[6] * (m[ 8] * m[13] - m[ 9] * m[12]);
1137 
1138  sum = m[0] * det[0] - m[1] * det[1] + m[2] * det[2] - m[3] * det[3];
1139 
1140  return sum;
1141 
1142 }
1143 
1144 
1145 int
1146 bn_mat_is_non_unif(const mat_t m)
1147 {
1148  double mag[3];
1149 
1150  mag[0] = MAGSQ(m);
1151  mag[1] = MAGSQ(&m[4]);
1152  mag[2] = MAGSQ(&m[8]);
1153 
1154  if (fabs(1.0 - (mag[1] / mag[0])) > .0005 ||
1155  fabs(1.0 - (mag[2] / mag[0])) > .0005) {
1156 
1157  return 1;
1158  }
1159 
1160  if (!ZERO(m[12]) || !ZERO(m[13]) || !ZERO(m[14])) {
1161  return 2;
1162  }
1163 
1164  return 0;
1165 }
1166 
1167 
1168 void
1169 bn_wrt_point_direc(mat_t out, const mat_t change, const mat_t in, const point_t point, const vect_t direc, const struct bn_tol *tol)
1170 {
1171  static mat_t t1;
1172  static mat_t pt_to_origin, origin_to_pt;
1173  static mat_t d_to_zaxis, zaxis_to_d;
1174  static vect_t zaxis;
1175 
1176  /* build "point to origin" matrix */
1177  MAT_IDN(pt_to_origin);
1178  MAT_DELTAS_VEC_NEG(pt_to_origin, point);
1179 
1180  /* build "origin to point" matrix */
1181  MAT_IDN(origin_to_pt);
1182  MAT_DELTAS_VEC_NEG(origin_to_pt, point);
1183 
1184  /* build "direc to zaxis" matrix */
1185  VSET(zaxis, 0.0, 0.0, 1.0);
1186  bn_mat_fromto(d_to_zaxis, direc, zaxis, tol);
1187 
1188  /* build "zaxis to direc" matrix */
1189  bn_mat_inv(zaxis_to_d, d_to_zaxis);
1190 
1191  /* apply change matrix...
1192  * t1 = change * d_to_zaxis * pt_to_origin * in
1193  */
1194  bn_mat_mul4(t1, change, d_to_zaxis, pt_to_origin, in);
1195 
1196  /* apply origin_to_pt matrix:
1197  * out = origin_to_pt * zaxis_to_d *
1198  * change * d_to_zaxis * pt_to_origin * in
1199  */
1200  bn_mat_mul3(out, origin_to_pt, zaxis_to_d, t1);
1201 }
1202 
1203 /*
1204  * Compute a perspective matrix for a right-handed coordinate system.
1205  * Reference: SGI Graphics Reference Appendix C
1206  * (Note: SGI is left-handed, but the fix is done in the Display Manger).
1207  */
1208 void
1209 persp_mat(mat_t m, fastf_t fovy, fastf_t aspect, fastf_t near1, fastf_t far1, fastf_t backoff)
1210 {
1211  mat_t m2, tra;
1212 
1213  fovy *= DEG2RAD;
1214 
1215  MAT_IDN(m2);
1216  m2[5] = cos(fovy/2.0) / sin(fovy/2.0);
1217  m2[0] = m2[5]/aspect;
1218  m2[10] = (far1+near1) / (far1-near1);
1219  m2[11] = 2*far1*near1 / (far1-near1); /* This should be negative */
1220 
1221  m2[14] = -1; /* XXX This should be positive */
1222  m2[15] = 0;
1223 
1224  /* Move eye to origin, then apply perspective */
1225  MAT_IDN(tra);
1226  tra[11] = -backoff;
1227  bn_mat_mul(m, m2, tra);
1228 }
1229 
1230 /*
1231  *
1232  * Create a perspective matrix that transforms the +/1 viewing cube,
1233  * with the actual eye position (not at Z=+1) specified in viewing coords,
1234  * into a related space where the eye has been sheared onto the Z axis
1235  * and repositioned at Z=(0, 0, 1), with the same perspective field of view
1236  * as before.
1237  *
1238  * The Zbuffer clips off stuff with negative Z values.
1239  *
1240  * pmat = persp * xlate * shear
1241  */
1242 void
1243 mike_persp_mat(fastf_t *pmat, const fastf_t *eye)
1244 {
1245  mat_t shear;
1246  mat_t persp;
1247  mat_t xlate;
1248  mat_t t1, t2;
1249  point_t sheared_eye;
1250 
1251  if (eye[Z] < SMALL) {
1252  VPRINT("mike_persp_mat(): ERROR, z<0, eye", eye);
1253  return;
1254  }
1255 
1256  /* Shear "eye" to +Z axis */
1257  MAT_IDN(shear);
1258  shear[2] = -eye[X]/eye[Z];
1259  shear[6] = -eye[Y]/eye[Z];
1260 
1261  MAT4X3VEC(sheared_eye, shear, eye);
1262  if (!NEAR_ZERO(sheared_eye[X], .01) || !NEAR_ZERO(sheared_eye[Y], .01)) {
1263  VPRINT("ERROR sheared_eye", sheared_eye);
1264  return;
1265  }
1266 
1267  /* Translate along +Z axis to put sheared_eye at (0, 0, 1). */
1268  MAT_IDN(xlate);
1269  /* XXX should I use MAT_DELTAS_VEC_NEG()? X and Y should be 0 now */
1270  MAT_DELTAS(xlate, 0, 0, 1-sheared_eye[Z]);
1271 
1272  /* Build perspective matrix in place, substituting fov=2*atan(1, Z) */
1273  MAT_IDN(persp);
1274  /* From page 492 of Graphics Gems */
1275  persp[0] = sheared_eye[Z]; /* scaling: fov aspect term */
1276  persp[5] = sheared_eye[Z]; /* scaling: determines fov */
1277 
1278  /* From page 158 of Rogers Mathematical Elements */
1279  /* Z center of projection at Z=+1, r=-1/1 */
1280  persp[14] = -1;
1281 
1282  bn_mat_mul(t1, xlate, shear);
1283  bn_mat_mul(t2, persp, t1);
1284  /* Now, move eye from Z=1 to Z=0, for clipping purposes */
1285  MAT_DELTAS(xlate, 0, 0, -1);
1286  bn_mat_mul(pmat, xlate, t2);
1287 }
1288 
1289 
1290 /*
1291  * Map "display plate coordinates" (which can just be the screen viewing cube),
1292  * into [-1, +1] coordinates, with perspective.
1293  * Per "High Resolution Virtual Reality" by Michael Deering,
1294  * Computer Graphics 26, 2, July 1992, pp 195-201.
1295  *
1296  * L is lower left corner of screen, H is upper right corner.
1297  * L[Z] is the front (near) clipping plane location.
1298  * H[Z] is the back (far) clipping plane location.
1299  *
1300  * This corresponds to the SGI "window()" routine, but taking into account
1301  * skew due to the eyepoint being offset parallel to the image plane.
1302  *
1303  * The gist of the algorithm is to translate the display plate to the
1304  * view center, shear the eye point to (0, 0, 1), translate back,
1305  * then apply an off-axis perspective projection.
1306  *
1307  * Another (partial) reference is "A comparison of stereoscopic cursors
1308  * for the interactive manipulation of B-splines" by Barham & McAllister,
1309  * SPIE Vol 1457 Stereoscopic Display & Applications, 1991, pg 19.
1310  */
1311 void
1312 deering_persp_mat(fastf_t *m, const fastf_t *l, const fastf_t *h, const fastf_t *eye)
1313  /* lower left corner of screen */
1314  /* upper right (high) corner of screen */
1315  /* eye location. Traditionally at (0, 0, 1) */
1316 {
1317  vect_t diff; /* H - L */
1318  vect_t sum; /* H + L */
1319 
1320  VSUB2(diff, h, l);
1321  VADD2(sum, h, l);
1322 
1323  m[0] = 2 * eye[Z] / diff[X];
1324  m[1] = 0;
1325  m[2] = (sum[X] - 2 * eye[X]) / diff[X];
1326  m[3] = -eye[Z] * sum[X] / diff[X];
1327 
1328  m[4] = 0;
1329  m[5] = 2 * eye[Z] / diff[Y];
1330  m[6] = (sum[Y] - 2 * eye[Y]) / diff[Y];
1331  m[7] = -eye[Z] * sum[Y] / diff[Y];
1332 
1333  /* Multiplied by -1, to do right-handed Z coords */
1334  m[8] = 0;
1335  m[9] = 0;
1336  m[10] = -(sum[Z] - 2 * eye[Z]) / diff[Z];
1337  m[11] = -(-eye[Z] + 2 * h[Z] * eye[Z]) / diff[Z];
1338 
1339  m[12] = 0;
1340  m[13] = 0;
1341  m[14] = -1;
1342  m[15] = eye[Z];
1343 
1344 /* XXX May need to flip Z ? (lefthand to righthand?) */
1345 }
1346 
1347 
1348 
1349 /** @} */
1350 /*
1351  * Local Variables:
1352  * mode: C
1353  * tab-width: 8
1354  * indent-tabs-mode: t
1355  * c-file-style: "stroustrup"
1356  * End:
1357  * ex: shiftwidth=4 tabstop=8
1358  */
void bn_ae_vec(fastf_t *azp, fastf_t *elp, const vect_t v)
Definition: mat.c:374
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
int bn_mat_is_identity(const mat_t m)
Definition: mat.c:980
fastf_t C[2 *MAX_CNT+1][2 *MAX_CNT+1]
Definition: dsp_brep.cpp:38
void bn_mat_lookat(mat_t rot, const vect_t dir, int yflip)
Definition: mat.c:785
fastf_t bn_mat_determinant(const mat_t m)
Definition: mat.c:1117
#define SMALL
Definition: defines.h:351
void bn_mat_inv(register mat_t output, const mat_t input)
Definition: mat.c:204
double dist
>= 0
Definition: tol.h:73
const mat_t bn_mat_identity
Matrix and vector functionality.
Definition: mat.c:46
if lu s
Definition: nmg_mod.c:3860
#define VSET(a, b, c, d)
Definition: color.c:53
#define VSETALL(a, s)
Definition: color.c:54
#define N
Definition: randmt.c:39
void bn_mat_ae(register fastf_t *m, double azimuth, double elev)
Definition: mat.c:340
#define SMALL_FASTF
Definition: defines.h:342
void bn_mat_angles(register fastf_t *mat, double alpha_in, double beta_in, double ggamma_in)
Definition: mat.c:458
Header file for the BRL-CAD common definitions.
void bn_mat_xrot(fastf_t *m, double sinx, double cosx)
Definition: mat.c:716
void bn_mat_mul4(mat_t o, const mat_t a, const mat_t b, const mat_t c, const mat_t d)
Definition: mat.c:166
fastf_t bn_mat_det3(const mat_t m)
Definition: mat.c:1104
void bn_vec_ortho(register vect_t out, register const vect_t in)
Definition: mat.c:840
int bn_mat_is_equal(const mat_t a, const mat_t b, const struct bn_tol *tol)
Definition: mat.c:925
Definition: color.c:49
void deering_persp_mat(fastf_t *m, const fastf_t *l, const fastf_t *h, const fastf_t *eye)
Definition: mat.c:1312
void bn_mat_angles_rad(register mat_t mat, double alpha, double beta, double ggamma)
Definition: mat.c:523
int bn_mat_is_non_unif(const mat_t m)
Definition: mat.c:1146
void bn_mat_print(const char *title, const mat_t m)
Definition: mat.c:81
int bn_mat_scale_about_pt(mat_t mat, const point_t pt, const double scale)
Definition: mat.c:884
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
void bn_vec_ae(vect_t vec, fastf_t az, fastf_t el)
Definition: mat.c:433
#define V3ARGS(a)
Definition: color.c:56
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
#define SQRT_SMALL_FASTF
Definition: defines.h:346
void bn_mat_fromto(mat_t m, const fastf_t *from, const fastf_t *to, const struct bn_tol *tol)
Definition: mat.c:639
void bn_eigen2x2(fastf_t *val1, fastf_t *val2, fastf_t *vec1, fastf_t *vec2, fastf_t a, fastf_t b, fastf_t c)
Definition: mat.c:565
void persp_mat(mat_t m, fastf_t fovy, fastf_t aspect, fastf_t near1, fastf_t far1, fastf_t backoff)
Definition: mat.c:1209
Coord * point
Definition: chull3d.cpp:52
matp_t bn_mat_dup(const mat_t in)
Definition: mat.c:1046
void bn_htov_move(register vect_t v, register const vect_t h)
Definition: mat.c:315
int bn_mat_ck(const char *title, const mat_t m)
Definition: mat.c:1058
int bn_lseg3_lseg3_parallel(const point_t sg1pt1, const point_t sg1pt2, const point_t sg2pt1, const point_t sg2pt2, const struct bn_tol *tol)
Definition: plane.c:2532
goto out
Definition: nmg_mod.c:3846
Support for uniform tolerances.
Definition: tol.h:71
#define BN_CK_TOL(_p)
Definition: tol.h:82
void bn_wrt_point_direc(mat_t out, const mat_t change, const mat_t in, const point_t point, const vect_t direc, const struct bn_tol *tol)
Definition: mat.c:1169
double bn_atan2(double x, double y)
Definition: mat.c:104
#define BU_DEBUG_MATH
Definition: debug.h:63
ustring alpha
#define BU_DEBUG_COREDUMP
Definition: debug.h:53
void bn_mat_trn(mat_t om, register const mat_t im)
Definition: mat.c:333
void mike_persp_mat(fastf_t *pmat, const fastf_t *eye)
Definition: mat.c:1243
void bn_vtoh_move(register vect_t h, register const vect_t v)
Definition: mat.c:307
void bn_aet_vec(fastf_t *az, fastf_t *el, fastf_t *twist, fastf_t *vec_ae, fastf_t *vec_twist, fastf_t accuracy)
Definition: mat.c:390
void bn_mat_print_vls(const char *title, const mat_t m, struct bu_vls *vls)
Definition: mat.c:91
double perp
nearly 0
Definition: tol.h:75
#define ZERO(val)
Definition: units.c:38
int bu_debug
Definition: globals.c:87
void bn_vec_aed(vect_t vec, fastf_t az, fastf_t el, fastf_t dist)
Definition: mat.c:446
void bn_mat_mul3(mat_t o, const mat_t a, const mat_t b, const mat_t c)
Definition: mat.c:156
void bn_mat_print_guts(const char *title, const mat_t m, char *buf, int buflen)
Definition: mat.c:50
#define R
Definition: msr.c:55
void bn_matXvec(register vect_t ov, register const mat_t im, register const vect_t iv)
Definition: mat.c:182
void bn_mat_mul2(register const mat_t i, register mat_t o)
Definition: mat.c:146
void bu_vls_printf(struct bu_vls *vls, const char *fmt,...) _BU_ATTR_PRINTF23
Definition: vls.c:694
#define A
Definition: msr.c:51
Definition: color.c:51
void bn_mat_yrot(fastf_t *m, double siny, double cosy)
Definition: mat.c:739
#define Q
Definition: msr.c:54
void bn_vec_perp(vect_t new_vec, const vect_t old_vec)
Definition: mat.c:616
void bn_mat_mul(register mat_t o, register const mat_t a, register const mat_t b)
Definition: mat.c:121
basis_s * b
Definition: chull3d.cpp:151
void bn_mat_xform_about_pt(mat_t mat, const mat_t xform, const point_t pt)
Definition: mat.c:909
Definition: vls.h:56
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
double fastf_t
Definition: defines.h:300
#define VPRINT(a, b)
Definition: raytrace.h:1881
void bn_mat_zrot(fastf_t *m, double sinz, double cosz)
Definition: mat.c:762
void bn_mat_arb_rot(mat_t m, const point_t pt, const vect_t dir, const fastf_t ang)
Definition: mat.c:987
Definition: color.c:50
#define bu_strlcat(dst, src, size)
Definition: str.h:50
int bn_mat_inverse(register mat_t output, const mat_t input)
Definition: mat.c:219
#define M
Definition: msr.c:52
#define UNLIKELY(expression)
Definition: common.h:282