/*****************************************************************************/
|
// 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.
|
|
}
|
|
/*****************************************************************************/
|
