1 /*****************************************************************************/ 2 // Copyright 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_lens_correction.cpp#1 $ */ 10 /* $DateTime: 2012/05/30 13:28:51 $ */ 11 /* $Change: 832332 $ */ 12 /* $Author: tknoll $ */ 13 14 /*****************************************************************************/ 15 16 #include <cfloat> 17 #include <limits.h> 18 19 #include "dng_1d_table.h" 20 #include "dng_assertions.h" 21 #include "dng_bottlenecks.h" 22 #include "dng_exceptions.h" 23 #include "dng_filter_task.h" 24 #include "dng_globals.h" 25 #include "dng_host.h" 26 #include "dng_image.h" 27 #include "dng_lens_correction.h" 28 #include "dng_negative.h" 29 #include "dng_safe_arithmetic.h" 30 #include "dng_sdk_limits.h" 31 #include "dng_tag_values.h" 32 #include "dng_utils.h" 33 34 /*****************************************************************************/ 35 36 dng_warp_params::dng_warp_params () 37 38 : fPlanes (1) 39 , fCenter (0.5, 0.5) 40 41 { 42 43 } 44 45 /*****************************************************************************/ 46 47 dng_warp_params::dng_warp_params (uint32 planes, 48 const dng_point_real64 ¢er) 49 50 : fPlanes (planes) 51 , fCenter (center) 52 53 { 54 55 DNG_ASSERT (planes >= 1, "Too few planes." ); 56 DNG_ASSERT (planes <= kMaxColorPlanes, "Too many planes."); 57 58 DNG_ASSERT (fCenter.h >= 0.0 && fCenter.h <= 1.0, 59 "Center (horizontal) out of range."); 60 61 DNG_ASSERT (fCenter.v >= 0.0 && fCenter.v <= 1.0, 62 "Center (vertical) out of range."); 63 64 } 65 66 /*****************************************************************************/ 67 68 dng_warp_params::~dng_warp_params () 69 { 70 71 } 72 73 /*****************************************************************************/ 74 75 bool dng_warp_params::IsNOPAll () const 76 { 77 78 return IsRadNOPAll () && 79 IsTanNOPAll (); 80 81 } 82 83 /*****************************************************************************/ 84 85 bool dng_warp_params::IsNOP (uint32 plane) const 86 { 87 88 return IsRadNOP (plane) && 89 IsTanNOP (plane); 90 91 } 92 93 /*****************************************************************************/ 94 95 bool dng_warp_params::IsRadNOPAll () const 96 { 97 98 for (uint32 plane = 0; plane < fPlanes; plane++) 99 { 100 101 if (!IsRadNOP (plane)) 102 { 103 return false; 104 } 105 106 } 107 108 return true; 109 110 } 111 112 /*****************************************************************************/ 113 114 bool dng_warp_params::IsRadNOP (uint32 /* plane */) const 115 { 116 117 return false; 118 119 } 120 121 /*****************************************************************************/ 122 123 bool dng_warp_params::IsTanNOPAll () const 124 { 125 126 for (uint32 plane = 0; plane < fPlanes; plane++) 127 { 128 129 if (!IsTanNOP (plane)) 130 { 131 return false; 132 } 133 134 } 135 136 return true; 137 138 } 139 140 /*****************************************************************************/ 141 142 bool dng_warp_params::IsTanNOP (uint32 /* plane */) const 143 { 144 145 return false; 146 147 } 148 149 /*****************************************************************************/ 150 151 bool dng_warp_params::IsValid () const 152 { 153 154 if (fPlanes < 1 || fPlanes > kMaxColorPlanes) 155 { 156 157 return false; 158 159 } 160 161 if (fCenter.h < 0.0 || 162 fCenter.h > 1.0 || 163 fCenter.v < 0.0 || 164 fCenter.v > 1.0) 165 { 166 167 return false; 168 169 } 170 171 return true; 172 173 } 174 175 /*****************************************************************************/ 176 177 bool dng_warp_params::IsValidForNegative (const dng_negative &negative) const 178 { 179 180 if (!IsValid ()) 181 { 182 return false; 183 } 184 185 if ((fPlanes != 1) && 186 (fPlanes != negative.ColorChannels ())) 187 { 188 return false; 189 } 190 191 return true; 192 193 } 194 195 /*****************************************************************************/ 196 197 real64 dng_warp_params::EvaluateInverse (uint32 plane, 198 real64 y) const 199 { 200 201 const uint32 kMaxIterations = 30; 202 const real64 kNearZero = 1.0e-10; 203 204 real64 x0 = 0.0; 205 real64 y0 = Evaluate (plane, 206 x0); 207 208 real64 x1 = 1.0; 209 real64 y1 = Evaluate (plane, 210 x1); 211 212 for (uint32 iteration = 0; iteration < kMaxIterations; iteration++) 213 { 214 215 if (Abs_real64 (y1 - y0) < kNearZero) 216 { 217 break; 218 } 219 220 const real64 x2 = Pin_real64 (0.0, 221 x1 + (y - y1) * (x1 - x0) / (y1 - y0), 222 1.0); 223 224 const real64 y2 = Evaluate (plane, 225 x2); 226 227 x0 = x1; 228 y0 = y1; 229 230 x1 = x2; 231 y1 = y2; 232 233 } 234 235 return x1; 236 237 } 238 239 /*****************************************************************************/ 240 241 dng_point_real64 dng_warp_params::EvaluateTangential2 (uint32 plane, 242 const dng_point_real64 &diff) const 243 { 244 245 const real64 dvdv = diff.v * diff.v; 246 const real64 dhdh = diff.h * diff.h; 247 248 const real64 rr = dvdv + dhdh; 249 250 dng_point_real64 diffSqr (dvdv, 251 dhdh); 252 253 return EvaluateTangential (plane, 254 rr, 255 diff, 256 diffSqr); 257 258 } 259 260 /*****************************************************************************/ 261 262 dng_point_real64 dng_warp_params::EvaluateTangential3 (uint32 plane, 263 real64 r2, 264 const dng_point_real64 &diff) const 265 { 266 267 dng_point_real64 diffSqr (diff.v * diff.v, 268 diff.h * diff.h); 269 270 return EvaluateTangential (plane, 271 r2, 272 diff, 273 diffSqr); 274 275 } 276 277 /*****************************************************************************/ 278 279 void dng_warp_params::Dump () const 280 { 281 282 #if qDNGValidate 283 284 printf ("Planes: %u\n", (unsigned) fPlanes); 285 286 printf (" Optical center:\n" 287 " h = %.6lf\n" 288 " v = %.6lf\n", 289 fCenter.h, 290 fCenter.v); 291 292 #endif 293 294 } 295 296 /*****************************************************************************/ 297 298 dng_warp_params_rectilinear::dng_warp_params_rectilinear () 299 300 : dng_warp_params () 301 302 { 303 304 for (uint32 plane = 0; plane < kMaxColorPlanes; plane++) 305 { 306 307 fRadParams [plane] = dng_vector (4); 308 fTanParams [plane] = dng_vector (2); 309 310 fRadParams [plane][0] = 1.0; 311 312 } 313 314 } 315 316 /*****************************************************************************/ 317 318 dng_warp_params_rectilinear::dng_warp_params_rectilinear (uint32 planes, 319 const dng_vector radParams [], 320 const dng_vector tanParams [], 321 const dng_point_real64 ¢er) 322 323 : dng_warp_params (planes, 324 center) 325 326 { 327 328 for (uint32 i = 0; i < fPlanes; i++) 329 { 330 fRadParams [i] = radParams [i]; 331 fTanParams [i] = tanParams [i]; 332 } 333 334 } 335 336 /*****************************************************************************/ 337 338 dng_warp_params_rectilinear::~dng_warp_params_rectilinear () 339 { 340 341 } 342 343 /*****************************************************************************/ 344 345 bool dng_warp_params_rectilinear::IsRadNOP (uint32 plane) const 346 { 347 348 DNG_ASSERT (plane < fPlanes, "plane out of range."); 349 350 const dng_vector &r = fRadParams [plane]; 351 352 return (r [0] == 1.0 && 353 r [1] == 0.0 && 354 r [2] == 0.0 && 355 r [3] == 0.0); 356 357 } 358 359 /*****************************************************************************/ 360 361 bool dng_warp_params_rectilinear::IsTanNOP (uint32 plane) const 362 { 363 364 DNG_ASSERT (plane < fPlanes, "plane out of range."); 365 366 const dng_vector &t = fTanParams [plane]; 367 368 return (t [0] == 0.0 && 369 t [1] == 0.0); 370 371 } 372 373 /*****************************************************************************/ 374 375 bool dng_warp_params_rectilinear::IsValid () const 376 { 377 378 for (uint32 plane = 0; plane < fPlanes; plane++) 379 { 380 381 if (fRadParams [plane].Count () != 4) 382 { 383 return false; 384 } 385 386 if (fTanParams [plane].Count () < 2) 387 { 388 return false; 389 } 390 391 } 392 393 return dng_warp_params::IsValid (); 394 395 } 396 397 /*****************************************************************************/ 398 399 void dng_warp_params_rectilinear::PropagateToAllPlanes (uint32 totalPlanes) 400 { 401 402 for (uint32 plane = fPlanes; plane < totalPlanes; plane++) 403 { 404 405 fRadParams [plane] = fRadParams [0]; 406 fTanParams [plane] = fTanParams [0]; 407 408 } 409 410 } 411 412 /*****************************************************************************/ 413 414 real64 dng_warp_params_rectilinear::Evaluate (uint32 plane, 415 real64 x) const 416 { 417 418 const dng_vector &K = fRadParams [plane]; // Coefficients. 419 420 const real64 x2 = x * x; 421 422 return x * (K [0] + x2 * (K [1] + x2 * (K [2] + x2 * K [3]))); 423 424 } 425 426 /*****************************************************************************/ 427 428 real64 dng_warp_params_rectilinear::EvaluateRatio (uint32 plane, 429 real64 r2) const 430 { 431 432 const dng_vector &K = fRadParams [plane]; // Coefficients. 433 434 return K [0] + r2 * (K [1] + r2 * (K [2] + r2 * K [3])); 435 436 } 437 438 /*****************************************************************************/ 439 440 dng_point_real64 dng_warp_params_rectilinear::EvaluateTangential (uint32 plane, 441 real64 r2, 442 const dng_point_real64 &diff, 443 const dng_point_real64 &diff2) const 444 { 445 446 const real64 kt0 = fTanParams [plane][0]; 447 const real64 kt1 = fTanParams [plane][1]; 448 449 const real64 dh = diff.h; 450 const real64 dv = diff.v; 451 452 const real64 dhdh = diff2.h; 453 const real64 dvdv = diff2.v; 454 455 return dng_point_real64 (kt0 * (r2 + 2.0 * dvdv) + (2.0 * kt1 * dh * dv), // v 456 kt1 * (r2 + 2.0 * dhdh) + (2.0 * kt0 * dh * dv)); // h 457 458 } 459 460 /*****************************************************************************/ 461 462 real64 dng_warp_params_rectilinear::MaxSrcRadiusGap (real64 maxDstGap) const 463 { 464 465 real64 maxSrcGap = 0.0; 466 467 for (uint32 plane = 0; plane < fPlanes; plane++) 468 { 469 470 const dng_vector &coefs = fRadParams [plane]; 471 472 const real64 k3 = coefs [1]; 473 const real64 k5 = coefs [2]; 474 const real64 k7 = coefs [3]; 475 476 // 477 // Let f (r) be the radius warp function. Consider the function 478 // 479 // gap (r) = f (r + maxDstGap) - f (r) 480 // 481 // We wish to maximize gap (r) over the domain [0, 1 - maxDstGap]. This will 482 // give us a reasonable upper bound on the src tile size, independent of 483 // dstBounds. 484 // 485 // As usual, we maximize gap (r) by examining its critical points, i.e., by 486 // considering the roots of its derivative which lie in the domain [0, 1 - 487 // maxDstGap]. gap' (r) is a 5th-order polynomial. One of its roots is 488 // -maxDstGap / 2, which is negative and hence lies outside the domain of 489 // interest. This leaves 4 other possible roots. We solve for these 490 // analytically below. 491 // 492 493 real64 roots [4]; 494 uint32 numRoots = 0; 495 496 if (k7 == 0.0) 497 { 498 499 if (k5 == 0.0) 500 { 501 502 // No roots in [0,1]. 503 504 } 505 506 else 507 { 508 509 // k7 is zero, but k5 is non-zero. At most two real roots. 510 511 const real64 discrim = 25.0 * (-6.0 * k3 * k5 - 5.0 * k5 * maxDstGap * maxDstGap); 512 513 if (discrim >= 0.0) 514 { 515 516 // Two real roots. 517 518 const real64 scale = 0.1 * k5; 519 const real64 offset = -5.0 * maxDstGap * k5; 520 const real64 sDiscrim = sqrt (discrim); 521 522 roots [numRoots++] = scale * (offset + sDiscrim); 523 roots [numRoots++] = scale * (offset - sDiscrim); 524 525 } 526 527 } 528 529 } 530 531 else 532 { 533 534 // k7 is non-zero. Up to 4 real roots. 535 536 const real64 d = maxDstGap; 537 const real64 d2 = d * d; 538 const real64 d4 = d2 * d2; 539 540 const real64 discrim = 25.0 * k5 * k5 541 - 63.0 * k3 * k7 542 + 35.0 * d2 * k5 * k7 543 + 49.0 * d4 * k7 * k7; 544 545 if (discrim >= 0.0) 546 { 547 548 const real64 sDiscrim = 4.0 * k7 * sqrt (discrim); 549 550 const real64 offset = -20.0 * k5 * k7 - 35.0 * d2 * k7 * k7; 551 552 const real64 discrim1 = offset - sDiscrim; 553 const real64 discrim2 = offset + sDiscrim; 554 555 const real64 scale = sqrt (21.0) / (42.0 * k7); 556 557 if (discrim1 >= 0.0) 558 { 559 560 const real64 offset1 = -d * 0.5; 561 const real64 sDiscrim1 = scale * sqrt (discrim1); 562 563 roots [numRoots++] = offset1 + sDiscrim1; 564 roots [numRoots++] = offset1 - sDiscrim1; 565 566 } 567 568 if (discrim2 >= 0.0) 569 { 570 571 const real64 offset2 = -d * 0.5; 572 const real64 sDiscrim2 = scale * sqrt (discrim2); 573 574 roots [numRoots++] = offset2 + sDiscrim2; 575 roots [numRoots++] = offset2 - sDiscrim2; 576 577 } 578 579 } 580 581 } 582 583 real64 planeMaxSrcGap = 0.0; 584 585 // Examine the endpoints. 586 587 { 588 589 // Check left endpoint: f (maxDstGap) - f (0). Remember that f (0) == 0. 590 591 const real64 gap1 = Evaluate (plane, maxDstGap); 592 593 planeMaxSrcGap = Max_real64 (planeMaxSrcGap, gap1); 594 595 // Check right endpoint: f (1) - f (1 - maxDstGap). 596 597 const real64 gap2 = Evaluate (plane, 1.0) 598 - Evaluate (plane, 1.0 - maxDstGap); 599 600 planeMaxSrcGap = Max_real64 (planeMaxSrcGap, gap2); 601 602 } 603 604 // Examine the roots we found, if any. 605 606 for (uint32 i = 0; i < numRoots; i++) 607 { 608 609 const real64 r = roots [i]; 610 611 if (r > 0.0 && r < 1.0 - maxDstGap) 612 { 613 614 const real64 gap = Evaluate (plane, r + maxDstGap) 615 - Evaluate (plane, r); 616 617 planeMaxSrcGap = Max_real64 (planeMaxSrcGap, gap); 618 619 } 620 621 } 622 623 maxSrcGap = Max_real64 (maxSrcGap, 624 planeMaxSrcGap); 625 626 } 627 628 return maxSrcGap; 629 630 } 631 632 /*****************************************************************************/ 633 634 dng_point_real64 dng_warp_params_rectilinear::MaxSrcTanGap (dng_point_real64 minDst, 635 dng_point_real64 maxDst) const 636 { 637 638 const real64 v [] = { minDst.v, maxDst.v, 0.0 }; 639 const real64 h [] = { minDst.h, maxDst.h, 0.0 }; 640 641 dng_point_real64 maxGap; 642 643 for (uint32 plane = 0; plane < fPlanes; plane++) 644 { 645 646 real64 hMin = +FLT_MAX; 647 real64 hMax = -FLT_MAX; 648 649 real64 vMin = +FLT_MAX; 650 real64 vMax = -FLT_MAX; 651 652 for (uint32 i = 0; i < 3; i++) 653 { 654 655 for (uint32 j = 0; j < 3; j++) 656 { 657 658 dng_point_real64 dstDiff (v [i], 659 h [j]); 660 661 dng_point_real64 srcDiff = EvaluateTangential2 (plane, 662 dstDiff); 663 664 hMin = Min_real64 (hMin, srcDiff.h); 665 hMax = Max_real64 (hMax, srcDiff.h); 666 667 vMin = Min_real64 (vMin, srcDiff.v); 668 vMax = Max_real64 (vMax, srcDiff.v); 669 670 } 671 672 } 673 674 const real64 hGap = hMax - hMin; 675 const real64 vGap = vMax - vMin; 676 677 maxGap.h = Max_real64 (maxGap.h, hGap); 678 maxGap.v = Max_real64 (maxGap.v, vGap); 679 680 } 681 682 return maxGap; 683 684 } 685 686 /*****************************************************************************/ 687 688 void dng_warp_params_rectilinear::Dump () const 689 { 690 691 #if qDNGValidate 692 693 dng_warp_params::Dump (); 694 695 for (uint32 plane = 0; plane < fPlanes; plane++) 696 { 697 698 printf (" Plane %u:\n", (unsigned) plane); 699 700 printf (" Radial params: %.6lf, %.6lf, %.6lf, %.6lf\n", 701 fRadParams [plane][0], 702 fRadParams [plane][1], 703 fRadParams [plane][2], 704 fRadParams [plane][3]); 705 706 printf (" Tangential params: %.6lf, %.6lf\n", 707 fTanParams [plane][0], 708 fTanParams [plane][1]); 709 710 } 711 712 #endif 713 714 } 715 716 /*****************************************************************************/ 717 718 dng_warp_params_fisheye::dng_warp_params_fisheye () 719 720 : dng_warp_params () 721 722 { 723 724 for (uint32 plane = 0; plane < kMaxColorPlanes; plane++) 725 { 726 727 fRadParams [plane] = dng_vector (4); 728 729 } 730 731 } 732 733 /*****************************************************************************/ 734 735 dng_warp_params_fisheye::dng_warp_params_fisheye (uint32 planes, 736 const dng_vector radParams [], 737 const dng_point_real64 ¢er) 738 739 : dng_warp_params (planes, center) 740 741 { 742 743 for (uint32 i = 0; i < fPlanes; i++) 744 { 745 746 fRadParams [i] = radParams [i]; 747 748 } 749 750 } 751 752 /*****************************************************************************/ 753 754 dng_warp_params_fisheye::~dng_warp_params_fisheye () 755 { 756 757 } 758 759 /*****************************************************************************/ 760 761 bool dng_warp_params_fisheye::IsRadNOP (uint32 /* plane */) const 762 { 763 764 return false; 765 766 } 767 768 /*****************************************************************************/ 769 770 bool dng_warp_params_fisheye::IsTanNOP (uint32 /* plane */) const 771 { 772 773 return true; 774 775 } 776 777 /*****************************************************************************/ 778 779 bool dng_warp_params_fisheye::IsValid () const 780 { 781 782 for (uint32 plane = 0; plane < fPlanes; plane++) 783 { 784 785 if (fRadParams [plane].Count () != 4) 786 { 787 return false; 788 } 789 790 } 791 792 return dng_warp_params::IsValid (); 793 794 } 795 796 /*****************************************************************************/ 797 798 void dng_warp_params_fisheye::PropagateToAllPlanes (uint32 totalPlanes) 799 { 800 801 for (uint32 plane = fPlanes; plane < totalPlanes; plane++) 802 { 803 804 fRadParams [plane] = fRadParams [0]; 805 806 } 807 808 } 809 810 /*****************************************************************************/ 811 812 real64 dng_warp_params_fisheye::Evaluate (uint32 plane, 813 real64 r) const 814 { 815 816 const real64 t = atan (r); 817 818 const dng_vector &K = fRadParams [plane]; 819 820 const real64 t2 = t * t; 821 822 return t * (K [0] + t2 * (K [1] + t2 * (K [2] + t2 * K [3]))); 823 824 } 825 826 /*****************************************************************************/ 827 828 real64 dng_warp_params_fisheye::EvaluateRatio (uint32 plane, 829 real64 rSqr) const 830 { 831 832 const real64 eps = 1.0e-12; 833 834 if (rSqr < eps) 835 { 836 837 // r is very close to zero. 838 839 return 1.0; 840 841 } 842 843 const real64 r = sqrt (rSqr); 844 845 return Evaluate (plane, r) / r; 846 847 } 848 849 /*****************************************************************************/ 850 851 dng_point_real64 dng_warp_params_fisheye::EvaluateTangential (uint32 /* plane */, 852 real64 /* r2 */, 853 const dng_point_real64 & /* diff */, 854 const dng_point_real64 & /* diff2 */) const 855 { 856 857 // This fisheye model does not support tangential warping. 858 859 ThrowProgramError (); 860 861 return dng_point_real64 (0.0, 0.0); 862 863 } 864 865 /*****************************************************************************/ 866 867 real64 dng_warp_params_fisheye::MaxSrcRadiusGap (real64 maxDstGap) const 868 { 869 870 // 871 // Let f (r) be the radius warp function. Consider the function 872 // 873 // gap (r) = f (r + maxDstGap) - f (r) 874 // 875 // We wish to maximize gap (r) over the domain [0, 1 - maxDstGap]. This will 876 // give us a reasonable upper bound on the src tile size, independent of 877 // dstBounds. 878 // 879 // Ideally, we'd like to maximize gap (r) by examining its critical points, 880 // i.e., by considering the roots of its derivative which lie in the domain [0, 881 // 1 - maxDstGap]. However, gap' (r) is too complex to find its roots 882 // analytically. 883 // 884 885 real64 maxSrcGap = 0.0; 886 887 DNG_REQUIRE (maxDstGap > 0.0, "maxDstGap must be positive."); 888 889 const real64 kMaxValue = 1.0 - maxDstGap; 890 891 const uint32 kSteps = 128; 892 893 const real64 kNorm = kMaxValue / real64 (kSteps - 1); 894 895 for (uint32 plane = 0; plane < fPlanes; plane++) 896 { 897 898 for (uint32 i = 0; i < kSteps; i++) 899 { 900 901 const real64 tt = i * kNorm; 902 903 const real64 gap = Evaluate (plane, tt + maxDstGap) 904 - Evaluate (plane, tt); 905 906 maxSrcGap = Max_real64 (maxSrcGap, 907 gap); 908 909 } 910 911 } 912 913 return maxSrcGap; 914 915 } 916 917 /*****************************************************************************/ 918 919 dng_point_real64 dng_warp_params_fisheye::MaxSrcTanGap (dng_point_real64 /* minDst */, 920 dng_point_real64 /* maxDst */) const 921 { 922 923 // This fisheye model does not support tangential distortion. 924 925 return dng_point_real64 (0.0, 0.0); 926 927 } 928 929 /*****************************************************************************/ 930 931 void dng_warp_params_fisheye::Dump () const 932 { 933 934 #if qDNGValidate 935 936 dng_warp_params::Dump (); 937 938 for (uint32 plane = 0; plane < fPlanes; plane++) 939 { 940 941 printf (" Plane %u:\n", (unsigned) plane); 942 943 printf (" Radial params: %.6lf, %.6lf, %.6lf, %.6lf\n", 944 fRadParams [plane][0], 945 fRadParams [plane][1], 946 fRadParams [plane][2], 947 fRadParams [plane][3]); 948 949 } 950 951 #endif 952 953 } 954 955 /*****************************************************************************/ 956 957 class dng_filter_warp: public dng_filter_task 958 { 959 960 protected: 961 962 AutoPtr<dng_warp_params> fParams; 963 964 dng_point_real64 fCenter; 965 966 dng_resample_weights_2d fWeights; 967 968 real64 fNormRadius; 969 real64 fInvNormRadius; 970 971 bool fIsRadNOP; 972 bool fIsTanNOP; 973 974 const real64 fPixelScaleV; 975 const real64 fPixelScaleVInv; 976 977 public: 978 979 dng_filter_warp (const dng_image &srcImage, 980 dng_image &dstImage, 981 const dng_negative &negative, 982 AutoPtr<dng_warp_params> ¶ms); 983 984 virtual void Initialize (dng_host &host); 985 986 virtual dng_rect SrcArea (const dng_rect &dstArea); 987 988 virtual dng_point SrcTileSize (const dng_point &dstTileSize); 989 990 virtual void ProcessArea (uint32 threadIndex, 991 dng_pixel_buffer &srcBuffer, 992 dng_pixel_buffer &dstBuffer); 993 994 virtual dng_point_real64 GetSrcPixelPosition (const dng_point_real64 &dst, 995 uint32 plane); 996 997 }; 998 999 /*****************************************************************************/ 1000 1001 dng_filter_warp::dng_filter_warp (const dng_image &srcImage, 1002 dng_image &dstImage, 1003 const dng_negative &negative, 1004 AutoPtr<dng_warp_params> ¶ms) 1005 1006 : dng_filter_task (srcImage, 1007 dstImage) 1008 1009 , fParams (params.Release ()) 1010 1011 , fCenter () 1012 1013 , fWeights () 1014 1015 , fNormRadius (1.0) 1016 , fInvNormRadius (1.0) 1017 1018 , fIsRadNOP (false) 1019 , fIsTanNOP (false) 1020 1021 , fPixelScaleV (1.0 / negative.PixelAspectRatio ()) 1022 , fPixelScaleVInv (1.0 / fPixelScaleV) 1023 1024 { 1025 1026 // Force float processing. 1027 1028 fSrcPixelType = ttFloat; 1029 fDstPixelType = ttFloat; 1030 1031 fIsRadNOP = fParams->IsRadNOPAll (); 1032 fIsTanNOP = fParams->IsTanNOPAll (); 1033 1034 const uint32 negPlanes = negative.ColorChannels (); 1035 1036 DNG_ASSERT (negPlanes >= 1, "Too few planes." ); 1037 DNG_ASSERT (negPlanes <= kMaxColorPlanes, "Too many planes."); 1038 1039 (void) negPlanes; 1040 1041 // At least one set of params must do something interesting. 1042 1043 if (fIsRadNOP && fIsTanNOP) 1044 { 1045 ThrowProgramError (); 1046 } 1047 1048 // Make sure the warp params are valid for this negative. 1049 1050 if (!fParams->IsValidForNegative (negative)) 1051 { 1052 ThrowBadFormat (); 1053 } 1054 1055 // Compute center. 1056 1057 const dng_rect bounds = srcImage.Bounds (); 1058 1059 fCenter.h = Lerp_real64 ((real64) bounds.l, 1060 (real64) bounds.r, 1061 fParams->fCenter.h); 1062 1063 fCenter.v = Lerp_real64 ((real64) bounds.t, 1064 (real64) bounds.b, 1065 fParams->fCenter.v); 1066 1067 // Compute max pixel radius and derive various normalized radius values. Note 1068 // that when computing the max pixel radius, we must take into account the pixel 1069 // aspect ratio. 1070 1071 { 1072 1073 dng_rect squareBounds (bounds); 1074 1075 squareBounds.b = squareBounds.t + 1076 Round_int32 (fPixelScaleV * (real64) squareBounds.H ()); 1077 1078 const dng_point_real64 squareCenter (Lerp_real64 ((real64) squareBounds.t, 1079 (real64) squareBounds.b, 1080 fParams->fCenter.v), 1081 1082 Lerp_real64 ((real64) squareBounds.l, 1083 (real64) squareBounds.r, 1084 fParams->fCenter.h)); 1085 1086 fNormRadius = MaxDistancePointToRect (squareCenter, 1087 squareBounds); 1088 1089 fInvNormRadius = 1.0 / fNormRadius; 1090 1091 } 1092 1093 // Propagate warp params to other planes. 1094 1095 fParams->PropagateToAllPlanes (fDstPlanes); 1096 1097 } 1098 1099 /*****************************************************************************/ 1100 1101 void dng_filter_warp::Initialize (dng_host &host) 1102 { 1103 1104 // Make resample weights. 1105 1106 const dng_resample_function &kernel = dng_resample_bicubic::Get (); 1107 1108 fWeights.Initialize (kernel, 1109 host.Allocator ()); 1110 1111 } 1112 1113 /*****************************************************************************/ 1114 1115 dng_rect dng_filter_warp::SrcArea (const dng_rect &dstArea) 1116 { 1117 1118 // Walk each pixel of the boundary of dstArea, map it to the uncorrected src 1119 // pixel position, and return the rectangle that contains all such src pixels. 1120 1121 int32 xMin = INT_MAX; 1122 int32 xMax = INT_MIN; 1123 int32 yMin = INT_MAX; 1124 int32 yMax = INT_MIN; 1125 1126 for (uint32 plane = 0; plane < fDstPlanes; plane++) 1127 { 1128 1129 // Top and bottom edges. 1130 1131 for (int32 c = dstArea.l; c < dstArea.r; c++) 1132 { 1133 1134 // Top edge. 1135 1136 { 1137 1138 const dng_point_real64 dst (dstArea.t, c); 1139 1140 const dng_point_real64 src = GetSrcPixelPosition (dst, plane); 1141 1142 const int32 y = ConvertDoubleToInt32 (floor (src.v)); 1143 1144 yMin = Min_int32 (yMin, y); 1145 1146 } 1147 1148 // Bottom edge. 1149 1150 { 1151 1152 const dng_point_real64 dst (dstArea.b - 1, c); 1153 1154 const dng_point_real64 src = GetSrcPixelPosition (dst, plane); 1155 1156 const int32 y = ConvertDoubleToInt32 (ceil (src.v)); 1157 1158 yMax = Max_int32 (yMax, y); 1159 1160 } 1161 1162 } 1163 1164 // Left and right edges. 1165 1166 for (int32 r = dstArea.t; r < dstArea.b; r++) 1167 { 1168 1169 // Left edge. 1170 1171 { 1172 1173 const dng_point_real64 dst (r, dstArea.l); 1174 1175 const dng_point_real64 src = GetSrcPixelPosition (dst, plane); 1176 1177 const int32 x = ConvertDoubleToInt32 (floor (src.h)); 1178 1179 xMin = Min_int32 (xMin, x); 1180 1181 } 1182 1183 // Right edge. 1184 1185 { 1186 1187 const dng_point_real64 dst (r, dstArea.r - 1); 1188 1189 const dng_point_real64 src = GetSrcPixelPosition (dst, plane); 1190 1191 const int32 x = ConvertDoubleToInt32 (ceil (src.h)); 1192 1193 xMax = Max_int32 (xMax, x); 1194 1195 } 1196 1197 } 1198 1199 } 1200 1201 // Pad each side by filter radius. 1202 1203 const int32 pad = ConvertUint32ToInt32(fWeights.Radius()); 1204 1205 xMin = SafeInt32Sub(xMin, pad); 1206 yMin = SafeInt32Sub(yMin, pad); 1207 xMax = SafeInt32Add(xMax, pad); 1208 yMax = SafeInt32Add(yMax, pad); 1209 1210 xMax = SafeInt32Add(xMax, 1); 1211 yMax = SafeInt32Add(yMax, 1); 1212 1213 const dng_rect srcArea (yMin, 1214 xMin, 1215 yMax, 1216 xMax); 1217 1218 return srcArea; 1219 1220 } 1221 1222 /*****************************************************************************/ 1223 1224 dng_point dng_filter_warp::SrcTileSize (const dng_point &dstTileSize) 1225 { 1226 1227 // Obtain an upper bound on the source tile size. We'll do this by considering 1228 // upper bounds on the radial and tangential warp components separately, then add 1229 // them together. This is not a tight bound but is good enough because the 1230 // tangential terms are usually quite small. 1231 1232 // Get upper bound on src tile size from radial warp. 1233 1234 DNG_REQUIRE (dstTileSize.v > 0, "Invalid tile height."); 1235 DNG_REQUIRE (dstTileSize.h > 0, "Invalid tile width."); 1236 1237 const real64 maxDstGap = fInvNormRadius * hypot ((real64) dstTileSize.h, 1238 (real64) dstTileSize.v); 1239 1240 dng_point srcTileSize; 1241 1242 if (maxDstGap >= 1.0) 1243 { 1244 1245 // The proposed tile size is unusually large, i.e., its hypotenuse is larger 1246 // than the maximum radius. Bite the bullet and just return a tile size big 1247 // enough to process the whole image. 1248 1249 srcTileSize = SrcArea (fDstImage.Bounds ()).Size (); 1250 1251 } 1252 1253 else 1254 { 1255 1256 // maxDstGap < 1.0. 1257 1258 const real64 maxSrcGap = fParams->MaxSrcRadiusGap (maxDstGap); 1259 1260 const int32 dim = ConvertDoubleToInt32 (ceil (maxSrcGap * fNormRadius)); 1261 1262 srcTileSize = dng_point (dim, dim); 1263 1264 } 1265 1266 srcTileSize.h += ConvertUint32ToInt32(fWeights.Width()); 1267 srcTileSize.v += ConvertUint32ToInt32(fWeights.Width()); 1268 1269 // Get upper bound on src tile size from tangential warp. 1270 1271 const dng_rect_real64 bounds (fSrcImage.Bounds ()); 1272 1273 const dng_point_real64 minDst ((bounds.t - fCenter.v) * fInvNormRadius, 1274 (bounds.l - fCenter.h) * fInvNormRadius); 1275 1276 const dng_point_real64 maxDst ((bounds.b - 1.0 - fCenter.v) * fInvNormRadius, 1277 (bounds.r - 1.0 - fCenter.h) * fInvNormRadius); 1278 1279 const dng_point_real64 srcTanGap = fParams->MaxSrcTanGap (minDst, 1280 maxDst); 1281 1282 // Add the two bounds together. 1283 1284 srcTileSize.v += ConvertDoubleToInt32 (ceil (srcTanGap.v * fNormRadius)); 1285 srcTileSize.h += ConvertDoubleToInt32 (ceil (srcTanGap.h * fNormRadius)); 1286 1287 return srcTileSize; 1288 1289 } 1290 1291 /*****************************************************************************/ 1292 1293 void dng_filter_warp::ProcessArea (uint32 /* threadIndex */, 1294 dng_pixel_buffer &srcBuffer, 1295 dng_pixel_buffer &dstBuffer) 1296 { 1297 1298 // Prepare resample constants. 1299 1300 const int32 wCount = fWeights.Width (); 1301 1302 const dng_point srcOffset (fWeights.Offset (), 1303 fWeights.Offset ()); 1304 1305 const real64 numSubsamples = (real64) kResampleSubsampleCount2D; 1306 1307 // Prepare area and step constants. 1308 1309 const dng_rect srcArea = srcBuffer.fArea; 1310 const dng_rect dstArea = dstBuffer.fArea; 1311 1312 const int32 srcRowStep = (int32) srcBuffer.RowStep (); 1313 1314 const int32 hMin = srcArea.l; 1315 const int32 hMax = SafeInt32Sub (SafeInt32Sub (srcArea.r, wCount), 1); 1316 1317 const int32 vMin = srcArea.t; 1318 const int32 vMax = SafeInt32Sub (SafeInt32Sub (srcArea.b, wCount), 1); 1319 1320 if (hMax < hMin || vMax < vMin) 1321 { 1322 1323 ThrowBadFormat ("Empty source area in dng_filter_warp."); 1324 1325 } 1326 1327 // Warp each plane. 1328 1329 for (uint32 plane = 0; plane < dstBuffer.fPlanes; plane++) 1330 { 1331 1332 real32 *dPtr = dstBuffer.DirtyPixel_real32 (dstArea.t, 1333 dstArea.l, 1334 plane); 1335 1336 for (int32 dstRow = dstArea.t; dstRow < dstArea.b; dstRow++) 1337 { 1338 1339 uint32 dstIndex = 0; 1340 1341 for (int32 dstCol = dstArea.l; dstCol < dstArea.r; dstCol++, dstIndex++) 1342 { 1343 1344 // Get destination (corrected) pixel position. 1345 1346 const dng_point_real64 dPos ((real64) dstRow, 1347 (real64) dstCol); 1348 1349 // Warp to source (uncorrected) pixel position. 1350 1351 const dng_point_real64 sPos = GetSrcPixelPosition (dPos, 1352 plane); 1353 1354 // Decompose into integer and fractional parts. 1355 1356 dng_point sInt (ConvertDoubleToInt32 (floor (sPos.v)), 1357 ConvertDoubleToInt32 (floor (sPos.h))); 1358 1359 dng_point sFct (ConvertDoubleToInt32 ((sPos.v - (real64) sInt.v) * numSubsamples), 1360 ConvertDoubleToInt32 ((sPos.h - (real64) sInt.h) * numSubsamples)); 1361 1362 // Add resample offset. 1363 1364 sInt = sInt + srcOffset; 1365 1366 // Clip. 1367 1368 if (sInt.h < hMin) 1369 { 1370 sInt.h = hMin; 1371 sFct.h = 0; 1372 } 1373 1374 else if (sInt.h > hMax) 1375 { 1376 sInt.h = hMax; 1377 sFct.h = 0; 1378 } 1379 1380 if (sInt.v < vMin) 1381 { 1382 sInt.v = vMin; 1383 sFct.v = 0; 1384 } 1385 1386 else if (sInt.v > vMax) 1387 { 1388 sInt.v = vMax; 1389 sFct.v = 0; 1390 } 1391 1392 // Perform 2D resample. 1393 1394 const real32 *w = fWeights.Weights32 (sFct); 1395 1396 const real32 *s = srcBuffer.ConstPixel_real32 (sInt.v, 1397 sInt.h, 1398 plane); 1399 1400 real32 total = 0.0f; 1401 1402 for (int32 i = 0; i < wCount; i++) 1403 { 1404 1405 for (int32 j = 0; j < wCount; j++) 1406 { 1407 1408 total += w [j] * s [j]; 1409 1410 } 1411 1412 w += wCount; 1413 s += srcRowStep; 1414 1415 } 1416 1417 // Store final pixel value. 1418 1419 dPtr [dstIndex] = Pin_real32 (total); 1420 1421 } 1422 1423 // Advance to next row. 1424 1425 dPtr += dstBuffer.RowStep (); 1426 1427 } 1428 1429 } 1430 1431 } 1432 1433 /*****************************************************************************/ 1434 1435 dng_point_real64 dng_filter_warp::GetSrcPixelPosition (const dng_point_real64 &dst, 1436 uint32 plane) 1437 { 1438 1439 const dng_point_real64 diff = dst - fCenter; 1440 1441 const dng_point_real64 diffNorm (diff.v * fInvNormRadius, 1442 diff.h * fInvNormRadius); 1443 1444 const dng_point_real64 diffNormScaled (diffNorm.v * fPixelScaleV, 1445 diffNorm.h); 1446 1447 const dng_point_real64 diffNormSqr (diffNormScaled.v * diffNormScaled.v, 1448 diffNormScaled.h * diffNormScaled.h); 1449 1450 const real64 rr = Min_real64 (diffNormSqr.v + diffNormSqr.h, 1451 1.0); 1452 1453 dng_point_real64 dSrc; 1454 1455 if (fIsTanNOP) 1456 { 1457 1458 // Radial only. 1459 1460 const real64 ratio = fParams->EvaluateRatio (plane, 1461 rr); 1462 1463 dSrc.h = diff.h * ratio; 1464 dSrc.v = diff.v * ratio; 1465 1466 } 1467 1468 else if (fIsRadNOP) 1469 { 1470 1471 // Tangential only. 1472 1473 const dng_point_real64 tan = fParams->EvaluateTangential (plane, 1474 rr, 1475 diffNormScaled, 1476 diffNormSqr); 1477 1478 dSrc.h = diff.h + (fNormRadius * tan.h); 1479 dSrc.v = diff.v + (fNormRadius * tan.v * fPixelScaleVInv); 1480 1481 } 1482 1483 else 1484 { 1485 1486 // Radial and tangential. 1487 1488 const real64 ratio = fParams->EvaluateRatio (plane, 1489 rr); 1490 1491 const dng_point_real64 tan = fParams->EvaluateTangential (plane, 1492 rr, 1493 diffNormScaled, 1494 diffNormSqr); 1495 1496 dSrc.h = fNormRadius * (diffNorm.h * ratio + tan.h); 1497 dSrc.v = fNormRadius * (diffNorm.v * ratio + tan.v * fPixelScaleVInv); 1498 1499 } 1500 1501 return fCenter + dSrc; 1502 1503 } 1504 1505 /*****************************************************************************/ 1506 1507 dng_opcode_WarpRectilinear::dng_opcode_WarpRectilinear (const dng_warp_params_rectilinear ¶ms, 1508 uint32 flags) 1509 1510 : dng_opcode (dngOpcode_WarpRectilinear, 1511 dngVersion_1_3_0_0, 1512 flags) 1513 1514 , fWarpParams (params) 1515 1516 { 1517 1518 if (!params.IsValid ()) 1519 { 1520 ThrowBadFormat (); 1521 } 1522 1523 } 1524 1525 /*****************************************************************************/ 1526 1527 dng_opcode_WarpRectilinear::dng_opcode_WarpRectilinear (dng_stream &stream) 1528 1529 : dng_opcode (dngOpcode_WarpRectilinear, 1530 stream, 1531 "WarpRectilinear") 1532 1533 , fWarpParams () 1534 1535 { 1536 1537 // Grab the size in bytes. 1538 1539 const uint32 bytes = stream.Get_uint32 (); 1540 1541 // Grab the number of planes to warp. 1542 1543 fWarpParams.fPlanes = stream.Get_uint32 (); 1544 1545 // Verify number of planes. 1546 1547 if (fWarpParams.fPlanes == 0 || 1548 fWarpParams.fPlanes > kMaxColorPlanes) 1549 { 1550 ThrowBadFormat (); 1551 } 1552 1553 // Verify the size. 1554 1555 if (bytes != ParamBytes (fWarpParams.fPlanes)) 1556 { 1557 ThrowBadFormat (); 1558 } 1559 1560 // Read warp parameters for each plane. 1561 1562 for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) 1563 { 1564 1565 fWarpParams.fRadParams [plane][0] = stream.Get_real64 (); 1566 fWarpParams.fRadParams [plane][1] = stream.Get_real64 (); 1567 fWarpParams.fRadParams [plane][2] = stream.Get_real64 (); 1568 fWarpParams.fRadParams [plane][3] = stream.Get_real64 (); 1569 1570 fWarpParams.fTanParams [plane][0] = stream.Get_real64 (); 1571 fWarpParams.fTanParams [plane][1] = stream.Get_real64 (); 1572 1573 } 1574 1575 // Read the image center. 1576 1577 fWarpParams.fCenter.h = stream.Get_real64 (); 1578 fWarpParams.fCenter.v = stream.Get_real64 (); 1579 1580 #if qDNGValidate 1581 1582 if (gVerbose) 1583 { 1584 1585 fWarpParams.Dump (); 1586 1587 } 1588 1589 #endif 1590 1591 if (!fWarpParams.IsValid ()) 1592 { 1593 ThrowBadFormat (); 1594 } 1595 1596 } 1597 1598 /*****************************************************************************/ 1599 1600 bool dng_opcode_WarpRectilinear::IsNOP () const 1601 { 1602 1603 return fWarpParams.IsNOPAll (); 1604 1605 } 1606 1607 /*****************************************************************************/ 1608 1609 bool dng_opcode_WarpRectilinear::IsValidForNegative (const dng_negative &negative) const 1610 { 1611 1612 return fWarpParams.IsValidForNegative (negative); 1613 1614 } 1615 1616 /*****************************************************************************/ 1617 1618 void dng_opcode_WarpRectilinear::PutData (dng_stream &stream) const 1619 { 1620 1621 const uint32 bytes = ParamBytes (fWarpParams.fPlanes); 1622 1623 stream.Put_uint32 (bytes); 1624 1625 stream.Put_uint32 (fWarpParams.fPlanes); 1626 1627 for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) 1628 { 1629 1630 stream.Put_real64 (fWarpParams.fRadParams [plane][0]); 1631 stream.Put_real64 (fWarpParams.fRadParams [plane][1]); 1632 stream.Put_real64 (fWarpParams.fRadParams [plane][2]); 1633 stream.Put_real64 (fWarpParams.fRadParams [plane][3]); 1634 1635 stream.Put_real64 (fWarpParams.fTanParams [plane][0]); 1636 stream.Put_real64 (fWarpParams.fTanParams [plane][1]); 1637 1638 } 1639 1640 stream.Put_real64 (fWarpParams.fCenter.h); 1641 stream.Put_real64 (fWarpParams.fCenter.v); 1642 1643 } 1644 1645 /*****************************************************************************/ 1646 1647 void dng_opcode_WarpRectilinear::Apply (dng_host &host, 1648 dng_negative &negative, 1649 AutoPtr<dng_image> &image) 1650 { 1651 1652 #if qDNGValidate 1653 1654 dng_timer timer ("WarpRectilinear time"); 1655 1656 #endif 1657 1658 // Allocate destination image. 1659 1660 AutoPtr<dng_image> dstImage (host.Make_dng_image (image->Bounds (), 1661 image->Planes (), 1662 image->PixelType ())); 1663 1664 // Warp the image. 1665 1666 AutoPtr<dng_warp_params> params (new dng_warp_params_rectilinear (fWarpParams)); 1667 1668 dng_filter_warp filter (*image, 1669 *dstImage, 1670 negative, 1671 params); 1672 1673 filter.Initialize (host); 1674 1675 host.PerformAreaTask (filter, 1676 image->Bounds ()); 1677 1678 // Return the new image. 1679 1680 image.Reset (dstImage.Release ()); 1681 1682 } 1683 1684 /*****************************************************************************/ 1685 1686 uint32 dng_opcode_WarpRectilinear::ParamBytes (uint32 planes) 1687 { 1688 1689 return (1 * (uint32) sizeof (uint32) ) + // Number of planes. 1690 (6 * (uint32) sizeof (real64) * planes) + // Warp coefficients. 1691 (2 * (uint32) sizeof (real64) ); // Optical center. 1692 1693 } 1694 1695 /*****************************************************************************/ 1696 1697 dng_opcode_WarpFisheye::dng_opcode_WarpFisheye (const dng_warp_params_fisheye ¶ms, 1698 uint32 flags) 1699 1700 : dng_opcode (dngOpcode_WarpFisheye, 1701 dngVersion_1_3_0_0, 1702 flags) 1703 1704 , fWarpParams (params) 1705 1706 { 1707 1708 if (!params.IsValid ()) 1709 { 1710 ThrowBadFormat (); 1711 } 1712 1713 } 1714 1715 /*****************************************************************************/ 1716 1717 dng_opcode_WarpFisheye::dng_opcode_WarpFisheye (dng_stream &stream) 1718 1719 : dng_opcode (dngOpcode_WarpFisheye, 1720 stream, 1721 "WarpFisheye") 1722 1723 , fWarpParams () 1724 1725 { 1726 1727 // Grab the size in bytes. 1728 1729 const uint32 bytes = stream.Get_uint32 (); 1730 1731 // Grab the number of planes to warp. 1732 1733 fWarpParams.fPlanes = stream.Get_uint32 (); 1734 1735 // Verify number of planes. 1736 1737 if (fWarpParams.fPlanes == 0 || 1738 fWarpParams.fPlanes > kMaxColorPlanes) 1739 { 1740 ThrowBadFormat (); 1741 } 1742 1743 // Verify the size. 1744 1745 if (bytes != ParamBytes (fWarpParams.fPlanes)) 1746 { 1747 ThrowBadFormat (); 1748 } 1749 1750 // Read warp parameters for each plane. 1751 1752 for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) 1753 { 1754 1755 fWarpParams.fRadParams [plane][0] = stream.Get_real64 (); 1756 fWarpParams.fRadParams [plane][1] = stream.Get_real64 (); 1757 fWarpParams.fRadParams [plane][2] = stream.Get_real64 (); 1758 fWarpParams.fRadParams [plane][3] = stream.Get_real64 (); 1759 1760 } 1761 1762 // Read the image center. 1763 1764 fWarpParams.fCenter.h = stream.Get_real64 (); 1765 fWarpParams.fCenter.v = stream.Get_real64 (); 1766 1767 #if qDNGValidate 1768 1769 if (gVerbose) 1770 { 1771 1772 fWarpParams.Dump (); 1773 1774 } 1775 1776 #endif 1777 1778 if (!fWarpParams.IsValid ()) 1779 { 1780 ThrowBadFormat (); 1781 } 1782 1783 } 1784 1785 /*****************************************************************************/ 1786 1787 bool dng_opcode_WarpFisheye::IsNOP () const 1788 { 1789 1790 return fWarpParams.IsNOPAll (); 1791 1792 } 1793 1794 /*****************************************************************************/ 1795 1796 bool dng_opcode_WarpFisheye::IsValidForNegative (const dng_negative &negative) const 1797 { 1798 1799 return fWarpParams.IsValidForNegative (negative); 1800 1801 } 1802 1803 /*****************************************************************************/ 1804 1805 void dng_opcode_WarpFisheye::PutData (dng_stream &stream) const 1806 { 1807 1808 const uint32 bytes = ParamBytes (fWarpParams.fPlanes); 1809 1810 stream.Put_uint32 (bytes); 1811 1812 // Write the number of planes. 1813 1814 stream.Put_uint32 (fWarpParams.fPlanes); 1815 1816 // Write the warp coefficients. 1817 1818 for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) 1819 { 1820 1821 stream.Put_real64 (fWarpParams.fRadParams [plane][0]); 1822 stream.Put_real64 (fWarpParams.fRadParams [plane][1]); 1823 stream.Put_real64 (fWarpParams.fRadParams [plane][2]); 1824 stream.Put_real64 (fWarpParams.fRadParams [plane][3]); 1825 1826 } 1827 1828 // Write the optical center. 1829 1830 stream.Put_real64 (fWarpParams.fCenter.h); 1831 stream.Put_real64 (fWarpParams.fCenter.v); 1832 1833 } 1834 1835 /*****************************************************************************/ 1836 1837 void dng_opcode_WarpFisheye::Apply (dng_host &host, 1838 dng_negative &negative, 1839 AutoPtr<dng_image> &image) 1840 { 1841 1842 #if qDNGValidate 1843 1844 dng_timer timer ("WarpFisheye time"); 1845 1846 #endif 1847 1848 // Allocate destination image. 1849 1850 AutoPtr<dng_image> dstImage (host.Make_dng_image (image->Bounds (), 1851 image->Planes (), 1852 image->PixelType ())); 1853 1854 // Warp the image. 1855 1856 AutoPtr<dng_warp_params> params (new dng_warp_params_fisheye (fWarpParams)); 1857 1858 dng_filter_warp filter (*image, 1859 *dstImage, 1860 negative, 1861 params); 1862 1863 filter.Initialize (host); 1864 1865 host.PerformAreaTask (filter, 1866 image->Bounds ()); 1867 1868 // Return the new image. 1869 1870 image.Reset (dstImage.Release ()); 1871 1872 } 1873 1874 /*****************************************************************************/ 1875 1876 uint32 dng_opcode_WarpFisheye::ParamBytes (uint32 planes) 1877 { 1878 1879 return (1 * (uint32) sizeof (uint32) ) + // Number of planes. 1880 (4 * (uint32) sizeof (real64) * planes) + // Warp coefficients. 1881 (2 * (uint32) sizeof (real64) ); // Optical center. 1882 1883 } 1884 1885 /*****************************************************************************/ 1886 1887 dng_vignette_radial_params::dng_vignette_radial_params () 1888 1889 : fParams (kNumTerms) 1890 , fCenter (0.5, 0.5) 1891 1892 { 1893 1894 } 1895 1896 /*****************************************************************************/ 1897 1898 dng_vignette_radial_params::dng_vignette_radial_params (const dng_std_vector<real64> ¶ms, 1899 const dng_point_real64 ¢er) 1900 1901 : fParams (params) 1902 , fCenter (center) 1903 1904 { 1905 1906 } 1907 1908 /*****************************************************************************/ 1909 1910 bool dng_vignette_radial_params::IsNOP () const 1911 { 1912 1913 for (uint32 i = 0; i < fParams.size (); i++) 1914 { 1915 1916 if (fParams [i] != 0.0) 1917 { 1918 return false; 1919 } 1920 1921 } 1922 1923 return true; 1924 1925 } 1926 1927 /*****************************************************************************/ 1928 1929 bool dng_vignette_radial_params::IsValid () const 1930 { 1931 1932 if (fParams.size () != kNumTerms) 1933 { 1934 return false; 1935 } 1936 1937 if (fCenter.h < 0.0 || 1938 fCenter.h > 1.0 || 1939 fCenter.v < 0.0 || 1940 fCenter.v > 1.0) 1941 { 1942 return false; 1943 } 1944 1945 return true; 1946 1947 } 1948 1949 /*****************************************************************************/ 1950 1951 void dng_vignette_radial_params::Dump () const 1952 { 1953 1954 #if qDNGValidate 1955 1956 printf (" Radial vignette params: "); 1957 1958 for (uint32 i = 0; i < fParams.size (); i++) 1959 { 1960 1961 printf ("%s%.6lf", 1962 (i == 0) ? "" : ", ", 1963 fParams [i]); 1964 1965 } 1966 1967 printf ("\n"); 1968 1969 printf (" Optical center:\n" 1970 " h = %.6lf\n" 1971 " v = %.6lf\n", 1972 fCenter.h, 1973 fCenter.v); 1974 1975 #endif 1976 1977 } 1978 1979 /*****************************************************************************/ 1980 1981 class dng_vignette_radial_function: public dng_1d_function 1982 { 1983 1984 protected: 1985 1986 const dng_vignette_radial_params fParams; 1987 1988 public: 1989 1990 explicit dng_vignette_radial_function (const dng_vignette_radial_params ¶ms) 1991 1992 : fParams (params) 1993 1994 { 1995 1996 } 1997 1998 // x represents r^2, where r is the normalized Euclidean distance from the 1999 // optical center to a pixel. r lies in [0,1], so r^2 (and hence x) also lies 2000 // in [0,1]. 2001 2002 virtual real64 Evaluate (real64 x) const 2003 { 2004 2005 DNG_REQUIRE (fParams.fParams.size () == dng_vignette_radial_params::kNumTerms, 2006 "Bad number of vignette opcode coefficients."); 2007 2008 real64 sum = 0.0; 2009 2010 const dng_std_vector<real64> &v = fParams.fParams; 2011 2012 for (dng_std_vector<real64>::const_reverse_iterator i = v.rbegin (); i != v.rend (); i++) 2013 { 2014 sum = x * ((*i) + sum); 2015 } 2016 2017 sum += 1.0; 2018 2019 return sum; 2020 2021 } 2022 2023 }; 2024 2025 /*****************************************************************************/ 2026 2027 dng_opcode_FixVignetteRadial::dng_opcode_FixVignetteRadial (const dng_vignette_radial_params ¶ms, 2028 uint32 flags) 2029 2030 : dng_inplace_opcode (dngOpcode_FixVignetteRadial, 2031 dngVersion_1_3_0_0, 2032 flags) 2033 2034 , fParams (params) 2035 2036 , fImagePlanes (1) 2037 2038 , fSrcOriginH (0) 2039 , fSrcOriginV (0) 2040 2041 , fSrcStepH (0) 2042 , fSrcStepV (0) 2043 2044 , fTableInputBits (0) 2045 , fTableOutputBits (0) 2046 2047 , fGainTable () 2048 2049 { 2050 2051 if (!params.IsValid ()) 2052 { 2053 ThrowBadFormat (); 2054 } 2055 2056 } 2057 2058 /*****************************************************************************/ 2059 2060 dng_opcode_FixVignetteRadial::dng_opcode_FixVignetteRadial (dng_stream &stream) 2061 2062 : dng_inplace_opcode (dngOpcode_FixVignetteRadial, 2063 stream, 2064 "FixVignetteRadial") 2065 2066 , fParams () 2067 2068 , fImagePlanes (1) 2069 2070 , fSrcOriginH (0) 2071 , fSrcOriginV (0) 2072 2073 , fSrcStepH (0) 2074 , fSrcStepV (0) 2075 2076 , fTableInputBits (0) 2077 , fTableOutputBits (0) 2078 2079 , fGainTable () 2080 2081 { 2082 2083 // Grab the size in bytes. 2084 2085 const uint32 bytes = stream.Get_uint32 (); 2086 2087 // Verify the size. 2088 2089 if (bytes != ParamBytes ()) 2090 { 2091 ThrowBadFormat (); 2092 } 2093 2094 // Read vignette coefficients. 2095 2096 fParams.fParams = dng_std_vector<real64> (dng_vignette_radial_params::kNumTerms); 2097 2098 for (uint32 i = 0; i < dng_vignette_radial_params::kNumTerms; i++) 2099 { 2100 fParams.fParams [i] = stream.Get_real64 (); 2101 } 2102 2103 // Read the image center. 2104 2105 fParams.fCenter.h = stream.Get_real64 (); 2106 fParams.fCenter.v = stream.Get_real64 (); 2107 2108 // Debug. 2109 2110 #if qDNGValidate 2111 2112 if (gVerbose) 2113 { 2114 2115 fParams.Dump (); 2116 2117 } 2118 2119 #endif 2120 2121 if (!fParams.IsValid ()) 2122 { 2123 ThrowBadFormat (); 2124 } 2125 2126 } 2127 2128 /*****************************************************************************/ 2129 2130 bool dng_opcode_FixVignetteRadial::IsNOP () const 2131 { 2132 2133 return fParams.IsNOP (); 2134 2135 } 2136 2137 /*****************************************************************************/ 2138 2139 bool dng_opcode_FixVignetteRadial::IsValidForNegative (const dng_negative & /* negative */) const 2140 { 2141 2142 return fParams.IsValid (); 2143 2144 } 2145 2146 /*****************************************************************************/ 2147 2148 void dng_opcode_FixVignetteRadial::PutData (dng_stream &stream) const 2149 { 2150 2151 const uint32 bytes = ParamBytes (); 2152 2153 stream.Put_uint32 (bytes); 2154 2155 DNG_REQUIRE (fParams.fParams.size () == dng_vignette_radial_params::kNumTerms, 2156 "Bad number of vignette opcode coefficients."); 2157 2158 for (uint32 i = 0; i < dng_vignette_radial_params::kNumTerms; i++) 2159 { 2160 stream.Put_real64 (fParams.fParams [i]); 2161 } 2162 2163 stream.Put_real64 (fParams.fCenter.h); 2164 stream.Put_real64 (fParams.fCenter.v); 2165 2166 } 2167 2168 /*****************************************************************************/ 2169 2170 void dng_opcode_FixVignetteRadial::Prepare (dng_negative &negative, 2171 uint32 threadCount, 2172 const dng_point &tileSize, 2173 const dng_rect &imageBounds, 2174 uint32 imagePlanes, 2175 uint32 bufferPixelType, 2176 dng_memory_allocator &allocator) 2177 { 2178 2179 // This opcode is restricted to 32-bit images. 2180 2181 if (bufferPixelType != ttFloat) 2182 { 2183 ThrowBadFormat (); 2184 } 2185 2186 // Sanity check number of planes. 2187 2188 DNG_ASSERT (imagePlanes >= 1 && imagePlanes <= kMaxColorPlanes, 2189 "Bad number of planes."); 2190 2191 if (imagePlanes < 1 || imagePlanes > kMaxColorPlanes) 2192 { 2193 ThrowProgramError (); 2194 } 2195 2196 fImagePlanes = imagePlanes; 2197 2198 // Set the vignette correction curve. 2199 2200 const dng_vignette_radial_function curve (fParams); 2201 2202 // Grab the destination image area. 2203 2204 const dng_rect_real64 bounds (imageBounds); 2205 2206 // Determine the optical center and maximum radius in pixel coordinates. 2207 2208 const dng_point_real64 centerPixel (Lerp_real64 (bounds.t, 2209 bounds.b, 2210 fParams.fCenter.v), 2211 2212 Lerp_real64 (bounds.l, 2213 bounds.r, 2214 fParams.fCenter.h)); 2215 2216 const real64 pixelScaleV = 1.0 / negative.PixelAspectRatio (); 2217 2218 const real64 maxRadius = hypot (Max_real64 (Abs_real64 (centerPixel.v - bounds.t), 2219 Abs_real64 (centerPixel.v - bounds.b)) * pixelScaleV, 2220 2221 Max_real64 (Abs_real64 (centerPixel.h - bounds.l), 2222 Abs_real64 (centerPixel.h - bounds.r))); 2223 2224 const dng_point_real64 radius (maxRadius, 2225 maxRadius); 2226 2227 // Find origin and scale. 2228 2229 const real64 pixelScaleH = 1.0; 2230 2231 fSrcOriginH = Real64ToFixed64 (-centerPixel.h * pixelScaleH / radius.h); 2232 fSrcOriginV = Real64ToFixed64 (-centerPixel.v * pixelScaleV / radius.v); 2233 2234 fSrcStepH = Real64ToFixed64 (pixelScaleH / radius.h); 2235 fSrcStepV = Real64ToFixed64 (pixelScaleV / radius.v); 2236 2237 // Adjust for pixel centers. 2238 2239 fSrcOriginH += fSrcStepH >> 1; 2240 fSrcOriginV += fSrcStepV >> 1; 2241 2242 // Evaluate 32-bit vignette correction table. 2243 2244 dng_1d_table table32; 2245 2246 table32.Initialize (allocator, 2247 curve, 2248 false); 2249 2250 // Find maximum scale factor. 2251 2252 const real64 maxScale = Max_real32 (table32.Interpolate (0.0f), 2253 table32.Interpolate (1.0f)); 2254 2255 // Find table input bits. 2256 2257 fTableInputBits = 16; 2258 2259 // Find table output bits. 2260 2261 fTableOutputBits = 15; 2262 2263 while ((1 << fTableOutputBits) * maxScale > 65535.0) 2264 { 2265 fTableOutputBits--; 2266 } 2267 2268 // Allocate 16-bit scale table. 2269 2270 const uint32 tableEntries = (1 << fTableInputBits) + 1; 2271 2272 fGainTable.Reset (allocator.Allocate (tableEntries * (uint32) sizeof (uint16))); 2273 2274 uint16 *table16 = fGainTable->Buffer_uint16 (); 2275 2276 // Interpolate 32-bit table into 16-bit table. 2277 2278 const real32 scale0 = 1.0f / (1 << fTableInputBits ); 2279 const real32 scale1 = 1.0f * (1 << fTableOutputBits); 2280 2281 for (uint32 index = 0; index < tableEntries; index++) 2282 { 2283 2284 real32 x = index * scale0; 2285 2286 real32 y = table32.Interpolate (x) * scale1; 2287 2288 table16 [index] = (uint16) Round_uint32 (y); 2289 2290 } 2291 2292 // Prepare vignette mask buffers. 2293 2294 { 2295 2296 const uint32 pixelType = ttShort; 2297 const uint32 bufferSize = ComputeBufferSize(pixelType, tileSize, 2298 imagePlanes, pad16Bytes); 2299 2300 for (uint32 threadIndex = 0; threadIndex < threadCount; threadIndex++) 2301 { 2302 2303 fMaskBuffers [threadIndex] . Reset (allocator.Allocate (bufferSize)); 2304 2305 } 2306 2307 } 2308 2309 } 2310 2311 /*****************************************************************************/ 2312 2313 void dng_opcode_FixVignetteRadial::ProcessArea (dng_negative & /* negative */, 2314 uint32 threadIndex, 2315 dng_pixel_buffer &buffer, 2316 const dng_rect &dstArea, 2317 const dng_rect & /* imageBounds */) 2318 { 2319 2320 // Setup mask pixel buffer. 2321 2322 dng_pixel_buffer maskPixelBuffer(dstArea, 0, fImagePlanes, ttShort, 2323 pcRowInterleavedAlign16, 2324 fMaskBuffers [threadIndex]->Buffer ()); 2325 2326 // Compute mask. 2327 2328 DoVignetteMask16 (maskPixelBuffer.DirtyPixel_uint16 (dstArea.t, dstArea.l), 2329 dstArea.H (), 2330 dstArea.W (), 2331 maskPixelBuffer.RowStep (), 2332 fSrcOriginH + fSrcStepH * dstArea.l, 2333 fSrcOriginV + fSrcStepV * dstArea.t, 2334 fSrcStepH, 2335 fSrcStepV, 2336 fTableInputBits, 2337 fGainTable->Buffer_uint16 ()); 2338 2339 // Apply mask. 2340 2341 DoVignette32 (buffer.DirtyPixel_real32 (dstArea.t, dstArea.l), 2342 maskPixelBuffer.ConstPixel_uint16 (dstArea.t, dstArea.l), 2343 dstArea.H (), 2344 dstArea.W (), 2345 fImagePlanes, 2346 buffer.RowStep (), 2347 buffer.PlaneStep (), 2348 maskPixelBuffer.RowStep (), 2349 fTableOutputBits); 2350 2351 } 2352 2353 /*****************************************************************************/ 2354 2355 uint32 dng_opcode_FixVignetteRadial::ParamBytes () 2356 { 2357 2358 const uint32 N = dng_vignette_radial_params::kNumTerms; 2359 2360 return ((N * sizeof (real64)) + // Vignette coefficients. 2361 (2 * sizeof (real64))); // Optical center. 2362 2363 } 2364 2365 /*****************************************************************************/ 2366
