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
23dng_matrix::dng_matrix ()
24
25 : fRows (0)
26 , fCols (0)
27
28 {
29
30 }
31
32/*****************************************************************************/
33
34dng_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
65dng_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
84void dng_matrix::Clear ()
85 {
86
87 fRows = 0;
88 fCols = 0;
89
90 }
91
92/*****************************************************************************/
93
94void 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
110bool 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
140bool 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
175real64 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
201real64 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
227void 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
242void 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
259void 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
292dng_matrix_3by3::dng_matrix_3by3 ()
293
294 : dng_matrix (3, 3)
295
296 {
297 }
298
299/*****************************************************************************/
300
301dng_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
319dng_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
344dng_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
358dng_matrix_4by3::dng_matrix_4by3 ()
359
360 : dng_matrix (4, 3)
361
362 {
363 }
364
365/*****************************************************************************/
366
367dng_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
385dng_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
415dng_vector::dng_vector ()
416
417 : fCount (0)
418
419 {
420
421 }
422
423/*****************************************************************************/
424
425dng_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
451dng_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
468void dng_vector::Clear ()
469 {
470
471 fCount = 0;
472
473 }
474
475/*****************************************************************************/
476
477void 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
493bool 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
521real64 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
546real64 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
571void 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
585void 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
601dng_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
619dng_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
637dng_vector_3::dng_vector_3 ()
638
639 : dng_vector (3)
640
641 {
642 }
643
644/******************************************************************************/
645
646dng_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
663dng_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
679dng_vector_4::dng_vector_4 ()
680
681 : dng_vector (4)
682
683 {
684 }
685
686/******************************************************************************/
687
688dng_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
705dng_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
723dng_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
761dng_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
798dng_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
812dng_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
826dng_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
854const 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
864static 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
921static 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
998dng_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
1017dng_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
1056dng_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