diff --git a/rtengine/lcp.cc b/rtengine/lcp.cc index ce1a6c67e..310361a72 100644 --- a/rtengine/lcp.cc +++ b/rtengine/lcp.cc @@ -16,103 +16,107 @@ * You should have received a copy of the GNU General Public License * along with RawTherapee. If not, see . */ + +#include #include #include "lcp.h" -#include "iccmatrices.h" -#include "iccstore.h" -#include "rawimagesource.h" -#include "improcfun.h" -#include "rt_math.h" +#include #ifdef WIN32 #include -// for GCC32 -#ifndef _WIN32_IE -#define _WIN32_IE 0x0600 -#endif #include #endif using namespace std; using namespace rtengine; -using namespace rtexif; -LCPModelCommon::LCPModelCommon() +LCPModelCommon::LCPModelCommon() : + foc_len_x(-1.0f), + foc_len_y(-1.0f), + img_center_x(0.5f), + img_center_y(0.5f), + param{{}}, + scale_factor(1.0f), + mean_error(0.0), + bad_error(false), + x0(0.0f), + y0(0.0f), + fx(0.0f), + fy(0.0f), + rfx(0.0f), + rfy(0.0f), + vign_param{{}} { - focLenX = focLenY = -1; - imgXCenter = imgYCenter = 0.5; - x0 = y0 = fx = fy = meanErr = 0; - badErr = false; - - for (int i = 0; i < 5; i++) { - param[i] = 0; - } - - scaleFac = 1; } bool LCPModelCommon::empty() const { - return param[0] == 0 && param[1] == 0 && param[2] == 0; + return param[0] == 0.0f && param[1] == 0.0f && param[2] == 0.0f; } void LCPModelCommon::print() const { - printf("focLen %g/%g; imgCenter %g/%g; scale %g; err %g\n", focLenX, focLenY, imgXCenter, imgYCenter, scaleFac, meanErr); + printf("focLen %g/%g; imgCenter %g/%g; scale %g; err %g\n", foc_len_x, foc_len_y, img_center_x, img_center_y, scale_factor, mean_error); printf("xy0 %g/%g fxy %g/%g\n", x0, y0, fx, fy); printf("param: %g/%g/%g/%g/%g\n", param[0], param[1], param[2], param[3], param[4]); } -// weightened merge two parameters +// weighted merge two parameters void LCPModelCommon::merge(const LCPModelCommon& a, const LCPModelCommon& b, float facA) { - float facB = 1 - facA; + const float facB = 1.0f - facA; - focLenX = facA * a.focLenX + facB * b.focLenX; - focLenY = facA * a.focLenY + facB * b.focLenY; - imgXCenter = facA * a.imgXCenter + facB * b.imgXCenter; - imgYCenter = facA * a.imgYCenter + facB * b.imgYCenter; - scaleFac = facA * a.scaleFac + facB * b.scaleFac; - meanErr = facA * a.meanErr + facB * b.meanErr; + foc_len_x = facA * a.foc_len_x + facB * b.foc_len_x; + foc_len_y = facA * a.foc_len_y + facB * b.foc_len_y; + img_center_x = facA * a.img_center_x + facB * b.img_center_x; + img_center_y = facA * a.img_center_y + facB * b.img_center_y; + scale_factor = facA * a.scale_factor + facB * b.scale_factor; + mean_error = facA * a.mean_error + facB * b.mean_error; for (int i = 0; i < 5; i++) { param[i] = facA * a.param[i] + facB * b.param[i]; } + + const float param0Sqr = param[0] * param[0]; + + vign_param[0] = -param[0]; + vign_param[1] = param0Sqr - param[1]; + vign_param[2] = param0Sqr * param[0] - 2.0f * param[0] * param[1] + param[2]; + vign_param[3] = param0Sqr * param0Sqr + param[1] * param[1] + 2.0f * param[0] * param[2] - 3.0f * param0Sqr * param[1]; + } void LCPModelCommon::prepareParams(int fullWidth, int fullHeight, float focalLength, float focalLength35mm, float sensorFormatFactor, bool swapXY, bool mirrorX, bool mirrorY) { // Mention that the Adobe technical paper has a bug here, the DMAX is handled differently for focLen and imgCenter - int Dmax = fullWidth; - - if (fullHeight > fullWidth) { - Dmax = fullHeight; - } + const int Dmax = std::max(fullWidth, fullHeight); // correct focLens - if (focLenX < 0) { // they may not be given + if (foc_len_x < 0.0f) { // they may not be given // and 35mm may not be given either - if (focalLength35mm < 1) { + if (focalLength35mm < 1.0f) { focalLength35mm = focalLength * sensorFormatFactor; } - focLenX = focLenY = focalLength / ( 35 * focalLength / focalLength35mm); // focLen must be calculated in pixels + foc_len_x = foc_len_y = focalLength / (35.0f * focalLength / focalLength35mm); // focLen must be calculated in pixels } if (swapXY) { - x0 = (mirrorX ? 1. - imgYCenter : imgYCenter) * fullWidth; - y0 = (mirrorY ? 1. - imgXCenter : imgXCenter) * fullHeight; - fx = focLenY * Dmax; - fy = focLenX * Dmax; + x0 = (mirrorX ? 1.0f - img_center_y : img_center_y) * fullWidth; + y0 = (mirrorY ? 1.0f - img_center_x : img_center_x) * fullHeight; + fx = foc_len_y * Dmax; + fy = foc_len_x * Dmax; } else { - x0 = (mirrorX ? 1. - imgXCenter : imgXCenter) * fullWidth; - y0 = (mirrorY ? 1. - imgYCenter : imgYCenter) * fullHeight; - fx = focLenX * Dmax; - fy = focLenY * Dmax; + x0 = (mirrorX ? 1.0f - img_center_x : img_center_x) * fullWidth; + y0 = (mirrorY ? 1.0f - img_center_y : img_center_y) * fullHeight; + fx = foc_len_x * Dmax; + fy = foc_len_y * Dmax; } + rfx = 1.0f / fx; + rfy = 1.0f / fy; //printf("FW %i /X0 %g FH %i /Y0 %g %g\n",fullWidth,x0,fullHeight,y0, imgYCenter); } @@ -125,9 +129,9 @@ LCPPersModel::LCPPersModel() // mode: 0=distortion, 1=vignette, 2=CA bool LCPPersModel::hasModeData(int mode) const { - return (mode == 0 && !vignette.empty() && !vignette.badErr) || (mode == 1 && !base.empty() && !base.badErr) + return (mode == 0 && !vignette.empty() && !vignette.bad_error) || (mode == 1 && !base.empty() && !base.bad_error) || (mode == 2 && !chromRG.empty() && !chromG.empty() && !chromBG.empty() && - !chromRG.badErr && !chromG.badErr && !chromBG.badErr); + !chromRG.bad_error && !chromG.bad_error && !chromBG.bad_error); } void LCPPersModel::print() const @@ -200,7 +204,7 @@ void LCPMapper::correctDistortion(double& x, double& y) const { double xd = (x - mc.x0) / mc.fx, yd = (y - mc.y0) / mc.fy; - const float* aDist = mc.param; + const LCPModelCommon::Param aDist = mc.param; double rsqr = xd * xd + yd * yd; double xfac = aDist[swapXY ? 3 : 4], yfac = aDist[swapXY ? 4 : 3]; @@ -228,7 +232,7 @@ void LCPMapper::correctCA(double& x, double& y, int channel) const // Green contains main distortion, just like base if (useCADist) { - const float* aDist = chrom[1].param; + const LCPModelCommon::Param aDist = chrom[1].param; double rsqr = xd * xd + yd * yd; double xfac = aDist[swapXY ? 3 : 4], yfac = aDist[swapXY ? 4 : 3]; @@ -252,30 +256,72 @@ void LCPMapper::correctCA(double& x, double& y, int channel) const yd = ygreen; rsqr = xd * xd + yd * yd; - const float* aCA = chrom[channel].param; + const LCPModelCommon::Param aCA = chrom[channel].param; double xfac = aCA[swapXY ? 3 : 4], yfac = aCA[swapXY ? 4 : 3]; double commonSum = 1. + rsqr * (aCA[0] + rsqr * (aCA[1] + aCA[2] * rsqr)) + 2. * (yfac * yd + xfac * xd); - x = (chrom[channel].scaleFac * ( xd * commonSum + xfac * rsqr )) * chrom[channel].fx + chrom[channel].x0; - y = (chrom[channel].scaleFac * ( yd * commonSum + yfac * rsqr )) * chrom[channel].fy + chrom[channel].y0; + x = (chrom[channel].scale_factor * ( xd * commonSum + xfac * rsqr )) * chrom[channel].fx + chrom[channel].x0; + y = (chrom[channel].scale_factor * ( yd * commonSum + yfac * rsqr )) * chrom[channel].fy + chrom[channel].y0; } } -float LCPMapper::calcVignetteFac(int x, int y) const +SSEFUNCTION void LCPMapper::processVignetteLine(int width, int y, float *line) const { // No need for swapXY, since vignette is in RAW and always before rotation - double xd = ((double)x - mc.x0) / mc.fx, yd = ((double)y - mc.y0) / mc.fy; + float yd = ((float)y - mc.y0) * mc.rfy; + yd *= yd; + int x = 0; +#ifdef __SSE2__ + const vfloat fourv = F2V(4.f); + const vfloat zerov = F2V(0.f); + const vfloat ydv = F2V(yd); + const vfloat p0 = F2V(mc.vign_param[0]); + const vfloat p1 = F2V(mc.vign_param[1]); + const vfloat p2 = F2V(mc.vign_param[2]); + const vfloat p3 = F2V(mc.vign_param[3]); + const vfloat x0v = F2V(mc.x0); + const vfloat rfxv = F2V(mc.rfx); - const float* aVig = mc.param; - double rsqr = xd * xd + yd * yd; - double param0Sqr = aVig[0] * aVig[0]; - - return 1. + rsqr * (-aVig[0] + rsqr * ((param0Sqr - aVig[1]) - - (param0Sqr * aVig[0] - 2.*aVig[0] * aVig[1] + aVig[2]) * rsqr - + (param0Sqr * param0Sqr + aVig[1] * aVig[1] - + 2.*aVig[0] * aVig[2] - 3.*param0Sqr * aVig[1]) * rsqr * rsqr)); + vfloat xv = _mm_setr_ps(0.f, 1.f, 2.f, 3.f); + for (; x < width-3; x+=4) { + vfloat xdv = (xv - x0v) * rfxv; + vfloat rsqr = xdv * xdv + ydv; + vfloat vignFactorv = rsqr * (p0 + rsqr * (p1 - p2 * rsqr + p3 * rsqr * rsqr)); + vfloat valv = LVFU(line[x]); + valv += valv * vselfzero(vmaskf_gt(valv, zerov), vignFactorv); + STVFU(line[x], valv); + xv += fourv; + } +#endif // __SSE2__ + for (; x < width; x++) { + if (line[x] > 0) { + float xd = ((float)x - mc.x0) * mc.rfx; + const LCPModelCommon::VignParam vignParam = mc.vign_param; + float rsqr = xd * xd + yd; + line[x] += line[x] * rsqr * (vignParam[0] + rsqr * ((vignParam[1]) - (vignParam[2]) * rsqr + (vignParam[3]) * rsqr * rsqr)); + } + } } +SSEFUNCTION void LCPMapper::processVignetteLine3Channels(int width, int y, float *line) const +{ + // No need for swapXY, since vignette is in RAW and always before rotation + float yd = ((float)y - mc.y0) * mc.rfy; + yd *= yd; + const LCPModelCommon::VignParam vignParam = mc.vign_param; + for (int x = 0; x < width; x++) { + float xd = ((float)x - mc.x0) * mc.rfx; + float rsqr = xd * xd + yd; + float vignetteFactor = rsqr * (vignParam[0] + rsqr * ((vignParam[1]) - (vignParam[2]) * rsqr + (vignParam[3]) * rsqr * rsqr)); + for(int c = 0;c < 3; ++c) { + if (line[3*x+c] > 0) { + line[3*x+c] += line[3*x+c] * vignetteFactor; + } + } + } +} + + LCPProfile::LCPProfile(const Glib::ustring &fname) { const int BufferSize = 8192; @@ -336,17 +382,17 @@ int LCPProfile::filterBadFrames(double maxAvgDevFac, int minFramesLeft) for (int pm = 0; pm < MaxPersModelCount && aPersModel[pm]; pm++) { if (aPersModel[pm]->hasModeData(0)) { - errVignette += aPersModel[pm]->vignette.meanErr; + errVignette += aPersModel[pm]->vignette.mean_error; vignetteCount++; } if (aPersModel[pm]->hasModeData(1)) { - errBase += aPersModel[pm]->base.meanErr; + errBase += aPersModel[pm]->base.mean_error; baseCount++; } if (aPersModel[pm]->hasModeData(2)) { - errChrom += std::max(std::max(aPersModel[pm]->chromRG.meanErr, aPersModel[pm]->chromG.meanErr), aPersModel[pm]->chromBG.meanErr); + errChrom += std::max(std::max(aPersModel[pm]->chromRG.mean_error, aPersModel[pm]->chromG.mean_error), aPersModel[pm]->chromBG.mean_error); chromCount++; } } @@ -369,20 +415,20 @@ int LCPProfile::filterBadFrames(double maxAvgDevFac, int minFramesLeft) // Now mark all the bad ones as bad, and hasModeData will return false; for (int pm = 0; pm < MaxPersModelCount && aPersModel[pm]; pm++) { - if (aPersModel[pm]->hasModeData(0) && aPersModel[pm]->vignette.meanErr > maxAvgDevFac * errVignette) { - aPersModel[pm]->vignette.badErr = true; + if (aPersModel[pm]->hasModeData(0) && aPersModel[pm]->vignette.mean_error > maxAvgDevFac * errVignette) { + aPersModel[pm]->vignette.bad_error = true; filtered++; } - if (aPersModel[pm]->hasModeData(1) && aPersModel[pm]->base.meanErr > maxAvgDevFac * errBase) { - aPersModel[pm]->base.badErr = true; + if (aPersModel[pm]->hasModeData(1) && aPersModel[pm]->base.mean_error > maxAvgDevFac * errBase) { + aPersModel[pm]->base.bad_error = true; filtered++; } if (aPersModel[pm]->hasModeData(2) && - (aPersModel[pm]->chromRG.meanErr > maxAvgDevFac * errChrom || aPersModel[pm]->chromG.meanErr > maxAvgDevFac * errChrom - || aPersModel[pm]->chromBG.meanErr > maxAvgDevFac * errChrom)) { - aPersModel[pm]->chromRG.badErr = aPersModel[pm]->chromG.badErr = aPersModel[pm]->chromBG.badErr = true; + (aPersModel[pm]->chromRG.mean_error > maxAvgDevFac * errChrom || aPersModel[pm]->chromG.mean_error > maxAvgDevFac * errChrom + || aPersModel[pm]->chromBG.mean_error > maxAvgDevFac * errChrom)) { + aPersModel[pm]->chromRG.bad_error = aPersModel[pm]->chromG.bad_error = aPersModel[pm]->chromBG.bad_error = true; filtered++; } } @@ -438,48 +484,48 @@ void LCPProfile::calcParams(int mode, float focalLength, float focusDist, float if (aPersModel[pm]->hasModeData(mode)) { if (mode == 0) { - meanErr = aPersModel[pm]->vignette.meanErr; + meanErr = aPersModel[pm]->vignette.mean_error; // by aperture (vignette), and max out focus distance // tests showed doing this by log(aperture) is not as advisable if (aPersModel[pm]->focLen == bestFocLenLow && ( - (aper == aperture && pLow->vignette.meanErr > meanErr) + (aper == aperture && pLow->vignette.mean_error > meanErr) || (aper >= aperture && aper < pLow->aperture && pLow->aperture > aperture) || (aper <= aperture && (pLow->aperture > aperture || fabs(aperture - aper) < fabs(aperture - pLow->aperture))))) { pLow = aPersModel[pm]; } if (aPersModel[pm]->focLen == bestFocLenHigh && ( - (aper == aperture && pHigh->vignette.meanErr > meanErr) + (aper == aperture && pHigh->vignette.mean_error > meanErr) || (aper <= aperture && aper > pHigh->aperture && pHigh->aperture < aperture) || (aper >= aperture && (pHigh->aperture < aperture || fabs(aperture - aper) < fabs(aperture - pHigh->aperture))))) { pHigh = aPersModel[pm]; } } else { - meanErr = (mode == 1 ? aPersModel[pm]->base.meanErr : aPersModel[pm]->chromG.meanErr); + meanErr = (mode == 1 ? aPersModel[pm]->base.mean_error : aPersModel[pm]->chromG.mean_error); if (focusDist > 0) { // by focus distance if (aPersModel[pm]->focLen == bestFocLenLow && ( - (focDist == focusDist && (mode == 1 ? pLow->base.meanErr : pLow->chromG.meanErr) > meanErr) + (focDist == focusDist && (mode == 1 ? pLow->base.mean_error : pLow->chromG.mean_error) > meanErr) || (focDist >= focusDist && focDist < pLow->focDist && pLow->focDist > focusDist) || (focDist <= focusDist && (pLow->focDist > focusDist || fabs(focusDistLog - focDistLog) < fabs(focusDistLog - (log(pLow->focDist) + euler)))))) { pLow = aPersModel[pm]; } if (aPersModel[pm]->focLen == bestFocLenHigh && ( - (focDist == focusDist && (mode == 1 ? pHigh->base.meanErr : pHigh->chromG.meanErr) > meanErr) + (focDist == focusDist && (mode == 1 ? pHigh->base.mean_error : pHigh->chromG.mean_error) > meanErr) || (focDist <= focusDist && focDist > pHigh->focDist && pHigh->focDist < focusDist) || (focDist >= focusDist && (pHigh->focDist < focusDist || fabs(focusDistLog - focDistLog) < fabs(focusDistLog - (log(pHigh->focDist) + euler)))))) { pHigh = aPersModel[pm]; } } else { // no focus distance available, just error - if (aPersModel[pm]->focLen == bestFocLenLow && (mode == 1 ? pLow->base.meanErr : pLow->chromG.meanErr) > meanErr) { + if (aPersModel[pm]->focLen == bestFocLenLow && (mode == 1 ? pLow->base.mean_error : pLow->chromG.mean_error) > meanErr) { pLow = aPersModel[pm]; } - if (aPersModel[pm]->focLen == bestFocLenHigh && (mode == 1 ? pHigh->base.meanErr : pHigh->chromG.meanErr) > meanErr) { + if (aPersModel[pm]->focLen == bestFocLenHigh && (mode == 1 ? pHigh->base.mean_error : pHigh->chromG.mean_error) > meanErr) { pHigh = aPersModel[pm]; } } @@ -715,17 +761,17 @@ void XMLCALL LCPProfile::XmlTextHandler(void *pLCPProfile, const XML_Char *s, in // Section depended if (!strcmp("FocalLengthX", tag)) { - pProf->pCurCommon->focLenX = atof(raw); + pProf->pCurCommon->foc_len_x = atof(raw); } else if (!strcmp("FocalLengthY", tag)) { - pProf->pCurCommon->focLenY = atof(raw); + pProf->pCurCommon->foc_len_y = atof(raw); } else if (!strcmp("ImageXCenter", tag)) { - pProf->pCurCommon->imgXCenter = atof(raw); + pProf->pCurCommon->img_center_x = atof(raw); } else if (!strcmp("ImageYCenter", tag)) { - pProf->pCurCommon->imgYCenter = atof(raw); + pProf->pCurCommon->img_center_y = atof(raw); } else if (!strcmp("ScaleFactor", tag)) { - pProf->pCurCommon->scaleFac = atof(raw); + pProf->pCurCommon->scale_factor = atof(raw); } else if (!strcmp("ResidualMeanError", tag)) { - pProf->pCurCommon->meanErr = atof(raw); + pProf->pCurCommon->mean_error = atof(raw); } else if (!strcmp("RadialDistortParam1", tag) || !strcmp("VignetteModelParam1", tag)) { pProf->pCurCommon->param[0] = atof(raw); } else if (!strcmp("RadialDistortParam2", tag) || !strcmp("VignetteModelParam2", tag)) { diff --git a/rtengine/lcp.h b/rtengine/lcp.h index 908bacf56..fa8cf226d 100644 --- a/rtengine/lcp.h +++ b/rtengine/lcp.h @@ -17,36 +17,52 @@ * along with RawTherapee. If not, see . */ -#ifndef _LCP_ -#define _LCP_ +#pragma once + +#include +#include +#include + +#include +#include #include "imagefloat.h" -#include "../rtgui/threadutils.h" -#include -#include -#include -#include -#include +#include "opthelper.h" namespace rtengine { + // Perspective model common data, also used for Vignette and Fisheye -class LCPModelCommon +class LCPModelCommon final { public: - float focLenX, focLenY, imgXCenter, imgYCenter; - float param[5]; // k1..k5, resp. alpha1..5 - float scaleFac; // alpha0 - double meanErr; - bool badErr; - - double x0, y0, fx, fy; // prepared params - LCPModelCommon(); bool empty() const; // is it empty void print() const; // printf all values void merge(const LCPModelCommon& a, const LCPModelCommon& b, float facA); void prepareParams(int fullWidth, int fullHeight, float focalLength, float focalLength35mm, float sensorFormatFactor, bool swapXY, bool mirrorX, bool mirrorY); + +//private: + using Param = std::array; + using VignParam = std::array; + + float foc_len_x; + float foc_len_y; + float img_center_x; + float img_center_y; + Param param; // k1..k5, resp. alpha1..5 + float scale_factor; // alpha0 + double mean_error; + bool bad_error; + + // prepared params + float x0; + float y0; + float fx; + float fy; + float rfx; + float rfy; + VignParam vign_param; }; class LCPPersModel @@ -130,9 +146,10 @@ 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 correctCA(double& x, double& y, int channel) const; - float calcVignetteFac (int x, int y) const; // MUST be in RAW + void correctDistortion(double& x, double& y) 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; }; + } -#endif diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index 9c8f9788b..9678fcaf5 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -1769,32 +1769,22 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le if (pLCPProf) { // don't check focal length to allow distortion correction for lenses without chip, also pass dummy focal length 1 in case of 0 LCPMapper map(pLCPProf, max(idata->getFocalLen(), 1.0), idata->getFocalLen35mm(), idata->getFocusDist(), idata->getFNumber(), true, false, W, H, coarse, -1); -if (ri->getSensorType() == ST_BAYER || ri->getSensorType() == ST_FUJI_XTRANS || ri->get_colors() == 1) { + if (ri->getSensorType() == ST_BAYER || ri->getSensorType() == ST_FUJI_XTRANS || ri->get_colors() == 1) { + #ifdef _OPENMP - #pragma omp parallel for + #pragma omp parallel for schedule(dynamic,16) #endif for (int y = 0; y < H; y++) { - for (int x = 0; x < W; x++) { - if (rawData[y][x] > 0) { - rawData[y][x] *= map.calcVignetteFac(x, y); - } - } + map.processVignetteLine(W, y, rawData[y]); } } else if(ri->get_colors() == 3) { #ifdef _OPENMP - #pragma omp parallel for + #pragma omp parallel for schedule(dynamic,16) #endif for (int y = 0; y < H; y++) { - for (int x = 0; x < W; x++) { - float vignFactor = map.calcVignetteFac(x, y); - for(int c = 0;c < 3; ++c) { - if (rawData[y][3 * x + c] > 0) { - rawData[y][3 * x + c] *= vignFactor; - } - } - } + map.processVignetteLine3Channels(W, y, rawData[y]); } } }