1 | /*****************************************************************************/ |
2 | // Copyright 2006-2008 Adobe Systems Incorporated |
3 | // All Rights Reserved. |
4 | // |
5 | // NOTICE: Adobe permits you to use, modify, and distribute this file in |
6 | // accordance with the terms of the Adobe license agreement accompanying it. |
7 | /*****************************************************************************/ |
8 | |
9 | /* $Id: //mondo/dng_sdk_1_4/dng_sdk/source/dng_matrix.cpp#1 $ */ |
10 | /* $DateTime: 2012/05/30 13:28:51 $ */ |
11 | /* $Change: 832332 $ */ |
12 | /* $Author: tknoll $ */ |
13 | |
14 | /*****************************************************************************/ |
15 | |
16 | #include "dng_matrix.h" |
17 | |
18 | #include "dng_exceptions.h" |
19 | #include "dng_utils.h" |
20 | |
21 | /*****************************************************************************/ |
22 | |
23 | dng_matrix::dng_matrix () |
24 | |
25 | : fRows (0) |
26 | , fCols (0) |
27 | |
28 | { |
29 | |
30 | } |
31 | |
32 | /*****************************************************************************/ |
33 | |
34 | dng_matrix::dng_matrix (uint32 rows, |
35 | uint32 cols) |
36 | |
37 | : fRows (0) |
38 | , fCols (0) |
39 | |
40 | { |
41 | |
42 | if (rows < 1 || rows > kMaxColorPlanes || |
43 | cols < 1 || cols > kMaxColorPlanes) |
44 | { |
45 | |
46 | ThrowProgramError (); |
47 | |
48 | } |
49 | |
50 | fRows = rows; |
51 | fCols = cols; |
52 | |
53 | for (uint32 row = 0; row < fRows; row++) |
54 | for (uint32 col = 0; col < fCols; col++) |
55 | { |
56 | |
57 | fData [row] [col] = 0.0; |
58 | |
59 | } |
60 | |
61 | } |
62 | |
63 | /*****************************************************************************/ |
64 | |
65 | dng_matrix::dng_matrix (const dng_matrix &m) |
66 | |
67 | : fRows (m.fRows) |
68 | , fCols (m.fCols) |
69 | |
70 | { |
71 | |
72 | for (uint32 row = 0; row < fRows; row++) |
73 | for (uint32 col = 0; col < fCols; col++) |
74 | { |
75 | |
76 | fData [row] [col] = m.fData [row] [col]; |
77 | |
78 | } |
79 | |
80 | } |
81 | |
82 | /*****************************************************************************/ |
83 | |
84 | void dng_matrix::Clear () |
85 | { |
86 | |
87 | fRows = 0; |
88 | fCols = 0; |
89 | |
90 | } |
91 | |
92 | /*****************************************************************************/ |
93 | |
94 | void dng_matrix::SetIdentity (uint32 count) |
95 | { |
96 | |
97 | *this = dng_matrix (count, count); |
98 | |
99 | for (uint32 j = 0; j < count; j++) |
100 | { |
101 | |
102 | fData [j] [j] = 1.0; |
103 | |
104 | } |
105 | |
106 | } |
107 | |
108 | /******************************************************************************/ |
109 | |
110 | bool dng_matrix::operator== (const dng_matrix &m) const |
111 | { |
112 | |
113 | if (Rows () != m.Rows () || |
114 | Cols () != m.Cols ()) |
115 | { |
116 | |
117 | return false; |
118 | |
119 | } |
120 | |
121 | for (uint32 j = 0; j < Rows (); j++) |
122 | for (uint32 k = 0; k < Cols (); k++) |
123 | { |
124 | |
125 | if (fData [j] [k] != m.fData [j] [k]) |
126 | { |
127 | |
128 | return false; |
129 | |
130 | } |
131 | |
132 | } |
133 | |
134 | return true; |
135 | |
136 | } |
137 | |
138 | /******************************************************************************/ |
139 | |
140 | bool dng_matrix::IsDiagonal () const |
141 | { |
142 | |
143 | if (IsEmpty ()) |
144 | { |
145 | return false; |
146 | } |
147 | |
148 | if (Rows () != Cols ()) |
149 | { |
150 | return false; |
151 | } |
152 | |
153 | for (uint32 j = 0; j < Rows (); j++) |
154 | for (uint32 k = 0; k < Cols (); k++) |
155 | { |
156 | |
157 | if (j != k) |
158 | { |
159 | |
160 | if (fData [j] [k] != 0.0) |
161 | { |
162 | return false; |
163 | } |
164 | |
165 | } |
166 | |
167 | } |
168 | |
169 | return true; |
170 | |
171 | } |
172 | |
173 | /******************************************************************************/ |
174 | |
175 | real64 dng_matrix::MaxEntry () const |
176 | { |
177 | |
178 | if (IsEmpty ()) |
179 | { |
180 | |
181 | return 0.0; |
182 | |
183 | } |
184 | |
185 | real64 m = fData [0] [0]; |
186 | |
187 | for (uint32 j = 0; j < Rows (); j++) |
188 | for (uint32 k = 0; k < Cols (); k++) |
189 | { |
190 | |
191 | m = Max_real64 (m, fData [j] [k]); |
192 | |
193 | } |
194 | |
195 | return m; |
196 | |
197 | } |
198 | |
199 | /******************************************************************************/ |
200 | |
201 | real64 dng_matrix::MinEntry () const |
202 | { |
203 | |
204 | if (IsEmpty ()) |
205 | { |
206 | |
207 | return 0.0; |
208 | |
209 | } |
210 | |
211 | real64 m = fData [0] [0]; |
212 | |
213 | for (uint32 j = 0; j < Rows (); j++) |
214 | for (uint32 k = 0; k < Cols (); k++) |
215 | { |
216 | |
217 | m = Min_real64 (m, fData [j] [k]); |
218 | |
219 | } |
220 | |
221 | return m; |
222 | |
223 | } |
224 | |
225 | /*****************************************************************************/ |
226 | |
227 | void dng_matrix::Scale (real64 factor) |
228 | { |
229 | |
230 | for (uint32 j = 0; j < Rows (); j++) |
231 | for (uint32 k = 0; k < Cols (); k++) |
232 | { |
233 | |
234 | fData [j] [k] *= factor; |
235 | |
236 | } |
237 | |
238 | } |
239 | |
240 | /*****************************************************************************/ |
241 | |
242 | void dng_matrix::Round (real64 factor) |
243 | { |
244 | |
245 | real64 invFactor = 1.0 / factor; |
246 | |
247 | for (uint32 j = 0; j < Rows (); j++) |
248 | for (uint32 k = 0; k < Cols (); k++) |
249 | { |
250 | |
251 | fData [j] [k] = Round_int32 (fData [j] [k] * factor) * invFactor; |
252 | |
253 | } |
254 | |
255 | } |
256 | |
257 | /*****************************************************************************/ |
258 | |
259 | void dng_matrix::SafeRound (real64 factor) |
260 | { |
261 | |
262 | real64 invFactor = 1.0 / factor; |
263 | |
264 | for (uint32 j = 0; j < Rows (); j++) |
265 | { |
266 | |
267 | // Round each row to the specified accuracy, but make sure the |
268 | // a rounding does not affect the total of the elements in a row |
269 | // more than necessary. |
270 | |
271 | real64 error = 0.0; |
272 | |
273 | for (uint32 k = 0; k < Cols (); k++) |
274 | { |
275 | |
276 | fData [j] [k] += error; |
277 | |
278 | real64 rounded = Round_int32 (fData [j] [k] * factor) * invFactor; |
279 | |
280 | error = fData [j] [k] - rounded; |
281 | |
282 | fData [j] [k] = rounded; |
283 | |
284 | } |
285 | |
286 | } |
287 | |
288 | } |
289 | |
290 | /*****************************************************************************/ |
291 | |
292 | dng_matrix_3by3::dng_matrix_3by3 () |
293 | |
294 | : dng_matrix (3, 3) |
295 | |
296 | { |
297 | } |
298 | |
299 | /*****************************************************************************/ |
300 | |
301 | dng_matrix_3by3::dng_matrix_3by3 (const dng_matrix &m) |
302 | |
303 | : dng_matrix (m) |
304 | |
305 | { |
306 | |
307 | if (Rows () != 3 || |
308 | Cols () != 3) |
309 | { |
310 | |
311 | ThrowMatrixMath (); |
312 | |
313 | } |
314 | |
315 | } |
316 | |
317 | /*****************************************************************************/ |
318 | |
319 | dng_matrix_3by3::dng_matrix_3by3 (real64 a00, real64 a01, real64 a02, |
320 | real64 a10, real64 a11, real64 a12, |
321 | real64 a20, real64 a21, real64 a22) |
322 | |
323 | |
324 | : dng_matrix (3, 3) |
325 | |
326 | { |
327 | |
328 | fData [0] [0] = a00; |
329 | fData [0] [1] = a01; |
330 | fData [0] [2] = a02; |
331 | |
332 | fData [1] [0] = a10; |
333 | fData [1] [1] = a11; |
334 | fData [1] [2] = a12; |
335 | |
336 | fData [2] [0] = a20; |
337 | fData [2] [1] = a21; |
338 | fData [2] [2] = a22; |
339 | |
340 | } |
341 | |
342 | /*****************************************************************************/ |
343 | |
344 | dng_matrix_3by3::dng_matrix_3by3 (real64 a00, real64 a11, real64 a22) |
345 | |
346 | : dng_matrix (3, 3) |
347 | |
348 | { |
349 | |
350 | fData [0] [0] = a00; |
351 | fData [1] [1] = a11; |
352 | fData [2] [2] = a22; |
353 | |
354 | } |
355 | |
356 | /*****************************************************************************/ |
357 | |
358 | dng_matrix_4by3::dng_matrix_4by3 () |
359 | |
360 | : dng_matrix (4, 3) |
361 | |
362 | { |
363 | } |
364 | |
365 | /*****************************************************************************/ |
366 | |
367 | dng_matrix_4by3::dng_matrix_4by3 (const dng_matrix &m) |
368 | |
369 | : dng_matrix (m) |
370 | |
371 | { |
372 | |
373 | if (Rows () != 4 || |
374 | Cols () != 3) |
375 | { |
376 | |
377 | ThrowMatrixMath (); |
378 | |
379 | } |
380 | |
381 | } |
382 | |
383 | /*****************************************************************************/ |
384 | |
385 | dng_matrix_4by3::dng_matrix_4by3 (real64 a00, real64 a01, real64 a02, |
386 | real64 a10, real64 a11, real64 a12, |
387 | real64 a20, real64 a21, real64 a22, |
388 | real64 a30, real64 a31, real64 a32) |
389 | |
390 | |
391 | : dng_matrix (4, 3) |
392 | |
393 | { |
394 | |
395 | fData [0] [0] = a00; |
396 | fData [0] [1] = a01; |
397 | fData [0] [2] = a02; |
398 | |
399 | fData [1] [0] = a10; |
400 | fData [1] [1] = a11; |
401 | fData [1] [2] = a12; |
402 | |
403 | fData [2] [0] = a20; |
404 | fData [2] [1] = a21; |
405 | fData [2] [2] = a22; |
406 | |
407 | fData [3] [0] = a30; |
408 | fData [3] [1] = a31; |
409 | fData [3] [2] = a32; |
410 | |
411 | } |
412 | |
413 | /*****************************************************************************/ |
414 | |
415 | dng_vector::dng_vector () |
416 | |
417 | : fCount (0) |
418 | |
419 | { |
420 | |
421 | } |
422 | |
423 | /*****************************************************************************/ |
424 | |
425 | dng_vector::dng_vector (uint32 count) |
426 | |
427 | : fCount (0) |
428 | |
429 | { |
430 | |
431 | if (count < 1 || count > kMaxColorPlanes) |
432 | { |
433 | |
434 | ThrowProgramError (); |
435 | |
436 | } |
437 | |
438 | fCount = count; |
439 | |
440 | for (uint32 index = 0; index < fCount; index++) |
441 | { |
442 | |
443 | fData [index] = 0.0; |
444 | |
445 | } |
446 | |
447 | } |
448 | |
449 | /*****************************************************************************/ |
450 | |
451 | dng_vector::dng_vector (const dng_vector &v) |
452 | |
453 | : fCount (v.fCount) |
454 | |
455 | { |
456 | |
457 | for (uint32 index = 0; index < fCount; index++) |
458 | { |
459 | |
460 | fData [index] = v.fData [index]; |
461 | |
462 | } |
463 | |
464 | } |
465 | |
466 | /*****************************************************************************/ |
467 | |
468 | void dng_vector::Clear () |
469 | { |
470 | |
471 | fCount = 0; |
472 | |
473 | } |
474 | |
475 | /*****************************************************************************/ |
476 | |
477 | void dng_vector::SetIdentity (uint32 count) |
478 | { |
479 | |
480 | *this = dng_vector (count); |
481 | |
482 | for (uint32 j = 0; j < count; j++) |
483 | { |
484 | |
485 | fData [j] = 1.0; |
486 | |
487 | } |
488 | |
489 | } |
490 | |
491 | /******************************************************************************/ |
492 | |
493 | bool dng_vector::operator== (const dng_vector &v) const |
494 | { |
495 | |
496 | if (Count () != v.Count ()) |
497 | { |
498 | |
499 | return false; |
500 | |
501 | } |
502 | |
503 | for (uint32 j = 0; j < Count (); j++) |
504 | { |
505 | |
506 | if (fData [j] != v.fData [j]) |
507 | { |
508 | |
509 | return false; |
510 | |
511 | } |
512 | |
513 | } |
514 | |
515 | return true; |
516 | |
517 | } |
518 | |
519 | /******************************************************************************/ |
520 | |
521 | real64 dng_vector::MaxEntry () const |
522 | { |
523 | |
524 | if (IsEmpty ()) |
525 | { |
526 | |
527 | return 0.0; |
528 | |
529 | } |
530 | |
531 | real64 m = fData [0]; |
532 | |
533 | for (uint32 j = 0; j < Count (); j++) |
534 | { |
535 | |
536 | m = Max_real64 (m, fData [j]); |
537 | |
538 | } |
539 | |
540 | return m; |
541 | |
542 | } |
543 | |
544 | /******************************************************************************/ |
545 | |
546 | real64 dng_vector::MinEntry () const |
547 | { |
548 | |
549 | if (IsEmpty ()) |
550 | { |
551 | |
552 | return 0.0; |
553 | |
554 | } |
555 | |
556 | real64 m = fData [0]; |
557 | |
558 | for (uint32 j = 0; j < Count (); j++) |
559 | { |
560 | |
561 | m = Min_real64 (m, fData [j]); |
562 | |
563 | } |
564 | |
565 | return m; |
566 | |
567 | } |
568 | |
569 | /*****************************************************************************/ |
570 | |
571 | void dng_vector::Scale (real64 factor) |
572 | { |
573 | |
574 | for (uint32 j = 0; j < Count (); j++) |
575 | { |
576 | |
577 | fData [j] *= factor; |
578 | |
579 | } |
580 | |
581 | } |
582 | |
583 | /*****************************************************************************/ |
584 | |
585 | void dng_vector::Round (real64 factor) |
586 | { |
587 | |
588 | real64 invFactor = 1.0 / factor; |
589 | |
590 | for (uint32 j = 0; j < Count (); j++) |
591 | { |
592 | |
593 | fData [j] = Round_int32 (fData [j] * factor) * invFactor; |
594 | |
595 | } |
596 | |
597 | } |
598 | |
599 | /*****************************************************************************/ |
600 | |
601 | dng_matrix dng_vector::AsDiagonal () const |
602 | { |
603 | |
604 | dng_matrix M (Count (), Count ()); |
605 | |
606 | for (uint32 j = 0; j < Count (); j++) |
607 | { |
608 | |
609 | M [j] [j] = fData [j]; |
610 | |
611 | } |
612 | |
613 | return M; |
614 | |
615 | } |
616 | |
617 | /*****************************************************************************/ |
618 | |
619 | dng_matrix dng_vector::AsColumn () const |
620 | { |
621 | |
622 | dng_matrix M (Count (), 1); |
623 | |
624 | for (uint32 j = 0; j < Count (); j++) |
625 | { |
626 | |
627 | M [j] [0] = fData [j]; |
628 | |
629 | } |
630 | |
631 | return M; |
632 | |
633 | } |
634 | |
635 | /******************************************************************************/ |
636 | |
637 | dng_vector_3::dng_vector_3 () |
638 | |
639 | : dng_vector (3) |
640 | |
641 | { |
642 | } |
643 | |
644 | /******************************************************************************/ |
645 | |
646 | dng_vector_3::dng_vector_3 (const dng_vector &v) |
647 | |
648 | : dng_vector (v) |
649 | |
650 | { |
651 | |
652 | if (Count () != 3) |
653 | { |
654 | |
655 | ThrowMatrixMath (); |
656 | |
657 | } |
658 | |
659 | } |
660 | |
661 | /******************************************************************************/ |
662 | |
663 | dng_vector_3::dng_vector_3 (real64 a0, |
664 | real64 a1, |
665 | real64 a2) |
666 | |
667 | : dng_vector (3) |
668 | |
669 | { |
670 | |
671 | fData [0] = a0; |
672 | fData [1] = a1; |
673 | fData [2] = a2; |
674 | |
675 | } |
676 | |
677 | /******************************************************************************/ |
678 | |
679 | dng_vector_4::dng_vector_4 () |
680 | |
681 | : dng_vector (4) |
682 | |
683 | { |
684 | } |
685 | |
686 | /******************************************************************************/ |
687 | |
688 | dng_vector_4::dng_vector_4 (const dng_vector &v) |
689 | |
690 | : dng_vector (v) |
691 | |
692 | { |
693 | |
694 | if (Count () != 4) |
695 | { |
696 | |
697 | ThrowMatrixMath (); |
698 | |
699 | } |
700 | |
701 | } |
702 | |
703 | /******************************************************************************/ |
704 | |
705 | dng_vector_4::dng_vector_4 (real64 a0, |
706 | real64 a1, |
707 | real64 a2, |
708 | real64 a3) |
709 | |
710 | : dng_vector (4) |
711 | |
712 | { |
713 | |
714 | fData [0] = a0; |
715 | fData [1] = a1; |
716 | fData [2] = a2; |
717 | fData [3] = a3; |
718 | |
719 | } |
720 | |
721 | /******************************************************************************/ |
722 | |
723 | dng_matrix operator* (const dng_matrix &A, |
724 | const dng_matrix &B) |
725 | { |
726 | |
727 | if (A.Cols () != B.Rows ()) |
728 | { |
729 | |
730 | ThrowMatrixMath (); |
731 | |
732 | } |
733 | |
734 | dng_matrix C (A.Rows (), B.Cols ()); |
735 | |
736 | for (uint32 j = 0; j < C.Rows (); j++) |
737 | for (uint32 k = 0; k < C.Cols (); k++) |
738 | { |
739 | |
740 | C [j] [k] = 0.0; |
741 | |
742 | for (uint32 m = 0; m < A.Cols (); m++) |
743 | { |
744 | |
745 | real64 aa = A [j] [m]; |
746 | |
747 | real64 bb = B [m] [k]; |
748 | |
749 | C [j] [k] += aa * bb; |
750 | |
751 | } |
752 | |
753 | } |
754 | |
755 | return C; |
756 | |
757 | } |
758 | |
759 | /******************************************************************************/ |
760 | |
761 | dng_vector operator* (const dng_matrix &A, |
762 | const dng_vector &B) |
763 | { |
764 | |
765 | if (A.Cols () != B.Count ()) |
766 | { |
767 | |
768 | ThrowMatrixMath (); |
769 | |
770 | } |
771 | |
772 | dng_vector C (A.Rows ()); |
773 | |
774 | for (uint32 j = 0; j < C.Count (); j++) |
775 | { |
776 | |
777 | C [j] = 0.0; |
778 | |
779 | for (uint32 m = 0; m < A.Cols (); m++) |
780 | { |
781 | |
782 | real64 aa = A [j] [m]; |
783 | |
784 | real64 bb = B [m]; |
785 | |
786 | C [j] += aa * bb; |
787 | |
788 | } |
789 | |
790 | } |
791 | |
792 | return C; |
793 | |
794 | } |
795 | |
796 | /******************************************************************************/ |
797 | |
798 | dng_matrix operator* (real64 scale, |
799 | const dng_matrix &A) |
800 | { |
801 | |
802 | dng_matrix B (A); |
803 | |
804 | B.Scale (scale); |
805 | |
806 | return B; |
807 | |
808 | } |
809 | |
810 | /******************************************************************************/ |
811 | |
812 | dng_vector operator* (real64 scale, |
813 | const dng_vector &A) |
814 | { |
815 | |
816 | dng_vector B (A); |
817 | |
818 | B.Scale (scale); |
819 | |
820 | return B; |
821 | |
822 | } |
823 | |
824 | /******************************************************************************/ |
825 | |
826 | dng_matrix operator+ (const dng_matrix &A, |
827 | const dng_matrix &B) |
828 | { |
829 | |
830 | if (A.Cols () != B.Cols () || |
831 | A.Rows () != B.Rows ()) |
832 | { |
833 | |
834 | ThrowMatrixMath (); |
835 | |
836 | } |
837 | |
838 | dng_matrix C (A); |
839 | |
840 | for (uint32 j = 0; j < C.Rows (); j++) |
841 | for (uint32 k = 0; k < C.Cols (); k++) |
842 | { |
843 | |
844 | C [j] [k] += B [j] [k]; |
845 | |
846 | } |
847 | |
848 | return C; |
849 | |
850 | } |
851 | |
852 | /******************************************************************************/ |
853 | |
854 | const real64 kNearZero = 1.0E-10; |
855 | |
856 | /******************************************************************************/ |
857 | |
858 | // Work around bug #1294195, which may be a hardware problem on a specific machine. |
859 | // This pragma turns on "improved" floating-point consistency. |
860 | #ifdef _MSC_VER |
861 | #pragma optimize ("p", on) |
862 | #endif |
863 | |
864 | static dng_matrix Invert3by3 (const dng_matrix &A) |
865 | { |
866 | |
867 | real64 a00 = A [0] [0]; |
868 | real64 a01 = A [0] [1]; |
869 | real64 a02 = A [0] [2]; |
870 | real64 a10 = A [1] [0]; |
871 | real64 a11 = A [1] [1]; |
872 | real64 a12 = A [1] [2]; |
873 | real64 a20 = A [2] [0]; |
874 | real64 a21 = A [2] [1]; |
875 | real64 a22 = A [2] [2]; |
876 | |
877 | real64 temp [3] [3]; |
878 | |
879 | temp [0] [0] = a11 * a22 - a21 * a12; |
880 | temp [0] [1] = a21 * a02 - a01 * a22; |
881 | temp [0] [2] = a01 * a12 - a11 * a02; |
882 | temp [1] [0] = a20 * a12 - a10 * a22; |
883 | temp [1] [1] = a00 * a22 - a20 * a02; |
884 | temp [1] [2] = a10 * a02 - a00 * a12; |
885 | temp [2] [0] = a10 * a21 - a20 * a11; |
886 | temp [2] [1] = a20 * a01 - a00 * a21; |
887 | temp [2] [2] = a00 * a11 - a10 * a01; |
888 | |
889 | real64 det = (a00 * temp [0] [0] + |
890 | a01 * temp [1] [0] + |
891 | a02 * temp [2] [0]); |
892 | |
893 | if (Abs_real64 (det) < kNearZero) |
894 | { |
895 | |
896 | ThrowMatrixMath (); |
897 | |
898 | } |
899 | |
900 | dng_matrix B (3, 3); |
901 | |
902 | for (uint32 j = 0; j < 3; j++) |
903 | for (uint32 k = 0; k < 3; k++) |
904 | { |
905 | |
906 | B [j] [k] = temp [j] [k] / det; |
907 | |
908 | } |
909 | |
910 | return B; |
911 | |
912 | } |
913 | |
914 | // Reset floating-point optimization. See comment above. |
915 | #ifdef _MSC_VER |
916 | #pragma optimize ("p", off) |
917 | #endif |
918 | |
919 | /******************************************************************************/ |
920 | |
921 | static dng_matrix InvertNbyN (const dng_matrix &A) |
922 | { |
923 | |
924 | uint32 i; |
925 | uint32 j; |
926 | uint32 k; |
927 | |
928 | uint32 n = A.Rows (); |
929 | |
930 | real64 temp [kMaxColorPlanes] [kMaxColorPlanes * 2]; |
931 | |
932 | for (i = 0; i < n; i++) |
933 | for (j = 0; j < n; j++) |
934 | { |
935 | |
936 | temp [i] [j ] = A [i] [j]; |
937 | |
938 | temp [i] [j + n] = (i == j ? 1.0 : 0.0); |
939 | |
940 | } |
941 | |
942 | for (i = 0; i < n; i++) |
943 | { |
944 | |
945 | real64 alpha = temp [i] [i]; |
946 | |
947 | if (Abs_real64 (alpha) < kNearZero) |
948 | { |
949 | |
950 | ThrowMatrixMath (); |
951 | |
952 | } |
953 | |
954 | for (j = 0; j < n * 2; j++) |
955 | { |
956 | |
957 | temp [i] [j] /= alpha; |
958 | |
959 | } |
960 | |
961 | for (k = 0; k < n; k++) |
962 | { |
963 | |
964 | if (i != k) |
965 | { |
966 | |
967 | real64 beta = temp [k] [i]; |
968 | |
969 | for (j = 0; j < n * 2; j++) |
970 | { |
971 | |
972 | temp [k] [j] -= beta * temp [i] [j]; |
973 | |
974 | } |
975 | |
976 | } |
977 | |
978 | } |
979 | |
980 | } |
981 | |
982 | dng_matrix B (n, n); |
983 | |
984 | for (i = 0; i < n; i++) |
985 | for (j = 0; j < n; j++) |
986 | { |
987 | |
988 | B [i] [j] = temp [i] [j + n]; |
989 | |
990 | } |
991 | |
992 | return B; |
993 | |
994 | } |
995 | |
996 | /******************************************************************************/ |
997 | |
998 | dng_matrix Transpose (const dng_matrix &A) |
999 | { |
1000 | |
1001 | dng_matrix B (A.Cols (), A.Rows ()); |
1002 | |
1003 | for (uint32 j = 0; j < B.Rows (); j++) |
1004 | for (uint32 k = 0; k < B.Cols (); k++) |
1005 | { |
1006 | |
1007 | B [j] [k] = A [k] [j]; |
1008 | |
1009 | } |
1010 | |
1011 | return B; |
1012 | |
1013 | } |
1014 | |
1015 | /******************************************************************************/ |
1016 | |
1017 | dng_matrix Invert (const dng_matrix &A) |
1018 | { |
1019 | |
1020 | if (A.Rows () < 2 || A.Cols () < 2) |
1021 | { |
1022 | |
1023 | ThrowMatrixMath (); |
1024 | |
1025 | } |
1026 | |
1027 | if (A.Rows () == A.Cols ()) |
1028 | { |
1029 | |
1030 | if (A.Rows () == 3) |
1031 | { |
1032 | |
1033 | return Invert3by3 (A); |
1034 | |
1035 | } |
1036 | |
1037 | return InvertNbyN (A); |
1038 | |
1039 | } |
1040 | |
1041 | else |
1042 | { |
1043 | |
1044 | // Compute the pseudo inverse. |
1045 | |
1046 | dng_matrix B = Transpose (A); |
1047 | |
1048 | return Invert (B * A) * B; |
1049 | |
1050 | } |
1051 | |
1052 | } |
1053 | |
1054 | /******************************************************************************/ |
1055 | |
1056 | dng_matrix Invert (const dng_matrix &A, |
1057 | const dng_matrix &hint) |
1058 | { |
1059 | |
1060 | if (A.Rows () == A .Cols () || |
1061 | A.Rows () != hint.Cols () || |
1062 | A.Cols () != hint.Rows ()) |
1063 | { |
1064 | |
1065 | return Invert (A); |
1066 | |
1067 | } |
1068 | |
1069 | else |
1070 | { |
1071 | |
1072 | // Use the specified hint matrix. |
1073 | |
1074 | return Invert (hint * A) * hint; |
1075 | |
1076 | } |
1077 | |
1078 | } |
1079 | |
1080 | /*****************************************************************************/ |
1081 | |