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
