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