From 871c75e4948322434e24982c432674217f7ca9ce Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Sat, 8 Apr 2017 10:13:53 +0200 Subject: [PATCH 1/2] 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; From 3464c0ec9221cb842af6b4f3bf8d9b19377d2418 Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Sat, 8 Apr 2017 19:58:32 +0200 Subject: [PATCH 2/2] fixed glitch in handling crop windows at the edges of the image when LCP distortion correction is enabled --- rtengine/dcrop.cc | 36 +++++++++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 9 deletions(-) diff --git a/rtengine/dcrop.cc b/rtengine/dcrop.cc index 18cb927aa..351ce1f08 100644 --- a/rtengine/dcrop.cc +++ b/rtengine/dcrop.cc @@ -1145,15 +1145,33 @@ bool Crop::setCropSizes (int rcx, int rcy, int rcw, int rch, int skip, bool inte parent->ipf.transCoord (parent->fw, parent->fh, bx1, by1, bw, bh, orx, ory, orw, orh); if (check_need_larger_crop_for_lcp_distortion(parent->fw, parent->fh, orx, ory, orw, orh, parent->params)) { - double dW = double(parent->fw) * 0.15 / skip; // TODO - this is hardcoded ATM! - double dH = double(parent->fh) * 0.15 / skip; // this is an estimate of the max - // distortion relative to the image - // size. BUT IS 15% REALLY ENOUGH? - // In fact, is there a better way?? - orx = max(int(orx - dW/2.0), 0); - ory = max(int(ory - dH/2.0), 0); - orw = min(int(orw + dW), parent->fw - orx); - orh = min(int(orh + dH), parent->fh - ory); + // TODO - this is an estimate of the max distortion relative to the image size. ATM it is hardcoded to be 15%, which seems enough. If not, need to revise + int dW = int(double(parent->fw) * 0.15 / (2 * skip)); + int dH = int(double(parent->fh) * 0.15 / (2 * skip)); + int x1 = orx - dW; + int x2 = orx + orw + dW; + int y1 = ory - dH; + int y2 = ory + orh + dH; + if (x1 < 0) { + x2 += -x1; + x1 = 0; + } + if (x2 > parent->fw) { + x1 -= x2 - parent->fw; + x2 = parent->fw; + } + if (y1 < 0) { + y2 += -y1; + y1 = 0; + } + if (y2 > parent->fh) { + y1 -= y2 - parent->fh; + y2 = parent->fh; + } + orx = max(x1, 0); + ory = max(y1, 0); + orw = min(x2 - x1, parent->fw - orx); + orh = min(y2 - y1, parent->fh - ory); }