/*****************************************************************************/ // Copyright 2008 Adobe Systems Incorporated // All Rights Reserved. // // NOTICE: Adobe permits you to use, modify, and distribute this file in // accordance with the terms of the Adobe license agreement accompanying it. /*****************************************************************************/ /* $Id: //mondo/dng_sdk_1_4/dng_sdk/source/dng_lens_correction.cpp#1 $ */ /* $DateTime: 2012/05/30 13:28:51 $ */ /* $Change: 832332 $ */ /* $Author: tknoll $ */ /*****************************************************************************/ #include <cfloat> #include <limits.h> #include "dng_1d_table.h" #include "dng_assertions.h" #include "dng_bottlenecks.h" #include "dng_exceptions.h" #include "dng_filter_task.h" #include "dng_globals.h" #include "dng_host.h" #include "dng_image.h" #include "dng_lens_correction.h" #include "dng_negative.h" #include "dng_safe_arithmetic.h" #include "dng_sdk_limits.h" #include "dng_tag_values.h" #include "dng_utils.h" /*****************************************************************************/ dng_warp_params::dng_warp_params () : fPlanes (1) , fCenter (0.5, 0.5) { } /*****************************************************************************/ dng_warp_params::dng_warp_params (uint32 planes, const dng_point_real64 ¢er) : fPlanes (planes) , fCenter (center) { DNG_ASSERT (planes >= 1, "Too few planes." ); DNG_ASSERT (planes <= kMaxColorPlanes, "Too many planes."); DNG_ASSERT (fCenter.h >= 0.0 && fCenter.h <= 1.0, "Center (horizontal) out of range."); DNG_ASSERT (fCenter.v >= 0.0 && fCenter.v <= 1.0, "Center (vertical) out of range."); } /*****************************************************************************/ dng_warp_params::~dng_warp_params () { } /*****************************************************************************/ bool dng_warp_params::IsNOPAll () const { return IsRadNOPAll () && IsTanNOPAll (); } /*****************************************************************************/ bool dng_warp_params::IsNOP (uint32 plane) const { return IsRadNOP (plane) && IsTanNOP (plane); } /*****************************************************************************/ bool dng_warp_params::IsRadNOPAll () const { for (uint32 plane = 0; plane < fPlanes; plane++) { if (!IsRadNOP (plane)) { return false; } } return true; } /*****************************************************************************/ bool dng_warp_params::IsRadNOP (uint32 /* plane */) const { return false; } /*****************************************************************************/ bool dng_warp_params::IsTanNOPAll () const { for (uint32 plane = 0; plane < fPlanes; plane++) { if (!IsTanNOP (plane)) { return false; } } return true; } /*****************************************************************************/ bool dng_warp_params::IsTanNOP (uint32 /* plane */) const { return false; } /*****************************************************************************/ bool dng_warp_params::IsValid () const { if (fPlanes < 1 || fPlanes > kMaxColorPlanes) { return false; } if (fCenter.h < 0.0 || fCenter.h > 1.0 || fCenter.v < 0.0 || fCenter.v > 1.0) { return false; } return true; } /*****************************************************************************/ bool dng_warp_params::IsValidForNegative (const dng_negative &negative) const { if (!IsValid ()) { return false; } if ((fPlanes != 1) && (fPlanes != negative.ColorChannels ())) { return false; } return true; } /*****************************************************************************/ real64 dng_warp_params::EvaluateInverse (uint32 plane, real64 y) const { const uint32 kMaxIterations = 30; const real64 kNearZero = 1.0e-10; real64 x0 = 0.0; real64 y0 = Evaluate (plane, x0); real64 x1 = 1.0; real64 y1 = Evaluate (plane, x1); for (uint32 iteration = 0; iteration < kMaxIterations; iteration++) { if (Abs_real64 (y1 - y0) < kNearZero) { break; } const real64 x2 = Pin_real64 (0.0, x1 + (y - y1) * (x1 - x0) / (y1 - y0), 1.0); const real64 y2 = Evaluate (plane, x2); x0 = x1; y0 = y1; x1 = x2; y1 = y2; } return x1; } /*****************************************************************************/ dng_point_real64 dng_warp_params::EvaluateTangential2 (uint32 plane, const dng_point_real64 &diff) const { const real64 dvdv = diff.v * diff.v; const real64 dhdh = diff.h * diff.h; const real64 rr = dvdv + dhdh; dng_point_real64 diffSqr (dvdv, dhdh); return EvaluateTangential (plane, rr, diff, diffSqr); } /*****************************************************************************/ dng_point_real64 dng_warp_params::EvaluateTangential3 (uint32 plane, real64 r2, const dng_point_real64 &diff) const { dng_point_real64 diffSqr (diff.v * diff.v, diff.h * diff.h); return EvaluateTangential (plane, r2, diff, diffSqr); } /*****************************************************************************/ void dng_warp_params::Dump () const { #if qDNGValidate printf ("Planes: %u\n", (unsigned) fPlanes); printf (" Optical center:\n" " h = %.6lf\n" " v = %.6lf\n", fCenter.h, fCenter.v); #endif } /*****************************************************************************/ dng_warp_params_rectilinear::dng_warp_params_rectilinear () : dng_warp_params () { for (uint32 plane = 0; plane < kMaxColorPlanes; plane++) { fRadParams [plane] = dng_vector (4); fTanParams [plane] = dng_vector (2); fRadParams [plane][0] = 1.0; } } /*****************************************************************************/ dng_warp_params_rectilinear::dng_warp_params_rectilinear (uint32 planes, const dng_vector radParams [], const dng_vector tanParams [], const dng_point_real64 ¢er) : dng_warp_params (planes, center) { for (uint32 i = 0; i < fPlanes; i++) { fRadParams [i] = radParams [i]; fTanParams [i] = tanParams [i]; } } /*****************************************************************************/ dng_warp_params_rectilinear::~dng_warp_params_rectilinear () { } /*****************************************************************************/ bool dng_warp_params_rectilinear::IsRadNOP (uint32 plane) const { DNG_ASSERT (plane < fPlanes, "plane out of range."); const dng_vector &r = fRadParams [plane]; return (r [0] == 1.0 && r [1] == 0.0 && r [2] == 0.0 && r [3] == 0.0); } /*****************************************************************************/ bool dng_warp_params_rectilinear::IsTanNOP (uint32 plane) const { DNG_ASSERT (plane < fPlanes, "plane out of range."); const dng_vector &t = fTanParams [plane]; return (t [0] == 0.0 && t [1] == 0.0); } /*****************************************************************************/ bool dng_warp_params_rectilinear::IsValid () const { for (uint32 plane = 0; plane < fPlanes; plane++) { if (fRadParams [plane].Count () != 4) { return false; } if (fTanParams [plane].Count () < 2) { return false; } } return dng_warp_params::IsValid (); } /*****************************************************************************/ void dng_warp_params_rectilinear::PropagateToAllPlanes (uint32 totalPlanes) { for (uint32 plane = fPlanes; plane < totalPlanes; plane++) { fRadParams [plane] = fRadParams [0]; fTanParams [plane] = fTanParams [0]; } } /*****************************************************************************/ real64 dng_warp_params_rectilinear::Evaluate (uint32 plane, real64 x) const { const dng_vector &K = fRadParams [plane]; // Coefficients. const real64 x2 = x * x; return x * (K [0] + x2 * (K [1] + x2 * (K [2] + x2 * K [3]))); } /*****************************************************************************/ real64 dng_warp_params_rectilinear::EvaluateRatio (uint32 plane, real64 r2) const { const dng_vector &K = fRadParams [plane]; // Coefficients. return K [0] + r2 * (K [1] + r2 * (K [2] + r2 * K [3])); } /*****************************************************************************/ dng_point_real64 dng_warp_params_rectilinear::EvaluateTangential (uint32 plane, real64 r2, const dng_point_real64 &diff, const dng_point_real64 &diff2) const { const real64 kt0 = fTanParams [plane][0]; const real64 kt1 = fTanParams [plane][1]; const real64 dh = diff.h; const real64 dv = diff.v; const real64 dhdh = diff2.h; const real64 dvdv = diff2.v; return dng_point_real64 (kt0 * (r2 + 2.0 * dvdv) + (2.0 * kt1 * dh * dv), // v kt1 * (r2 + 2.0 * dhdh) + (2.0 * kt0 * dh * dv)); // h } /*****************************************************************************/ real64 dng_warp_params_rectilinear::MaxSrcRadiusGap (real64 maxDstGap) const { real64 maxSrcGap = 0.0; for (uint32 plane = 0; plane < fPlanes; plane++) { const dng_vector &coefs = fRadParams [plane]; const real64 k3 = coefs [1]; const real64 k5 = coefs [2]; const real64 k7 = coefs [3]; // // Let f (r) be the radius warp function. Consider the function // // gap (r) = f (r + maxDstGap) - f (r) // // We wish to maximize gap (r) over the domain [0, 1 - maxDstGap]. This will // give us a reasonable upper bound on the src tile size, independent of // dstBounds. // // As usual, we maximize gap (r) by examining its critical points, i.e., by // considering the roots of its derivative which lie in the domain [0, 1 - // maxDstGap]. gap' (r) is a 5th-order polynomial. One of its roots is // -maxDstGap / 2, which is negative and hence lies outside the domain of // interest. This leaves 4 other possible roots. We solve for these // analytically below. // real64 roots [4]; uint32 numRoots = 0; if (k7 == 0.0) { if (k5 == 0.0) { // No roots in [0,1]. } else { // k7 is zero, but k5 is non-zero. At most two real roots. const real64 discrim = 25.0 * (-6.0 * k3 * k5 - 5.0 * k5 * maxDstGap * maxDstGap); if (discrim >= 0.0) { // Two real roots. const real64 scale = 0.1 * k5; const real64 offset = -5.0 * maxDstGap * k5; const real64 sDiscrim = sqrt (discrim); roots [numRoots++] = scale * (offset + sDiscrim); roots [numRoots++] = scale * (offset - sDiscrim); } } } else { // k7 is non-zero. Up to 4 real roots. const real64 d = maxDstGap; const real64 d2 = d * d; const real64 d4 = d2 * d2; const real64 discrim = 25.0 * k5 * k5 - 63.0 * k3 * k7 + 35.0 * d2 * k5 * k7 + 49.0 * d4 * k7 * k7; if (discrim >= 0.0) { const real64 sDiscrim = 4.0 * k7 * sqrt (discrim); const real64 offset = -20.0 * k5 * k7 - 35.0 * d2 * k7 * k7; const real64 discrim1 = offset - sDiscrim; const real64 discrim2 = offset + sDiscrim; const real64 scale = sqrt (21.0) / (42.0 * k7); if (discrim1 >= 0.0) { const real64 offset1 = -d * 0.5; const real64 sDiscrim1 = scale * sqrt (discrim1); roots [numRoots++] = offset1 + sDiscrim1; roots [numRoots++] = offset1 - sDiscrim1; } if (discrim2 >= 0.0) { const real64 offset2 = -d * 0.5; const real64 sDiscrim2 = scale * sqrt (discrim2); roots [numRoots++] = offset2 + sDiscrim2; roots [numRoots++] = offset2 - sDiscrim2; } } } real64 planeMaxSrcGap = 0.0; // Examine the endpoints. { // Check left endpoint: f (maxDstGap) - f (0). Remember that f (0) == 0. const real64 gap1 = Evaluate (plane, maxDstGap); planeMaxSrcGap = Max_real64 (planeMaxSrcGap, gap1); // Check right endpoint: f (1) - f (1 - maxDstGap). const real64 gap2 = Evaluate (plane, 1.0) - Evaluate (plane, 1.0 - maxDstGap); planeMaxSrcGap = Max_real64 (planeMaxSrcGap, gap2); } // Examine the roots we found, if any. for (uint32 i = 0; i < numRoots; i++) { const real64 r = roots [i]; if (r > 0.0 && r < 1.0 - maxDstGap) { const real64 gap = Evaluate (plane, r + maxDstGap) - Evaluate (plane, r); planeMaxSrcGap = Max_real64 (planeMaxSrcGap, gap); } } maxSrcGap = Max_real64 (maxSrcGap, planeMaxSrcGap); } return maxSrcGap; } /*****************************************************************************/ dng_point_real64 dng_warp_params_rectilinear::MaxSrcTanGap (dng_point_real64 minDst, dng_point_real64 maxDst) const { const real64 v [] = { minDst.v, maxDst.v, 0.0 }; const real64 h [] = { minDst.h, maxDst.h, 0.0 }; dng_point_real64 maxGap; for (uint32 plane = 0; plane < fPlanes; plane++) { real64 hMin = +FLT_MAX; real64 hMax = -FLT_MAX; real64 vMin = +FLT_MAX; real64 vMax = -FLT_MAX; for (uint32 i = 0; i < 3; i++) { for (uint32 j = 0; j < 3; j++) { dng_point_real64 dstDiff (v [i], h [j]); dng_point_real64 srcDiff = EvaluateTangential2 (plane, dstDiff); hMin = Min_real64 (hMin, srcDiff.h); hMax = Max_real64 (hMax, srcDiff.h); vMin = Min_real64 (vMin, srcDiff.v); vMax = Max_real64 (vMax, srcDiff.v); } } const real64 hGap = hMax - hMin; const real64 vGap = vMax - vMin; maxGap.h = Max_real64 (maxGap.h, hGap); maxGap.v = Max_real64 (maxGap.v, vGap); } return maxGap; } /*****************************************************************************/ void dng_warp_params_rectilinear::Dump () const { #if qDNGValidate dng_warp_params::Dump (); for (uint32 plane = 0; plane < fPlanes; plane++) { printf (" Plane %u:\n", (unsigned) plane); printf (" Radial params: %.6lf, %.6lf, %.6lf, %.6lf\n", fRadParams [plane][0], fRadParams [plane][1], fRadParams [plane][2], fRadParams [plane][3]); printf (" Tangential params: %.6lf, %.6lf\n", fTanParams [plane][0], fTanParams [plane][1]); } #endif } /*****************************************************************************/ dng_warp_params_fisheye::dng_warp_params_fisheye () : dng_warp_params () { for (uint32 plane = 0; plane < kMaxColorPlanes; plane++) { fRadParams [plane] = dng_vector (4); } } /*****************************************************************************/ dng_warp_params_fisheye::dng_warp_params_fisheye (uint32 planes, const dng_vector radParams [], const dng_point_real64 ¢er) : dng_warp_params (planes, center) { for (uint32 i = 0; i < fPlanes; i++) { fRadParams [i] = radParams [i]; } } /*****************************************************************************/ dng_warp_params_fisheye::~dng_warp_params_fisheye () { } /*****************************************************************************/ bool dng_warp_params_fisheye::IsRadNOP (uint32 /* plane */) const { return false; } /*****************************************************************************/ bool dng_warp_params_fisheye::IsTanNOP (uint32 /* plane */) const { return true; } /*****************************************************************************/ bool dng_warp_params_fisheye::IsValid () const { for (uint32 plane = 0; plane < fPlanes; plane++) { if (fRadParams [plane].Count () != 4) { return false; } } return dng_warp_params::IsValid (); } /*****************************************************************************/ void dng_warp_params_fisheye::PropagateToAllPlanes (uint32 totalPlanes) { for (uint32 plane = fPlanes; plane < totalPlanes; plane++) { fRadParams [plane] = fRadParams [0]; } } /*****************************************************************************/ real64 dng_warp_params_fisheye::Evaluate (uint32 plane, real64 r) const { const real64 t = atan (r); const dng_vector &K = fRadParams [plane]; const real64 t2 = t * t; return t * (K [0] + t2 * (K [1] + t2 * (K [2] + t2 * K [3]))); } /*****************************************************************************/ real64 dng_warp_params_fisheye::EvaluateRatio (uint32 plane, real64 rSqr) const { const real64 eps = 1.0e-12; if (rSqr < eps) { // r is very close to zero. return 1.0; } const real64 r = sqrt (rSqr); return Evaluate (plane, r) / r; } /*****************************************************************************/ dng_point_real64 dng_warp_params_fisheye::EvaluateTangential (uint32 /* plane */, real64 /* r2 */, const dng_point_real64 & /* diff */, const dng_point_real64 & /* diff2 */) const { // This fisheye model does not support tangential warping. ThrowProgramError (); return dng_point_real64 (0.0, 0.0); } /*****************************************************************************/ real64 dng_warp_params_fisheye::MaxSrcRadiusGap (real64 maxDstGap) const { // // Let f (r) be the radius warp function. Consider the function // // gap (r) = f (r + maxDstGap) - f (r) // // We wish to maximize gap (r) over the domain [0, 1 - maxDstGap]. This will // give us a reasonable upper bound on the src tile size, independent of // dstBounds. // // Ideally, we'd like to maximize gap (r) by examining its critical points, // i.e., by considering the roots of its derivative which lie in the domain [0, // 1 - maxDstGap]. However, gap' (r) is too complex to find its roots // analytically. // real64 maxSrcGap = 0.0; DNG_REQUIRE (maxDstGap > 0.0, "maxDstGap must be positive."); const real64 kMaxValue = 1.0 - maxDstGap; const uint32 kSteps = 128; const real64 kNorm = kMaxValue / real64 (kSteps - 1); for (uint32 plane = 0; plane < fPlanes; plane++) { for (uint32 i = 0; i < kSteps; i++) { const real64 tt = i * kNorm; const real64 gap = Evaluate (plane, tt + maxDstGap) - Evaluate (plane, tt); maxSrcGap = Max_real64 (maxSrcGap, gap); } } return maxSrcGap; } /*****************************************************************************/ dng_point_real64 dng_warp_params_fisheye::MaxSrcTanGap (dng_point_real64 /* minDst */, dng_point_real64 /* maxDst */) const { // This fisheye model does not support tangential distortion. return dng_point_real64 (0.0, 0.0); } /*****************************************************************************/ void dng_warp_params_fisheye::Dump () const { #if qDNGValidate dng_warp_params::Dump (); for (uint32 plane = 0; plane < fPlanes; plane++) { printf (" Plane %u:\n", (unsigned) plane); printf (" Radial params: %.6lf, %.6lf, %.6lf, %.6lf\n", fRadParams [plane][0], fRadParams [plane][1], fRadParams [plane][2], fRadParams [plane][3]); } #endif } /*****************************************************************************/ class dng_filter_warp: public dng_filter_task { protected: AutoPtr<dng_warp_params> fParams; dng_point_real64 fCenter; dng_resample_weights_2d fWeights; real64 fNormRadius; real64 fInvNormRadius; bool fIsRadNOP; bool fIsTanNOP; const real64 fPixelScaleV; const real64 fPixelScaleVInv; public: dng_filter_warp (const dng_image &srcImage, dng_image &dstImage, const dng_negative &negative, AutoPtr<dng_warp_params> ¶ms); virtual void Initialize (dng_host &host); virtual dng_rect SrcArea (const dng_rect &dstArea); virtual dng_point SrcTileSize (const dng_point &dstTileSize); virtual void ProcessArea (uint32 threadIndex, dng_pixel_buffer &srcBuffer, dng_pixel_buffer &dstBuffer); virtual dng_point_real64 GetSrcPixelPosition (const dng_point_real64 &dst, uint32 plane); }; /*****************************************************************************/ dng_filter_warp::dng_filter_warp (const dng_image &srcImage, dng_image &dstImage, const dng_negative &negative, AutoPtr<dng_warp_params> ¶ms) : dng_filter_task (srcImage, dstImage) , fParams (params.Release ()) , fCenter () , fWeights () , fNormRadius (1.0) , fInvNormRadius (1.0) , fIsRadNOP (false) , fIsTanNOP (false) , fPixelScaleV (1.0 / negative.PixelAspectRatio ()) , fPixelScaleVInv (1.0 / fPixelScaleV) { // Force float processing. fSrcPixelType = ttFloat; fDstPixelType = ttFloat; fIsRadNOP = fParams->IsRadNOPAll (); fIsTanNOP = fParams->IsTanNOPAll (); const uint32 negPlanes = negative.ColorChannels (); DNG_ASSERT (negPlanes >= 1, "Too few planes." ); DNG_ASSERT (negPlanes <= kMaxColorPlanes, "Too many planes."); (void) negPlanes; // At least one set of params must do something interesting. if (fIsRadNOP && fIsTanNOP) { ThrowProgramError (); } // Make sure the warp params are valid for this negative. if (!fParams->IsValidForNegative (negative)) { ThrowBadFormat (); } // Compute center. const dng_rect bounds = srcImage.Bounds (); fCenter.h = Lerp_real64 ((real64) bounds.l, (real64) bounds.r, fParams->fCenter.h); fCenter.v = Lerp_real64 ((real64) bounds.t, (real64) bounds.b, fParams->fCenter.v); // Compute max pixel radius and derive various normalized radius values. Note // that when computing the max pixel radius, we must take into account the pixel // aspect ratio. { dng_rect squareBounds (bounds); squareBounds.b = squareBounds.t + Round_int32 (fPixelScaleV * (real64) squareBounds.H ()); const dng_point_real64 squareCenter (Lerp_real64 ((real64) squareBounds.t, (real64) squareBounds.b, fParams->fCenter.v), Lerp_real64 ((real64) squareBounds.l, (real64) squareBounds.r, fParams->fCenter.h)); fNormRadius = MaxDistancePointToRect (squareCenter, squareBounds); fInvNormRadius = 1.0 / fNormRadius; } // Propagate warp params to other planes. fParams->PropagateToAllPlanes (fDstPlanes); } /*****************************************************************************/ void dng_filter_warp::Initialize (dng_host &host) { // Make resample weights. const dng_resample_function &kernel = dng_resample_bicubic::Get (); fWeights.Initialize (kernel, host.Allocator ()); } /*****************************************************************************/ dng_rect dng_filter_warp::SrcArea (const dng_rect &dstArea) { // Walk each pixel of the boundary of dstArea, map it to the uncorrected src // pixel position, and return the rectangle that contains all such src pixels. int32 xMin = INT_MAX; int32 xMax = INT_MIN; int32 yMin = INT_MAX; int32 yMax = INT_MIN; for (uint32 plane = 0; plane < fDstPlanes; plane++) { // Top and bottom edges. for (int32 c = dstArea.l; c < dstArea.r; c++) { // Top edge. { const dng_point_real64 dst (dstArea.t, c); const dng_point_real64 src = GetSrcPixelPosition (dst, plane); const int32 y = ConvertDoubleToInt32 (floor (src.v)); yMin = Min_int32 (yMin, y); } // Bottom edge. { const dng_point_real64 dst (dstArea.b - 1, c); const dng_point_real64 src = GetSrcPixelPosition (dst, plane); const int32 y = ConvertDoubleToInt32 (ceil (src.v)); yMax = Max_int32 (yMax, y); } } // Left and right edges. for (int32 r = dstArea.t; r < dstArea.b; r++) { // Left edge. { const dng_point_real64 dst (r, dstArea.l); const dng_point_real64 src = GetSrcPixelPosition (dst, plane); const int32 x = ConvertDoubleToInt32 (floor (src.h)); xMin = Min_int32 (xMin, x); } // Right edge. { const dng_point_real64 dst (r, dstArea.r - 1); const dng_point_real64 src = GetSrcPixelPosition (dst, plane); const int32 x = ConvertDoubleToInt32 (ceil (src.h)); xMax = Max_int32 (xMax, x); } } } // Pad each side by filter radius. const int32 pad = ConvertUint32ToInt32(fWeights.Radius()); xMin = SafeInt32Sub(xMin, pad); yMin = SafeInt32Sub(yMin, pad); xMax = SafeInt32Add(xMax, pad); yMax = SafeInt32Add(yMax, pad); xMax = SafeInt32Add(xMax, 1); yMax = SafeInt32Add(yMax, 1); const dng_rect srcArea (yMin, xMin, yMax, xMax); return srcArea; } /*****************************************************************************/ dng_point dng_filter_warp::SrcTileSize (const dng_point &dstTileSize) { // Obtain an upper bound on the source tile size. We'll do this by considering // upper bounds on the radial and tangential warp components separately, then add // them together. This is not a tight bound but is good enough because the // tangential terms are usually quite small. // Get upper bound on src tile size from radial warp. DNG_REQUIRE (dstTileSize.v > 0, "Invalid tile height."); DNG_REQUIRE (dstTileSize.h > 0, "Invalid tile width."); const real64 maxDstGap = fInvNormRadius * hypot ((real64) dstTileSize.h, (real64) dstTileSize.v); dng_point srcTileSize; if (maxDstGap >= 1.0) { // The proposed tile size is unusually large, i.e., its hypotenuse is larger // than the maximum radius. Bite the bullet and just return a tile size big // enough to process the whole image. srcTileSize = SrcArea (fDstImage.Bounds ()).Size (); } else { // maxDstGap < 1.0. const real64 maxSrcGap = fParams->MaxSrcRadiusGap (maxDstGap); const int32 dim = ConvertDoubleToInt32 (ceil (maxSrcGap * fNormRadius)); srcTileSize = dng_point (dim, dim); } srcTileSize.h += ConvertUint32ToInt32(fWeights.Width()); srcTileSize.v += ConvertUint32ToInt32(fWeights.Width()); // Get upper bound on src tile size from tangential warp. const dng_rect_real64 bounds (fSrcImage.Bounds ()); const dng_point_real64 minDst ((bounds.t - fCenter.v) * fInvNormRadius, (bounds.l - fCenter.h) * fInvNormRadius); const dng_point_real64 maxDst ((bounds.b - 1.0 - fCenter.v) * fInvNormRadius, (bounds.r - 1.0 - fCenter.h) * fInvNormRadius); const dng_point_real64 srcTanGap = fParams->MaxSrcTanGap (minDst, maxDst); // Add the two bounds together. srcTileSize.v += ConvertDoubleToInt32 (ceil (srcTanGap.v * fNormRadius)); srcTileSize.h += ConvertDoubleToInt32 (ceil (srcTanGap.h * fNormRadius)); return srcTileSize; } /*****************************************************************************/ void dng_filter_warp::ProcessArea (uint32 /* threadIndex */, dng_pixel_buffer &srcBuffer, dng_pixel_buffer &dstBuffer) { // Prepare resample constants. const int32 wCount = fWeights.Width (); const dng_point srcOffset (fWeights.Offset (), fWeights.Offset ()); const real64 numSubsamples = (real64) kResampleSubsampleCount2D; // Prepare area and step constants. const dng_rect srcArea = srcBuffer.fArea; const dng_rect dstArea = dstBuffer.fArea; const int32 srcRowStep = (int32) srcBuffer.RowStep (); const int32 hMin = srcArea.l; const int32 hMax = SafeInt32Sub (SafeInt32Sub (srcArea.r, wCount), 1); const int32 vMin = srcArea.t; const int32 vMax = SafeInt32Sub (SafeInt32Sub (srcArea.b, wCount), 1); if (hMax < hMin || vMax < vMin) { ThrowBadFormat ("Empty source area in dng_filter_warp."); } // Warp each plane. for (uint32 plane = 0; plane < dstBuffer.fPlanes; plane++) { real32 *dPtr = dstBuffer.DirtyPixel_real32 (dstArea.t, dstArea.l, plane); for (int32 dstRow = dstArea.t; dstRow < dstArea.b; dstRow++) { uint32 dstIndex = 0; for (int32 dstCol = dstArea.l; dstCol < dstArea.r; dstCol++, dstIndex++) { // Get destination (corrected) pixel position. const dng_point_real64 dPos ((real64) dstRow, (real64) dstCol); // Warp to source (uncorrected) pixel position. const dng_point_real64 sPos = GetSrcPixelPosition (dPos, plane); // Decompose into integer and fractional parts. dng_point sInt (ConvertDoubleToInt32 (floor (sPos.v)), ConvertDoubleToInt32 (floor (sPos.h))); dng_point sFct (ConvertDoubleToInt32 ((sPos.v - (real64) sInt.v) * numSubsamples), ConvertDoubleToInt32 ((sPos.h - (real64) sInt.h) * numSubsamples)); // Add resample offset. sInt = sInt + srcOffset; // Clip. if (sInt.h < hMin) { sInt.h = hMin; sFct.h = 0; } else if (sInt.h > hMax) { sInt.h = hMax; sFct.h = 0; } if (sInt.v < vMin) { sInt.v = vMin; sFct.v = 0; } else if (sInt.v > vMax) { sInt.v = vMax; sFct.v = 0; } // Perform 2D resample. const real32 *w = fWeights.Weights32 (sFct); const real32 *s = srcBuffer.ConstPixel_real32 (sInt.v, sInt.h, plane); real32 total = 0.0f; for (int32 i = 0; i < wCount; i++) { for (int32 j = 0; j < wCount; j++) { total += w [j] * s [j]; } w += wCount; s += srcRowStep; } // Store final pixel value. dPtr [dstIndex] = Pin_real32 (total); } // Advance to next row. dPtr += dstBuffer.RowStep (); } } } /*****************************************************************************/ dng_point_real64 dng_filter_warp::GetSrcPixelPosition (const dng_point_real64 &dst, uint32 plane) { const dng_point_real64 diff = dst - fCenter; const dng_point_real64 diffNorm (diff.v * fInvNormRadius, diff.h * fInvNormRadius); const dng_point_real64 diffNormScaled (diffNorm.v * fPixelScaleV, diffNorm.h); const dng_point_real64 diffNormSqr (diffNormScaled.v * diffNormScaled.v, diffNormScaled.h * diffNormScaled.h); const real64 rr = Min_real64 (diffNormSqr.v + diffNormSqr.h, 1.0); dng_point_real64 dSrc; if (fIsTanNOP) { // Radial only. const real64 ratio = fParams->EvaluateRatio (plane, rr); dSrc.h = diff.h * ratio; dSrc.v = diff.v * ratio; } else if (fIsRadNOP) { // Tangential only. const dng_point_real64 tan = fParams->EvaluateTangential (plane, rr, diffNormScaled, diffNormSqr); dSrc.h = diff.h + (fNormRadius * tan.h); dSrc.v = diff.v + (fNormRadius * tan.v * fPixelScaleVInv); } else { // Radial and tangential. const real64 ratio = fParams->EvaluateRatio (plane, rr); const dng_point_real64 tan = fParams->EvaluateTangential (plane, rr, diffNormScaled, diffNormSqr); dSrc.h = fNormRadius * (diffNorm.h * ratio + tan.h); dSrc.v = fNormRadius * (diffNorm.v * ratio + tan.v * fPixelScaleVInv); } return fCenter + dSrc; } /*****************************************************************************/ dng_opcode_WarpRectilinear::dng_opcode_WarpRectilinear (const dng_warp_params_rectilinear ¶ms, uint32 flags) : dng_opcode (dngOpcode_WarpRectilinear, dngVersion_1_3_0_0, flags) , fWarpParams (params) { if (!params.IsValid ()) { ThrowBadFormat (); } } /*****************************************************************************/ dng_opcode_WarpRectilinear::dng_opcode_WarpRectilinear (dng_stream &stream) : dng_opcode (dngOpcode_WarpRectilinear, stream, "WarpRectilinear") , fWarpParams () { // Grab the size in bytes. const uint32 bytes = stream.Get_uint32 (); // Grab the number of planes to warp. fWarpParams.fPlanes = stream.Get_uint32 (); // Verify number of planes. if (fWarpParams.fPlanes == 0 || fWarpParams.fPlanes > kMaxColorPlanes) { ThrowBadFormat (); } // Verify the size. if (bytes != ParamBytes (fWarpParams.fPlanes)) { ThrowBadFormat (); } // Read warp parameters for each plane. for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) { fWarpParams.fRadParams [plane][0] = stream.Get_real64 (); fWarpParams.fRadParams [plane][1] = stream.Get_real64 (); fWarpParams.fRadParams [plane][2] = stream.Get_real64 (); fWarpParams.fRadParams [plane][3] = stream.Get_real64 (); fWarpParams.fTanParams [plane][0] = stream.Get_real64 (); fWarpParams.fTanParams [plane][1] = stream.Get_real64 (); } // Read the image center. fWarpParams.fCenter.h = stream.Get_real64 (); fWarpParams.fCenter.v = stream.Get_real64 (); #if qDNGValidate if (gVerbose) { fWarpParams.Dump (); } #endif if (!fWarpParams.IsValid ()) { ThrowBadFormat (); } } /*****************************************************************************/ bool dng_opcode_WarpRectilinear::IsNOP () const { return fWarpParams.IsNOPAll (); } /*****************************************************************************/ bool dng_opcode_WarpRectilinear::IsValidForNegative (const dng_negative &negative) const { return fWarpParams.IsValidForNegative (negative); } /*****************************************************************************/ void dng_opcode_WarpRectilinear::PutData (dng_stream &stream) const { const uint32 bytes = ParamBytes (fWarpParams.fPlanes); stream.Put_uint32 (bytes); stream.Put_uint32 (fWarpParams.fPlanes); for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) { stream.Put_real64 (fWarpParams.fRadParams [plane][0]); stream.Put_real64 (fWarpParams.fRadParams [plane][1]); stream.Put_real64 (fWarpParams.fRadParams [plane][2]); stream.Put_real64 (fWarpParams.fRadParams [plane][3]); stream.Put_real64 (fWarpParams.fTanParams [plane][0]); stream.Put_real64 (fWarpParams.fTanParams [plane][1]); } stream.Put_real64 (fWarpParams.fCenter.h); stream.Put_real64 (fWarpParams.fCenter.v); } /*****************************************************************************/ void dng_opcode_WarpRectilinear::Apply (dng_host &host, dng_negative &negative, AutoPtr<dng_image> &image) { #if qDNGValidate dng_timer timer ("WarpRectilinear time"); #endif // Allocate destination image. AutoPtr<dng_image> dstImage (host.Make_dng_image (image->Bounds (), image->Planes (), image->PixelType ())); // Warp the image. AutoPtr<dng_warp_params> params (new dng_warp_params_rectilinear (fWarpParams)); dng_filter_warp filter (*image, *dstImage, negative, params); filter.Initialize (host); host.PerformAreaTask (filter, image->Bounds ()); // Return the new image. image.Reset (dstImage.Release ()); } /*****************************************************************************/ uint32 dng_opcode_WarpRectilinear::ParamBytes (uint32 planes) { return (1 * (uint32) sizeof (uint32) ) + // Number of planes. (6 * (uint32) sizeof (real64) * planes) + // Warp coefficients. (2 * (uint32) sizeof (real64) ); // Optical center. } /*****************************************************************************/ dng_opcode_WarpFisheye::dng_opcode_WarpFisheye (const dng_warp_params_fisheye ¶ms, uint32 flags) : dng_opcode (dngOpcode_WarpFisheye, dngVersion_1_3_0_0, flags) , fWarpParams (params) { if (!params.IsValid ()) { ThrowBadFormat (); } } /*****************************************************************************/ dng_opcode_WarpFisheye::dng_opcode_WarpFisheye (dng_stream &stream) : dng_opcode (dngOpcode_WarpFisheye, stream, "WarpFisheye") , fWarpParams () { // Grab the size in bytes. const uint32 bytes = stream.Get_uint32 (); // Grab the number of planes to warp. fWarpParams.fPlanes = stream.Get_uint32 (); // Verify number of planes. if (fWarpParams.fPlanes == 0 || fWarpParams.fPlanes > kMaxColorPlanes) { ThrowBadFormat (); } // Verify the size. if (bytes != ParamBytes (fWarpParams.fPlanes)) { ThrowBadFormat (); } // Read warp parameters for each plane. for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) { fWarpParams.fRadParams [plane][0] = stream.Get_real64 (); fWarpParams.fRadParams [plane][1] = stream.Get_real64 (); fWarpParams.fRadParams [plane][2] = stream.Get_real64 (); fWarpParams.fRadParams [plane][3] = stream.Get_real64 (); } // Read the image center. fWarpParams.fCenter.h = stream.Get_real64 (); fWarpParams.fCenter.v = stream.Get_real64 (); #if qDNGValidate if (gVerbose) { fWarpParams.Dump (); } #endif if (!fWarpParams.IsValid ()) { ThrowBadFormat (); } } /*****************************************************************************/ bool dng_opcode_WarpFisheye::IsNOP () const { return fWarpParams.IsNOPAll (); } /*****************************************************************************/ bool dng_opcode_WarpFisheye::IsValidForNegative (const dng_negative &negative) const { return fWarpParams.IsValidForNegative (negative); } /*****************************************************************************/ void dng_opcode_WarpFisheye::PutData (dng_stream &stream) const { const uint32 bytes = ParamBytes (fWarpParams.fPlanes); stream.Put_uint32 (bytes); // Write the number of planes. stream.Put_uint32 (fWarpParams.fPlanes); // Write the warp coefficients. for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) { stream.Put_real64 (fWarpParams.fRadParams [plane][0]); stream.Put_real64 (fWarpParams.fRadParams [plane][1]); stream.Put_real64 (fWarpParams.fRadParams [plane][2]); stream.Put_real64 (fWarpParams.fRadParams [plane][3]); } // Write the optical center. stream.Put_real64 (fWarpParams.fCenter.h); stream.Put_real64 (fWarpParams.fCenter.v); } /*****************************************************************************/ void dng_opcode_WarpFisheye::Apply (dng_host &host, dng_negative &negative, AutoPtr<dng_image> &image) { #if qDNGValidate dng_timer timer ("WarpFisheye time"); #endif // Allocate destination image. AutoPtr<dng_image> dstImage (host.Make_dng_image (image->Bounds (), image->Planes (), image->PixelType ())); // Warp the image. AutoPtr<dng_warp_params> params (new dng_warp_params_fisheye (fWarpParams)); dng_filter_warp filter (*image, *dstImage, negative, params); filter.Initialize (host); host.PerformAreaTask (filter, image->Bounds ()); // Return the new image. image.Reset (dstImage.Release ()); } /*****************************************************************************/ uint32 dng_opcode_WarpFisheye::ParamBytes (uint32 planes) { return (1 * (uint32) sizeof (uint32) ) + // Number of planes. (4 * (uint32) sizeof (real64) * planes) + // Warp coefficients. (2 * (uint32) sizeof (real64) ); // Optical center. } /*****************************************************************************/ dng_vignette_radial_params::dng_vignette_radial_params () : fParams (kNumTerms) , fCenter (0.5, 0.5) { } /*****************************************************************************/ dng_vignette_radial_params::dng_vignette_radial_params (const dng_std_vector<real64> ¶ms, const dng_point_real64 ¢er) : fParams (params) , fCenter (center) { } /*****************************************************************************/ bool dng_vignette_radial_params::IsNOP () const { for (uint32 i = 0; i < fParams.size (); i++) { if (fParams [i] != 0.0) { return false; } } return true; } /*****************************************************************************/ bool dng_vignette_radial_params::IsValid () const { if (fParams.size () != kNumTerms) { return false; } if (fCenter.h < 0.0 || fCenter.h > 1.0 || fCenter.v < 0.0 || fCenter.v > 1.0) { return false; } return true; } /*****************************************************************************/ void dng_vignette_radial_params::Dump () const { #if qDNGValidate printf (" Radial vignette params: "); for (uint32 i = 0; i < fParams.size (); i++) { printf ("%s%.6lf", (i == 0) ? "" : ", ", fParams [i]); } printf ("\n"); printf (" Optical center:\n" " h = %.6lf\n" " v = %.6lf\n", fCenter.h, fCenter.v); #endif } /*****************************************************************************/ class dng_vignette_radial_function: public dng_1d_function { protected: const dng_vignette_radial_params fParams; public: explicit dng_vignette_radial_function (const dng_vignette_radial_params ¶ms) : fParams (params) { } // x represents r^2, where r is the normalized Euclidean distance from the // optical center to a pixel. r lies in [0,1], so r^2 (and hence x) also lies // in [0,1]. virtual real64 Evaluate (real64 x) const { DNG_REQUIRE (fParams.fParams.size () == dng_vignette_radial_params::kNumTerms, "Bad number of vignette opcode coefficients."); real64 sum = 0.0; const dng_std_vector<real64> &v = fParams.fParams; for (dng_std_vector<real64>::const_reverse_iterator i = v.rbegin (); i != v.rend (); i++) { sum = x * ((*i) + sum); } sum += 1.0; return sum; } }; /*****************************************************************************/ dng_opcode_FixVignetteRadial::dng_opcode_FixVignetteRadial (const dng_vignette_radial_params ¶ms, uint32 flags) : dng_inplace_opcode (dngOpcode_FixVignetteRadial, dngVersion_1_3_0_0, flags) , fParams (params) , fImagePlanes (1) , fSrcOriginH (0) , fSrcOriginV (0) , fSrcStepH (0) , fSrcStepV (0) , fTableInputBits (0) , fTableOutputBits (0) , fGainTable () { if (!params.IsValid ()) { ThrowBadFormat (); } } /*****************************************************************************/ dng_opcode_FixVignetteRadial::dng_opcode_FixVignetteRadial (dng_stream &stream) : dng_inplace_opcode (dngOpcode_FixVignetteRadial, stream, "FixVignetteRadial") , fParams () , fImagePlanes (1) , fSrcOriginH (0) , fSrcOriginV (0) , fSrcStepH (0) , fSrcStepV (0) , fTableInputBits (0) , fTableOutputBits (0) , fGainTable () { // Grab the size in bytes. const uint32 bytes = stream.Get_uint32 (); // Verify the size. if (bytes != ParamBytes ()) { ThrowBadFormat (); } // Read vignette coefficients. fParams.fParams = dng_std_vector<real64> (dng_vignette_radial_params::kNumTerms); for (uint32 i = 0; i < dng_vignette_radial_params::kNumTerms; i++) { fParams.fParams [i] = stream.Get_real64 (); } // Read the image center. fParams.fCenter.h = stream.Get_real64 (); fParams.fCenter.v = stream.Get_real64 (); // Debug. #if qDNGValidate if (gVerbose) { fParams.Dump (); } #endif if (!fParams.IsValid ()) { ThrowBadFormat (); } } /*****************************************************************************/ bool dng_opcode_FixVignetteRadial::IsNOP () const { return fParams.IsNOP (); } /*****************************************************************************/ bool dng_opcode_FixVignetteRadial::IsValidForNegative (const dng_negative & /* negative */) const { return fParams.IsValid (); } /*****************************************************************************/ void dng_opcode_FixVignetteRadial::PutData (dng_stream &stream) const { const uint32 bytes = ParamBytes (); stream.Put_uint32 (bytes); DNG_REQUIRE (fParams.fParams.size () == dng_vignette_radial_params::kNumTerms, "Bad number of vignette opcode coefficients."); for (uint32 i = 0; i < dng_vignette_radial_params::kNumTerms; i++) { stream.Put_real64 (fParams.fParams [i]); } stream.Put_real64 (fParams.fCenter.h); stream.Put_real64 (fParams.fCenter.v); } /*****************************************************************************/ void dng_opcode_FixVignetteRadial::Prepare (dng_negative &negative, uint32 threadCount, const dng_point &tileSize, const dng_rect &imageBounds, uint32 imagePlanes, uint32 bufferPixelType, dng_memory_allocator &allocator) { // This opcode is restricted to 32-bit images. if (bufferPixelType != ttFloat) { ThrowBadFormat (); } // Sanity check number of planes. DNG_ASSERT (imagePlanes >= 1 && imagePlanes <= kMaxColorPlanes, "Bad number of planes."); if (imagePlanes < 1 || imagePlanes > kMaxColorPlanes) { ThrowProgramError (); } fImagePlanes = imagePlanes; // Set the vignette correction curve. const dng_vignette_radial_function curve (fParams); // Grab the destination image area. const dng_rect_real64 bounds (imageBounds); // Determine the optical center and maximum radius in pixel coordinates. const dng_point_real64 centerPixel (Lerp_real64 (bounds.t, bounds.b, fParams.fCenter.v), Lerp_real64 (bounds.l, bounds.r, fParams.fCenter.h)); const real64 pixelScaleV = 1.0 / negative.PixelAspectRatio (); const real64 maxRadius = hypot (Max_real64 (Abs_real64 (centerPixel.v - bounds.t), Abs_real64 (centerPixel.v - bounds.b)) * pixelScaleV, Max_real64 (Abs_real64 (centerPixel.h - bounds.l), Abs_real64 (centerPixel.h - bounds.r))); const dng_point_real64 radius (maxRadius, maxRadius); // Find origin and scale. const real64 pixelScaleH = 1.0; fSrcOriginH = Real64ToFixed64 (-centerPixel.h * pixelScaleH / radius.h); fSrcOriginV = Real64ToFixed64 (-centerPixel.v * pixelScaleV / radius.v); fSrcStepH = Real64ToFixed64 (pixelScaleH / radius.h); fSrcStepV = Real64ToFixed64 (pixelScaleV / radius.v); // Adjust for pixel centers. fSrcOriginH += fSrcStepH >> 1; fSrcOriginV += fSrcStepV >> 1; // Evaluate 32-bit vignette correction table. dng_1d_table table32; table32.Initialize (allocator, curve, false); // Find maximum scale factor. const real64 maxScale = Max_real32 (table32.Interpolate (0.0f), table32.Interpolate (1.0f)); // Find table input bits. fTableInputBits = 16; // Find table output bits. fTableOutputBits = 15; while ((1 << fTableOutputBits) * maxScale > 65535.0) { fTableOutputBits--; } // Allocate 16-bit scale table. const uint32 tableEntries = (1 << fTableInputBits) + 1; fGainTable.Reset (allocator.Allocate (tableEntries * (uint32) sizeof (uint16))); uint16 *table16 = fGainTable->Buffer_uint16 (); // Interpolate 32-bit table into 16-bit table. const real32 scale0 = 1.0f / (1 << fTableInputBits ); const real32 scale1 = 1.0f * (1 << fTableOutputBits); for (uint32 index = 0; index < tableEntries; index++) { real32 x = index * scale0; real32 y = table32.Interpolate (x) * scale1; table16 [index] = (uint16) Round_uint32 (y); } // Prepare vignette mask buffers. { const uint32 pixelType = ttShort; const uint32 bufferSize = ComputeBufferSize(pixelType, tileSize, imagePlanes, pad16Bytes); for (uint32 threadIndex = 0; threadIndex < threadCount; threadIndex++) { fMaskBuffers [threadIndex] . Reset (allocator.Allocate (bufferSize)); } } } /*****************************************************************************/ void dng_opcode_FixVignetteRadial::ProcessArea (dng_negative & /* negative */, uint32 threadIndex, dng_pixel_buffer &buffer, const dng_rect &dstArea, const dng_rect & /* imageBounds */) { // Setup mask pixel buffer. dng_pixel_buffer maskPixelBuffer(dstArea, 0, fImagePlanes, ttShort, pcRowInterleavedAlign16, fMaskBuffers [threadIndex]->Buffer ()); // Compute mask. DoVignetteMask16 (maskPixelBuffer.DirtyPixel_uint16 (dstArea.t, dstArea.l), dstArea.H (), dstArea.W (), maskPixelBuffer.RowStep (), fSrcOriginH + fSrcStepH * dstArea.l, fSrcOriginV + fSrcStepV * dstArea.t, fSrcStepH, fSrcStepV, fTableInputBits, fGainTable->Buffer_uint16 ()); // Apply mask. DoVignette32 (buffer.DirtyPixel_real32 (dstArea.t, dstArea.l), maskPixelBuffer.ConstPixel_uint16 (dstArea.t, dstArea.l), dstArea.H (), dstArea.W (), fImagePlanes, buffer.RowStep (), buffer.PlaneStep (), maskPixelBuffer.RowStep (), fTableOutputBits); } /*****************************************************************************/ uint32 dng_opcode_FixVignetteRadial::ParamBytes () { const uint32 N = dng_vignette_radial_params::kNumTerms; return ((N * sizeof (real64)) + // Vignette coefficients. (2 * sizeof (real64))); // Optical center. } /*****************************************************************************/