liz/rawTherapee
liz/rawTherapee
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #3713 from Beep6581/lcpvignette_speedup

Browse Source 
Speedup for lcp vignette correction
This commit is contained in:
Ingo Weyrich 2017-02-23 21:41:00 +01:00 committed by GitHub
commit bc3003bc0d
3 changed files with 178 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

222
rtengine/lcp.cc
View File

@ -16,103 +16,107 @@
* You should have received a copy of the GNU General Public License
* along with RawTherapee. If not, see <http://www.gnu.org/licenses/>.
*/
#include <algorithm>
#include <cstring>
#include "lcp.h"
#include "iccmatrices.h"
#include "iccstore.h"
#include "rawimagesource.h"
#include "improcfun.h"
#include "rt_math.h"
#include <glib/gstdio.h>
#ifdef WIN32
#include <windows.h>
// for GCC32
#ifndef _WIN32_IE
#define _WIN32_IE 0x0600
#endif
#include <shlobj.h>
#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,29 +256,71 @@ 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)
{
@ -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)) {

55
rtengine/lcp.h
View File

@ -17,36 +17,52 @@
* along with RawTherapee. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _LCP_
#define _LCP_
#pragma once
#include <array>
#include <map>
#include <string>
#include <glibmm.h>
#include <expat.h>
#include "imagefloat.h"
#include "../rtgui/threadutils.h"
#include <glibmm.h>
#include <map>
#include <list>
#include <string>
#include <expat.h>
#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<float, 5>;
using VignParam = std::array<float, 4>;
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
@ -132,7 +148,8 @@ public:
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 processVignetteLine(int width, int y, float *line) const;
void processVignetteLine3Channels(int width, int y, float *line) const;
};
}
#endif

20
rtengine/rawimagesource.cc
View File

@ -1770,31 +1770,21 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le
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) {
#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]);
}
}
}