From 871c75e4948322434e24982c432674217f7ca9ce Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Sat, 8 Apr 2017 10:13:53 +0200 Subject: [PATCH] Correctly apply LCP "distortion" correction for Fisheye lenses Note: as discussed e.g. at http://lensfun.sourceforge.net/manual/corrections.html, this is really a combination of distortion correction and change of projection (from fisheye to rectilinear), but "distortion correction" is how the Adobe camera model calls it, and this is how it appears in the RT gui --- rtengine/iptransform.cc | 38 +++++++++++++++++++------------ rtengine/lcp.cc | 50 ++++++++++++++++++++++++++++++++--------- rtengine/lcp.h | 3 ++- 3 files changed, 65 insertions(+), 26 deletions(-) diff --git a/rtengine/iptransform.cc b/rtengine/iptransform.cc index 04ba4966d..399711c2e 100644 --- a/rtengine/iptransform.cc +++ b/rtengine/iptransform.cc @@ -85,16 +85,17 @@ float normn (float a, float b, int n) } } -void correct_distortion (const rtengine::LCPMapper *lcp, double &x, double &y, - int cx, int cy) + +void correct_distortion(const rtengine::LCPMapper *lcp, double &x, double &y, + int cx, int cy, double scale) { assert (lcp); x += cx; y += cy; - lcp->correctDistortion (x, y); - x -= cx; - y -= cy; + lcp->correctDistortion(x, y, scale); + x -= (cx * scale); + y -= (cy * scale); } } @@ -156,11 +157,14 @@ bool ImProcFunctions::transCoord (int W, int H, const std::vector &src, double x_d = src[i].x, y_d = src[i].y; if (pLCPMap && params->lensProf.useDist) { - correct_distortion (pLCPMap, x_d, y_d, 0, 0); + correct_distortion (pLCPMap, x_d, y_d, 0, 0, ascale); + } else { + x_d *= ascale; + y_d *= ascale; } - y_d = ascale * (y_d - h2); - x_d = ascale * (x_d - w2); + x_d += ascale * (0 - w2); // centering x coord & scale + y_d += ascale * (0 - h2); // centering y coord & scale if (needsPerspective()) { // horizontal perspective transformation @@ -799,11 +803,14 @@ void ImProcFunctions::transformHighQuality (Imagefloat* original, Imagefloat* tr double x_d = x, y_d = y; if (enableLCPDist) { - correct_distortion (pLCPMap, x_d, y_d, cx, cy); // must be first transform + correct_distortion(pLCPMap, x_d, y_d, cx, cy, ascale); // must be first transform + } else { + x_d *= ascale; + y_d *= ascale; } - x_d = ascale * (x_d + cx - w2); // centering x coord & scale - y_d = ascale * (y_d + cy - h2); // centering y coord & scale + x_d += ascale * (cx - w2); // centering x coord & scale + y_d += ascale * (cy - h2); // centering y coord & scale double vig_x_d = 0., vig_y_d = 0.; @@ -976,11 +983,14 @@ void ImProcFunctions::transformPreview (Imagefloat* original, Imagefloat* transf double x_d = x, y_d = y; if (pLCPMap && params->lensProf.useDist) { - correct_distortion (pLCPMap, x_d, y_d, cx, cy); // must be first transform + correct_distortion(pLCPMap, x_d, y_d, cx, cy, ascale); // must be first transform + } else { + x_d *= ascale; + y_d *= ascale; } - y_d = ascale * (y_d + cy - h2); // centering y coord & scale - x_d = ascale * (x_d + cx - w2); // centering x coord & scale + x_d += ascale * (cx - w2); // centering x coord & scale + y_d += ascale * (cy - h2); // centering y coord & scale double vig_x_d = 0., vig_y_d = 0.; diff --git a/rtengine/lcp.cc b/rtengine/lcp.cc index fdf46270f..d2d18b7d8 100644 --- a/rtengine/lcp.cc +++ b/rtengine/lcp.cc @@ -198,24 +198,52 @@ LCPMapper::LCPMapper(LCPProfile* pProf, float focalLength, float focalLength35mm } enableCA = !vignette && focusDist > 0; + isFisheye = pProf->isFisheye; } -void LCPMapper::correctDistortion(double& x, double& y) const +void LCPMapper::correctDistortion(double& x, double& y, double scale) const { - double xd = (x - mc.x0) / mc.fx, yd = (y - mc.y0) / mc.fy; + if (isFisheye) { + double u = x * scale; + double v = y * scale; + double u0 = mc.x0 * scale; + double v0 = mc.y0 * scale; + double du = (u - u0); + double dv = (v - v0); + double fx = mc.fx; + double fy = mc.fy; + double k1 = mc.param[0]; + double k2 = mc.param[1]; + double r = sqrt(du * du + dv * dv); + double f = sqrt(fx*fy / (scale * scale)); + double th = atan2(r, f); + double th2 = th * th; + double cfact = (((k2 * th2 + k1) * th2 + 1) * th) / r; + double ud = cfact * fx * du + u0; + double vd = cfact * fy * dv + v0; - const LCPModelCommon::Param aDist = mc.param; - double rsqr = xd * xd + yd * yd; - double xfac = aDist[swapXY ? 3 : 4], yfac = aDist[swapXY ? 4 : 3]; + x = ud; + y = vd; + } else { + x *= scale; + y *= scale; + double x0 = mc.x0 * scale; + double y0 = mc.y0 * scale; + double xd = (x - x0) / mc.fx, yd = (y - y0) / mc.fy; - double commonFac = (((aDist[2] * rsqr + aDist[1]) * rsqr + aDist[0]) * rsqr + 1.) - + 2. * (yfac * yd + xfac * xd); + const LCPModelCommon::Param aDist = mc.param; + double rsqr = xd * xd + yd * yd; + double xfac = aDist[swapXY ? 3 : 4], yfac = aDist[swapXY ? 4 : 3]; - double xnew = xd * commonFac + xfac * rsqr; - double ynew = yd * commonFac + yfac * rsqr; + double commonFac = (((aDist[2] * rsqr + aDist[1]) * rsqr + aDist[0]) * rsqr + 1.) + + 2. * (yfac * yd + xfac * xd); - x = xnew * mc.fx + mc.x0; - y = ynew * mc.fy + mc.y0; + double xnew = xd * commonFac + xfac * rsqr; + double ynew = yd * commonFac + yfac * rsqr; + + x = xnew * mc.fx + x0; + y = ynew * mc.fy + y0; + } } void LCPMapper::correctCA(double& x, double& y, int channel) const diff --git a/rtengine/lcp.h b/rtengine/lcp.h index e108e61b7..8cfcd5694 100644 --- a/rtengine/lcp.h +++ b/rtengine/lcp.h @@ -142,6 +142,7 @@ class LCPMapper bool swapXY; LCPModelCommon mc; LCPModelCommon chrom[3]; // in order RedGreen/Green/BlueGreen + bool isFisheye; public: bool enableCA; // is the mapper capable if CA correction? @@ -150,7 +151,7 @@ public: LCPMapper(LCPProfile* pProf, float focalLength, float focalLength35mm, float focusDist, float aperture, bool vignette, bool useCADistP, int fullWidth, int fullHeight, const CoarseTransformParams& coarse, int rawRotationDeg); - void correctDistortion(double& x, double& y) const; // MUST be the first stage + void correctDistortion(double& x, double& y, double scale) const; // MUST be the first stage void correctCA(double& x, double& y, int channel) const; void processVignetteLine(int width, int y, float *line) const; void processVignetteLine3Channels(int width, int y, float *line) const;