1 /**************************************************************************\
2 * Copyright (c) Kongsberg Oil & Gas Technologies AS
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are
7 * met:
8 *
9 * Redistributions of source code must retain the above copyright notice,
10 * this list of conditions and the following disclaimer.
11 *
12 * Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 *
16 * Neither the name of the copyright holder nor the names of its
17 * contributors may be used to endorse or promote products derived from
18 * this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 \**************************************************************************/
32
33 /*!
34 \class SbDPMatrix SbDPMatrix.h Inventor/SbDPMatrix.h
35 \brief The SbDPMatrix class is a 4x4 dimensional representation of a double-precision matrix.
36
37 \ingroup base
38
39 This class is like the SbMatrix class, but uses double-precision
40 floating point values for its elements. For more class
41 documentation, see SbMatrix.
42
43 \since Coin 2.0.
44 */
45
46 // FIXME: we should _really_ have double-precision classes compatible
47 // with those in TGS' API, for several good reasons. So either rename
48 // this, make a typedef (if that is sufficient), or write a "wrapper
49 // class" around this with inline functions, using with TGS' name for
50 // it. 20020225 mortene.
51
52 // FIXME:
53 //
54 // * The SbDPMatrix::factor() function has not been implemented yet.
55 //
56 // * The element access methods should be inlined.
57 //
58 // * Optimizations are not done yet, so there's a lot of room for
59 // improvements.
60
61
62 #include <Inventor/SbDPMatrix.h>
63
64 #include <cassert>
65 #include <cstring>
66 #include <cfloat>
67
68 #include <Inventor/SbDPRotation.h>
69 #include <Inventor/SbDPLine.h>
70 #include <Inventor/SbMatrix.h>
71
72 #include "coindefs.h" // COIN_STUB()
73
74 #if COIN_DEBUG
75 #include <Inventor/errors/SoDebugError.h>
76 #endif // COIN_DEBUG
77
78 #ifndef COIN_WORKAROUND_NO_USING_STD_FUNCS
79 using std::memmove;
80 using std::memcmp;
81 using std::memcpy;
82 #endif // !COIN_WORKAROUND_NO_USING_STD_FUNCS
83
84 // FIXME: should merge all the PD code we're using from GGIV into
85 // SbDPMatrix, SbDPRotation and SbVec3d proper (for two reasons: 1)
86 // there's a lot of duplicated code here (like for instance the
87 // matrix->quaternion decomposition, which also exists in
88 // SbDPRotation::setValue(SbDPMatrix&)), and 2) the remaining code
89 // snippets look generally useful outside the purpose of breaking down
90 // a matrix into it's transformation components). 20010114 mortene.
91
92 /*
93 * declarations for polar_decomp algorithm from Graphics Gems IV,
94 * by Ken Shoemake <shoemake@graphics.cis.upenn.edu>
95 */
96 enum QuatPart {X, Y, Z, W};
97 typedef double HMatrix[4][4]; /* Right-handed, for column vectors */
98 typedef struct {
99 SbVec4d t; /* Translation components */
100 SbDPRotation q; /* Essential rotation */
101 SbDPRotation u; /* Stretch rotation */
102 SbVec4d k; /* Stretch factors */
103 double f; /* Sign of determinant */
104 } AffineParts;
105 static double polar_decomp(HMatrix M, HMatrix Q, HMatrix S);
106 static SbVec4d spect_decomp(HMatrix S, HMatrix U);
107 static SbDPRotation snuggle(SbDPRotation q, SbVec4d & k);
108 static void decomp_affine(HMatrix A, AffineParts * parts);
109
110 static const SbDPMat IDENTITYMATRIX = {
111 { 1.0, 0.0, 0.0, 0.0 },
112 { 0.0, 1.0, 0.0, 0.0 },
113 { 0.0, 0.0, 1.0, 0.0 },
114 { 0.0, 0.0, 0.0, 1.0 }
115 };
116
117 static inline
SbDPMatrix_isIdentity(const double fm[][4])118 SbBool SbDPMatrix_isIdentity(const double fm[][4])
119 {
120 #if 0 // I would assume that the memcmp() version is faster..? Should run some profile checks.
121 return ((fm[0][0] == 1.0) && (fm[0][1] == 0.0) && (fm[0][2] == 0.0) && (fm[0][3] == 0.0) &&
122 (fm[1][0] == 0.0) && (fm[1][1] == 1.0) && (fm[1][2] == 0.0) && (fm[1][3] == 0.0) &&
123 (fm[2][0] == 0.0) && (fm[2][1] == 0.0) && (fm[2][2] == 1.0) && (fm[2][3] == 0.0) &&
124 (fm[3][0] == 0.0) && (fm[3][1] == 0.0) && (fm[3][2] == 0.0) && (fm[3][3] == 1.0));
125 #else
126 // Note: as far as I know, memcmp() only compares bytes until
127 // there's a mismatch (and does *not* run over the full array and
128 // adds up a total, as it sometimes seems from documentation). So
129 // this should be very quick for non-identity matrices.
130 //
131 // Also, we check the first value on it's own, to avoid the function
132 // call for the most common case.
133 return (fm[0][0]==1.0) && memcmp(&fm[0][1], &IDENTITYMATRIX[0][1], (4 * 3 + 3) * sizeof(double)) == 0;
134 #endif
135 }
136
137 /*!
138 The default constructor does nothing. The matrix will be uninitialized.
139 */
SbDPMatrix(void)140 SbDPMatrix::SbDPMatrix(void)
141 {
142 }
143
144
145 /*!
146 Constructs a matrix instance with the given initial elements.
147 */
SbDPMatrix(const double a11,const double a12,const double a13,const double a14,const double a21,const double a22,const double a23,const double a24,const double a31,const double a32,const double a33,const double a34,const double a41,const double a42,const double a43,const double a44)148 SbDPMatrix::SbDPMatrix(const double a11, const double a12,
149 const double a13, const double a14,
150 const double a21, const double a22,
151 const double a23, const double a24,
152 const double a31, const double a32,
153 const double a33, const double a34,
154 const double a41, const double a42,
155 const double a43, const double a44)
156 {
157 const SbDPMat m = { { a11, a12, a13, a14 },
158 { a21, a22, a23, a24 },
159 { a31, a32, a33, a34 },
160 { a41, a42, a43, a44 } };
161 this->setValue(m);
162 }
163
164 /*!
165 Constructs a matrix instance with the initial elements from the
166 \a matrix argument.
167 */
SbDPMatrix(const SbDPMat & matrixref)168 SbDPMatrix::SbDPMatrix(const SbDPMat & matrixref)
169 {
170 this->setValue(matrixref);
171 }
172
173 /*!
174 This constructor is courtesy of the Microsoft Visual C++ compiler.
175 */
SbDPMatrix(const SbDPMat * matrixptr)176 SbDPMatrix::SbDPMatrix(const SbDPMat * matrixptr)
177 {
178 this->setValue(*matrixptr);
179 }
180
181 /*!
182 This constructor converts a single-precision matrix to a double-precision matrix.
183 */
SbDPMatrix(const SbMatrix & matrixref)184 SbDPMatrix::SbDPMatrix(const SbMatrix & matrixref)
185 {
186
187 const SbDPMat m = { { matrixref[0][0], matrixref[0][1], matrixref[0][2], matrixref[0][3] },
188 { matrixref[1][0], matrixref[1][1], matrixref[1][2], matrixref[1][3] },
189 { matrixref[2][0], matrixref[2][1], matrixref[2][2], matrixref[2][3] },
190 { matrixref[3][0], matrixref[3][1], matrixref[3][2], matrixref[3][3] } };
191 this->setValue(m);
192 }
193
194 /*!
195 Default destructor does nothing.
196 */
~SbDPMatrix(void)197 SbDPMatrix::~SbDPMatrix(void)
198 {
199 }
200
201 /*!
202 Returns a pointer to the 2 dimensional double array with the matrix
203 elements.
204
205 \sa setValue().
206 */
207 const SbDPMat &
getValue(void) const208 SbDPMatrix::getValue(void) const
209 {
210 return this->matrix;
211 }
212
213 /*!
214 Copies the elements from \a m into the matrix.
215
216 \sa getValue().
217 */
218 void
setValue(const SbDPMat & m)219 SbDPMatrix::setValue(const SbDPMat & m)
220 {
221 (void)memmove(this->matrix, m, sizeof(double)*4*4);
222 }
223
224 /*!
225 Copies the elements from \a m into the matrix.
226
227 \sa getValue().
228 */
229 void
setValue(const double * m)230 SbDPMatrix::setValue(const double * m)
231 {
232 (void)memmove(this->matrix, m, sizeof(double)*4*4);
233 }
234
235 /*!
236 Copies the elements from \a m into the matrix.
237
238 \sa getValue().
239 */
240 void
setValue(const SbMatrix & m)241 SbDPMatrix::setValue(const SbMatrix & m)
242 {
243 const SbDPMat dmat = { { m[0][0], m[0][1], m[0][2], m[0][3] },
244 { m[1][0], m[1][1], m[1][2], m[1][3] },
245 { m[2][0], m[2][1], m[2][2], m[2][3] },
246 { m[3][0], m[3][1], m[3][2], m[3][3] } };
247 this->setValue(dmat);
248 }
249
250
251 /*!
252 Assignment operator. Copies the elements from \a m to the matrix.
253 */
254 SbDPMatrix &
operator =(const SbDPMat & m)255 SbDPMatrix::operator=(const SbDPMat & m)
256 {
257 this->setValue(m);
258 return *this;
259 }
260
261 /*!
262 Assignment operator. Copies the elements from \a m to the matrix.
263 */
264 SbDPMatrix &
operator =(const SbDPMatrix & m)265 SbDPMatrix::operator=(const SbDPMatrix & m)
266 {
267 this->setValue(m.matrix);
268 return *this;
269 }
270
271 /*!
272 Set the matrix to be the identity matrix.
273
274 \sa identity().
275 */
276 void
makeIdentity(void)277 SbDPMatrix::makeIdentity(void)
278 {
279 this->matrix[0][0]=this->matrix[1][1]=
280 this->matrix[2][2]=this->matrix[3][3] = 1.0;
281
282 this->matrix[0][1]=this->matrix[0][2]=this->matrix[0][3]=
283 this->matrix[1][0]=this->matrix[1][2]=this->matrix[1][3]=
284 this->matrix[2][0]=this->matrix[2][1]=this->matrix[2][3]=
285 this->matrix[3][0]=this->matrix[3][1]=this->matrix[3][2] = 0.0;
286 }
287
288 /*!
289 Set matrix to be a rotation matrix with the given rotation.
290
291 \sa setTranslate(), setScale().
292 */
293 void
setRotate(const SbDPRotation & q)294 SbDPMatrix::setRotate(const SbDPRotation & q)
295 {
296 q.getValue(*this);
297 }
298
299 /*!
300 Multiply all element values in the matrix with \a v.
301 */
302 void
operator *=(const double v)303 SbDPMatrix::operator*=(const double v)
304 {
305 for (int i=0; i < 4; i++) {
306 for (int j=0; j < 4; j++) {
307 this->matrix[i][j] *= v;
308 }
309 }
310 }
311
312 /*!
313 Divide all element values in the matrix on \a v.
314 */
315 void
operator /=(const double v)316 SbDPMatrix::operator/=(const double v)
317 {
318 #if COIN_DEBUG
319 if (v == 0.0)
320 SoDebugError::postWarning("SbDPMatrix::operator/=",
321 "Division by zero.");
322 #endif // COIN_DEBUG
323
324 this->operator*=(1.0/v);
325 }
326
327 /*!
328 Returns the determinant of the 3x3 submatrix specified by the row and
329 column indices.
330 */
331 double
det3(int r1,int r2,int r3,int c1,int c2,int c3) const332 SbDPMatrix::det3(int r1, int r2, int r3,
333 int c1, int c2, int c3) const
334 {
335 #if COIN_EXTRA_DEBUG
336 // Check indices.
337 if (r1<0 || r1>3 || r2<0 || r2>3 || r3<0 || r3>3 ||
338 c1<0 || c1>3 || c2<0 || c2>3 || c3<0 || c3>3) {
339 SoDebugError::post("SbDPMatrix::det3",
340 "At least one idx out of bounds [0, 3]. ");
341 }
342 if (r1==r2 || r1==r3 || r2==r3 ||
343 c1==c2 || c1==c3 || c2==c3)
344 SoDebugError::post("SbDPMatrix::det3", "Indices should be distinct.");
345 #endif // COIN_EXTRA_DEBUG
346
347 // More or less directly from "Advanced Engineering Mathematics"
348 // (E. Kreyszig), 6th edition.
349
350 double a11 = this->matrix[r1][c1];
351 double a12 = this->matrix[r1][c2];
352 double a13 = this->matrix[r1][c3];
353 double a21 = this->matrix[r2][c1];
354 double a22 = this->matrix[r2][c2];
355 double a23 = this->matrix[r2][c3];
356 double a31 = this->matrix[r3][c1];
357 double a32 = this->matrix[r3][c2];
358 double a33 = this->matrix[r3][c3];
359
360 double M11 = a22 * a33 - a32 * a23;
361 double M21 = -(a12 * a33 - a32 * a13);
362 double M31 = a12 * a23 - a22 * a13;
363
364 return (a11 * M11 + a21 * M21 + a31 * M31);
365 }
366
367 /*!
368 Returns the determinant of the upper left 3x3 submatrix.
369 */
370 double
det3(void) const371 SbDPMatrix::det3(void) const
372 {
373 return this->det3(0, 1, 2, 0, 1, 2);
374 }
375
376 /*!
377 Returns the determinant of the matrix.
378 */
379 double
det4(void) const380 SbDPMatrix::det4(void) const
381 {
382 double det = 0.0;
383 det += this->matrix[0][0] * det3(1, 2, 3, 1, 2, 3);
384 det -= this->matrix[1][0] * det3(0, 2, 3, 1, 2, 3);
385 det += this->matrix[2][0] * det3(0, 1, 3, 1, 2, 3);
386 det -= this->matrix[3][0] * det3(0, 1, 2, 1, 2, 3);
387 return det;
388 }
389
390 /*!
391 Return a new matrix which is the inverse matrix of this.
392
393 The user is responsible for checking that this is a valid operation
394 to execute, by first making sure that the result of SbDPMatrix::det4()
395 is not equal to zero.
396 */
397 SbDPMatrix
inverse(void) const398 SbDPMatrix::inverse(void) const
399 {
400 // check for identity matrix
401 if (SbDPMatrix_isIdentity(this->matrix)) { return SbMatrix::identity(); }
402
403 SbDPMatrix result;
404
405 // use local pointers for speed
406 double (*dst)[4];
407 dst = reinterpret_cast<double (*)[4]>(result.matrix[0]);
408 const double (*src)[4];
409 src = reinterpret_cast<const double (*)[4]>(this->matrix[0]);
410
411 // check for affine matrix (common case)
412 if (src[0][3] == 0.0 && src[1][3] == 0.0 &&
413 src[2][3] == 0.0 && src[3][3] == 1.0) {
414
415 // More or less directly from:
416 // Kevin Wu, "Fast Matrix Inversion", Graphics Gems II
417 double det_1;
418 double pos, neg, temp;
419
420 #define ACCUMULATE \
421 if (temp >= 0.0) \
422 pos += temp; \
423 else \
424 neg += temp
425
426 /*
427 * Calculate the determinant of submatrix A and determine if the
428 * the matrix is singular as limited by floating-point data
429 * representation.
430 */
431 pos = neg = 0.0;
432 temp = src[0][0] * src[1][1] * src[2][2];
433 ACCUMULATE;
434 temp = src[0][1] * src[1][2] * src[2][0];
435 ACCUMULATE;
436 temp = src[0][2] * src[1][0] * src[2][1];
437 ACCUMULATE;
438 temp = -src[0][2] * src[1][1] * src[2][0];
439 ACCUMULATE;
440 temp = -src[0][1] * src[1][0] * src[2][2];
441 ACCUMULATE;
442 temp = -src[0][0] * src[1][2] * src[2][1];
443 ACCUMULATE;
444 det_1 = pos + neg;
445
446 #undef ACCUMULATE
447
448 /* Is the submatrix A singular? */
449 if ((det_1 == 0.0) || (SbAbs(det_1 / (pos - neg)) < DBL_EPSILON)) {
450 /* Matrix M has no inverse */
451 #if COIN_DEBUG
452 SoDebugError::postWarning("SbMatrix::inverse",
453 "Matrix is singular.");
454 #endif // COIN_DEBUG
455 return *this;
456 }
457 else {
458 /* Calculate inverse(A) = adj(A) / det(A) */
459 det_1 = 1.0 / det_1;
460 dst[0][0] = (src[1][1] * src[2][2] -
461 src[1][2] * src[2][1]) * det_1;
462 dst[1][0] = - (src[1][0] * src[2][2] -
463 src[1][2] * src[2][0]) * det_1;
464 dst[2][0] = (src[1][0] * src[2][1] -
465 src[1][1] * src[2][0]) * det_1;
466 dst[0][1] = - (src[0][1] * src[2][2] -
467 src[0][2] * src[2][1]) * det_1;
468 dst[1][1] = (src[0][0] * src[2][2] -
469 src[0][2] * src[2][0]) * det_1;
470 dst[2][1] = - (src[0][0] * src[2][1] -
471 src[0][1] * src[2][0]) * det_1;
472 dst[0][2] = (src[0][1] * src[1][2] -
473 src[0][2] * src[1][1]) * det_1;
474 dst[1][2] = - (src[0][0] * src[1][2] -
475 src[0][2] * src[1][0]) * det_1;
476 dst[2][2] = (src[0][0] * src[1][1] -
477 src[0][1] * src[1][0]) * det_1;
478
479 /* Calculate -C * inverse(A) */
480 dst[3][0] = - (src[3][0] * dst[0][0] +
481 src[3][1] * dst[1][0] +
482 src[3][2] * dst[2][0]);
483 dst[3][1] = - (src[3][0] * dst[0][1] +
484 src[3][1] * dst[1][1] +
485 src[3][2] * dst[2][1]);
486 dst[3][2] = - (src[3][0] * dst[0][2] +
487 src[3][1] * dst[1][2] +
488 src[3][2] * dst[2][2]);
489
490 /* Fill in last column */
491 dst[0][3] = dst[1][3] = dst[2][3] = 0.0;
492 dst[3][3] = 1.0;
493 }
494 }
495 else { // non-affine matrix
496 double max, sum, tmp, inv_pivot;
497 int p[4];
498 int i, j, k;
499
500 // algorithm from: Schwarz, "Numerische Mathematik"
501 result = *this;
502
503 for (k = 0; k < 4; k++) {
504 max = 0.0;
505 p[k] = 0;
506
507 for (i = k; i < 4; i++) {
508 sum = 0.0;
509 for (j = k; j < 4; j++)
510 sum += SbAbs(dst[i][j]);
511 if (sum > 0.0) {
512 tmp = SbAbs(dst[i][k]) / sum;
513 if (tmp > max) {
514 max = tmp;
515 p[k] = i;
516 }
517 }
518 }
519
520 if (max == 0.0) {
521 #if COIN_DEBUG
522 SoDebugError::postWarning("SbMatrix::inverse",
523 "Matrix is singular.");
524 #endif // COIN_DEBUG
525 return *this;
526 }
527
528 if (p[k] != k) {
529 for (j = 0; j < 4; j++) {
530 tmp = dst[k][j];
531 dst[k][j] = dst[p[k]][j];
532 dst[p[k]][j] = tmp;
533 }
534 }
535
536 inv_pivot = 1.0 / dst[k][k];
537 for (j = 0; j < 4; j++) {
538 if (j != k) {
539 dst[k][j] = - dst[k][j] * inv_pivot;
540 for (i = 0; i < 4; i++) {
541 if (i != k) dst[i][j] += dst[i][k] * dst[k][j];
542 }
543 }
544 }
545
546 for (i = 0; i < 4; i++) dst[i][k] *= inv_pivot;
547 dst[k][k] = inv_pivot;
548 }
549
550 for (k = 2; k >= 0; k--) {
551 if (p[k] != k) {
552 for (i = 0; i < 4; i++) {
553 tmp = dst[i][k];
554 dst[i][k] = dst[i][p[k]];
555 dst[i][p[k]] = tmp;
556 }
557 }
558 }
559 }
560 return result;
561 }
562
563 /*!
564 Check if the \a m matrix is equal to this one, within the given tolerance
565 value. The tolerance value is applied in the comparison on a component by
566 component basis.
567 */
568 SbBool
equals(const SbDPMatrix & m,double tolerance) const569 SbDPMatrix::equals(const SbDPMatrix & m, double tolerance) const
570 {
571 #if COIN_DEBUG
572 if (tolerance < 0.0)
573 SoDebugError::postWarning("SbDPMatrix::equals",
574 "tolerance should be >=0.0.");
575 #endif // COIN_DEBUG
576
577 for (int i=0; i < 4; i++) {
578 for (int j=0; j< 4; j++) {
579 if (fabs(this->matrix[i][j] - m.matrix[i][j]) > tolerance) return FALSE;
580 }
581 }
582
583 return TRUE;
584 }
585
586
587 /*!
588 Return pointer to the matrix' 4x4 double array.
589 */
operator double*(void)590 SbDPMatrix::operator double*(void)
591 {
592 return &(this->matrix[0][0]);
593 }
594
595
596 /*!
597 Return pointer to the matrix' 4x4 double array.
598 */
operator SbDPMat&(void)599 SbDPMatrix::operator SbDPMat&(void)
600 {
601 return this->matrix;
602 }
603
604 /*!
605 Returns pointer to the 4 element array representing a matrix row.
606 \a i should be within [0, 3].
607
608 \sa getValue(), setValue().
609 */
610 double *
operator [](int i)611 SbDPMatrix::operator [](int i)
612 {
613 #if COIN_EXTRA_DEBUG
614 if (i<0 || i>3) {
615 SoDebugError::post("SbDPMatrix::operator[]", "Index out of bounds. ");
616 }
617 #endif // COIN_EXTRA_DEBUG
618
619 return this->matrix[i];
620 }
621
622 /*!
623 Returns pointer to the 4 element array representing a matrix row.
624 \a i should be within [0, 3].
625
626 \sa getValue(), setValue().
627 */
628 const double *
operator [](int i) const629 SbDPMatrix::operator [](int i) const
630 {
631 #if COIN_EXTRA_DEBUG
632 if (i<0 || i>3) {
633 SoDebugError::postWarning("SbDPMatrix::operator[]", "Index out of bounds. ");
634 }
635 #endif // COIN_EXTRA_DEBUG
636
637 return this->matrix[i];
638 }
639
640 /*!
641 Set matrix to be a rotation matrix with the given rotation.
642
643 \sa setRotate().
644 */
645 SbDPMatrix&
operator =(const SbDPRotation & q)646 SbDPMatrix::operator =(const SbDPRotation & q)
647 {
648 this->setRotate(q);
649 return *this;
650 }
651
652 /*!
653 Right-multiply with the \a m matrix.
654
655 \sa multRight().
656 */
657 SbDPMatrix&
operator *=(const SbDPMatrix & m)658 SbDPMatrix::operator *=(const SbDPMatrix & m)
659 {
660 return this->multRight(m);
661 }
662
663 /*!
664 \relates SbDPMatrix
665
666 Multiplies matrix \a m1 with matrix \a m2 and returns the resultant
667 matrix.
668 */
669 SbDPMatrix
operator *(const SbDPMatrix & m1,const SbDPMatrix & m2)670 operator *(const SbDPMatrix & m1, const SbDPMatrix & m2)
671 {
672 SbDPMatrix result = m1;
673 result *= m2;
674 return result;
675 }
676
677 /*!
678 \relates SbDPMatrix
679
680 Compare matrices to see if they are equal. For two matrices to be equal,
681 all their individual elements must be equal.
682
683 \sa equals().
684 */
685 int
operator ==(const SbDPMatrix & m1,const SbDPMatrix & m2)686 operator ==(const SbDPMatrix & m1, const SbDPMatrix & m2)
687 {
688 for (int i=0; i < 4; i++) {
689 for (int j=0; j < 4; j++) {
690 if (m1.matrix[i][j] != m2.matrix[i][j]) return FALSE;
691 }
692 }
693
694 return TRUE;
695 }
696
697 /*!
698 \relates SbDPMatrix
699
700 Compare matrices to see if they are not equal. For two matrices to not be
701 equal, it is enough that at least one of their elements are unequal.
702
703 \sa equals().
704 */
705 int
operator !=(const SbDPMatrix & m1,const SbDPMatrix & m2)706 operator !=(const SbDPMatrix & m1, const SbDPMatrix & m2)
707 {
708 return !(m1 == m2);
709 }
710
711 /*!
712 Return matrix components in the SbDPMat structure.
713
714 \sa setValue().
715 */
716 void
getValue(SbDPMat & m) const717 SbDPMatrix::getValue(SbDPMat & m) const
718 {
719 (void)memmove(&m[0][0], &(this->matrix[0][0]), sizeof(double)*4*4);
720 }
721
722 /*!
723 Return the identity matrix.
724
725 \sa makeIdentity().
726 */
727 SbDPMatrix
identity(void)728 SbDPMatrix::identity(void)
729 {
730 return SbDPMatrix(&IDENTITYMATRIX);
731 }
732
733 /*!
734 Set matrix to be a pure scaling matrix. Scale factors are specified
735 by \a s.
736
737 \sa setRotate(), setTranslate().
738 */
739 void
setScale(const double s)740 SbDPMatrix::setScale(const double s)
741 {
742 this->makeIdentity();
743 this->matrix[0][0] = s;
744 this->matrix[1][1] = s;
745 this->matrix[2][2] = s;
746 }
747
748 /*!
749 Set matrix to be a pure scaling matrix. Scale factors in x, y and z
750 is specified by the \a s vector.
751
752 \sa setRotate(), setTranslate().
753 */
754 void
setScale(const SbVec3d & s)755 SbDPMatrix::setScale(const SbVec3d & s)
756 {
757 this->makeIdentity();
758 this->matrix[0][0] = s[0];
759 this->matrix[1][1] = s[1];
760 this->matrix[2][2] = s[2];
761 }
762
763 /*!
764 Make this matrix into a pure translation matrix (no scale or rotation
765 components) with the given vector \a t as the translation.
766
767 \sa setRotate(), setScale().
768 */
769 void
setTranslate(const SbVec3d & t)770 SbDPMatrix::setTranslate(const SbVec3d & t)
771 {
772 this->makeIdentity();
773 this->matrix[3][0] = t[0];
774 this->matrix[3][1] = t[1];
775 this->matrix[3][2] = t[2];
776 }
777
778 /*!
779 Set translation, rotation and scaling all at once. The resulting
780 matrix gets calculated like this:
781
782 \code
783 M = S * R * T
784 \endcode
785
786 where \a S, \a R and \a T is scaling, rotation and translation
787 matrices.
788
789 \sa setTranslate(), setRotate(), setScale() and getTransform().
790 */
791 void
setTransform(const SbVec3d & t,const SbDPRotation & r,const SbVec3d & s)792 SbDPMatrix::setTransform(const SbVec3d & t, const SbDPRotation & r, const SbVec3d & s)
793 {
794 SbDPMatrix tmp;
795
796 this->setScale(s);
797
798 tmp.setRotate(r);
799 this->multRight(tmp);
800
801 tmp.setTranslate(t);
802 this->multRight(tmp);
803 }
804
805 /*!
806 Set translation, rotation and scaling all at once with a specified
807 scale orientation. The resulting matrix gets calculated like this:
808
809 \code
810 M = Ro-� * S * Ro * R * T
811 \endcode
812
813 where \a Ro is the scale orientation, and \a S, \a R
814 and \a T is scaling, rotation and translation.
815
816 \sa setTranslate(), setRotate(), setScale() and getTransform().
817 */
818 void
setTransform(const SbVec3d & t,const SbDPRotation & r,const SbVec3d & s,const SbDPRotation & so)819 SbDPMatrix::setTransform(const SbVec3d & t, const SbDPRotation & r,
820 const SbVec3d & s, const SbDPRotation & so)
821 {
822 SbDPMatrix tmp;
823
824 this->setRotate(so.inverse());
825
826 tmp.setScale(s);
827 this->multRight(tmp);
828
829 tmp.setRotate(so);
830 this->multRight(tmp);
831
832 tmp.setRotate(r);
833 this->multRight(tmp);
834
835 tmp.setTranslate(t);
836 this->multRight(tmp);
837 }
838
839 /*!
840 Set translation, rotation and scaling all at once with a specified
841 scale orientation and center point. The resulting matrix gets
842 calculated like this:
843
844 \code
845 M = -Tc * Ro-� * S * Ro * R * T * Tc
846 \endcode
847
848 where \a Tc is the center point, \a Ro the scale orientation, \a S,
849 \a R and \a T is scaling, rotation and translation.
850
851 \sa setTranslate(), setRotate(), setScale() and getTransform().
852 */
853 void
setTransform(const SbVec3d & translation,const SbDPRotation & rotation,const SbVec3d & scaleFactor,const SbDPRotation & scaleOrientation,const SbVec3d & center)854 SbDPMatrix::setTransform(const SbVec3d & translation,
855 const SbDPRotation & rotation,
856 const SbVec3d & scaleFactor,
857 const SbDPRotation & scaleOrientation,
858 const SbVec3d & center)
859 {
860 SbDPMatrix tmp;
861
862 this->setTranslate(-center);
863
864 tmp.setRotate(scaleOrientation.inverse());
865 this->multRight(tmp);
866
867 tmp.setScale(scaleFactor);
868 this->multRight(tmp);
869
870 tmp.setRotate(scaleOrientation);
871 this->multRight(tmp);
872
873 tmp.setRotate(rotation);
874 this->multRight(tmp);
875
876 tmp.setTranslate(translation);
877 this->multRight(tmp);
878
879 tmp.setTranslate(center);
880 this->multRight(tmp);
881 }
882
883 /*!
884 Factor the matrix back into its translation, rotation, scale and
885 scaleorientation components.
886
887 \sa factor()
888 */
889 void
getTransform(SbVec3d & t,SbDPRotation & r,SbVec3d & s,SbDPRotation & so) const890 SbDPMatrix::getTransform(SbVec3d & t, SbDPRotation & r, SbVec3d & s,
891 SbDPRotation & so) const
892 {
893 // FIXME: test if this code works with non-affine matrices.
894 // pederb, 2000-01-17
895
896 AffineParts parts;
897 HMatrix hmatrix;
898
899 // transpose-copy
900 hmatrix[0][0] = this->matrix[0][0];
901 hmatrix[0][1] = this->matrix[1][0];
902 hmatrix[0][2] = this->matrix[2][0];
903 hmatrix[0][3] = this->matrix[3][0];
904
905 hmatrix[1][0] = this->matrix[0][1];
906 hmatrix[1][1] = this->matrix[1][1];
907 hmatrix[1][2] = this->matrix[2][1];
908 hmatrix[1][3] = this->matrix[3][1];
909
910 hmatrix[2][0] = this->matrix[0][2];
911 hmatrix[2][1] = this->matrix[1][2];
912 hmatrix[2][2] = this->matrix[2][2];
913 hmatrix[2][3] = this->matrix[3][2];
914
915 hmatrix[3][0] = this->matrix[0][3];
916 hmatrix[3][1] = this->matrix[1][3];
917 hmatrix[3][2] = this->matrix[2][3];
918 hmatrix[3][3] = this->matrix[3][3];
919
920 decomp_affine(hmatrix, &parts);
921
922 double mul = 1.0;
923 if (parts.t[W] != 0.0) mul = 1.0 / parts.t[W];
924 t[0] = parts.t[X] * mul;
925 t[1] = parts.t[Y] * mul;
926 t[2] = parts.t[Z] * mul;
927
928 r = parts.q;
929 mul = 1.0;
930 if (parts.k[W] != 0.0) mul = 1.0 / parts.k[W];
931 // mul be sign of determinant to support negative scales.
932 mul *= parts.f;
933 s[0] = parts.k[X] * mul;
934 s[1] = parts.k[Y] * mul;
935 s[2] = parts.k[Z] * mul;
936
937 so = parts.u;
938 }
939
940 /*!
941 Factor the matrix back into its \a translation, \a rotation,
942 \a scaleFactor and \a scaleorientation components. Will eliminate
943 the \a center variable from the matrix.
944
945 \sa factor()
946 */
947 void
getTransform(SbVec3d & translation,SbDPRotation & rotation,SbVec3d & scaleFactor,SbDPRotation & scaleOrientation,const SbVec3d & center) const948 SbDPMatrix::getTransform(SbVec3d & translation,
949 SbDPRotation & rotation,
950 SbVec3d & scaleFactor,
951 SbDPRotation & scaleOrientation,
952 const SbVec3d & center) const
953 {
954 SbDPMatrix m2 = *this;
955 SbDPMatrix trans;
956 trans.setTranslate(center);
957 m2.multLeft(trans);
958 trans.setTranslate(-center);
959 m2.multRight(trans);
960
961 m2.getTransform(translation, rotation, scaleFactor, scaleOrientation);
962 }
963
964 /*!
965 This function is not implemented in Coin.
966
967 \sa getTransform()
968 */
969 SbBool
factor(SbDPMatrix & COIN_UNUSED_ARG (r),SbVec3d & COIN_UNUSED_ARG (s),SbDPMatrix & COIN_UNUSED_ARG (u),SbVec3d & COIN_UNUSED_ARG (t),SbDPMatrix & COIN_UNUSED_ARG (proj))970 SbDPMatrix::factor(SbDPMatrix & COIN_UNUSED_ARG(r), SbVec3d & COIN_UNUSED_ARG(s), SbDPMatrix & COIN_UNUSED_ARG(u), SbVec3d & COIN_UNUSED_ARG(t),
971 SbDPMatrix & COIN_UNUSED_ARG(proj))
972 {
973 // FIXME: not implemented, not documented. 1998MMDD mortene.
974 COIN_STUB();
975 return FALSE;
976 }
977
978 /*!
979 This function produces a permuted LU decomposition of the matrix. It
980 uses the common single-row-pivoting strategy.
981
982 \a FALSE is returned if the matrix is singular, which it never is, because
983 of small adjustment values inserted if a singularity is found (as Open
984 Inventor does too).
985
986 The parity argument is always set to 1.0 or -1.0. Don't really know what
987 it's for, so it's not checked for correctness.
988
989 The index[] argument returns the permutation that was done on the matrix
990 to LU-decompose it. index[i] is the row that row i was swapped with at
991 step i in the decomposition, so index[] is not the actual permutation of
992 the row indexes!
993
994 BUGS: The function does not produce results that are numerically identical
995 with those produced by Open Inventor for the same matrices, because the
996 pivoting strategy in OI was never fully understood.
997
998 \sa SbDPMatrix::LUBackSubstitution
999 */
1000
1001
1002 SbBool
LUDecomposition(int index[4],double & d)1003 SbDPMatrix::LUDecomposition(int index[4], double & d)
1004 {
1005 int i;
1006 for (i = 0; i < 4; i++) index[i] = i;
1007 d = 1.0;
1008
1009 const double MINIMUM_PIVOT = 1e-6f; // Inventor fix...
1010
1011 for (int row = 1; row < 4; row++) {
1012 int swap_row = row;
1013 double max_pivot = 0.0;
1014 for (int test_row = row; test_row < 4; test_row++) {
1015 const double test_pivot = SbAbs(matrix[test_row][row]);
1016 if (test_pivot > max_pivot) {
1017 swap_row = test_row;
1018 max_pivot = test_pivot;
1019 }
1020 }
1021
1022 // swap rows
1023 if (swap_row != row) {
1024 d = -d;
1025 index[row] = swap_row;
1026 for (i = 0; i < 4; i++)
1027 SbSwap(matrix[row][i], matrix[swap_row][i]);
1028 }
1029
1030 double pivot = matrix[row][row];
1031 if (matrix[row][row] == 0.0) {
1032 // return FALSE;
1033 // instead of returning FALSE on singulars...
1034 matrix[row][row] = pivot = MINIMUM_PIVOT;
1035 }
1036
1037 // the factorization
1038 for (i = row + 1; i < 4; i++) {
1039 const double factor = (matrix[i][row] /= pivot);
1040 for (int j = row + 1; j < 4; j++)
1041 matrix[i][j] -= factor * matrix[row][j];
1042 }
1043 }
1044 return TRUE;
1045 }
1046
1047 /*!
1048 This function does a solve on the "Ax = b" system, given that the matrix
1049 is LU-decomposed in advance. First, a forward substitution is done on the
1050 lower system (Ly = b), and then a backwards substitution is done on the
1051 upper triangular system (Ux = y).
1052
1053 The index[] argument is the one returned from
1054 SbDPMatrix::LUDecomposition(), so see that function for an explanation.
1055
1056 The b[] argument must contain the b vector in "Ax = b" when calling the
1057 function. After the function has solved the system, the b[] vector contains
1058 the x vector.
1059
1060 BUGS: As is done by Open Inventor, unsolvable x values will not return
1061 NaN but 0.
1062 */
1063
1064 void
LUBackSubstitution(int index[4],double b[4]) const1065 SbDPMatrix::LUBackSubstitution(int index[4], double b[4]) const
1066 {
1067 int i;
1068
1069 // permute b[] the way matrix[][] is permuted
1070 for (i = 0; i < 4; i++)
1071 if (i != index[i])
1072 SbSwap(b[i], b[index[i]]);
1073
1074 // forward substitution on L (Ly = b)
1075 double y[4];
1076 for (i = 0; i < 4; i++) {
1077 double sum = 0.0;
1078 for (int j = 0; j < i; j++)
1079 sum += matrix[i][j] * y[j];
1080 y[i] = b[i] - sum;
1081 }
1082
1083 // backwards substitution on U (Ux = y)
1084 double x[4];
1085 for (i = 3; i >= 0; i--) {
1086 double sum = 0.0;
1087 for (int j = i + 1; j < 4; j++)
1088 sum += matrix[i][j] * x[j];
1089 if (matrix[i][i] != 0.0)
1090 x[i] = (y[i] - sum) / matrix[i][i];
1091 else
1092 x[i] = 0.0;
1093 }
1094
1095 // de-permute x[]? doesn't look like it
1096 // for (i = 3; i >= 0; i--)
1097 // if (i != index[i])
1098 // SbSwap(x[i], x[index[i]]);
1099
1100 // copy x[] into b[] for "return to sender"
1101 for (i = 0; i < 4; i++) b[i] = x[i];
1102 }
1103
1104 /*!
1105 Returns the transpose of this matrix.
1106 */
1107
1108 SbDPMatrix
transpose(void) const1109 SbDPMatrix::transpose(void) const
1110 {
1111 SbDPMatrix trans = (*this);
1112
1113 for (int i=0; i < 3; i++) {
1114 for (int j=i+1; j < 4; j++) {
1115 SbSwap(trans[i][j], trans[j][i]);
1116 }
1117 }
1118
1119 return trans;
1120 }
1121
1122 /*!
1123 Let this matrix be right-multiplied by \a m. Returns reference to
1124 self.
1125
1126 \sa multLeft()
1127 */
1128 SbDPMatrix &
multRight(const SbDPMatrix & m)1129 SbDPMatrix::multRight(const SbDPMatrix & m)
1130 {
1131 // Checks if one or the other matrix is equal to the identity matrix
1132 // before multiplying them. We do this because it's a major
1133 // optimization if one of them _is_, and the isIdentity() check
1134 // should be very quick in the common case where a matrix is not the
1135 // identity matrix.
1136 const SbDPMat & mfm = m.matrix;
1137 if (SbDPMatrix_isIdentity(mfm)) { return *this; }
1138 SbDPMat & tfm = this->matrix;
1139 if (SbDPMatrix_isIdentity(tfm)) { *this = m; return *this; }
1140
1141 SbDPMat tmp;
1142 (void)memcpy(tmp, tfm, 4*4*sizeof(double));
1143
1144 for (int i=0; i < 4; i++) {
1145 for (int j=0; j < 4; j++) {
1146 tfm[i][j] =
1147 tmp[i][0] * mfm[0][j] +
1148 tmp[i][1] * mfm[1][j] +
1149 tmp[i][2] * mfm[2][j] +
1150 tmp[i][3] * mfm[3][j];
1151 }
1152 }
1153
1154 return *this;
1155 }
1156
1157 /*!
1158 Let this matrix be left-multiplied by \a m. Returns reference to
1159 self.
1160
1161 \sa multRight()
1162 */
1163 SbDPMatrix&
multLeft(const SbDPMatrix & m)1164 SbDPMatrix::multLeft(const SbDPMatrix & m)
1165 {
1166 // Checks if one or the other matrix is equal to the identity
1167 // matrix. See also code comments at the start of
1168 // SbDPMatrix::multRight().
1169 const SbDPMat & mfm = m.matrix;
1170 if (SbDPMatrix_isIdentity(mfm)) { return *this; }
1171 SbDPMat & tfm = this->matrix;
1172 if (SbDPMatrix_isIdentity(tfm)) { *this = m; return *this; }
1173
1174 SbDPMat tmp;
1175 (void)memcpy(tmp, tfm, 4*4*sizeof(double));
1176
1177 for (int i=0; i < 4; i++) {
1178 for (int j=0; j < 4; j++) {
1179 tfm[i][j] =
1180 tmp[0][j] * mfm[i][0] +
1181 tmp[1][j] * mfm[i][1] +
1182 tmp[2][j] * mfm[i][2] +
1183 tmp[3][j] * mfm[i][3];
1184 }
1185 }
1186 return *this;
1187 }
1188
1189 /*!
1190 Multiply \a src vector with this matrix and return the result in \a dst.
1191 Multiplication is done with the vector on the right side of the
1192 expression, i.e. dst = M * src.
1193
1194 \sa multVecMatrix(), multDirMatrix() and multLineMatrix().
1195 */
1196 void
multMatrixVec(const SbVec3d & src,SbVec3d & dst) const1197 SbDPMatrix::multMatrixVec(const SbVec3d & src, SbVec3d & dst) const
1198 {
1199 // Checks if the "this" matrix is equal to the identity matrix. See
1200 // also code comments at the start of SbDPMatrix::multRight().
1201 if (SbDPMatrix_isIdentity(this->matrix)) { dst = src; return; }
1202
1203 const double * t0 = (*this)[0];
1204 const double * t1 = (*this)[1];
1205 const double * t2 = (*this)[2];
1206 const double * t3 = (*this)[3];
1207 // Copy the src vector, just in case src and dst is the same vector.
1208 SbVec3d s = src;
1209
1210 double W = s[0]*t3[0] + s[1]*t3[1] + s[2]*t3[2] + t3[3];
1211
1212 dst[0] = (s[0]*t0[0] + s[1]*t0[1] + s[2]*t0[2] + t0[3])/W;
1213 dst[1] = (s[0]*t1[0] + s[1]*t1[1] + s[2]*t1[2] + t1[3])/W;
1214 dst[2] = (s[0]*t2[0] + s[1]*t2[1] + s[2]*t2[2] + t2[3])/W;
1215 }
1216
1217 /*!
1218 Multiply \a src vector with this matrix and return the result in \a dst.
1219 Multiplication is done with the vector on the left side of the
1220 expression, i.e. dst = src * M.
1221
1222 It is safe to let \a src and \a dst be the same SbVec3d instance.
1223
1224 \sa multMatrixVec(), multDirMatrix() and multLineMatrix().
1225 */
1226 void
multVecMatrix(const SbVec3d & src,SbVec3d & dst) const1227 SbDPMatrix::multVecMatrix(const SbVec3d & src, SbVec3d & dst) const
1228 {
1229 // Checks if the "this" matrix is equal to the identity matrix. See
1230 // also code comments at the start of SbDPMatrix::multRight().
1231 if (SbDPMatrix_isIdentity(this->matrix)) { dst = src; return; }
1232
1233 const double * t0 = this->matrix[0];
1234 const double * t1 = this->matrix[1];
1235 const double * t2 = this->matrix[2];
1236 const double * t3 = this->matrix[3];
1237 // Copy the src vector, just in case src and dst is the same vector.
1238 SbVec3d s = src;
1239
1240 double W = s[0]*t0[3] + s[1]*t1[3] + s[2]*t2[3] + t3[3];
1241
1242 dst[0] = (s[0]*t0[0] + s[1]*t1[0] + s[2]*t2[0] + t3[0])/W;
1243 dst[1] = (s[0]*t0[1] + s[1]*t1[1] + s[2]*t2[1] + t3[1])/W;
1244 dst[2] = (s[0]*t0[2] + s[1]*t1[2] + s[2]*t2[2] + t3[2])/W;
1245 }
1246
1247 /*!
1248 \overload
1249 */
1250 void
multVecMatrix(const SbVec4d & src,SbVec4d & dst) const1251 SbDPMatrix::multVecMatrix(const SbVec4d & src, SbVec4d & dst) const
1252 {
1253 // Checks if the "this" matrix is equal to the identity matrix. See
1254 // also code comments at the start of SbDPMatrix::multRight().
1255 if (SbDPMatrix_isIdentity(this->matrix)) { dst = src; return; }
1256
1257 const double * t0 = (*this)[0];
1258 const double * t1 = (*this)[1];
1259 const double * t2 = (*this)[2];
1260 const double * t3 = (*this)[3];
1261
1262 SbVec4d s = src;
1263
1264 dst[0] = (s[0]*t0[0] + s[1]*t1[0] + s[2]*t2[0] + s[3]*t3[0]);
1265 dst[1] = (s[0]*t0[1] + s[1]*t1[1] + s[2]*t2[1] + s[3]*t3[1]);
1266 dst[2] = (s[0]*t0[2] + s[1]*t1[2] + s[2]*t2[2] + s[3]*t3[2]);
1267 dst[3] = (s[0]*t0[3] + s[1]*t1[3] + s[2]*t2[3] + s[3]*t3[3]);
1268 }
1269
1270 /*!
1271 Multiplies \a src by the matrix. \a src is assumed to be a direction
1272 vector, and the translation components of the matrix are therefore
1273 ignored.
1274
1275 Multiplication is done with the vector on the left side of the
1276 expression, i.e. dst = src * M.
1277
1278 \sa multVecMatrix(), multMatrixVec() and multLineMatrix().
1279 */
1280 void
multDirMatrix(const SbVec3d & src,SbVec3d & dst) const1281 SbDPMatrix::multDirMatrix(const SbVec3d & src, SbVec3d & dst) const
1282 {
1283 // Checks if the "this" matrix is equal to the identity matrix. See
1284 // also code comments at the start of SbDPMatrix::multRight().
1285 if (SbDPMatrix_isIdentity(this->matrix)) { dst = src; return; }
1286
1287 const double * t0 = (*this)[0];
1288 const double * t1 = (*this)[1];
1289 const double * t2 = (*this)[2];
1290 // Copy the src vector, just in case src and dst is the same vector.
1291 SbVec3d s = src;
1292
1293 dst[0] = s[0]*t0[0] + s[1]*t1[0] + s[2]*t2[0];
1294 dst[1] = s[0]*t0[1] + s[1]*t1[1] + s[2]*t2[1];
1295 dst[2] = s[0]*t0[2] + s[1]*t1[2] + s[2]*t2[2];
1296 }
1297
1298 /*!
1299 Multiplies line point with the full matrix and multiplies the
1300 line direction with the matrix without the translation components.
1301
1302 \sa multVecMatrix(), multMatrixVec() and multDirMatrix().
1303 */
1304 void
multLineMatrix(const SbDPLine & src,SbDPLine & dst) const1305 SbDPMatrix::multLineMatrix(const SbDPLine & src, SbDPLine & dst) const
1306 {
1307 SbVec3d newpt, newdir;
1308 this->multVecMatrix(src.getPosition(), newpt);
1309 this->multDirMatrix(src.getDirection(), newdir);
1310
1311 dst.setPosDir(newpt, newdir);
1312 }
1313
1314 /*!
1315 Write out the matrix contents to the given file.
1316 */
1317 void
print(FILE * fp) const1318 SbDPMatrix::print(FILE * fp) const
1319 {
1320 for (int i=0; i < 4; i++) {
1321 fprintf(fp, "%10.5g\t%10.5g\t%10.5g\t%10.5g\n",
1322 this->matrix[i][0], this->matrix[i][1],
1323 this->matrix[i][2], this->matrix[i][3]);
1324 }
1325 }
1326
1327 /***********************************************************************
1328 below is the polar_decomp implementation by Ken Shoemake
1329 <shoemake@graphics.cis.upenn.edu>. It was part of the
1330 Graphics Gems IV archive.
1331 ************************************************************************/
1332
1333 // FIXME: should merge all the PD code we're using from GGIV into
1334 // SbDPMatrix, SbDPRotation and SbVec3d proper (for two reasons: 1)
1335 // there's a lot of duplicated code here (like for instance the
1336 // matrix->quaternion decomposition, which also exists in
1337 // SbDPRotation::setValue(SbDPMatrix&)), and 2) the remaining code
1338 // snippets look generally useful outside the purpose of breaking down
1339 // a matrix into it's transformation components). 20010114 mortene.
1340
1341
1342 /**** Decompose.c ****/
1343 /* Ken Shoemake, 1993 */
1344
1345 /******* Matrix Preliminaries *******/
1346
1347 /** Fill out 3x3 matrix to 4x4 **/
1348 #define mat_pad(A) (A[W][X]=A[X][W]=A[W][Y]=A[Y][W]=A[W][Z]=A[Z][W]=0, A[W][W]=1)
1349
1350 /** Copy nxn matrix A to C using "gets" for assignment **/
1351 #define mat_copy(C, gets, A, n) {int i, j; for (i=0;i<n;i++) for (j=0;j<n;j++)\
1352 C[i][j] gets (A[i][j]);}
1353
1354 /** Copy transpose of nxn matrix A to C using "gets" for assignment **/
1355 #define mat_tpose(AT, gets, A, n) {int i, j; for (i=0;i<n;i++) for (j=0;j<n;j++)\
1356 AT[i][j] gets (A[j][i]);}
1357
1358 /** Assign nxn matrix C the element-wise combination of A and B using "op" **/
1359 #define mat_binop(C, gets, A, op, B, n) {int i, j; for (i=0;i<n;i++) for (j=0;j<n;j++)\
1360 C[i][j] gets (A[i][j]) op (B[i][j]);}
1361
1362 /** Multiply the upper left 3x3 parts of A and B to get AB **/
1363 static void
mat_mult(HMatrix A,HMatrix B,HMatrix AB)1364 mat_mult(HMatrix A, HMatrix B, HMatrix AB)
1365 {
1366 int i, j;
1367 for (i=0; i<3; i++) for (j=0; j<3; j++)
1368 AB[i][j] = A[i][0]*B[0][j] + A[i][1]*B[1][j] + A[i][2]*B[2][j];
1369 }
1370
1371 /** Return dot product of length 3 vectors va and vb **/
1372 static double
vdot(double * va,double * vb)1373 vdot(double * va, double * vb)
1374 {
1375 return (va[0]*vb[0] + va[1]*vb[1] + va[2]*vb[2]);
1376 }
1377
1378 /** Set v to cross product of length 3 vectors va and vb **/
1379 static void
vcross(double * va,double * vb,double * v)1380 vcross(double * va, double * vb, double * v)
1381 {
1382 v[0] = va[1]*vb[2] - va[2]*vb[1];
1383 v[1] = va[2]*vb[0] - va[0]*vb[2];
1384 v[2] = va[0]*vb[1] - va[1]*vb[0];
1385 }
1386
1387 /** Set MadjT to transpose of inverse of M times determinant of M **/
1388 static void
adjoint_transpose(HMatrix M,HMatrix MadjT)1389 adjoint_transpose(HMatrix M, HMatrix MadjT)
1390 {
1391 vcross(M[1], M[2], MadjT[0]);
1392 vcross(M[2], M[0], MadjT[1]);
1393 vcross(M[0], M[1], MadjT[2]);
1394 }
1395
1396 /******* Decomp Auxiliaries *******/
1397
1398 static HMatrix mat_id = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
1399
1400 /** Compute either the 1 or infinity norm of M, depending on tpose **/
1401 static double
mat_norm(HMatrix M,int tpose)1402 mat_norm(HMatrix M, int tpose)
1403 {
1404 int i;
1405 double sum, max;
1406 max = 0.0;
1407 for (i=0; i<3; i++) {
1408 if (tpose) sum = static_cast<double>(fabs(M[0][i])+fabs(M[1][i])+fabs(M[2][i]));
1409 else sum = static_cast<double>(fabs(M[i][0])+fabs(M[i][1])+fabs(M[i][2]));
1410 if (max<sum) max = sum;
1411 }
1412 return max;
1413 }
1414
norm_inf(HMatrix M)1415 static double norm_inf(HMatrix M) {return mat_norm(M, 0);}
norm_one(HMatrix M)1416 static double norm_one(HMatrix M) {return mat_norm(M, 1);}
1417
1418 /** Return index of column of M containing maximum abs entry, or -1 if M=0 **/
1419 static int
find_max_col(HMatrix M)1420 find_max_col(HMatrix M)
1421 {
1422 double abs, max;
1423 int i, j, col;
1424 max = 0.0; col = -1;
1425 for (i=0; i<3; i++) for (j=0; j<3; j++) {
1426 abs = M[i][j]; if (abs<0.0) abs = -abs;
1427 if (abs>max) {max = abs; col = j;}
1428 }
1429 return col;
1430 }
1431
1432 /** Setup u for Household reflection to zero all v components but first **/
1433 static void
make_reflector(double * v,double * u)1434 make_reflector(double * v, double * u)
1435 {
1436 double s = static_cast<double>(sqrt(vdot(v, v)));
1437 u[0] = v[0]; u[1] = v[1];
1438 u[2] = v[2] + ((v[2]<0.0) ? -s : s);
1439 s = static_cast<double>(sqrt(2.0/vdot(u, u)));
1440 u[0] = u[0]*s; u[1] = u[1]*s; u[2] = u[2]*s;
1441 }
1442
1443 /** Apply Householder reflection represented by u to column vectors of M **/
1444 static void
reflect_cols(HMatrix M,double * u)1445 reflect_cols(HMatrix M, double * u)
1446 {
1447 int i, j;
1448 for (i=0; i<3; i++) {
1449 double s = u[0]*M[0][i] + u[1]*M[1][i] + u[2]*M[2][i];
1450 for (j=0; j<3; j++) M[j][i] -= u[j]*s;
1451 }
1452 }
1453 /** Apply Householder reflection represented by u to row vectors of M **/
1454 static void
reflect_rows(HMatrix M,double * u)1455 reflect_rows(HMatrix M, double * u)
1456 {
1457 int i, j;
1458 for (i=0; i<3; i++) {
1459 double s = vdot(u, M[i]);
1460 for (j=0; j<3; j++) M[i][j] -= u[j]*s;
1461 }
1462 }
1463
1464 /** Find orthogonal factor Q of rank 1 (or less) M **/
1465 static void
do_rank1(HMatrix M,HMatrix Q)1466 do_rank1(HMatrix M, HMatrix Q)
1467 {
1468 double v1[3], v2[3], s;
1469 int col;
1470 mat_copy(Q, =, mat_id, 4);
1471 /* If rank(M) is 1, we should find a non-zero column in M */
1472 col = find_max_col(M);
1473 if (col<0) return; /* Rank is 0 */
1474 v1[0] = M[0][col]; v1[1] = M[1][col]; v1[2] = M[2][col];
1475 make_reflector(v1, v1); reflect_cols(M, v1);
1476 v2[0] = M[2][0]; v2[1] = M[2][1]; v2[2] = M[2][2];
1477 make_reflector(v2, v2); reflect_rows(M, v2);
1478 s = M[2][2];
1479 if (s<0.0) Q[2][2] = -1.0;
1480 reflect_cols(Q, v1); reflect_rows(Q, v2);
1481 }
1482
1483 /** Find orthogonal factor Q of rank 2 (or less) M using adjoint transpose **/
1484 static void
do_rank2(HMatrix M,HMatrix MadjT,HMatrix Q)1485 do_rank2(HMatrix M, HMatrix MadjT, HMatrix Q)
1486 {
1487 double v1[3], v2[3];
1488 double w, x, y, z, c, s, d;
1489 int col;
1490 /* If rank(M) is 2, we should find a non-zero column in MadjT */
1491 col = find_max_col(MadjT);
1492 if (col<0) {do_rank1(M, Q); return;} /* Rank<2 */
1493 v1[0] = MadjT[0][col]; v1[1] = MadjT[1][col]; v1[2] = MadjT[2][col];
1494 make_reflector(v1, v1); reflect_cols(M, v1);
1495 vcross(M[0], M[1], v2);
1496 make_reflector(v2, v2); reflect_rows(M, v2);
1497 w = M[0][0]; x = M[0][1]; y = M[1][0]; z = M[1][1];
1498 if (w*z>x*y) {
1499 c = z+w; s = y-x; d = static_cast<double>(sqrt(c*c+s*s)); c = c/d; s = s/d;
1500 Q[0][0] = Q[1][1] = c; Q[0][1] = -(Q[1][0] = s);
1501 }
1502 else {
1503 c = z-w; s = y+x; d = static_cast<double>(sqrt(c*c+s*s)); c = c/d; s = s/d;
1504 Q[0][0] = -(Q[1][1] = c); Q[0][1] = Q[1][0] = s;
1505 }
1506 Q[0][2] = Q[2][0] = Q[1][2] = Q[2][1] = 0.0; Q[2][2] = 1.0;
1507 reflect_cols(Q, v1); reflect_rows(Q, v2);
1508 }
1509
1510
1511 /******* Polar Decomposition *******/
1512
1513 /* Polar Decomposition of 3x3 matrix in 4x4,
1514 * M = QS. See Nicholas Higham and Robert S. Schreiber,
1515 * Fast Polar Decomposition of An Arbitrary Matrix,
1516 * Technical Report 88-942, October 1988,
1517 * Department of Computer Science, Cornell University.
1518 */
1519 static double
polar_decomp(HMatrix M,HMatrix Q,HMatrix S)1520 polar_decomp(HMatrix M, HMatrix Q, HMatrix S)
1521 {
1522 #define TOL 1.0e-6
1523 HMatrix Mk, MadjTk, Ek;
1524 double det, M_one, M_inf, MadjT_one, MadjT_inf, E_one, gamma, g1, g2;
1525 int i, j;
1526 mat_tpose(Mk, =, M, 3);
1527 M_one = norm_one(Mk); M_inf = norm_inf(Mk);
1528 do {
1529 adjoint_transpose(Mk, MadjTk);
1530 det = vdot(Mk[0], MadjTk[0]);
1531 if (det==0.0) {do_rank2(Mk, MadjTk, Mk); break;}
1532 MadjT_one = norm_one(MadjTk); MadjT_inf = norm_inf(MadjTk);
1533 gamma = static_cast<double>(sqrt(sqrt((MadjT_one*MadjT_inf)/(M_one*M_inf))/fabs(det)));
1534 g1 = gamma*0.5;
1535 g2 = 0.5/(gamma*det);
1536 mat_copy(Ek, =, Mk, 3);
1537 mat_binop(Mk, =, g1*Mk, +, g2*MadjTk, 3);
1538 mat_copy(Ek, -=, Mk, 3);
1539 E_one = norm_one(Ek);
1540 M_one = norm_one(Mk); M_inf = norm_inf(Mk);
1541 } while (E_one>(M_one*TOL));
1542 mat_tpose(Q, =, Mk, 3); mat_pad(Q);
1543 mat_mult(Mk, M, S); mat_pad(S);
1544 for (i=0; i<3; i++) for (j=i; j<3; j++)
1545 S[i][j] = S[j][i] = 0.5*(S[i][j]+S[j][i]);
1546 return (det);
1547 }
1548
1549 /******* Spectral Decomposition *******/
1550
1551 /* Compute the spectral decomposition of symmetric positive semi-definite S.
1552 * Returns rotation in U and scale factors in result, so that if K is a diagonal
1553 * matrix of the scale factors, then S = U K (U transpose). Uses Jacobi method.
1554 * See Gene H. Golub and Charles F. Van Loan. Matrix Computations. Hopkins 1983.
1555 */
1556 static SbVec4d
spect_decomp(HMatrix S,HMatrix U)1557 spect_decomp(HMatrix S, HMatrix U)
1558 {
1559 SbVec4d kv;
1560 double Diag[3], OffD[3]; /* OffD is off-diag (by omitted index) */
1561 double g, h, fabsh, fabsOffDi, t, theta, c, s, tau, ta, OffDq, a, b;
1562 static char nxt[] = {Y, Z, X};
1563 int sweep, i, j;
1564 mat_copy(U, =, mat_id, 4);
1565 Diag[X] = S[X][X]; Diag[Y] = S[Y][Y]; Diag[Z] = S[Z][Z];
1566 OffD[X] = S[Y][Z]; OffD[Y] = S[Z][X]; OffD[Z] = S[X][Y];
1567 for (sweep=20; sweep>0; sweep--) {
1568 double sm = static_cast<double>(fabs(OffD[X])+fabs(OffD[Y])+fabs(OffD[Z]));
1569 if (sm==0.0) break;
1570 for (i=Z; i>=X; i--) {
1571 int p = nxt[i]; int q = nxt[p];
1572 fabsOffDi = fabs(OffD[i]);
1573 g = 100.0*fabsOffDi;
1574 if (fabsOffDi>0.0) {
1575 h = Diag[q] - Diag[p];
1576 fabsh = fabs(h);
1577 if (fabsh+g==fabsh) {
1578 t = OffD[i]/h;
1579 }
1580 else {
1581 theta = 0.5*h/OffD[i];
1582 t = 1.0/(fabs(theta)+sqrt(theta*theta+1.0));
1583 if (theta<0.0) t = -t;
1584 }
1585 c = 1.0/sqrt(t*t+1.0); s = t*c;
1586 tau = s/(c+1.0);
1587 ta = t*OffD[i]; OffD[i] = 0.0;
1588 Diag[p] -= ta; Diag[q] += ta;
1589 OffDq = OffD[q];
1590 OffD[q] -= s*(OffD[p] + tau*OffD[q]);
1591 OffD[p] += s*(OffDq - tau*OffD[p]);
1592 for (j=Z; j>=X; j--) {
1593 a = U[j][p]; b = U[j][q];
1594 U[j][p] -= static_cast<double>(s*(b + tau*a));
1595 U[j][q] += static_cast<double>(s*(a - tau*b));
1596 }
1597 }
1598 }
1599 }
1600 kv[X] = static_cast<double>(Diag[X]);
1601 kv[Y] = static_cast<double>(Diag[Y]);
1602 kv[Z] = static_cast<double>(Diag[Z]);
1603 kv[W] = 1.0;
1604 return (kv);
1605 }
1606
1607 /* Helper function for the snuggle() function below. */
1608 static inline void
cycle(double * a,SbBool flip)1609 cycle(double * a, SbBool flip)
1610 {
1611 if (flip) {
1612 a[3]=a[0]; a[0]=a[1]; a[1]=a[2]; a[2]=a[3];
1613 }
1614 else {
1615 a[3]=a[2]; a[2]=a[1]; a[1]=a[0]; a[0]=a[3];
1616 }
1617 }
1618
1619 /******* Spectral Axis Adjustment *******/
1620
1621 /* Given a unit quaternion, q, and a scale vector, k, find a unit quaternion, p,
1622 * which permutes the axes and turns freely in the plane of duplicate scale
1623 * factors, such that q p has the largest possible w component, i.e. the
1624 * smallest possible angle. Permutes k's components to go with q p instead of q.
1625 * See Ken Shoemake and Tom Duff. Matrix Animation and Polar Decomposition.
1626 * Proceedings of Graphics Interface 1992. Details on p. 262-263.
1627 */
1628 static SbDPRotation
snuggle(SbDPRotation q,SbVec4d & k)1629 snuggle(SbDPRotation q, SbVec4d & k)
1630 {
1631 #define SQRTHALF (0.7071067811865475244f)
1632 #define sgn(n, v) ((n)?-(v):(v))
1633 #define swap(a, i, j) {a[3]=a[i]; a[i]=a[j]; a[j]=a[3];}
1634
1635 SbDPRotation p;
1636 double ka[4];
1637 int i, turn = -1;
1638 ka[X] = k[X]; ka[Y] = k[Y]; ka[Z] = k[Z];
1639 if (ka[X]==ka[Y]) {if (ka[X]==ka[Z]) turn = W; else turn = Z;}
1640 else {if (ka[X]==ka[Z]) turn = Y; else if (ka[Y]==ka[Z]) turn = X;}
1641 if (turn>=0) {
1642 SbDPRotation qtoz, qp;
1643 unsigned neg[3], win;
1644 double mag[3], t;
1645 static SbDPRotation qxtoz(0.0, SQRTHALF, 0.0, SQRTHALF);
1646 static SbDPRotation qytoz(SQRTHALF, 0.0, 0.0, SQRTHALF);
1647 static SbDPRotation qppmm(0.5, 0.5, -0.5, -0.5);
1648 static SbDPRotation qpppp(0.5, 0.5, 0.5, 0.5);
1649 static SbDPRotation qmpmm(-0.5, 0.5, -0.5, -0.5);
1650 static SbDPRotation qpppm(0.5, 0.5, 0.5, -0.5);
1651 static SbDPRotation q0001(0.0, 0.0, 0.0, 1.0);
1652 static SbDPRotation q1000(1.0, 0.0, 0.0, 0.0);
1653 switch (turn) {
1654 default: return SbDPRotation(q).invert();
1655 case X: q = (qtoz = qxtoz) * q; swap(ka, X, Z) break;
1656 case Y: q = (qtoz = qytoz) * q; swap(ka, Y, Z) break;
1657 case Z: qtoz = q0001; break;
1658 }
1659 q.invert();
1660 mag[0] = static_cast<double>(q.getValue()[Z]*q.getValue()[Z])+static_cast<double>(q.getValue()[W]*q.getValue()[W]-0.5);
1661 mag[1] = static_cast<double>(q.getValue()[X]*q.getValue()[Z])-static_cast<double>(q.getValue()[Y]*q.getValue()[W]);
1662 mag[2] = static_cast<double>(q.getValue()[Y]*q.getValue()[Z]+static_cast<double>(q.getValue()[X]*q.getValue()[W]));
1663 for (i=0; i<3; i++) if ((neg[i] = (mag[i] < 0.0))) mag[i] = -mag[i];
1664 if (mag[0]>mag[1]) {if (mag[0]>mag[2]) win = 0; else win = 2;}
1665 else {if (mag[1]>mag[2]) win = 1; else win = 2;}
1666 switch (win) {
1667 case 0: if (neg[0]) p = q1000; else p = q0001; break;
1668 case 1: if (neg[1]) p = qppmm; else p = qpppp; cycle(ka, FALSE); break;
1669 case 2: if (neg[2]) p = qmpmm; else p = qpppm; cycle(ka, TRUE); break;
1670 }
1671 qp = p * q;
1672 t = sqrt(mag[win]+0.5);
1673 p = SbDPRotation(0.0, 0.0, -qp.getValue()[Z]/static_cast<double>(t), qp.getValue()[W]/static_cast<double>(t)) * p;
1674 p = SbDPRotation(p).invert() * qtoz;
1675 }
1676 else {
1677 double qa[4], pa[4];
1678 unsigned lo, hi, neg[4], par = 0;
1679 double all, big, two;
1680 qa[0] = q.getValue()[X]; qa[1] = q.getValue()[Y]; qa[2] = q.getValue()[Z]; qa[3] = q.getValue()[W];
1681 for (i=0; i<4; i++) {
1682 pa[i] = 0.0;
1683 if ((neg[i] = (qa[i]<0.0))) qa[i] = -qa[i];
1684 par ^= neg[i];
1685 }
1686 /* Find two largest components, indices in hi and lo */
1687 if (qa[0]>qa[1]) lo = 0; else lo = 1;
1688 if (qa[2]>qa[3]) hi = 2; else hi = 3;
1689 if (qa[lo]>qa[hi]) {
1690 if (qa[lo^1]>qa[hi]) {hi = lo; lo ^= 1;}
1691 else {hi ^= lo; lo ^= hi; hi ^= lo;}
1692 } else {if (qa[hi^1]>qa[lo]) lo = hi^1;}
1693 all = (qa[0]+qa[1]+qa[2]+qa[3])*0.5;
1694 two = (qa[hi]+qa[lo])*SQRTHALF;
1695 big = qa[hi];
1696 if (all>two) {
1697 if (all>big) {/*all*/
1698 {int i; for (i=0; i<4; i++) pa[i] = sgn(neg[i], 0.5);}
1699 cycle(ka, par);
1700 }
1701 else {/*big*/ pa[hi] = sgn(neg[hi], 1.0);}
1702 }
1703 else {
1704 if (two>big) {/*two*/
1705 pa[hi] = sgn(neg[hi], SQRTHALF); pa[lo] = sgn(neg[lo], SQRTHALF);
1706 if (lo>hi) {hi ^= lo; lo ^= hi; hi ^= lo;}
1707 if (hi==W) {hi = "\001\002\000"[lo]; lo = 3-hi-lo;}
1708 swap(ka, hi, lo)
1709 } else {/*big*/ pa[hi] = sgn(neg[hi], 1.0);}
1710 }
1711 // FIXME: p = conjugate(pa)? 20010114 mortene.
1712 p.setValue(-pa[0], -pa[1], -pa[2], pa[3]);
1713 }
1714 k[X] = ka[X]; k[Y] = ka[Y]; k[Z] = ka[Z];
1715 return (p);
1716 }
1717
1718 /******* Decompose Affine Matrix *******/
1719
1720 /* Decompose 4x4 affine matrix A as TFRUK(U transpose), where t contains the
1721 * translation components, q contains the rotation R, u contains U, k contains
1722 * scale factors, and f contains the sign of the determinant.
1723 * Assumes A transforms column vectors in right-handed coordinates.
1724 * See Ken Shoemake and Tom Duff. Matrix Animation and Polar Decomposition.
1725 * Proceedings of Graphics Interface 1992.
1726 */
1727 static void
decomp_affine(HMatrix A,AffineParts * parts)1728 decomp_affine(HMatrix A, AffineParts * parts)
1729 {
1730 HMatrix Q, S, U;
1731 SbDPRotation p;
1732 parts->t = SbVec4d(A[X][W], A[Y][W], A[Z][W], 0);
1733 double det = polar_decomp(A, Q, S);
1734 if (det<0.0) {
1735 mat_copy(Q, =, -Q, 3);
1736 parts->f = -1;
1737 }
1738 else parts->f = 1;
1739
1740 // Transpose for our code (we use OpenGL's convention for numbering
1741 // rows and columns).
1742
1743 // reinterpret_cast to humor MSVC6, this should have no effect on
1744 // other compilers.
1745 SbDPMatrix TQ(reinterpret_cast<const SbDPMat *>(&Q));
1746 parts->q = SbDPRotation(TQ.transpose());
1747 parts->k = spect_decomp(S, U);
1748 // Transpose for our code (we use OpenGL's convention for numbering
1749 // rows and columns).
1750
1751 // reinterpret_cast to humor MSVC6, this should have no effect on
1752 // other compilers.
1753 SbDPMatrix TU(reinterpret_cast<const SbDPMat *>(&U));
1754 parts->u = SbDPRotation(TU.transpose());
1755 p = snuggle(parts->u, parts->k);
1756 parts->u = p * parts->u;
1757 }
1758
1759 #undef mat_pad
1760 #undef mat_copy
1761 #undef mat_tpose
1762 #undef mat_binop
1763 #undef TOL
1764 #undef SQRTHALF
1765 #undef sgn
1766 #undef swap
1767
1768 #ifdef COIN_TEST_SUITE
1769 #include <Inventor/SbMatrix.h>
1770
BOOST_AUTO_TEST_CASE(constructFromSbMatrix)1771 BOOST_AUTO_TEST_CASE(constructFromSbMatrix) {
1772 SbMatrix a(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15);
1773 SbMatrixd b;
1774 double c[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};
1775 b.setValue(c);
1776 SbDPMatrix d = SbMatrixd(a);
1777 BOOST_CHECK_MESSAGE(b == d,
1778 "Equality comparrison failed!");
1779 }
1780 #endif //COIN_TEST_SUITE
1781