BRL-CAD
mat.h
Go to the documentation of this file.
1 /* M A T . H
2  * BRL-CAD
3  *
4  * Copyright (c) 2004-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 /*----------------------------------------------------------------------*/
22 /* @file mat.h */
23 /** @addtogroup mat */
24 /** @{ */
25 
26 /**
27  * @brief Matrix and vector functionality
28  */
29 
30 #ifndef BN_MAT_H
31 #define BN_MAT_H
32 
33 #include "common.h"
34 #include <stdio.h> /* For FILE */
35 #include "vmath.h"
36 #include "bu/vls.h"
37 #include "bn/defines.h"
38 #include "bn/tol.h"
39 
41 
42 /*
43  * 4x4 Matrix math
44  */
45 
46 BN_EXPORT extern const mat_t bn_mat_identity;
47 
48 BN_EXPORT extern void bn_mat_print(const char *title,
49  const mat_t m);
50 
51 BN_EXPORT extern void bn_mat_print_guts(const char *title,
52  const mat_t m,
53  char *buf,
54  int buflen);
55 
56 BN_EXPORT extern void bn_mat_print_vls(const char *title,
57  const mat_t m,
58  struct bu_vls *vls);
59 
60 /**
61  * A wrapper for the system atan2(). On the Silicon Graphics, and
62  * perhaps on others, x==0 incorrectly returns infinity.
63  */
64 BN_EXPORT extern double bn_atan2(double x, double y);
65 
66 #define bn_mat_zero(_m) { \
67  bu_log("%s:%d bn_mat_zero() is deprecated, use MAT_ZERO()\n", \
68  __FILE__, __LINE__); \
69  (_m)[0] = (_m)[1] = (_m)[2] = (_m)[3] = \
70  (_m)[4] = (_m)[5] = (_m)[6] = (_m)[7] = \
71  (_m)[8] = (_m)[9] = (_m)[10] = (_m)[11] = \
72  (_m)[12] = (_m)[13] = (_m)[14] = (_m)[15] = 0.0; }
73 /*
74  #define bn_mat_zero(_m) (void)memset((void *)_m, 0, sizeof(mat_t))
75 */
76 #define bn_mat_idn(_m) { \
77  bu_log("%s:%d bn_mat_idn() is deprecated, use MAT_IDN()\n", \
78  __FILE__, __LINE__); \
79  (_m)[1] = (_m)[2] = (_m)[3] = (_m)[4] = \
80  (_m)[6] = (_m)[7] = (_m)[8] = (_m)[9] = \
81  (_m)[11] = (_m)[12] = (_m)[13] = (_m)[14] = 0.0; \
82  (_m)[0] = (_m)[5] = (_m)[10] = (_m)[15] = 1.0; }
83 /*
84  #define bn_mat_idn(_m) (void)memcpy((void *)_m, (const void *)bn_mat_identity, sizeof(mat_t))
85 */
86 
87 #define bn_mat_copy(_d, _s) { \
88  bu_log("%s:%d bn_mat_copy() is deprecated, use MAT_COPY()\n", \
89  __FILE__, __LINE__); \
90  (_d)[0] = (_s)[0];\
91  (_d)[1] = (_s)[1];\
92  (_d)[2] = (_s)[2];\
93  (_d)[3] = (_s)[3];\
94  (_d)[4] = (_s)[4];\
95  (_d)[5] = (_s)[5];\
96  (_d)[6] = (_s)[6];\
97  (_d)[7] = (_s)[7];\
98  (_d)[8] = (_s)[8];\
99  (_d)[9] = (_s)[9];\
100  (_d)[10] = (_s)[10];\
101  (_d)[11] = (_s)[11];\
102  (_d)[12] = (_s)[12];\
103  (_d)[13] = (_s)[13];\
104  (_d)[14] = (_s)[14];\
105  (_d)[15] = (_s)[15]; }
106 /*
107  #define bn_mat_copy(_d, _s) (void)memcpy((void *)_d, (const void *)(_s), sizeof(mat_t))
108 */
109 
110 
111 /**
112  * Multiply matrix "a" by "b" and store the result in "o".
113  *
114  * This is different from multiplying "b" by "a" (most of the time!)
115  * Also, "o" must not be the same as either of the inputs.
116  */
117 BN_EXPORT extern void bn_mat_mul(mat_t o,
118  const mat_t a,
119  const mat_t b);
120 
121 /**
122  * o = i * o
123  *
124  * A convenience wrapper for bn_mat_mul() to update a matrix in place.
125  * The argument ordering is confusing either way.
126  */
127 BN_EXPORT extern void bn_mat_mul2(const mat_t i,
128  mat_t o);
129 
130 /**
131  * o = a * b * c
132  *
133  * The output matrix may be the same as 'b' or 'c', but may not be
134  * 'a'.
135  */
136 BN_EXPORT extern void bn_mat_mul3(mat_t o,
137  const mat_t a,
138  const mat_t b,
139  const mat_t c);
140 
141 /**
142  * o = a * b * c * d
143  *
144  * The output matrix may be the same as any input matrix.
145  */
146 BN_EXPORT extern void bn_mat_mul4(mat_t o,
147  const mat_t a,
148  const mat_t b,
149  const mat_t c,
150  const mat_t d);
151 
152 /**
153  * Multiply the matrix "im" by the vector "iv" and store the result in
154  * the vector "ov". Note this is post-multiply, and operates on
155  * 4-tuples. Use MAT4X3VEC() to operate on 3-tuples.
156  */
157 BN_EXPORT extern void bn_matXvec(hvect_t ov,
158  const mat_t im,
159  const hvect_t iv);
160 
161 /**
162  * The matrix pointed at by "input" is inverted and stored in the area
163  * pointed at by "output".
164  *
165  * Calls bu_bomb if matrix is singular.
166  */
167 BN_EXPORT extern void bn_mat_inv(mat_t output,
168  const mat_t input);
169 
170 /**
171  * The matrix pointed at by "input" is inverted and stored in the area
172  * pointed at by "output".
173  *
174  * Invert a 4-by-4 matrix using Algorithm 120 from ACM. This is a
175  * modified Gauss-Jordan algorithm. Note: Inversion is done in place,
176  * with 3 work vectors.
177  *
178  * @return 1 if OK.
179  * @return 0 if matrix is singular.
180  */
181 BN_EXPORT extern int bn_mat_inverse(mat_t output,
182  const mat_t input);
183 
184 /**
185  * Takes a pointer to a [x, y, z] vector, and a pointer to space for a
186  * homogeneous vector [x, y, z, w], and builds [x, y, z, 1].
187  */
188 BN_EXPORT extern void bn_vtoh_move(vect_t h,
189  const vect_t v);
190 
191 /**
192  * Takes a pointer to [x, y, z, w], and converts it to an ordinary
193  * vector [x/w, y/w, z/w]. Optimization for the case of w==1 is
194  * performed.
195  *
196  * FIXME: make tolerance configurable
197  */
198 BN_EXPORT extern void bn_htov_move(vect_t v,
199  const vect_t h);
200 
201 BN_EXPORT extern void bn_mat_trn(mat_t om,
202  const mat_t im);
203 
204 /**
205  * Compute a 4x4 rotation matrix given Azimuth and Elevation.
206  *
207  * Azimuth is +X, Elevation is +Z, both in degrees.
208  *
209  * Formula due to Doug Gwyn, BRL.
210  */
211 BN_EXPORT extern void bn_mat_ae(mat_t m,
212  double azimuth,
213  double elev);
214 
215 /**
216  * Find the azimuth and elevation angles that correspond to the
217  * direction (not including twist) given by a direction vector.
218  */
219 BN_EXPORT extern void bn_ae_vec(fastf_t *azp,
220  fastf_t *elp,
221  const vect_t v);
222 
223 /**
224  * Find the azimuth, elevation, and twist from two vectors. Vec_ae is
225  * in the direction of view (+z in mged view) and vec_twist points to
226  * the viewers right (+x in mged view). Accuracy (degrees) is used to
227  * stabilize flutter between equivalent extremes of atan2(), and to
228  * snap twist to zero when elevation is near +/- 90
229  */
230 BN_EXPORT extern void bn_aet_vec(fastf_t *az,
231  fastf_t *el,
232  fastf_t *twist,
233  vect_t vec_ae,
234  vect_t vec_twist,
235  fastf_t accuracy);
236 
237 /**
238  * Find a unit vector from the origin given azimuth and elevation.
239  */
240 BN_EXPORT extern void bn_vec_ae(vect_t vec,
241  fastf_t az,
242  fastf_t el);
243 
244 /**
245  * Find a vector from the origin given azimuth, elevation, and distance.
246  */
247 BN_EXPORT extern void bn_vec_aed(vect_t vec,
248  fastf_t az,
249  fastf_t el,
250  fastf_t dist);
251 
252 /**
253  * This routine builds a Homogeneous rotation matrix, given alpha,
254  * beta, and gamma as angles of rotation, in degrees.
255  *
256  * Alpha is angle of rotation about the X axis, and is done third.
257  * Beta is angle of rotation about the Y axis, and is done second.
258  * Gamma is angle of rotation about Z axis, and is done first.
259  *
260  * FIXME: make tolerance configurable
261  */
262 BN_EXPORT extern void bn_mat_angles(mat_t mat,
263  double alpha,
264  double beta, double ggamma);
265 
266 /**
267  * This routine builds a Homogeneous rotation matrix, given alpha,
268  * beta, and gamma as angles of rotation, in radians.
269  *
270  * Alpha is angle of rotation about the X axis, and is done third.
271  * Beta is angle of rotation about the Y axis, and is done second.
272  * Gamma is angle of rotation about Z axis, and is done first.
273  *
274  * FIXME: make tolerance configurable
275  */
276 BN_EXPORT extern void bn_mat_angles_rad(mat_t mat,
277  double alpha,
278  double beta,
279  double ggamma);
280 
281 /**
282  * Find the eigenvalues and eigenvectors of a symmetric 2x2 matrix.
283  * (a b)
284  * (b c)
285  *
286  * The eigenvalue with the smallest absolute value is returned in
287  * val1, with its eigenvector in vec1.
288  */
289 BN_EXPORT extern void bn_eigen2x2(fastf_t *val1,
290  fastf_t *val2,
291  vect_t vec1,
292  vect_t vec2,
293  fastf_t a,
294  fastf_t b,
295  fastf_t c);
296 
297 /**
298  * Given a vector, create another vector which is perpendicular to it.
299  * The output vector will have unit length only if the input vector
300  * did.
301  *
302  * FIXME: make tolerance configurable
303  */
304 BN_EXPORT extern void bn_vec_perp(vect_t new_vec,
305  const vect_t old_vec);
306 
307 /**
308  * Given two vectors, compute a rotation matrix that will transform
309  * space by the angle between the two. There are many candidate
310  * matrices.
311  *
312  * The input 'from' and 'to' vectors need not be unit length.
313  * MAT4X3VEC(to, m, from) is the identity that is created.
314  *
315  */
316 BN_EXPORT extern void bn_mat_fromto(mat_t m,
317  const fastf_t *from,
318  const fastf_t *to,
319  const struct bn_tol *tol);
320 
321 /**
322  * Given the sin and cos of an X rotation angle, produce the rotation
323  * matrix.
324  */
325 BN_EXPORT extern void bn_mat_xrot(mat_t m,
326  double sinx,
327  double cosx);
328 
329 /**
330  * Given the sin and cos of a Y rotation angle, produce the rotation
331  * matrix.
332  */
333 BN_EXPORT extern void bn_mat_yrot(mat_t m,
334  double siny,
335  double cosy);
336 
337 /**
338  * Given the sin and cos of a Z rotation angle, produce the rotation
339  * matrix.
340  */
341 BN_EXPORT extern void bn_mat_zrot(mat_t m,
342  double sinz,
343  double cosz);
344 
345 /**
346  * Given a direction vector D of unit length, product a matrix which
347  * rotates that vector D onto the -Z axis. This matrix will be
348  * suitable for use as a "model2view" matrix.
349  *
350  * XXX This routine will fail if the vector is already more or less
351  * aligned with the Z axis.
352  *
353  * This is done in several steps.
354  *
355  @code
356  1) Rotate D about Z to match +X axis. Azimuth adjustment.
357  2) Rotate D about Y to match -Y axis. Elevation adjustment.
358  3) Rotate D about Z to make projection of X axis again point
359  in the +X direction. Twist adjustment.
360  4) Optionally, flip sign on Y axis if original Z becomes inverted.
361  This can be nice for static frames, but is astonishing when
362  used in animation.
363  @endcode
364 */
365 BN_EXPORT extern void bn_mat_lookat(mat_t rot,
366  const vect_t dir,
367  int yflip);
368 
369 /**
370  * Given a vector, create another vector which is perpendicular to it,
371  * and with unit length. This algorithm taken from Gift's arvec.f; a
372  * faster algorithm may be possible.
373  *
374  * FIXME: make tolerance configurable
375  */
376 BN_EXPORT extern void bn_vec_ortho(vect_t out,
377  const vect_t in);
378 
379 /**
380  * Build a matrix to scale uniformly around a given point.
381  *
382  * @return -1 if scale is too small.
383  * @return 0 if OK.
384  *
385  * FIXME: make tolerance configurable
386  */
387 BN_EXPORT extern int bn_mat_scale_about_pt(mat_t mat,
388  const point_t pt,
389  const double scale);
390 
391 /**
392  * Build a matrix to apply arbitrary 4x4 transformation around a given
393  * point.
394  */
395 BN_EXPORT extern void bn_mat_xform_about_pt(mat_t mat,
396  const mat_t xform,
397  const point_t pt);
398 
399 /**
400  * @return 0 When matrices are not equal
401  * @return 1 When matrices are equal
402  */
403 BN_EXPORT extern int bn_mat_is_equal(const mat_t a,
404  const mat_t b,
405  const struct bn_tol *tol);
406 
407 /**
408  * This routine is intended for detecting identity matrices read in
409  * from ascii or binary files, where the numbers are pure ones or
410  * zeros. This routine is *not* intended for tolerance-based
411  * "near-zero" comparisons; as such, it shouldn't be used on matrices
412  * which are the result of calculation.
413  *
414  *
415  * @return 0 non-identity
416  * @return 1 a perfect identity matrix
417  */
418 BN_EXPORT extern int bn_mat_is_identity(const mat_t m);
419 
420 /**
421  * Construct a transformation matrix for rotation about an arbitrary
422  * axis. The axis is defined by a point (pt) and a unit direction
423  * vector (dir). The angle of rotation is "ang"
424  *
425  * FIXME: make tolerance configurable
426  */
427 BN_EXPORT extern void bn_mat_arb_rot(mat_t m,
428  const point_t pt,
429  const vect_t dir,
430  const fastf_t ang);
431 
432 /**
433  * Return a pointer to a copy of the matrix in dynamically allocated
434  * memory.
435  */
436 BN_EXPORT extern matp_t bn_mat_dup(const mat_t in);
437 
438 /**
439  * Check to ensure that a rotation matrix preserves axis
440  * perpendicularity. Note that not all matrices are rotation
441  * matrices.
442  *
443  *
444  * @return -1 FAIL
445  * @return 0 OK
446  */
447 BN_EXPORT extern int bn_mat_ck(const char *title,
448  const mat_t m);
449 
450 /**
451  * Calculates the determinant of the 3X3 "rotation" part of the passed
452  * matrix.
453  */
454 BN_EXPORT extern fastf_t bn_mat_det3(const mat_t m);
455 
456 /**
457  * Calculates the determinant of the 4X4 matrix
458  */
459 BN_EXPORT extern fastf_t bn_mat_determinant(const mat_t m);
460 
461 /**
462  * FIXME: make tolerance configurable
463  */
464 BN_EXPORT extern int bn_mat_is_non_unif(const mat_t m);
465 
466 /**
467  * Given a model-space transformation matrix "change",
468  * return a matrix which applies the change with-respect-to
469  * given "point" and "direc".
470  */
471 BN_EXPORT extern void bn_wrt_point_direc(mat_t out,
472  const mat_t change,
473  const mat_t in,
474  const point_t point,
475  const vect_t direc,
476  const struct bn_tol *tol);
477 
478 
479 /* Matrix stuff from MGED's dozoom.c - somewhat unfinished,
480  * judging by the comments. Can probably be simplified/combined */
481 BN_EXPORT extern void persp_mat(mat_t m, fastf_t fovy, fastf_t aspect, fastf_t near1, fastf_t far1, fastf_t backoff);
482 BN_EXPORT extern void mike_persp_mat(fastf_t *pmat, const fastf_t *eye);
483 BN_EXPORT extern void deering_persp_mat(fastf_t *m, const fastf_t *l, const fastf_t *h, const fastf_t *eye);
484 
486 
487 #endif /* BN_MAT_H */
488 /** @} */
489 /*
490  * Local Variables:
491  * mode: C
492  * tab-width: 8
493  * indent-tabs-mode: t
494  * c-file-style: "stroustrup"
495  * End:
496  * ex: shiftwidth=4 tabstop=8
497  */
void bn_ae_vec(fastf_t *azp, fastf_t *elp, const vect_t v)
Definition: mat.c:374
void bn_mat_ae(mat_t m, double azimuth, double elev)
int bn_mat_is_identity(const mat_t m)
Definition: mat.c:980
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
const mat_t bn_mat_identity
Matrix and vector functionality.
Definition: mat.c:46
void bn_htov_move(vect_t v, const vect_t h)
void bn_mat_xrot(mat_t m, double sinx, double cosx)
Header file for the BRL-CAD common definitions.
int bn_mat_inverse(mat_t output, const mat_t input)
void bn_eigen2x2(fastf_t *val1, fastf_t *val2, vect_t vec1, vect_t vec2, fastf_t a, fastf_t b, fastf_t c)
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_mat_mul2(const mat_t i, mat_t o)
int bn_mat_is_equal(const mat_t a, const mat_t b, const struct bn_tol *tol)
Definition: mat.c:925
void deering_persp_mat(fastf_t *m, const fastf_t *l, const fastf_t *h, const fastf_t *eye)
Definition: mat.c:1312
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
void bn_matXvec(hvect_t ov, const mat_t im, const hvect_t iv)
#define __BEGIN_DECLS
Definition: common.h:73
void bn_vec_ae(vect_t vec, fastf_t az, fastf_t el)
Definition: mat.c:433
void bn_mat_inv(mat_t output, const mat_t input)
void bn_vtoh_move(vect_t h, const vect_t v)
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 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_mat_trn(mat_t om, const mat_t im)
int bn_mat_ck(const char *title, const mat_t m)
Definition: mat.c:1058
goto out
Definition: nmg_mod.c:3846
void bn_mat_mul(mat_t o, const mat_t a, const mat_t b)
Support for uniform tolerances.
Definition: tol.h:71
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
void bn_mat_angles(mat_t mat, double alpha, double beta, double ggamma)
void bn_mat_angles_rad(mat_t mat, double alpha, double beta, double ggamma)
ustring alpha
void mike_persp_mat(fastf_t *pmat, const fastf_t *eye)
Definition: mat.c:1243
void bn_mat_print_vls(const char *title, const mat_t m, struct bu_vls *vls)
Definition: mat.c:91
void bn_vec_aed(vect_t vec, fastf_t az, fastf_t el, fastf_t dist)
Definition: mat.c:446
void bn_mat_yrot(mat_t m, double siny, double cosy)
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
void bn_vec_ortho(vect_t out, const vect_t in)
void bn_aet_vec(fastf_t *az, fastf_t *el, fastf_t *twist, vect_t vec_ae, vect_t vec_twist, fastf_t accuracy)
#define __END_DECLS
Definition: common.h:74
void bn_vec_perp(vect_t new_vec, const vect_t old_vec)
Definition: mat.c:616
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
double fastf_t
Definition: defines.h:300
void bn_mat_zrot(mat_t m, double sinz, double cosz)
void bn_mat_arb_rot(mat_t m, const point_t pt, const vect_t dir, const fastf_t ang)
Definition: mat.c:987