1 /*****************************************************************************/
2 // Copyright 2006-2019 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 #include "dng_matrix.h"
10
11 #include "dng_assertions.h"
12 #include "dng_exceptions.h"
13 #include "dng_utils.h"
14
15 /*****************************************************************************/
16
dng_matrix()17 dng_matrix::dng_matrix ()
18
19 : fRows (0)
20 , fCols (0)
21
22 {
23
24 }
25
26 /*****************************************************************************/
27
dng_matrix(uint32 rows,uint32 cols)28 dng_matrix::dng_matrix (uint32 rows,
29 uint32 cols)
30
31 : fRows (0)
32 , fCols (0)
33
34 {
35
36 if (rows < 1 || rows > kMaxColorPlanes ||
37 cols < 1 || cols > kMaxColorPlanes)
38 {
39
40 ThrowProgramError ();
41
42 }
43
44 fRows = rows;
45 fCols = cols;
46
47 for (uint32 row = 0; row < fRows; row++)
48 for (uint32 col = 0; col < fCols; col++)
49 {
50
51 fData [row] [col] = 0.0;
52
53 }
54
55 }
56
57 /*****************************************************************************/
58
dng_matrix(const dng_matrix & m)59 dng_matrix::dng_matrix (const dng_matrix &m)
60
61 : fRows (m.fRows)
62 , fCols (m.fCols)
63
64 {
65
66 for (uint32 row = 0; row < fRows; row++)
67 for (uint32 col = 0; col < fCols; col++)
68 {
69
70 fData [row] [col] = m.fData [row] [col];
71
72 }
73
74 }
75
76 /*****************************************************************************/
77
Clear()78 void dng_matrix::Clear ()
79 {
80
81 fRows = 0;
82 fCols = 0;
83
84 }
85
86 /*****************************************************************************/
87
SetIdentity(uint32 count)88 void dng_matrix::SetIdentity (uint32 count)
89 {
90
91 *this = dng_matrix (count, count);
92
93 for (uint32 j = 0; j < count; j++)
94 {
95
96 fData [j] [j] = 1.0;
97
98 }
99
100 }
101
102 /******************************************************************************/
103
operator ==(const dng_matrix & m) const104 bool dng_matrix::operator== (const dng_matrix &m) const
105 {
106
107 if (Rows () != m.Rows () ||
108 Cols () != m.Cols ())
109 {
110
111 return false;
112
113 }
114
115 for (uint32 j = 0; j < Rows (); j++)
116 for (uint32 k = 0; k < Cols (); k++)
117 {
118
119 if (fData [j] [k] != m.fData [j] [k])
120 {
121
122 return false;
123
124 }
125
126 }
127
128 return true;
129
130 }
131
132 /******************************************************************************/
133
IsDiagonal() const134 bool dng_matrix::IsDiagonal () const
135 {
136
137 if (IsEmpty ())
138 {
139 return false;
140 }
141
142 if (Rows () != Cols ())
143 {
144 return false;
145 }
146
147 for (uint32 j = 0; j < Rows (); j++)
148 for (uint32 k = 0; k < Cols (); k++)
149 {
150
151 if (j != k)
152 {
153
154 if (fData [j] [k] != 0.0)
155 {
156 return false;
157 }
158
159 }
160
161 }
162
163 return true;
164
165 }
166
167 /******************************************************************************/
168
IsIdentity() const169 bool dng_matrix::IsIdentity () const
170 {
171
172 if (IsDiagonal ())
173 {
174
175 for (uint32 j = 0; j < Rows (); j++)
176 {
177
178 if (fData [j] [j] != 1.0)
179 {
180 return false;
181 }
182
183 }
184
185 return true;
186
187 }
188
189 return false;
190
191 }
192
193 /******************************************************************************/
194
MaxEntry() const195 real64 dng_matrix::MaxEntry () const
196 {
197
198 if (IsEmpty ())
199 {
200
201 return 0.0;
202
203 }
204
205 real64 m = fData [0] [0];
206
207 for (uint32 j = 0; j < Rows (); j++)
208 for (uint32 k = 0; k < Cols (); k++)
209 {
210
211 m = Max_real64 (m, fData [j] [k]);
212
213 }
214
215 return m;
216
217 }
218
219 /******************************************************************************/
220
MinEntry() const221 real64 dng_matrix::MinEntry () const
222 {
223
224 if (IsEmpty ())
225 {
226
227 return 0.0;
228
229 }
230
231 real64 m = fData [0] [0];
232
233 for (uint32 j = 0; j < Rows (); j++)
234 for (uint32 k = 0; k < Cols (); k++)
235 {
236
237 m = Min_real64 (m, fData [j] [k]);
238
239 }
240
241 return m;
242
243 }
244
245 /*****************************************************************************/
246
Scale(real64 factor)247 void dng_matrix::Scale (real64 factor)
248 {
249
250 for (uint32 j = 0; j < Rows (); j++)
251 for (uint32 k = 0; k < Cols (); k++)
252 {
253
254 fData [j] [k] *= factor;
255
256 }
257
258 }
259
260 /*****************************************************************************/
261
Round(real64 factor)262 void dng_matrix::Round (real64 factor)
263 {
264
265 real64 invFactor = 1.0 / factor;
266
267 for (uint32 j = 0; j < Rows (); j++)
268 for (uint32 k = 0; k < Cols (); k++)
269 {
270
271 fData [j] [k] = Round_int32 (fData [j] [k] * factor) * invFactor;
272
273 }
274
275 }
276
277 /*****************************************************************************/
278
SafeRound(real64 factor)279 void dng_matrix::SafeRound (real64 factor)
280 {
281
282 real64 invFactor = 1.0 / factor;
283
284 for (uint32 j = 0; j < Rows (); j++)
285 {
286
287 // Round each row to the specified accuracy, but make sure the
288 // a rounding does not affect the total of the elements in a row
289 // more than necessary.
290
291 real64 error = 0.0;
292
293 for (uint32 k = 0; k < Cols (); k++)
294 {
295
296 fData [j] [k] += error;
297
298 real64 rounded = Round_int32 (fData [j] [k] * factor) * invFactor;
299
300 error = fData [j] [k] - rounded;
301
302 fData [j] [k] = rounded;
303
304 }
305
306 }
307
308 }
309
310 /*****************************************************************************/
311
AlmostEqual(const dng_matrix & m,real64 slop) const312 bool dng_matrix::AlmostEqual (const dng_matrix &m,
313 real64 slop) const
314 {
315
316 if (Rows () != m.Rows () ||
317 Cols () != m.Cols ())
318 {
319 return false;
320 }
321
322 for (uint32 j = 0; j < Rows (); j++)
323 {
324
325 for (uint32 k = 0; k < Cols (); k++)
326 {
327
328 if (Abs_real64 (fData [j] [k] - m [j] [k]) > slop)
329 {
330 return false;
331 }
332
333 }
334
335 }
336
337 return true;
338
339 }
340
341 /*****************************************************************************/
342
AlmostIdentity(real64 slop) const343 bool dng_matrix::AlmostIdentity (real64 slop) const
344 {
345
346 dng_matrix m;
347
348 m.SetIdentity (Rows ());
349
350 return AlmostEqual (m, slop);
351
352 }
353
354 /*****************************************************************************/
355
dng_matrix_3by3()356 dng_matrix_3by3::dng_matrix_3by3 ()
357
358 : dng_matrix (3, 3)
359
360 {
361 }
362
363 /*****************************************************************************/
364
dng_matrix_3by3(const dng_matrix & m)365 dng_matrix_3by3::dng_matrix_3by3 (const dng_matrix &m)
366
367 : dng_matrix (m)
368
369 {
370
371 if (Rows () != 3 ||
372 Cols () != 3)
373 {
374
375 ThrowMatrixMath ();
376
377 }
378
379 }
380
381 /*****************************************************************************/
382
dng_matrix_3by3(real64 a00,real64 a01,real64 a02,real64 a10,real64 a11,real64 a12,real64 a20,real64 a21,real64 a22)383 dng_matrix_3by3::dng_matrix_3by3 (real64 a00, real64 a01, real64 a02,
384 real64 a10, real64 a11, real64 a12,
385 real64 a20, real64 a21, real64 a22)
386
387
388 : dng_matrix (3, 3)
389
390 {
391
392 fData [0] [0] = a00;
393 fData [0] [1] = a01;
394 fData [0] [2] = a02;
395
396 fData [1] [0] = a10;
397 fData [1] [1] = a11;
398 fData [1] [2] = a12;
399
400 fData [2] [0] = a20;
401 fData [2] [1] = a21;
402 fData [2] [2] = a22;
403
404 }
405
406 /*****************************************************************************/
407
dng_matrix_3by3(real64 a00,real64 a11,real64 a22)408 dng_matrix_3by3::dng_matrix_3by3 (real64 a00, real64 a11, real64 a22)
409
410 : dng_matrix (3, 3)
411
412 {
413
414 fData [0] [0] = a00;
415 fData [1] [1] = a11;
416 fData [2] [2] = a22;
417
418 }
419
420 /*****************************************************************************/
421
dng_matrix_4by3()422 dng_matrix_4by3::dng_matrix_4by3 ()
423
424 : dng_matrix (4, 3)
425
426 {
427 }
428
429 /*****************************************************************************/
430
dng_matrix_4by3(const dng_matrix & m)431 dng_matrix_4by3::dng_matrix_4by3 (const dng_matrix &m)
432
433 : dng_matrix (m)
434
435 {
436
437 if (Rows () != 4 ||
438 Cols () != 3)
439 {
440
441 ThrowMatrixMath ();
442
443 }
444
445 }
446
447 /*****************************************************************************/
448
dng_matrix_4by3(real64 a00,real64 a01,real64 a02,real64 a10,real64 a11,real64 a12,real64 a20,real64 a21,real64 a22,real64 a30,real64 a31,real64 a32)449 dng_matrix_4by3::dng_matrix_4by3 (real64 a00, real64 a01, real64 a02,
450 real64 a10, real64 a11, real64 a12,
451 real64 a20, real64 a21, real64 a22,
452 real64 a30, real64 a31, real64 a32)
453
454
455 : dng_matrix (4, 3)
456
457 {
458
459 fData [0] [0] = a00;
460 fData [0] [1] = a01;
461 fData [0] [2] = a02;
462
463 fData [1] [0] = a10;
464 fData [1] [1] = a11;
465 fData [1] [2] = a12;
466
467 fData [2] [0] = a20;
468 fData [2] [1] = a21;
469 fData [2] [2] = a22;
470
471 fData [3] [0] = a30;
472 fData [3] [1] = a31;
473 fData [3] [2] = a32;
474
475 }
476
477 /*****************************************************************************/
478
dng_matrix_4by4()479 dng_matrix_4by4::dng_matrix_4by4 ()
480
481 : dng_matrix (4, 4)
482
483 {
484
485 }
486
487 /*****************************************************************************/
488
dng_matrix_4by4(const dng_matrix & m)489 dng_matrix_4by4::dng_matrix_4by4 (const dng_matrix &m)
490
491 : dng_matrix (m)
492
493 {
494
495 // Input must be either 3x3 or 4x4.
496
497 const bool is3by3 = (m.Rows () == 3 &&
498 m.Cols () == 3);
499
500 const bool is4by4 = (m.Rows () == 4 &&
501 m.Cols () == 4);
502
503 if (!is3by3 && !is4by4)
504 {
505
506 ThrowMatrixMath ();
507
508 }
509
510 // For 3x3 case, pad to 4x4 (equivalent 4x4 matrix).
511
512 if (is3by3)
513 {
514
515 fRows = 4;
516 fCols = 4;
517
518 fData [0] [3] = 0.0;
519 fData [1] [3] = 0.0;
520 fData [2] [3] = 0.0;
521
522 fData [3] [0] = 0.0;
523 fData [3] [1] = 0.0;
524 fData [3] [2] = 0.0;
525
526 fData [3] [3] = 1.0;
527
528 }
529
530 }
531
532 /*****************************************************************************/
533
dng_matrix_4by4(real64 a00,real64 a01,real64 a02,real64 a03,real64 a10,real64 a11,real64 a12,real64 a13,real64 a20,real64 a21,real64 a22,real64 a23,real64 a30,real64 a31,real64 a32,real64 a33)534 dng_matrix_4by4::dng_matrix_4by4 (real64 a00, real64 a01, real64 a02, real64 a03,
535 real64 a10, real64 a11, real64 a12, real64 a13,
536 real64 a20, real64 a21, real64 a22, real64 a23,
537 real64 a30, real64 a31, real64 a32, real64 a33)
538
539 : dng_matrix (4, 4)
540
541 {
542
543 fData [0] [0] = a00;
544 fData [0] [1] = a01;
545 fData [0] [2] = a02;
546 fData [0] [3] = a03;
547
548 fData [1] [0] = a10;
549 fData [1] [1] = a11;
550 fData [1] [2] = a12;
551 fData [1] [3] = a13;
552
553 fData [2] [0] = a20;
554 fData [2] [1] = a21;
555 fData [2] [2] = a22;
556 fData [2] [3] = a23;
557
558 fData [3] [0] = a30;
559 fData [3] [1] = a31;
560 fData [3] [2] = a32;
561 fData [3] [3] = a33;
562
563 }
564
565 /*****************************************************************************/
566
dng_matrix_4by4(real64 a00,real64 a11,real64 a22,real64 a33)567 dng_matrix_4by4::dng_matrix_4by4 (real64 a00,
568 real64 a11,
569 real64 a22,
570 real64 a33)
571
572 : dng_matrix (4, 4)
573
574 {
575
576 fData [0] [0] = a00;
577 fData [1] [1] = a11;
578 fData [2] [2] = a22;
579 fData [3] [3] = a33;
580
581 }
582
583 /*****************************************************************************/
584
dng_vector()585 dng_vector::dng_vector ()
586
587 : fCount (0)
588
589 {
590
591 }
592
593 /*****************************************************************************/
594
dng_vector(uint32 count)595 dng_vector::dng_vector (uint32 count)
596
597 : fCount (0)
598
599 {
600
601 if (count < 1 || count > kMaxColorPlanes)
602 {
603
604 ThrowProgramError ();
605
606 }
607
608 fCount = count;
609
610 for (uint32 index = 0; index < fCount; index++)
611 {
612
613 fData [index] = 0.0;
614
615 }
616
617 }
618
619 /*****************************************************************************/
620
dng_vector(const dng_vector & v)621 dng_vector::dng_vector (const dng_vector &v)
622
623 : fCount (v.fCount)
624
625 {
626
627 for (uint32 index = 0; index < fCount; index++)
628 {
629
630 fData [index] = v.fData [index];
631
632 }
633
634 }
635
636 /*****************************************************************************/
637
Clear()638 void dng_vector::Clear ()
639 {
640
641 fCount = 0;
642
643 }
644
645 /*****************************************************************************/
646
SetIdentity(uint32 count)647 void dng_vector::SetIdentity (uint32 count)
648 {
649
650 *this = dng_vector (count);
651
652 for (uint32 j = 0; j < count; j++)
653 {
654
655 fData [j] = 1.0;
656
657 }
658
659 }
660
661 /******************************************************************************/
662
operator ==(const dng_vector & v) const663 bool dng_vector::operator== (const dng_vector &v) const
664 {
665
666 if (Count () != v.Count ())
667 {
668
669 return false;
670
671 }
672
673 for (uint32 j = 0; j < Count (); j++)
674 {
675
676 if (fData [j] != v.fData [j])
677 {
678
679 return false;
680
681 }
682
683 }
684
685 return true;
686
687 }
688
689 /******************************************************************************/
690
MaxEntry() const691 real64 dng_vector::MaxEntry () const
692 {
693
694 if (IsEmpty ())
695 {
696
697 return 0.0;
698
699 }
700
701 real64 m = fData [0];
702
703 for (uint32 j = 0; j < Count (); j++)
704 {
705
706 m = Max_real64 (m, fData [j]);
707
708 }
709
710 return m;
711
712 }
713
714 /******************************************************************************/
715
MinEntry() const716 real64 dng_vector::MinEntry () const
717 {
718
719 if (IsEmpty ())
720 {
721
722 return 0.0;
723
724 }
725
726 real64 m = fData [0];
727
728 for (uint32 j = 0; j < Count (); j++)
729 {
730
731 m = Min_real64 (m, fData [j]);
732
733 }
734
735 return m;
736
737 }
738
739 /*****************************************************************************/
740
Scale(real64 factor)741 void dng_vector::Scale (real64 factor)
742 {
743
744 for (uint32 j = 0; j < Count (); j++)
745 {
746
747 fData [j] *= factor;
748
749 }
750
751 }
752
753 /*****************************************************************************/
754
Round(real64 factor)755 void dng_vector::Round (real64 factor)
756 {
757
758 real64 invFactor = 1.0 / factor;
759
760 for (uint32 j = 0; j < Count (); j++)
761 {
762
763 fData [j] = Round_int32 (fData [j] * factor) * invFactor;
764
765 }
766
767 }
768
769 /*****************************************************************************/
770
AsDiagonal() const771 dng_matrix dng_vector::AsDiagonal () const
772 {
773
774 dng_matrix M (Count (), Count ());
775
776 for (uint32 j = 0; j < Count (); j++)
777 {
778
779 M [j] [j] = fData [j];
780
781 }
782
783 return M;
784
785 }
786
787 /*****************************************************************************/
788
AsColumn() const789 dng_matrix dng_vector::AsColumn () const
790 {
791
792 dng_matrix M (Count (), 1);
793
794 for (uint32 j = 0; j < Count (); j++)
795 {
796
797 M [j] [0] = fData [j];
798
799 }
800
801 return M;
802
803 }
804
805 /******************************************************************************/
806
dng_vector_3()807 dng_vector_3::dng_vector_3 ()
808
809 : dng_vector (3)
810
811 {
812 }
813
814 /******************************************************************************/
815
dng_vector_3(const dng_vector & v)816 dng_vector_3::dng_vector_3 (const dng_vector &v)
817
818 : dng_vector (v)
819
820 {
821
822 if (Count () != 3)
823 {
824
825 ThrowMatrixMath ();
826
827 }
828
829 }
830
831 /******************************************************************************/
832
dng_vector_3(real64 a0,real64 a1,real64 a2)833 dng_vector_3::dng_vector_3 (real64 a0,
834 real64 a1,
835 real64 a2)
836
837 : dng_vector (3)
838
839 {
840
841 fData [0] = a0;
842 fData [1] = a1;
843 fData [2] = a2;
844
845 }
846
847 /******************************************************************************/
848
dng_vector_4()849 dng_vector_4::dng_vector_4 ()
850
851 : dng_vector (4)
852
853 {
854 }
855
856 /******************************************************************************/
857
dng_vector_4(const dng_vector & v)858 dng_vector_4::dng_vector_4 (const dng_vector &v)
859
860 : dng_vector (v)
861
862 {
863
864 if (Count () != 4)
865 {
866
867 ThrowMatrixMath ();
868
869 }
870
871 }
872
873 /******************************************************************************/
874
dng_vector_4(real64 a0,real64 a1,real64 a2,real64 a3)875 dng_vector_4::dng_vector_4 (real64 a0,
876 real64 a1,
877 real64 a2,
878 real64 a3)
879
880 : dng_vector (4)
881
882 {
883
884 fData [0] = a0;
885 fData [1] = a1;
886 fData [2] = a2;
887 fData [3] = a3;
888
889 }
890
891 /******************************************************************************/
892
operator *(const dng_matrix & A,const dng_matrix & B)893 dng_matrix operator* (const dng_matrix &A,
894 const dng_matrix &B)
895 {
896
897 if (A.Cols () != B.Rows ())
898 {
899
900 ThrowMatrixMath ();
901
902 }
903
904 dng_matrix C (A.Rows (), B.Cols ());
905
906 for (uint32 j = 0; j < C.Rows (); j++)
907 for (uint32 k = 0; k < C.Cols (); k++)
908 {
909
910 C [j] [k] = 0.0;
911
912 for (uint32 m = 0; m < A.Cols (); m++)
913 {
914
915 real64 aa = A [j] [m];
916
917 real64 bb = B [m] [k];
918
919 C [j] [k] += aa * bb;
920
921 }
922
923 }
924
925 return C;
926
927 }
928
929 /******************************************************************************/
930
operator *(const dng_matrix & A,const dng_vector & B)931 dng_vector operator* (const dng_matrix &A,
932 const dng_vector &B)
933 {
934
935 if (A.Cols () != B.Count ())
936 {
937
938 ThrowMatrixMath ();
939
940 }
941
942 dng_vector C (A.Rows ());
943
944 for (uint32 j = 0; j < C.Count (); j++)
945 {
946
947 C [j] = 0.0;
948
949 for (uint32 m = 0; m < A.Cols (); m++)
950 {
951
952 real64 aa = A [j] [m];
953
954 real64 bb = B [m];
955
956 C [j] += aa * bb;
957
958 }
959
960 }
961
962 return C;
963
964 }
965
966 /******************************************************************************/
967
operator *(real64 scale,const dng_matrix & A)968 dng_matrix operator* (real64 scale,
969 const dng_matrix &A)
970 {
971
972 dng_matrix B (A);
973
974 B.Scale (scale);
975
976 return B;
977
978 }
979
980 /******************************************************************************/
981
operator *(real64 scale,const dng_vector & A)982 dng_vector operator* (real64 scale,
983 const dng_vector &A)
984 {
985
986 dng_vector B (A);
987
988 B.Scale (scale);
989
990 return B;
991
992 }
993
994 /******************************************************************************/
995
operator +(const dng_matrix & A,const dng_matrix & B)996 dng_matrix operator+ (const dng_matrix &A,
997 const dng_matrix &B)
998 {
999
1000 if (A.Cols () != B.Cols () ||
1001 A.Rows () != B.Rows ())
1002 {
1003
1004 ThrowMatrixMath ();
1005
1006 }
1007
1008 dng_matrix C (A);
1009
1010 for (uint32 j = 0; j < C.Rows (); j++)
1011 for (uint32 k = 0; k < C.Cols (); k++)
1012 {
1013
1014 C [j] [k] += B [j] [k];
1015
1016 }
1017
1018 return C;
1019
1020 }
1021
1022 /******************************************************************************/
1023
operator -(const dng_vector & a,const dng_vector & b)1024 dng_vector operator- (const dng_vector &a,
1025 const dng_vector &b)
1026 {
1027
1028 uint32 count = a.Count ();
1029
1030 DNG_REQUIRE (count == b.Count (),
1031 "Mismatch count in Dot");
1032
1033 if (!count)
1034 {
1035 return dng_vector ();
1036 }
1037
1038 dng_vector result (count);
1039
1040 for (uint32 i = 0; i < count; i++)
1041 {
1042
1043 result [i] = a [i] - b [i];
1044
1045 }
1046
1047 return result;
1048
1049 }
1050
1051 /******************************************************************************/
1052
1053 const real64 kNearZero = 1.0E-10;
1054
1055 /******************************************************************************/
1056
1057 // Work around bug #1294195, which may be a hardware problem on a specific machine.
1058 // This pragma turns on "improved" floating-point consistency.
1059 #ifdef _MSC_VER
1060 #pragma optimize ("p", on)
1061 #endif
1062
Invert3by3(const dng_matrix & A)1063 static dng_matrix Invert3by3 (const dng_matrix &A)
1064 {
1065
1066 real64 a00 = A [0] [0];
1067 real64 a01 = A [0] [1];
1068 real64 a02 = A [0] [2];
1069 real64 a10 = A [1] [0];
1070 real64 a11 = A [1] [1];
1071 real64 a12 = A [1] [2];
1072 real64 a20 = A [2] [0];
1073 real64 a21 = A [2] [1];
1074 real64 a22 = A [2] [2];
1075
1076 real64 temp [3] [3];
1077
1078 temp [0] [0] = a11 * a22 - a21 * a12;
1079 temp [0] [1] = a21 * a02 - a01 * a22;
1080 temp [0] [2] = a01 * a12 - a11 * a02;
1081 temp [1] [0] = a20 * a12 - a10 * a22;
1082 temp [1] [1] = a00 * a22 - a20 * a02;
1083 temp [1] [2] = a10 * a02 - a00 * a12;
1084 temp [2] [0] = a10 * a21 - a20 * a11;
1085 temp [2] [1] = a20 * a01 - a00 * a21;
1086 temp [2] [2] = a00 * a11 - a10 * a01;
1087
1088 real64 det = (a00 * temp [0] [0] +
1089 a01 * temp [1] [0] +
1090 a02 * temp [2] [0]);
1091
1092 if (Abs_real64 (det) < kNearZero)
1093 {
1094
1095 ThrowMatrixMath ();
1096
1097 }
1098
1099 dng_matrix B (3, 3);
1100
1101 for (uint32 j = 0; j < 3; j++)
1102 for (uint32 k = 0; k < 3; k++)
1103 {
1104
1105 B [j] [k] = temp [j] [k] / det;
1106
1107 }
1108
1109 return B;
1110
1111 }
1112
1113 // Reset floating-point optimization. See comment above.
1114 #ifdef _MSC_VER
1115 #pragma optimize ("p", off)
1116 #endif
1117
1118 /******************************************************************************/
1119
InvertNbyN(const dng_matrix & A)1120 static dng_matrix InvertNbyN (const dng_matrix &A)
1121 {
1122
1123 uint32 i;
1124 uint32 j;
1125 uint32 k;
1126
1127 const uint32 n = A.Rows ();
1128
1129 const uint32 augmented_cols = 2 * n;
1130
1131 real64 temp [kMaxColorPlanes] [kMaxColorPlanes * 2];
1132
1133 memset (temp, 0, sizeof (temp));
1134
1135 for (i = 0; i < n; i++)
1136 for (j = 0; j < n; j++)
1137 {
1138
1139 temp [i] [j ] = A [i] [j];
1140
1141 temp [i] [j + n] = (i == j ? 1.0 : 0.0);
1142
1143 }
1144
1145 for (i = 0; i < n; i++)
1146 {
1147
1148 // Find row iMax with largest absolute entry in column i.
1149
1150 uint32 iMax = i;
1151 real64 vMax = -1.0;
1152
1153 for (k = i; k < n; k++)
1154 {
1155
1156 real64 v = Abs_real64 (A [k] [i]);
1157
1158 if (v > vMax)
1159 {
1160 vMax = v;
1161 iMax = k;
1162 }
1163
1164 }
1165
1166 real64 alpha = temp [iMax] [i];
1167
1168 if (Abs_real64 (alpha) < kNearZero)
1169 {
1170
1171 ThrowMatrixMath ();
1172
1173 }
1174
1175 // Swap rows i and iMax, column by column.
1176
1177 if (i != iMax)
1178 {
1179
1180 for (j = 0; j < augmented_cols; j++)
1181 {
1182
1183 std::swap (temp [i ] [j],
1184 temp [iMax] [j]);
1185
1186 }
1187
1188 }
1189
1190 for (j = 0; j < augmented_cols; j++)
1191 {
1192
1193 temp [i] [j] /= alpha;
1194
1195 }
1196
1197 for (k = 0; k < n; k++)
1198 {
1199
1200 if (i != k)
1201 {
1202
1203 real64 beta = temp [k] [i];
1204
1205 for (j = 0; j < augmented_cols; j++)
1206 {
1207
1208 temp [k] [j] -= beta * temp [i] [j];
1209
1210 }
1211
1212 }
1213
1214 }
1215
1216 }
1217
1218 dng_matrix B (n, n);
1219
1220 for (i = 0; i < n; i++)
1221 for (j = 0; j < n; j++)
1222 {
1223
1224 B [i] [j] = temp [i] [j + n];
1225
1226 }
1227
1228 return B;
1229
1230 }
1231
1232 /******************************************************************************/
1233
Transpose(const dng_matrix & A)1234 dng_matrix Transpose (const dng_matrix &A)
1235 {
1236
1237 dng_matrix B (A.Cols (), A.Rows ());
1238
1239 for (uint32 j = 0; j < B.Rows (); j++)
1240 for (uint32 k = 0; k < B.Cols (); k++)
1241 {
1242
1243 B [j] [k] = A [k] [j];
1244
1245 }
1246
1247 return B;
1248
1249 }
1250
1251 /******************************************************************************/
1252
Invert(const dng_matrix & A)1253 dng_matrix Invert (const dng_matrix &A)
1254 {
1255
1256 if (A.Rows () < 2 || A.Cols () < 2)
1257 {
1258
1259 ThrowMatrixMath ();
1260
1261 }
1262
1263 if (A.Rows () == A.Cols ())
1264 {
1265
1266 if (A.Rows () == 3)
1267 {
1268
1269 return Invert3by3 (A);
1270
1271 }
1272
1273 return InvertNbyN (A);
1274
1275 }
1276
1277 else
1278 {
1279
1280 // Compute the pseudo inverse.
1281
1282 dng_matrix B = Transpose (A);
1283
1284 return Invert (B * A) * B;
1285
1286 }
1287
1288 }
1289
1290 /******************************************************************************/
1291
Invert(const dng_matrix & A,const dng_matrix & hint)1292 dng_matrix Invert (const dng_matrix &A,
1293 const dng_matrix &hint)
1294 {
1295
1296 if (A.Rows () == A .Cols () ||
1297 A.Rows () != hint.Cols () ||
1298 A.Cols () != hint.Rows ())
1299 {
1300
1301 return Invert (A);
1302
1303 }
1304
1305 else
1306 {
1307
1308 // Use the specified hint matrix.
1309
1310 return Invert (hint * A) * hint;
1311
1312 }
1313
1314 }
1315
1316 /*****************************************************************************/
1317
Dot(const dng_vector & a,const dng_vector & b)1318 real64 Dot (const dng_vector &a,
1319 const dng_vector &b)
1320 {
1321
1322 DNG_REQUIRE (a.Count () == b.Count (),
1323 "Cannot take dot product between vectors of different size.");
1324
1325 // DNG_REQUIRE (a.Count () > 0,
1326 // "Cannot take dot product with an empty vector.");
1327
1328 real64 sum = 0.0;
1329
1330 for (uint32 j = 0; j < a.Count (); j++)
1331 {
1332
1333 sum += (a [j] * b [j]);
1334
1335 }
1336
1337 return sum;
1338
1339 }
1340
1341 /*****************************************************************************/
1342
Distance(const dng_vector & a,const dng_vector & b)1343 real64 Distance (const dng_vector &a,
1344 const dng_vector &b)
1345 {
1346
1347 dng_vector c = a - b;
1348
1349 return sqrt (Dot (c, c));
1350
1351 }
1352
1353 /*****************************************************************************/
1354