From 41c1f21c7687612fbd215322e26fb347b9482ba7 Mon Sep 17 00:00:00 2001 From: heckflosse Date: Fri, 1 Jun 2018 19:13:59 +0200 Subject: [PATCH] Pixelshift: replaced checkbox to use lmmse by a combobox to allow further demosaicers for parts with motion, #4579 --- rtdata/languages/default | 1 + rtengine/amaze_demosaic_RT.cc | 2 +- rtengine/amaze_vng4_demosaic_RT.cc | 12 +- rtengine/demosaic_algos.cc | 2 +- rtengine/demosaic_algos.cc.save-failed | 3861 ++++++++++++++++++++++++ rtengine/pixelshift.cc | 13 +- rtengine/procparams.cc | 43 +- rtengine/procparams.h | 11 +- rtengine/rawimagesource.cc | 2 +- rtengine/rawimagesource.h | 6 +- rtengine/simpleprocess.cc | 2 - rtgui/bayerprocess.cc | 66 +- rtgui/bayerprocess.h | 4 +- rtgui/paramsedited.cc | 14 +- rtgui/paramsedited.h | 2 +- rtgui/partialpastedlg.cc | 2 +- rtgui/ppversion.h | 4 +- 17 files changed, 3997 insertions(+), 50 deletions(-) create mode 100644 rtengine/demosaic_algos.cc.save-failed diff --git a/rtdata/languages/default b/rtdata/languages/default index 5f16404c2..8231932fb 100644 --- a/rtdata/languages/default +++ b/rtdata/languages/default @@ -1790,6 +1790,7 @@ TP_RAW_MONO;Mono TP_RAW_NONE;None (Shows sensor pattern) TP_RAW_PIXELSHIFT;Pixel Shift TP_RAW_PIXELSHIFTBLUR;Blur motion mask +TP_RAW_PIXELSHIFTDMETHOD;Demosaic method for motion TP_RAW_PIXELSHIFTEPERISO;Sensitivity TP_RAW_PIXELSHIFTEPERISO_TOOLTIP;The default value of 0 should work fine for base ISO.\nHigher values increase sensitivity of motion detection.\nChange in small steps and watch the motion mask while changing.\nIncrease sensitivity for underexposed or high ISO images. TP_RAW_PIXELSHIFTEQUALBRIGHT;Equalize brightness of frames diff --git a/rtengine/amaze_demosaic_RT.cc b/rtengine/amaze_demosaic_RT.cc index 4fe0bee69..1e1b97768 100644 --- a/rtengine/amaze_demosaic_RT.cc +++ b/rtengine/amaze_demosaic_RT.cc @@ -38,7 +38,7 @@ namespace rtengine { -void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, array2D &rawData, array2D &red, array2D &green, array2D &blue) +void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, const array2D &rawData, array2D &red, array2D &green, array2D &blue) { BENCHFUN diff --git a/rtengine/amaze_vng4_demosaic_RT.cc b/rtengine/amaze_vng4_demosaic_RT.cc index 0dae2c49b..841abf9d3 100644 --- a/rtengine/amaze_vng4_demosaic_RT.cc +++ b/rtengine/amaze_vng4_demosaic_RT.cc @@ -35,22 +35,24 @@ using namespace std; namespace rtengine { -void RawImageSource::amaze_vng4_demosaic_RT(int winw, int winh, array2D &rawData, array2D &red, array2D &green, array2D &blue, double contrast) +void RawImageSource::amaze_vng4_demosaic_RT(int winw, int winh, const array2D &rawData, array2D &red, array2D &green, array2D &blue, double contrast) { BENCHFUN if (contrast == 0.0) { // contrast == 0.0 means only AMaZE will be used - amaze_demosaic_RT (0, 0, winw, winh, rawData, red, green, blue); + amaze_demosaic_RT(0, 0, winw, winh, rawData, red, green, blue); return; } - vng4_demosaic (); + vng4_demosaic(rawData, red, green, blue); + array2D redTmp(winw, winh); array2D greenTmp(winw, winh); array2D blueTmp(winw, winh); array2D L(winw, winh); - amaze_demosaic_RT (0, 0, winw, winh, rawData, redTmp, greenTmp, blueTmp); + + amaze_demosaic_RT(0, 0, winw, winh, rawData, redTmp, greenTmp, blueTmp); const float xyz_rgb[3][3] = { // XYZ from RGB { 0.412453, 0.357580, 0.180423 }, { 0.212671, 0.715160, 0.072169 }, @@ -63,7 +65,7 @@ void RawImageSource::amaze_vng4_demosaic_RT(int winw, int winh, array2D & Color::RGB2L(redTmp[i], greenTmp[i], blueTmp[i], L[i], xyz_rgb, winw); } } - // calculate contrast based blend factors to reduce sharpening in regions with low contrast + // calculate contrast based blend factors to use vng4 in regions with low contrast JaggedArray blend(winw, winh); buildBlendMask(L, blend, winw, winh, contrast / 100.f); diff --git a/rtengine/demosaic_algos.cc b/rtengine/demosaic_algos.cc index 582b00629..d0f4eebd8 100644 --- a/rtengine/demosaic_algos.cc +++ b/rtengine/demosaic_algos.cc @@ -562,7 +562,7 @@ void RawImageSource::hphd_demosaic () //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #define fc(row,col) (prefilters >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3) typedef unsigned short ushort; -void RawImageSource::vng4_demosaic () +void RawImageSource::vng4_demosaic (const array2D &rawData, array2D &red, array2D &green, array2D &blue) { const signed short int *cp, terms[] = { -2, -2, +0, -1, 0, 0x01, -2, -2, +0, +0, 1, 0x01, -2, -1, -1, +0, 0, 0x01, diff --git a/rtengine/demosaic_algos.cc.save-failed b/rtengine/demosaic_algos.cc.save-failed new file mode 100644 index 000000000..29ee67332 --- /dev/null +++ b/rtengine/demosaic_algos.cc.save-failed @@ -0,0 +1,3861 @@ +/* + * This file is part of RawTherapee. + * + * Copyright (c) 2004-2010 Gabor Horvath + * + * RawTherapee is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * RawTherapee is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with RawTherapee. If not, see . + */ +#include +#include + +#include "rawimagesource.h" +#include "rawimagesource_i.h" +#include "jaggedarray.h" +#include "rawimage.h" +#include "mytime.h" +#include "iccmatrices.h" +#include "iccstore.h" +#include "image8.h" +#include "curves.h" +#include "dfmanager.h" +#include "slicer.h" +#include "rt_math.h" +#include "color.h" +#include "../rtgui/multilangmgr.h" +#include "procparams.h" +#include "sleef.c" +#include "opthelper.h" +#include "median.h" +#include "StopWatch.h" +#ifdef _OPENMP +#include +#endif + +using namespace std; + +namespace rtengine +{ +#undef ABS +#undef DIST + +#define ABS(a) ((a)<0?-(a):(a)) +#define DIST(a,b) (ABS(a-b)) +#define CLIREF(x) LIM(x,-200000.0f,200000.0f) // avoid overflow : do not act directly on image[] or pix[] +#define x1125(a) (a + xdivf(a, 3)) +#define x0875(a) (a - xdivf(a, 3)) +#define x0250(a) xdivf(a, 2) +#define x00625(a) xdivf(a, 4) +#define x0125(a) xdivf(a, 3) + +extern const Settings* settings; +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +void RawImageSource::eahd_demosaic () +{ + if (plistener) { + plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::EAHD))); + plistener->setProgress (0.0); + } + + // prepare cache and constants for cielab conversion + //TODO: revisit after conversion to D50 illuminant + lc00 = (0.412453 * imatrices.rgb_cam[0][0] + 0.357580 * imatrices.rgb_cam[0][1] + 0.180423 * imatrices.rgb_cam[0][2]) ;// / 0.950456; + lc01 = (0.412453 * imatrices.rgb_cam[1][0] + 0.357580 * imatrices.rgb_cam[1][1] + 0.180423 * imatrices.rgb_cam[1][2]) ;// / 0.950456; + lc02 = (0.412453 * imatrices.rgb_cam[2][0] + 0.357580 * imatrices.rgb_cam[2][1] + 0.180423 * imatrices.rgb_cam[2][2]) ;// / 0.950456; + + lc10 = 0.212671 * imatrices.rgb_cam[0][0] + 0.715160 * imatrices.rgb_cam[0][1] + 0.072169 * imatrices.rgb_cam[0][2]; + lc11 = 0.212671 * imatrices.rgb_cam[1][0] + 0.715160 * imatrices.rgb_cam[1][1] + 0.072169 * imatrices.rgb_cam[1][2]; + lc12 = 0.212671 * imatrices.rgb_cam[2][0] + 0.715160 * imatrices.rgb_cam[2][1] + 0.072169 * imatrices.rgb_cam[2][2]; + + lc20 = (0.019334 * imatrices.rgb_cam[0][0] + 0.119193 * imatrices.rgb_cam[0][1] + 0.950227 * imatrices.rgb_cam[0][2]) ;// / 1.088754; + lc21 = (0.019334 * imatrices.rgb_cam[1][0] + 0.119193 * imatrices.rgb_cam[1][1] + 0.950227 * imatrices.rgb_cam[1][2]) ;// / 1.088754; + lc22 = (0.019334 * imatrices.rgb_cam[2][0] + 0.119193 * imatrices.rgb_cam[2][1] + 0.950227 * imatrices.rgb_cam[2][2]) ;// / 1.088754; + + int maxindex = 3 * 65536; //2*65536 3 = avoid crash 3/2013 J.Desmis + cache = new double[maxindex]; + threshold = (int)(0.008856 * MAXVALD); + + for (int i = 0; i < maxindex; i++) { + cache[i] = std::cbrt(double(i) / MAXVALD); + } + + // end of cielab preparation + + JaggedArray + rh (W, 3), gh (W, 4), bh (W, 3), + rv (W, 3), gv (W, 4), bv (W, 3), + lLh (W, 3), lah (W, 3), lbh (W, 3), + lLv (W, 3), lav (W, 3), lbv (W, 3), + homh (W, 3), homv (W, 3); + + // interpolate first two lines + interpolate_row_g (gh[0], gv[0], 0); + interpolate_row_g (gh[1], gv[1], 1); + interpolate_row_g (gh[2], gv[2], 2); + interpolate_row_rb (rh[0], bh[0], nullptr, gh[0], gh[1], 0); + interpolate_row_rb (rv[0], bv[0], nullptr, gv[0], gv[1], 0); + interpolate_row_rb (rh[1], bh[1], gh[0], gh[1], gh[2], 1); + interpolate_row_rb (rv[1], bv[1], gv[0], gv[1], gv[2], 1); + + convert_to_cielab_row (rh[0], gh[0], bh[0], lLh[0], lah[0], lbh[0]); + convert_to_cielab_row (rv[0], gv[0], bv[0], lLv[0], lav[0], lbv[0]); + convert_to_cielab_row (rh[1], gh[1], bh[1], lLh[1], lah[1], lbh[1]); + convert_to_cielab_row (rv[1], gv[1], bv[1], lLv[1], lav[1], lbv[1]); + + for (int j = 0; j < W; j++) { + homh[0][j] = 0; + homv[0][j] = 0; + homh[1][j] = 0; + homv[1][j] = 0; + } + + int dLmaph[9]; + int dLmapv[9]; + int dCamaph[9]; + int dCamapv[9]; + int dCbmaph[9]; + int dCbmapv[9]; + + for (int i = 1; i < H - 1; i++) { + int ix = i % 3; + int imx = (i - 1) % 3; + int ipx = (i + 1) % 3; + + if (i < H - 2) { + interpolate_row_g (gh[(i + 2) % 4], gv[(i + 2) % 4], i + 2); + interpolate_row_rb (rh[(i + 1) % 3], bh[(i + 1) % 3], gh[i % 4], gh[(i + 1) % 4], gh[(i + 2) % 4], i + 1); + interpolate_row_rb (rv[(i + 1) % 3], bv[(i + 1) % 3], gv[i % 4], gv[(i + 1) % 4], gv[(i + 2) % 4], i + 1); + } else { + interpolate_row_rb (rh[(i + 1) % 3], bh[(i + 1) % 3], gh[i % 4], gh[(i + 1) % 4], nullptr, i + 1); + interpolate_row_rb (rv[(i + 1) % 3], bv[(i + 1) % 3], gv[i % 4], gv[(i + 1) % 4], nullptr, i + 1); + } + + convert_to_cielab_row (rh[(i + 1) % 3], gh[(i + 1) % 4], bh[(i + 1) % 3], lLh[(i + 1) % 3], lah[(i + 1) % 3], lbh[(i + 1) % 3]); + convert_to_cielab_row (rv[(i + 1) % 3], gv[(i + 1) % 4], bv[(i + 1) % 3], lLv[(i + 1) % 3], lav[(i + 1) % 3], lbv[(i + 1) % 3]); + + for (int j = 0; j < W; j++) { + homh[ipx][j] = 0; + homv[ipx][j] = 0; + } + + int sh, sv, idx; + + for (int j = 1; j < W - 1; j++) { + int dmi = 0; + + for (int x = -1; x <= 1; x++) { + idx = (i + x) % 3; + + for (int y = -1; y <= 1; y++) { + // compute distance in a, b, and L + if (dmi < 4) { + sh = homh[idx][j + y]; + sv = homv[idx][j + y]; + + if (sh > sv) { // fixate horizontal pixel + dLmaph[dmi] = DIST(lLh[ix][j], lLh[idx][j + y]); + dCamaph[dmi] = DIST(lah[ix][j], lah[idx][j + y]); + dCbmaph[dmi] = DIST(lbh[ix][j], lbh[idx][j + y]); + dLmapv[dmi] = DIST(lLv[ix][j], lLh[idx][j + y]); + dCamapv[dmi] = DIST(lav[ix][j], lah[idx][j + y]); + dCbmapv[dmi] = DIST(lbv[ix][j], lbh[idx][j + y]); + } else if (sh < sv) { + dLmaph[dmi] = DIST(lLh[ix][j], lLv[idx][j + y]); + dCamaph[dmi] = DIST(lah[ix][j], lav[idx][j + y]); + dCbmaph[dmi] = DIST(lbh[ix][j], lbv[idx][j + y]); + dLmapv[dmi] = DIST(lLv[ix][j], lLv[idx][j + y]); + dCamapv[dmi] = DIST(lav[ix][j], lav[idx][j + y]); + dCbmapv[dmi] = DIST(lbv[ix][j], lbv[idx][j + y]); + } else { + dLmaph[dmi] = DIST(lLh[ix][j], lLh[idx][j + y]); + dCamaph[dmi] = DIST(lah[ix][j], lah[idx][j + y]); + dCbmaph[dmi] = DIST(lbh[ix][j], lbh[idx][j + y]); + dLmapv[dmi] = DIST(lLv[ix][j], lLv[idx][j + y]); + dCamapv[dmi] = DIST(lav[ix][j], lav[idx][j + y]); + dCbmapv[dmi] = DIST(lbv[ix][j], lbv[idx][j + y]); + } + } else { + dLmaph[dmi] = DIST(lLh[ix][j], lLh[idx][j + y]); + dCamaph[dmi] = DIST(lah[ix][j], lah[idx][j + y]); + dCbmaph[dmi] = DIST(lbh[ix][j], lbh[idx][j + y]); + dLmapv[dmi] = DIST(lLv[ix][j], lLv[idx][j + y]); + dCamapv[dmi] = DIST(lav[ix][j], lav[idx][j + y]); + dCbmapv[dmi] = DIST(lbv[ix][j], lbv[idx][j + y]); + } + + dmi++; + } + } + + // compute eL & eC + int eL = min(max(dLmaph[3], dLmaph[5]), max(dLmapv[1], dLmapv[7])); + int eCa = min(max(dCamaph[3], dCamaph[5]), max(dCamapv[1], dCamapv[7])); + int eCb = min(max(dCbmaph[3], dCbmaph[5]), max(dCbmapv[1], dCbmapv[7])); + + int wh = 0; + + for (int dmi = 0; dmi < 9; dmi++) + if (dLmaph[dmi] <= eL && dCamaph[dmi] <= eCa && dCbmaph[dmi] <= eCb) { + wh++; + } + + int wv = 0; + + for (int dmi = 0; dmi < 9; dmi++) + if (dLmapv[dmi] <= eL && dCamapv[dmi] <= eCa && dCbmapv[dmi] <= eCb) { + wv++; + } + + homh[imx][j - 1] += wh; + homh[imx][j] += wh; + homh[imx][j + 1] += wh; + homh[ix][j - 1] += wh; + homh[ix][j] += wh; + homh[ix][j + 1] += wh; + homh[ipx][j - 1] += wh; + homh[ipx][j] += wh; + homh[ipx][j + 1] += wh; + + homv[imx][j - 1] += wv; + homv[imx][j] += wv; + homv[imx][j + 1] += wv; + homv[ix][j - 1] += wv; + homv[ix][j] += wv; + homv[ix][j + 1] += wv; + homv[ipx][j - 1] += wv; + homv[ipx][j] += wv; + homv[ipx][j + 1] += wv; + } + +//} + // finalize image + int hc, vc; + + for (int j = 0; j < W; j++) { + if (ri->ISGREEN(i - 1, j)) { + green[i - 1][j] = rawData[i - 1][j]; + } else { + hc = homh[imx][j]; + vc = homv[imx][j]; + + if (hc > vc) { + green[i - 1][j] = gh[(i - 1) % 4][j]; + } else if (hc < vc) { + green[i - 1][j] = gv[(i - 1) % 4][j]; + } else { + green[i - 1][j] = (gh[(i - 1) % 4][j] + gv[(i - 1) % 4][j]) / 2; + } + } + } + + if (!(i % 20) && plistener) { + plistener->setProgress ((double)i / (H - 2)); + } + } + + // finish H-2th and H-1th row, homogenity value is still valailable + int hc, vc; + + for (int i = H - 1; i < H + 1; i++) + for (int j = 0; j < W; j++) { + hc = homh[(i - 1) % 3][j]; + vc = homv[(i - 1) % 3][j]; + + if (hc > vc) { + green[i - 1][j] = gh[(i - 1) % 4][j]; + } else if (hc < vc) { + green[i - 1][j] = gv[(i - 1) % 4][j]; + } else { + green[i - 1][j] = (gh[(i - 1) % 4][j] + gv[(i - 1) % 4][j]) / 2; + } + } + + // Interpolate R and B + for (int i = 0; i < H; i++) { + if (i == 0) { + interpolate_row_rb_mul_pp (red[i], blue[i], nullptr, green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1); + } else if (i == H - 1) { + interpolate_row_rb_mul_pp (red[i], blue[i], green[i - 1], green[i], nullptr, i, 1.0, 1.0, 1.0, 0, W, 1); + } else { + interpolate_row_rb_mul_pp (red[i], blue[i], green[i - 1], green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1); + } + } +} + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +void RawImageSource::hphd_vertical (float** hpmap, int col_from, int col_to) +{ + float* temp = new float[max(W, H)]; + float* avg = new float[max(W, H)]; + float* dev = new float[max(W, H)]; + + memset (temp, 0, max(W, H)*sizeof(float)); + memset (avg, 0, max(W, H)*sizeof(float)); + memset (dev, 0, max(W, H)*sizeof(float)); + + for (int k = col_from; k < col_to; k++) { + for (int i = 5; i < H - 5; i++) { + temp[i] = (rawData[i - 5][k] - 8 * rawData[i - 4][k] + 27 * rawData[i - 3][k] - 48 * rawData[i - 2][k] + 42 * rawData[i - 1][k] - + (rawData[i + 5][k] - 8 * rawData[i + 4][k] + 27 * rawData[i + 3][k] - 48 * rawData[i + 2][k] + 42 * rawData[i + 1][k])) / 100.0; + temp[i] = ABS(temp[i]); + } + + for (int j = 4; j < H - 4; j++) { + float avgL = (temp[j - 4] + temp[j - 3] + temp[j - 2] + temp[j - 1] + temp[j] + temp[j + 1] + temp[j + 2] + temp[j + 3] + temp[j + 4]) / 9.0; + avg[j] = avgL; + float devL = ((temp[j - 4] - avgL) * (temp[j - 4] - avgL) + (temp[j - 3] - avgL) * (temp[j - 3] - avgL) + (temp[j - 2] - avgL) * (temp[j - 2] - avgL) + (temp[j - 1] - avgL) * (temp[j - 1] - avgL) + (temp[j] - avgL) * (temp[j] - avgL) + (temp[j + 1] - avgL) * (temp[j + 1] - avgL) + (temp[j + 2] - avgL) * (temp[j + 2] - avgL) + (temp[j + 3] - avgL) * (temp[j + 3] - avgL) + (temp[j + 4] - avgL) * (temp[j + 4] - avgL)) / 9.0; + + if (devL < 0.001) { + devL = 0.001; + } + + dev[j] = devL; + } + + for (int j = 5; j < H - 5; j++) { + float avgL = avg[j - 1]; + float avgR = avg[j + 1]; + float devL = dev[j - 1]; + float devR = dev[j + 1]; + hpmap[j][k] = avgL + (avgR - avgL) * devL / (devL + devR); + } + } + + delete [] temp; + delete [] avg; + delete [] dev; +} + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +void RawImageSource::hphd_horizontal (float** hpmap, int row_from, int row_to) +{ + float* temp = new float[max(W, H)]; + float* avg = new float[max(W, H)]; + float* dev = new float[max(W, H)]; + + memset (temp, 0, max(W, H)*sizeof(float)); + memset (avg, 0, max(W, H)*sizeof(float)); + memset (dev, 0, max(W, H)*sizeof(float)); + + for (int i = row_from; i < row_to; i++) { + for (int j = 5; j < W - 5; j++) { + temp[j] = (rawData[i][j - 5] - 8 * rawData[i][j - 4] + 27 * rawData[i][j - 3] - 48 * rawData[i][j - 2] + 42 * rawData[i][j - 1] - + (rawData[i][j + 5] - 8 * rawData[i][j + 4] + 27 * rawData[i][j + 3] - 48 * rawData[i][j + 2] + 42 * rawData[i][j + 1])) / 100; + temp[j] = ABS(temp[j]); + } + + for (int j = 4; j < W - 4; j++) { + float avgL = (temp[j - 4] + temp[j - 3] + temp[j - 2] + temp[j - 1] + temp[j] + temp[j + 1] + temp[j + 2] + temp[j + 3] + temp[j + 4]) / 9.0; + avg[j] = avgL; + float devL = ((temp[j - 4] - avgL) * (temp[j - 4] - avgL) + (temp[j - 3] - avgL) * (temp[j - 3] - avgL) + (temp[j - 2] - avgL) * (temp[j - 2] - avgL) + (temp[j - 1] - avgL) * (temp[j - 1] - avgL) + (temp[j] - avgL) * (temp[j] - avgL) + (temp[j + 1] - avgL) * (temp[j + 1] - avgL) + (temp[j + 2] - avgL) * (temp[j + 2] - avgL) + (temp[j + 3] - avgL) * (temp[j + 3] - avgL) + (temp[j + 4] - avgL) * (temp[j + 4] - avgL)) / 9.0; + + if (devL < 0.001) { + devL = 0.001; + } + + dev[j] = devL; + } + + for (int j = 5; j < W - 5; j++) { + float avgL = avg[j - 1]; + float avgR = avg[j + 1]; + float devL = dev[j - 1]; + float devR = dev[j + 1]; + float hpv = avgL + (avgR - avgL) * devL / (devL + devR); + + if (hpmap[i][j] < 0.8 * hpv) { + hpmap[i][j] = 2; + } else if (hpv < 0.8 * hpmap[i][j]) { + hpmap[i][j] = 1; + } else { + hpmap[i][j] = 0; + } + } + } + + delete [] temp; + delete [] avg; + delete [] dev; +} + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +void RawImageSource::hphd_green (float** hpmap) +{ +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int i = 3; i < H - 3; i++) { + for (int j = 3; j < W - 3; j++) { + if (ri->ISGREEN(i, j)) { + green[i][j] = rawData[i][j]; + } else { + if (hpmap[i][j] == 1) { + int g2 = rawData[i][j + 1] + ((rawData[i][j] - rawData[i][j + 2]) / 2); + int g4 = rawData[i][j - 1] + ((rawData[i][j] - rawData[i][j - 2]) / 2); + + int dx = rawData[i][j + 1] - rawData[i][j - 1]; + int d1 = rawData[i][j + 3] - rawData[i][j + 1]; + int d2 = rawData[i][j + 2] - rawData[i][j]; + int d3 = (rawData[i - 1][j + 2] - rawData[i - 1][j]) / 2; + int d4 = (rawData[i + 1][j + 2] - rawData[i + 1][j]) / 2; + + double e2 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + d1 = rawData[i][j - 3] - rawData[i][j - 1]; + d2 = rawData[i][j - 2] - rawData[i][j]; + d3 = (rawData[i - 1][j - 2] - rawData[i - 1][j]) / 2; + d4 = (rawData[i + 1][j - 2] - rawData[i + 1][j]) / 2; + + double e4 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + green[i][j] = (e2 * g2 + e4 * g4) / (e2 + e4); + } else if (hpmap[i][j] == 2) { + int g1 = rawData[i - 1][j] + ((rawData[i][j] - rawData[i - 2][j]) / 2); + int g3 = rawData[i + 1][j] + ((rawData[i][j] - rawData[i + 2][j]) / 2); + + int dy = rawData[i + 1][j] - rawData[i - 1][j]; + int d1 = rawData[i - 1][j] - rawData[i - 3][j]; + int d2 = rawData[i][j] - rawData[i - 2][j]; + int d3 = (rawData[i][j - 1] - rawData[i - 2][j - 1]) / 2; + int d4 = (rawData[i][j + 1] - rawData[i - 2][j + 1]) / 2; + + double e1 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + d1 = rawData[i + 1][j] - rawData[i + 3][j]; + d2 = rawData[i][j] - rawData[i + 2][j]; + d3 = (rawData[i][j - 1] - rawData[i + 2][j - 1]) / 2; + d4 = (rawData[i][j + 1] - rawData[i + 2][j + 1]) / 2; + + double e3 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + green[i][j] = (e1 * g1 + e3 * g3) / (e1 + e3); + } else { + int g1 = rawData[i - 1][j] + ((rawData[i][j] - rawData[i - 2][j]) / 2); + int g2 = rawData[i][j + 1] + ((rawData[i][j] - rawData[i][j + 2]) / 2); + int g3 = rawData[i + 1][j] + ((rawData[i][j] - rawData[i + 2][j]) / 2); + int g4 = rawData[i][j - 1] + ((rawData[i][j] - rawData[i][j - 2]) / 2); + + int dx = rawData[i][j + 1] - rawData[i][j - 1]; + int dy = rawData[i + 1][j] - rawData[i - 1][j]; + + int d1 = rawData[i - 1][j] - rawData[i - 3][j]; + int d2 = rawData[i][j] - rawData[i - 2][j]; + int d3 = (rawData[i][j - 1] - rawData[i - 2][j - 1]) / 2; + int d4 = (rawData[i][j + 1] - rawData[i - 2][j + 1]) / 2; + + double e1 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + d1 = rawData[i][j + 3] - rawData[i][j + 1]; + d2 = rawData[i][j + 2] - rawData[i][j]; + d3 = (rawData[i - 1][j + 2] - rawData[i - 1][j]) / 2; + d4 = (rawData[i + 1][j + 2] - rawData[i + 1][j]) / 2; + + double e2 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + d1 = rawData[i + 1][j] - rawData[i + 3][j]; + d2 = rawData[i][j] - rawData[i + 2][j]; + d3 = (rawData[i][j - 1] - rawData[i + 2][j - 1]) / 2; + d4 = (rawData[i][j + 1] - rawData[i + 2][j + 1]) / 2; + + double e3 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + d1 = rawData[i][j - 3] - rawData[i][j - 1]; + d2 = rawData[i][j - 2] - rawData[i][j]; + d3 = (rawData[i - 1][j - 2] - rawData[i - 1][j]) / 2; + d4 = (rawData[i + 1][j - 2] - rawData[i + 1][j]) / 2; + + double e4 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + green[i][j] = (e1 * g1 + e2 * g2 + e3 * g3 + e4 * g4) / (e1 + e2 + e3 + e4); + } + } + } + } +} + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +void RawImageSource::hphd_demosaic () +{ + if (plistener) { + plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::HPHD))); + plistener->setProgress (0.0); + } + + JaggedArray hpmap (W, H, true); + +#ifdef _OPENMP + #pragma omp parallel + { + int tid = omp_get_thread_num(); + int nthreads = omp_get_num_threads(); + int blk = W / nthreads; + + if (tid < nthreads - 1) { + hphd_vertical (hpmap, tid * blk, (tid + 1)*blk); + } else { + hphd_vertical (hpmap, tid * blk, W); + } + } +#else + hphd_vertical (hpmap, 0, W); +#endif + + if (plistener) { + plistener->setProgress (0.33); + } + +#ifdef _OPENMP + #pragma omp parallel + { + int tid = omp_get_thread_num(); + int nthreads = omp_get_num_threads(); + int blk = H / nthreads; + + if (tid < nthreads - 1) { + hphd_horizontal (hpmap, tid * blk, (tid + 1)*blk); + } else { + hphd_horizontal (hpmap, tid * blk, H); + } + } +#else + hphd_horizontal (hpmap, 0, H); +#endif + + hphd_green (hpmap); + + if (plistener) { + plistener->setProgress (0.66); + } + + for (int i = 0; i < H; i++) { + if (i == 0) { + interpolate_row_rb_mul_pp (red[i], blue[i], nullptr, green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1); + } else if (i == H - 1) { + interpolate_row_rb_mul_pp (red[i], blue[i], green[i - 1], green[i], nullptr, i, 1.0, 1.0, 1.0, 0, W, 1); + } else { + interpolate_row_rb_mul_pp (red[i], blue[i], green[i - 1], green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1); + } + } + + if (plistener) { + plistener->setProgress (1.0); + } +} + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +#define fc(row,col) (prefilters >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3) +typedef unsigned short ushort; +void RawImageSource::vng4_demosaic () +{ + const signed short int *cp, terms[] = { + -2, -2, +0, -1, 0, 0x01, -2, -2, +0, +0, 1, 0x01, -2, -1, -1, +0, 0, 0x01, + -2, -1, +0, -1, 0, 0x02, -2, -1, +0, +0, 0, 0x03, -2, -1, +0, +1, 1, 0x01, + -2, +0, +0, -1, 0, 0x06, -2, +0, +0, +0, 1, 0x02, -2, +0, +0, +1, 0, 0x03, + -2, +1, -1, +0, 0, 0x04, -2, +1, +0, -1, 1, 0x04, -2, +1, +0, +0, 0, 0x06, + -2, +1, +0, +1, 0, 0x02, -2, +2, +0, +0, 1, 0x04, -2, +2, +0, +1, 0, 0x04, + -1, -2, -1, +0, 0, 0x80, -1, -2, +0, -1, 0, 0x01, -1, -2, +1, -1, 0, 0x01, + -1, -2, +1, +0, 1, 0x01, -1, -1, -1, +1, 0, 0x88, -1, -1, +1, -2, 0, 0x40, + -1, -1, +1, -1, 0, 0x22, -1, -1, +1, +0, 0, 0x33, -1, -1, +1, +1, 1, 0x11, + -1, +0, -1, +2, 0, 0x08, -1, +0, +0, -1, 0, 0x44, -1, +0, +0, +1, 0, 0x11, + -1, +0, +1, -2, 1, 0x40, -1, +0, +1, -1, 0, 0x66, -1, +0, +1, +0, 1, 0x22, + -1, +0, +1, +1, 0, 0x33, -1, +0, +1, +2, 1, 0x10, -1, +1, +1, -1, 1, 0x44, + -1, +1, +1, +0, 0, 0x66, -1, +1, +1, +1, 0, 0x22, -1, +1, +1, +2, 0, 0x10, + -1, +2, +0, +1, 0, 0x04, -1, +2, +1, +0, 1, 0x04, -1, +2, +1, +1, 0, 0x04, + +0, -2, +0, +0, 1, 0x80, +0, -1, +0, +1, 1, 0x88, +0, -1, +1, -2, 0, 0x40, + +0, -1, +1, +0, 0, 0x11, +0, -1, +2, -2, 0, 0x40, +0, -1, +2, -1, 0, 0x20, + +0, -1, +2, +0, 0, 0x30, +0, -1, +2, +1, 1, 0x10, +0, +0, +0, +2, 1, 0x08, + +0, +0, +2, -2, 1, 0x40, +0, +0, +2, -1, 0, 0x60, +0, +0, +2, +0, 1, 0x20, + +0, +0, +2, +1, 0, 0x30, +0, +0, +2, +2, 1, 0x10, +0, +1, +1, +0, 0, 0x44, + +0, +1, +1, +2, 0, 0x10, +0, +1, +2, -1, 1, 0x40, +0, +1, +2, +0, 0, 0x60, + +0, +1, +2, +1, 0, 0x20, +0, +1, +2, +2, 0, 0x10, +1, -2, +1, +0, 0, 0x80, + +1, -1, +1, +1, 0, 0x88, +1, +0, +1, +2, 0, 0x08, +1, +0, +2, -1, 0, 0x40, + +1, +0, +2, +1, 0, 0x10 + }, + chood[] = { -1, -1, -1, 0, -1, +1, 0, +1, +1, +1, +1, 0, +1, -1, 0, -1 }; + + double progress = 0.0; + const bool plistenerActive = plistener; + + if (plistenerActive) { + plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::VNG4))); + plistener->setProgress (progress); + } + + const unsigned prefilters = ri->prefilters; + const int width = W, height = H; + constexpr unsigned int colors = 4; + float (*image)[4]; + + image = (float (*)[4]) calloc (height * width, sizeof * image); + +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int ii = 0; ii < H; ii++) + for (int jj = 0; jj < W; jj++) { + image[ii * W + jj][fc(ii, jj)] = rawData[ii][jj]; + } + + { + int lcode[16][16][32]; + float mul[16][16][8]; + float csum[16][16][4]; + +// first linear interpolation + for (int row = 0; row < 16; row++) + for (int col = 0; col < 16; col++) { + int * ip = lcode[row][col]; + int mulcount = 0; + float sum[4]; + memset (sum, 0, sizeof sum); + + for (int y = -1; y <= 1; y++) + for (int x = -1; x <= 1; x++) { + int shift = (y == 0) + (x == 0); + + if (shift == 2) { + continue; + } + + int color = fc(row + y, col + x); + *ip++ = (width * y + x) * 4 + color; + + mul[row][col][mulcount] = (1 << shift); + *ip++ = color; + sum[color] += (1 << shift); + mulcount++; + } + + int colcount = 0; + + for (unsigned int c = 0; c < colors; c++) + if (c != fc(row, col)) { + *ip++ = c; + csum[row][col][colcount] = sum[c]; + colcount ++; + } + } + +StopWatch Stop1("loop 1"); +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int row = 1; row < height - 1; row++) { + for (int col = 1; col < width - 1; col++) { + float * pix = image[row * width + col]; + int * ip = lcode[row & 15][col & 15]; + float sum[4]; + memset (sum, 0, sizeof sum); + + for (int i = 0; i < 8; i++, ip += 2) { + sum[ip[1]] += pix[ip[0]] * mul[row & 15][col & 15][i]; + } + + for (unsigned int i = 0; i < colors - 1; i++, ip++) { + pix[ip[0]] = sum[ip[0]] / csum[row & 15][col & 15][i]; + } + } + } +Stop1.stop(); + } + const int prow = 7, pcol = 1; + int *code[8][2]; + int t, g; + int * ip = (int *) calloc ((prow + 1) * (pcol + 1), 1280); + + for (int row = 0; row <= prow; row++) /* Precalculate for VNG */ + for (int col = 0; col <= pcol; col++) { + code[row][col] = ip; + + for (cp = terms, t = 0; t < 64; t++) { + int y1 = *cp++; + int x1 = *cp++; + int y2 = *cp++; + int x2 = *cp++; + int weight = *cp++; + int grads = *cp++; + unsigned int color = fc(row + y1, col + x1); + + if (fc(row + y2, col + x2) != color) { + continue; + } + + int diag = (fc(row, col + 1) == color && fc(row + 1, col) == color) ? 2 : 1; + + if (abs(y1 - y2) == diag && abs(x1 - x2) == diag) { + continue; + } + + *ip++ = (y1 * width + x1) * 4 + color; + *ip++ = (y2 * width + x2) * 4 + color; + *ip++ = weight; + + for (g = 0; g < 8; g++) + if (grads & (1 << g)) { + *ip++ = g; + } + + *ip++ = -1; + } + + *ip++ = INT_MAX; + + for (cp = chood, g = 0; g < 8; g++) { + int y = *cp++; + int x = *cp++; + *ip++ = (y * width + x) * 4; + unsigned int color = fc(row, col); + + if (fc(row + y, col + x) != color && fc(row + y * 2, col + x * 2) == color) { + *ip++ = (y * width + x) * 8 + color; + } else { + *ip++ = 0; + } + } + } + + if(plistenerActive) { + progress = 0.1; + plistener->setProgress (progress); + } + +StopWatch Stop2("loop 2"); + +#ifdef _OPENMP + #pragma omp parallel +#endif + { + float thold; + int g; + constexpr int progressStep = 64; + const double progressInc = (0.98 - progress) / ((height - 2) / progressStep); +#ifdef _OPENMP + #pragma omp for schedule(dynamic,16) nowait +#endif + + for (int row = 2; row < height - 2; row++) { /* Do VNG interpolation */ + for (int col = 2; col < width - 2; col++) { + float * pix = image[row * width + col]; + int * ip = code[row & prow][col & pcol]; + float gval[8] = {}; +// gval[0] = gval[1] = gval[2] = gval[3] = gval[4] = gval[5] = gval[6] = gval[7] = 0.f; +// memset (gval, 0, sizeof gval); + + while ((g = ip[0]) != INT_MAX) { /* Calculate gradients */ + float diff = fabsf(pix[g] - pix[ip[1]]) * (1 << ip[2]); + gval[ip[3]] += diff; + ip += 4; + + while ((g = *ip++) != -1) { + gval[g] += diff; + } + } + + ip++; + { + float gmin, gmax; + gmin = gmax = gval[0]; /* Choose a threshold */ + + for (g = 1; g < 8; g++) { + gmin = rtengine::min(gmin, gval[g]); +// if (gmin > gval[g]) { +// gmin = gval[g]; +// } + + gmax = rtengine::max(gmax, gval[g]); +// if (gmax < gval[g]) { +// gmax = gval[g]; +// } + } + + thold = gmin + (gmax / 2); + } + float sum[3] = {}; +// sum[0] = sum[1] = sum[2] = 0.f; +// memset (sum, 0, sizeof sum); + float t1, t2; + int color = fc(row, col); + t1 = t2 = pix[color]; + + if(color & 1) { + int num = 0; + + for (g = 0; g < 8; g++, ip += 2) { /* Average the neighbors */ + if (gval[g] <= thold) { + if(ip[1]) { + sum[0] += (t1 + pix[ip[1]]) * 0.5f; + } + + sum[1] += pix[ip[0] + (color ^ 2)]; + num++; + } + } + + t1 += (sum[1] - sum[0]) / num; + } else { + int num = 0; + + for (g = 0; g < 8; g++, ip += 2) { /* Average the neighbors */ + if (gval[g] <= thold) { + sum[1] += pix[ip[0] + 1]; + sum[2] += pix[ip[0] + 3]; + + if(ip[1]) { + sum[0] += (t1 + pix[ip[1]]) * 0.5f; + } + + num++; + } + } + + t1 += (sum[1] - sum[0]) / num; + t2 += (sum[2] - sum[0]) / num; + } + + green[row][col] = 0.5f * (t1 + t2); + } + + if(plistenerActive) { + if((row % progressStep) == 0) +#ifdef _OPENMP + #pragma omp critical (updateprogress) +#endif + { + progress += progressInc; + plistener->setProgress (progress); + } + } + } + + } +Stop2.stop(); + free (code[0][0]); + free (image); + + if(plistenerActive) { + plistener->setProgress (0.98); + } +StopWatch Stop3("loop 3"); + + // Interpolate R and B +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int i = 0; i < H; i++) { + if (i == 0) + // rm, gm, bm must be recovered + //interpolate_row_rb_mul_pp (red, blue, NULL, green[i], green[i+1], i, rm, gm, bm, 0, W, 1); + { + interpolate_row_rb_mul_pp (red[i], blue[i], nullptr, green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1); + } else if (i == H - 1) { + interpolate_row_rb_mul_pp (red[i], blue[i], green[i - 1], green[i], nullptr, i, 1.0, 1.0, 1.0, 0, W, 1); + } else { + interpolate_row_rb_mul_pp (red[i], blue[i], green[i - 1], green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1); + } + } +Stop3.stop(); + if(plistenerActive) { + plistener->setProgress (1.0); + } +} + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +#undef fc +#define fc(row,col) \ + (ri->get_filters() >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3) + +#define FORCC for (unsigned int c=0; c < colors; c++) + +/* + Patterned Pixel Grouping Interpolation by Alain Desbiolles +*/ +void RawImageSource::ppg_demosaic() +{ + int width = W, height = H; + int dir[5] = { 1, width, -1, -width, 1 }; + int row, col, diff[2] = {}, guess[2], c, d, i; + float (*pix)[4]; + + float (*image)[4]; + + if (plistener) { + // looks like ppg isn't supported anymore + //plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::ppg))); + plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), "xxx")); + plistener->setProgress (0.0); + } + + image = (float (*)[4]) calloc (H * W, sizeof * image); + + for (int ii = 0; ii < H; ii++) + for (int jj = 0; jj < W; jj++) { + image[ii * W + jj][fc(ii, jj)] = rawData[ii][jj]; + } + + border_interpolate(3, image); + + /* Fill in the green layer with gradients and pattern recognition: */ + for (row = 3; row < height - 3; row++) { + for (col = 3 + (FC(row, 3) & 1), c = FC(row, col); col < width - 3; col += 2) { + pix = image + row * width + col; + + for (i = 0; (d = dir[i]) > 0; i++) { + guess[i] = (pix[-d][1] + pix[0][c] + pix[d][1]) * 2 + - pix[-2 * d][c] - pix[2 * d][c]; + diff[i] = ( ABS(pix[-2 * d][c] - pix[ 0][c]) + + ABS(pix[ 2 * d][c] - pix[ 0][c]) + + ABS(pix[ -d][1] - pix[ d][1]) ) * 3 + + ( ABS(pix[ 3 * d][1] - pix[ d][1]) + + ABS(pix[-3 * d][1] - pix[-d][1]) ) * 2; + } + + d = dir[i = diff[0] > diff[1]]; + pix[0][1] = median(static_cast(guess[i] >> 2), pix[d][1], pix[-d][1]); + } + + if(plistener) { + plistener->setProgress(0.33 * row / (height - 3)); + } + } + + /* Calculate red and blue for each green pixel: */ + for (row = 1; row < height - 1; row++) { + for (col = 1 + (FC(row, 2) & 1), c = FC(row, col + 1); col < width - 1; col += 2) { + pix = image + row * width + col; + + for (i = 0; (d = dir[i]) > 0; c = 2 - c, i++) + pix[0][c] = CLIP(0.5 * (pix[-d][c] + pix[d][c] + 2 * pix[0][1] + - pix[-d][1] - pix[d][1]) ); + } + + if(plistener) { + plistener->setProgress(0.33 + 0.33 * row / (height - 1)); + } + } + + /* Calculate blue for red pixels and vice versa: */ + for (row = 1; row < height - 1; row++) { + for (col = 1 + (FC(row, 1) & 1), c = 2 - FC(row, col); col < width - 1; col += 2) { + pix = image + row * width + col; + + for (i = 0; (d = dir[i] + dir[i + 1]) > 0; i++) { + diff[i] = ABS(pix[-d][c] - pix[d][c]) + + ABS(pix[-d][1] - pix[0][1]) + + ABS(pix[ d][1] - pix[0][1]); + guess[i] = pix[-d][c] + pix[d][c] + 2 * pix[0][1] + - pix[-d][1] - pix[d][1]; + } + + if (diff[0] != diff[1]) { + pix[0][c] = CLIP(guess[diff[0] > diff[1]] / 2); + } else { + pix[0][c] = CLIP((guess[0] + guess[1]) / 4); + } + } + + if(plistener) { + plistener->setProgress(0.67 + 0.33 * row / (height - 1)); + } + } + + red(W, H); + + for (int i = 0; i < H; i++) + for (int j = 0; j < W; j++) { + red[i][j] = image[i * W + j][0]; + } + + green(W, H); + + for (int i = 0; i < H; i++) + for (int j = 0; j < W; j++) { + green[i][j] = image[i * W + j][1]; + } + + blue(W, H); + + for (int i = 0; i < H; i++) + for (int j = 0; j < W; j++) { + blue[i][j] = image[i * W + j][2]; + } + + free (image); +} + +void RawImageSource::border_interpolate(unsigned int border, float (*image)[4], unsigned int start, unsigned int end) +{ + unsigned row, col, y, x, f, sum[8]; + unsigned int width = W, height = H; + unsigned int colors = 3; + + if (end == 0 ) { + end = H; + } + + for (row = start; row < end; row++) + for (col = 0; col < width; col++) { + if (col == border && row >= border && row < height - border) { + col = width - border; + } + + memset (sum, 0, sizeof sum); + + for (y = row - 1; y != row + 2; y++) + for (x = col - 1; x != col + 2; x++) + if (y < height && x < width) { + f = fc(y, x); + sum[f] += image[y * width + x][f]; + sum[f + 4]++; + } + + f = fc(row, col); + + FORCC if (c != f && sum[c + 4]) { + image[row * width + col][c] = sum[c] / sum[c + 4]; + } + } +} + +void RawImageSource::border_interpolate2( int winw, int winh, int lborders) +{ + int bord = lborders; + int width = winw; + int height = winh; + + for (int i = 0; i < height; i++) { + + float sum[6]; + + for (int j = 0; j < bord; j++) { //first few columns + for (int c = 0; c < 6; c++) { + sum[c] = 0; + } + + for (int i1 = i - 1; i1 < i + 2; i1++) + for (int j1 = j - 1; j1 < j + 2; j1++) { + if ((i1 > -1) && (i1 < height) && (j1 > -1)) { + int c = FC(i1, j1); + sum[c] += rawData[i1][j1]; + sum[c + 3]++; + } + } + + int c = FC(i, j); + + if (c == 1) { + red[i][j] = sum[0] / sum[3]; + green[i][j] = rawData[i][j]; + blue[i][j] = sum[2] / sum[5]; + } else { + green[i][j] = sum[1] / sum[4]; + + if (c == 0) { + red[i][j] = rawData[i][j]; + blue[i][j] = sum[2] / sum[5]; + } else { + red[i][j] = sum[0] / sum[3]; + blue[i][j] = rawData[i][j]; + } + } + }//j + + for (int j = width - bord; j < width; j++) { //last few columns + for (int c = 0; c < 6; c++) { + sum[c] = 0; + } + + for (int i1 = i - 1; i1 < i + 2; i1++) + for (int j1 = j - 1; j1 < j + 2; j1++) { + if ((i1 > -1) && (i1 < height ) && (j1 < width)) { + int c = FC(i1, j1); + sum[c] += rawData[i1][j1]; + sum[c + 3]++; + } + } + + int c = FC(i, j); + + if (c == 1) { + red[i][j] = sum[0] / sum[3]; + green[i][j] = rawData[i][j]; + blue[i][j] = sum[2] / sum[5]; + } else { + green[i][j] = sum[1] / sum[4]; + + if (c == 0) { + red[i][j] = rawData[i][j]; + blue[i][j] = sum[2] / sum[5]; + } else { + red[i][j] = sum[0] / sum[3]; + blue[i][j] = rawData[i][j]; + } + } + }//j + }//i + + for (int i = 0; i < bord; i++) { + + float sum[6]; + + for (int j = bord; j < width - bord; j++) { //first few rows + for (int c = 0; c < 6; c++) { + sum[c] = 0; + } + + for (int i1 = i - 1; i1 < i + 2; i1++) + for (int j1 = j - 1; j1 < j + 2; j1++) { + if ((i1 > -1) && (i1 < height) && (j1 > -1)) { + int c = FC(i1, j1); + sum[c] += rawData[i1][j1]; + sum[c + 3]++; + } + } + + int c = FC(i, j); + + if (c == 1) { + red[i][j] = sum[0] / sum[3]; + green[i][j] = rawData[i][j]; + blue[i][j] = sum[2] / sum[5]; + } else { + green[i][j] = sum[1] / sum[4]; + + if (c == 0) { + red[i][j] = rawData[i][j]; + blue[i][j] = sum[2] / sum[5]; + } else { + red[i][j] = sum[0] / sum[3]; + blue[i][j] = rawData[i][j]; + } + } + }//j + } + + for (int i = height - bord; i < height; i++) { + + float sum[6]; + + for (int j = bord; j < width - bord; j++) { //last few rows + for (int c = 0; c < 6; c++) { + sum[c] = 0; + } + + for (int i1 = i - 1; i1 < i + 2; i1++) + for (int j1 = j - 1; j1 < j + 2; j1++) { + if ((i1 > -1) && (i1 < height) && (j1 < width)) { + int c = FC(i1, j1); + sum[c] += rawData[i1][j1]; + sum[c + 3]++; + } + } + + int c = FC(i, j); + + if (c == 1) { + red[i][j] = sum[0] / sum[3]; + green[i][j] = rawData[i][j]; + blue[i][j] = sum[2] / sum[5]; + } else { + green[i][j] = sum[1] / sum[4]; + + if (c == 0) { + red[i][j] = rawData[i][j]; + blue[i][j] = sum[2] / sum[5]; + } else { + red[i][j] = sum[0] / sum[3]; + blue[i][j] = rawData[i][j]; + } + } + }//j + } + +} + +// Joint Demosaicing and Denoising using High Order Interpolation Techniques +// Revision 0.9.1a - 09/02/2010 - Contact info: luis.sanz.rodriguez@gmail.com +// Copyright Luis Sanz Rodriguez 2010 +// Adapted to RawTherapee by Jacques Desmis 3/2013 + +void RawImageSource::jdl_interpolate_omp() // from "Lassus" +{ + int width = W, height = H; + int row, col, c, d, i, u = width, v = 2 * u, w = 3 * u, x = 4 * u, y = 5 * u, z = 6 * u, indx, (*dif)[2], (*chr)[2]; + float f[4], g[4]; + float (*image)[4]; + image = (float (*)[4]) calloc (width * height, sizeof * image); + dif = (int (*)[2]) calloc(width * height, sizeof * dif); + chr = (int (*)[2]) calloc(width * height, sizeof * chr); + + if (plistener) { + // this function seems to be unused + //plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::jdl))); + plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), "xxx")); + plistener->setProgress (0.0); + } + +#ifdef _OPENMP + #pragma omp parallel default(none) shared(image,width,height,u,w,v,y,x,z,dif,chr) private(row,col,f,g,indx,c,d,i) +#endif + { +#ifdef _OPENMP + #pragma omp for +#endif + + for (int ii = 0; ii < height; ii++) + for (int jj = 0; jj < width; jj++) { + image[ii * width + jj][fc(ii, jj)] = rawData[ii][jj]; + } + + border_interpolate(6, image); + +#ifdef _OPENMP + #pragma omp for +#endif + + for (row = 5; row < height - 5; row++) + for (col = 5 + (FC(row, 1) & 1), indx = row * width + col, c = FC(row, col); col < u - 5; col += 2, indx += 2) { + f[0] = 1.f + abs(image[indx - u][1] - image[indx - w][1]) + abs(image[indx - u][1] - image[indx + u][1]) + abs(image[indx][c] - image[indx - v][c]) + abs(image[indx - v][c] - image[indx - x][c]); + f[1] = 1.f + abs(image[indx + 1][1] - image[indx + 3][1]) + abs(image[indx + 1][1] - image[indx - 1][1]) + abs(image[indx][c] - image[indx + 2][c]) + abs(image[indx + 2][c] - image[indx + 4][c]); + f[2] = 1.f + abs(image[indx - 1][1] - image[indx - 3][1]) + abs(image[indx - 1][1] - image[indx + 1][1]) + abs(image[indx][c] - image[indx - 2][c]) + abs(image[indx - 2][c] - image[indx - 4][c]); + f[3] = 1.f + abs(image[indx + u][1] - image[indx + w][1]) + abs(image[indx + u][1] - image[indx - u][1]) + abs(image[indx][c] - image[indx + v][c]) + abs(image[indx + v][c] - image[indx + x][c]); + g[0] = CLIP((22.f * image[indx - u][1] + 22.f * image[indx - w][1] + 2.f * image[indx - y][1] + 2.f * image[indx + u][1] + 40.f * image[indx][c] - 32.f * image[indx - v][c] - 8.f * image[indx - x][c]) / 48.f); + g[1] = CLIP((22.f * image[indx + 1][1] + 22.f * image[indx + 3][1] + 2.f * image[indx + 5][1] + 2.f * image[indx - 1][1] + 40.f * image[indx][c] - 32.f * image[indx + 2][c] - 8.f * image[indx + 4][c]) / 48.f); + g[2] = CLIP((22.f * image[indx - 1][1] + 22.f * image[indx - 3][1] + 2.f * image[indx - 5][1] + 2.f * image[indx + 1][1] + 40.f * image[indx][c] - 32.f * image[indx - 2][c] - 8.f * image[indx - 4][c]) / 48.f); + g[3] = CLIP((22.f * image[indx + u][1] + 22.f * image[indx + w][1] + 2.f * image[indx + y][1] + 2.f * image[indx - u][1] + 40.f * image[indx][c] - 32.f * image[indx + v][c] - 8.f * image[indx + x][c]) / 48.f); + dif[indx][0] = CLIP((f[3] * g[0] + f[0] * g[3]) / (f[0] + f[3])) - image[indx][c]; + dif[indx][1] = CLIP((f[2] * g[1] + f[1] * g[2]) / (f[1] + f[2])) - image[indx][c]; + } + +#ifdef _OPENMP + #pragma omp for +#endif + + for (row = 6; row < height - 6; row++) + for (col = 6 + (FC(row, 2) & 1), indx = row * width + col, c = FC(row, col) / 2; col < u - 6; col += 2, indx += 2) { + f[0] = 1.f + 78.f * SQR((float)dif[indx][0]) + 69.f * (SQR((float) dif[indx - v][0]) + SQR((float)dif[indx + v][0])) + 51.f * (SQR((float)dif[indx - x][0]) + SQR((float)dif[indx + x][0])) + 21.f * (SQR((float)dif[indx - z][0]) + SQR((float)dif[indx + z][0])) - 6.f * SQR((float)dif[indx - v][0] + dif[indx][0] + dif[indx + v][0]) - 10.f * (SQR((float)dif[indx - x][0] + dif[indx - v][0] + dif[indx][0]) + SQR((float)dif[indx][0] + dif[indx + v][0] + dif[indx + x][0])) - 7.f * (SQR((float)dif[indx - z][0] + dif[indx - x][0] + dif[indx - v][0]) + SQR((float)dif[indx + v][0] + dif[indx + x][0] + dif[indx + z][0])); + f[1] = 1.f + 78.f * SQR((float)dif[indx][1]) + 69.f * (SQR((float)dif[indx - 2][1]) + SQR((float)dif[indx + 2][1])) + 51.f * (SQR((float)dif[indx - 4][1]) + SQR((float)dif[indx + 4][1])) + 21.f * (SQR((float)dif[indx - 6][1]) + SQR((float)dif[indx + 6][1])) - 6.f * SQR((float)dif[indx - 2][1] + dif[indx][1] + dif[indx + 2][1]) - 10.f * (SQR((float)dif[indx - 4][1] + dif[indx - 2][1] + dif[indx][1]) + SQR((float)dif[indx][1] + dif[indx + 2][1] + dif[indx + 4][1])) - 7.f * (SQR((float)dif[indx - 6][1] + dif[indx - 4][1] + dif[indx - 2][1]) + SQR((float)dif[indx + 2][1] + dif[indx + 4][1] + dif[indx + 6][1])); + g[0] = median(0.725f * dif[indx][0] + 0.1375f * dif[indx - v][0] + 0.1375f * dif[indx + v][0], static_cast(dif[indx - v][0]), static_cast(dif[indx + v][0])); + g[1] = median(0.725f * dif[indx][1] + 0.1375f * dif[indx - 2][1] + 0.1375f * dif[indx + 2][1], static_cast(dif[indx - 2][1]), static_cast(dif[indx + 2][1])); + chr[indx][c] = (f[1] * g[0] + f[0] * g[1]) / (f[0] + f[1]); + } + +#ifdef _OPENMP + #pragma omp for +#endif + + for (row = 6; row < height - 6; row++) + for (col = 6 + (FC(row, 2) & 1), indx = row * width + col, c = 1 - FC(row, col) / 2, d = 2 * c; col < u - 6; col += 2, indx += 2) { + f[0] = 1.f / (float)(1.f + fabs((float)chr[indx - u - 1][c] - chr[indx + u + 1][c]) + fabs((float)chr[indx - u - 1][c] - chr[indx - w - 3][c]) + fabs((float)chr[indx + u + 1][c] - chr[indx - w - 3][c])); + f[1] = 1.f / (float)(1.f + fabs((float)chr[indx - u + 1][c] - chr[indx + u - 1][c]) + fabs((float)chr[indx - u + 1][c] - chr[indx - w + 3][c]) + fabs((float)chr[indx + u - 1][c] - chr[indx - w + 3][c])); + f[2] = 1.f / (float)(1.f + fabs((float)chr[indx + u - 1][c] - chr[indx - u + 1][c]) + fabs((float)chr[indx + u - 1][c] - chr[indx + w + 3][c]) + fabs((float)chr[indx - u + 1][c] - chr[indx + w - 3][c])); + f[3] = 1.f / (float)(1.f + fabs((float)chr[indx + u + 1][c] - chr[indx - u - 1][c]) + fabs((float)chr[indx + u + 1][c] - chr[indx + w - 3][c]) + fabs((float)chr[indx - u - 1][c] - chr[indx + w + 3][c])); + g[0] = median(chr[indx - u - 1][c], chr[indx - w - 1][c], chr[indx - u - 3][c]); + g[1] = median(chr[indx - u + 1][c], chr[indx - w + 1][c], chr[indx - u + 3][c]); + g[2] = median(chr[indx + u - 1][c], chr[indx + w - 1][c], chr[indx + u - 3][c]); + g[3] = median(chr[indx + u + 1][c], chr[indx + w + 1][c], chr[indx + u + 3][c]); + chr[indx][c] = (f[0] * g[0] + f[1] * g[1] + f[2] * g[2] + f[3] * g[3]) / (f[0] + f[1] + f[2] + f[3]); + image[indx][1] = CLIP(image[indx][2 - d] + chr[indx][1 - c]); + image[indx][d] = CLIP(image[indx][1] - chr[indx][c]); + } + +#ifdef _OPENMP + #pragma omp for +#endif + + for (row = 6; row < height - 6; row++) + for (col = 6 + (FC(row, 1) & 1), indx = row * width + col, c = FC(row, col + 1) / 2, d = 2 * c; col < u - 6; col += 2, indx += 2) + for(i = 0; i <= 1; c = 1 - c, d = 2 * c, i++) { + f[0] = 1.f / (float)(1.f + fabs((float)chr[indx - u][c] - chr[indx + u][c]) + fabs((float)chr[indx - u][c] - chr[indx - w][c]) + fabs((float)chr[indx + u][c] - chr[indx - w][c])); + f[1] = 1.f / (float)(1.f + fabs((float)chr[indx + 1][c] - chr[indx - 1][c]) + fabs((float)chr[indx + 1][c] - chr[indx + 3][c]) + fabs((float)chr[indx - 1][c] - chr[indx + 3][c])); + f[2] = 1.f / (float)(1.f + fabs((float)chr[indx - 1][c] - chr[indx + 1][c]) + fabs((float)chr[indx - 1][c] - chr[indx - 3][c]) + fabs((float)chr[indx + 1][c] - chr[indx - 3][c])); + f[3] = 1.f / (float)(1.f + fabs((float)chr[indx + u][c] - chr[indx - u][c]) + fabs((float)chr[indx + u][c] - chr[indx + w][c]) + fabs((float)chr[indx - u][c] - chr[indx + w][c])); + g[0] = 0.875f * chr[indx - u][c] + 0.125f * chr[indx - w][c]; + g[1] = 0.875f * chr[indx + 1][c] + 0.125f * chr[indx + 3][c]; + g[2] = 0.875f * chr[indx - 1][c] + 0.125f * chr[indx - 3][c]; + g[3] = 0.875f * chr[indx + u][c] + 0.125f * chr[indx + w][c]; + image[indx][d] = CLIP(image[indx][1] - (f[0] * g[0] + f[1] * g[1] + f[2] * g[2] + f[3] * g[3]) / (f[0] + f[1] + f[2] + f[3])); + } + +#ifdef _OPENMP + #pragma omp for +#endif + + for (int ii = 0; ii < height; ii++) { + for (int jj = 0; jj < width; jj++) { + red[ii][jj] = CLIP(image[ii * width + jj][0]); + green[ii][jj] = CLIP(image[ii * width + jj][1]); + blue[ii][jj] = CLIP(image[ii * width + jj][2]); + } + } + } // End of parallelization + free (image); + free(dif); + free(chr); + //RawImageSource::refinement_lassus(); +} + +// LSMME demosaicing algorithm +// L. Zhang and X. Wu, +// Color demozaicing via directional Linear Minimum Mean Square-error Estimation, +// IEEE Trans. on Image Processing, vol. 14, pp. 2167-2178, +// Dec. 2005. +// Adapted to RawTherapee by Jacques Desmis 3/2013 +// Improved speed and reduced memory consumption by Ingo Weyrich 2/2015 +//TODO Tiles to reduce memory consumption +void RawImageSource::lmmse_interpolate_omp(int winw, int winh, array2D &rawData, array2D &red, array2D &green, array2D &blue, int iterations) +{ + const int width = winw, height = winh; + const int ba = 10; + const int rr1 = height + 2 * ba; + const int cc1 = width + 2 * ba; + const int w1 = cc1; + const int w2 = 2 * w1; + const int w3 = 3 * w1; + const int w4 = 4 * w1; + float h0, h1, h2, h3, h4, hs; + h0 = 1.0f; + h1 = exp( -1.0f / 8.0f); + h2 = exp( -4.0f / 8.0f); + h3 = exp( -9.0f / 8.0f); + h4 = exp(-16.0f / 8.0f); + hs = h0 + 2.0f * (h1 + h2 + h3 + h4); + h0 /= hs; + h1 /= hs; + h2 /= hs; + h3 /= hs; + h4 /= hs; + int passref = 0; + int iter = 0; + + if(iterations <= 4) { + iter = iterations - 1; + passref = 0; + } else if (iterations <= 6) { + iter = 3; + passref = iterations - 4; + } else if (iterations <= 8) { + iter = 3; + passref = iterations - 6; + } + + bool applyGamma = true; + + if(iterations == 0) { + applyGamma = false; + iter = 0; + } else { + applyGamma = true; + } + + float *rix[5]; + float *qix[5]; + float *buffer = (float *)calloc(rr1 * cc1 * 5 * sizeof(float), 1); + + if(buffer == nullptr) { // allocation of big block of memory failed, try to get 5 smaller ones + printf("lmmse_interpolate_omp: allocation of big memory block failed, try to get 5 smaller ones now...\n"); + bool allocationFailed = false; + + for(int i = 0; i < 5; i++) { + qix[i] = (float *)calloc(rr1 * cc1 * sizeof(float), 1); + + if(!qix[i]) { // allocation of at least one small block failed + allocationFailed = true; + } + } + + if(allocationFailed) { // fall back to igv_interpolate + printf("lmmse_interpolate_omp: allocation of 5 small memory blocks failed, falling back to igv_interpolate...\n"); + + for(int i = 0; i < 5; i++) { // free the already allocated buffers + if(qix[i]) { + free(qix[i]); + } + } + + igv_interpolate(winw, winh); + return; + } + } else { + qix[0] = buffer; + + for(int i = 1; i < 5; i++) { + qix[i] = qix[i - 1] + rr1 * cc1; + } + } + + if (plistener) { + plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::LMMSE))); + plistener->setProgress (0.0); + } + + + LUTf *gamtab; + + if(applyGamma) { + gamtab = &(Color::gammatab_24_17a); + } else { + gamtab = new LUTf(65536, LUT_CLIP_ABOVE | LUT_CLIP_BELOW); + gamtab->makeIdentity(65535.f); + } + + +#ifdef _OPENMP + #pragma omp parallel private(rix) +#endif + { +#ifdef _OPENMP + #pragma omp for +#endif + + for (int rrr = ba; rrr < rr1 - ba; rrr++) { + for (int ccc = ba, row = rrr - ba; ccc < cc1 - ba; ccc++) { + int col = ccc - ba; + float *rix = qix[4] + rrr * cc1 + ccc; + rix[0] = (*gamtab)[rawData[row][col]]; + } + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.1); + } + } + + // G-R(B) +#ifdef _OPENMP + #pragma omp for schedule(dynamic,16) +#endif + + for (int rr = 2; rr < rr1 - 2; rr++) { + // G-R(B) at R(B) location + for (int cc = 2 + (FC(rr, 2) & 1); cc < cc1 - 2; cc += 2) { + rix[4] = qix[4] + rr * cc1 + cc; + float v0 = x00625(rix[4][-w1 - 1] + rix[4][-w1 + 1] + rix[4][w1 - 1] + rix[4][w1 + 1]) + x0250(rix[4][0]); + // horizontal + rix[0] = qix[0] + rr * cc1 + cc; + rix[0][0] = - x0250(rix[4][ -2] + rix[4][ 2]) + xdiv2f(rix[4][ -1] + rix[4][0] + rix[4][ 1]); + float Y = v0 + xdiv2f(rix[0][0]); + + if (rix[4][0] > 1.75f * Y) { + rix[0][0] = median(rix[0][0], rix[4][ -1], rix[4][ 1]); + } else { + rix[0][0] = LIM(rix[0][0], 0.0f, 1.0f); + } + + rix[0][0] -= rix[4][0]; + // vertical + rix[1] = qix[1] + rr * cc1 + cc; + rix[1][0] = -x0250(rix[4][-w2] + rix[4][w2]) + xdiv2f(rix[4][-w1] + rix[4][0] + rix[4][w1]); + Y = v0 + xdiv2f(rix[1][0]); + + if (rix[4][0] > 1.75f * Y) { + rix[1][0] = median(rix[1][0], rix[4][-w1], rix[4][w1]); + } else { + rix[1][0] = LIM(rix[1][0], 0.0f, 1.0f); + } + + rix[1][0] -= rix[4][0]; + } + + // G-R(B) at G location + for (int ccc = 2 + (FC(rr, 3) & 1); ccc < cc1 - 2; ccc += 2) { + rix[0] = qix[0] + rr * cc1 + ccc; + rix[1] = qix[1] + rr * cc1 + ccc; + rix[4] = qix[4] + rr * cc1 + ccc; + rix[0][0] = x0250(rix[4][ -2] + rix[4][ 2]) - xdiv2f(rix[4][ -1] + rix[4][0] + rix[4][ 1]); + rix[1][0] = x0250(rix[4][-w2] + rix[4][w2]) - xdiv2f(rix[4][-w1] + rix[4][0] + rix[4][w1]); + rix[0][0] = LIM(rix[0][0], -1.0f, 0.0f) + rix[4][0]; + rix[1][0] = LIM(rix[1][0], -1.0f, 0.0f) + rix[4][0]; + } + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.2); + } + } + + + // apply low pass filter on differential colors +#ifdef _OPENMP + #pragma omp for +#endif + + for (int rr = 4; rr < rr1 - 4; rr++) + for (int cc = 4; cc < cc1 - 4; cc++) { + rix[0] = qix[0] + rr * cc1 + cc; + rix[2] = qix[2] + rr * cc1 + cc; + rix[2][0] = h0 * rix[0][0] + h1 * (rix[0][ -1] + rix[0][ 1]) + h2 * (rix[0][ -2] + rix[0][ 2]) + h3 * (rix[0][ -3] + rix[0][ 3]) + h4 * (rix[0][ -4] + rix[0][ 4]); + rix[1] = qix[1] + rr * cc1 + cc; + rix[3] = qix[3] + rr * cc1 + cc; + rix[3][0] = h0 * rix[1][0] + h1 * (rix[1][-w1] + rix[1][w1]) + h2 * (rix[1][-w2] + rix[1][w2]) + h3 * (rix[1][-w3] + rix[1][w3]) + h4 * (rix[1][-w4] + rix[1][w4]); + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.3); + } + } + + // interpolate G-R(B) at R(B) +#ifdef _OPENMP + #pragma omp for +#endif + + for (int rr = 4; rr < rr1 - 4; rr++) { + int cc = 4 + (FC(rr, 4) & 1); +#ifdef __SSE2__ + __m128 p1v, p2v, p3v, p4v, p5v, p6v, p7v, p8v, p9v, muv, vxv, vnv, xhv, vhv, xvv, vvv; + __m128 epsv = _mm_set1_ps(1e-7); + __m128 ninev = _mm_set1_ps(9.f); + + for (; cc < cc1 - 10; cc += 8) { + rix[0] = qix[0] + rr * cc1 + cc; + rix[1] = qix[1] + rr * cc1 + cc; + rix[2] = qix[2] + rr * cc1 + cc; + rix[3] = qix[3] + rr * cc1 + cc; + rix[4] = qix[4] + rr * cc1 + cc; + // horizontal + p1v = LC2VFU(rix[2][-4]); + p2v = LC2VFU(rix[2][-3]); + p3v = LC2VFU(rix[2][-2]); + p4v = LC2VFU(rix[2][-1]); + p5v = LC2VFU(rix[2][ 0]); + p6v = LC2VFU(rix[2][ 1]); + p7v = LC2VFU(rix[2][ 2]); + p8v = LC2VFU(rix[2][ 3]); + p9v = LC2VFU(rix[2][ 4]); + muv = (p1v + p2v + p3v + p4v + p5v + p6v + p7v + p8v + p9v) / ninev; + vxv = epsv + SQRV(p1v - muv) + SQRV(p2v - muv) + SQRV(p3v - muv) + SQRV(p4v - muv) + SQRV(p5v - muv) + SQRV(p6v - muv) + SQRV(p7v - muv) + SQRV(p8v - muv) + SQRV(p9v - muv); + p1v -= LC2VFU(rix[0][-4]); + p2v -= LC2VFU(rix[0][-3]); + p3v -= LC2VFU(rix[0][-2]); + p4v -= LC2VFU(rix[0][-1]); + p5v -= LC2VFU(rix[0][ 0]); + p6v -= LC2VFU(rix[0][ 1]); + p7v -= LC2VFU(rix[0][ 2]); + p8v -= LC2VFU(rix[0][ 3]); + p9v -= LC2VFU(rix[0][ 4]); + vnv = epsv + SQRV(p1v) + SQRV(p2v) + SQRV(p3v) + SQRV(p4v) + SQRV(p5v) + SQRV(p6v) + SQRV(p7v) + SQRV(p8v) + SQRV(p9v); + xhv = (LC2VFU(rix[0][0]) * vxv + LC2VFU(rix[2][0]) * vnv) / (vxv + vnv); + vhv = vxv * vnv / (vxv + vnv); + + // vertical + p1v = LC2VFU(rix[3][-w4]); + p2v = LC2VFU(rix[3][-w3]); + p3v = LC2VFU(rix[3][-w2]); + p4v = LC2VFU(rix[3][-w1]); + p5v = LC2VFU(rix[3][ 0]); + p6v = LC2VFU(rix[3][ w1]); + p7v = LC2VFU(rix[3][ w2]); + p8v = LC2VFU(rix[3][ w3]); + p9v = LC2VFU(rix[3][ w4]); + muv = (p1v + p2v + p3v + p4v + p5v + p6v + p7v + p8v + p9v) / ninev; + vxv = epsv + SQRV(p1v - muv) + SQRV(p2v - muv) + SQRV(p3v - muv) + SQRV(p4v - muv) + SQRV(p5v - muv) + SQRV(p6v - muv) + SQRV(p7v - muv) + SQRV(p8v - muv) + SQRV(p9v - muv); + p1v -= LC2VFU(rix[1][-w4]); + p2v -= LC2VFU(rix[1][-w3]); + p3v -= LC2VFU(rix[1][-w2]); + p4v -= LC2VFU(rix[1][-w1]); + p5v -= LC2VFU(rix[1][ 0]); + p6v -= LC2VFU(rix[1][ w1]); + p7v -= LC2VFU(rix[1][ w2]); + p8v -= LC2VFU(rix[1][ w3]); + p9v -= LC2VFU(rix[1][ w4]); + vnv = epsv + SQRV(p1v) + SQRV(p2v) + SQRV(p3v) + SQRV(p4v) + SQRV(p5v) + SQRV(p6v) + SQRV(p7v) + SQRV(p8v) + SQRV(p9v); + xvv = (LC2VFU(rix[1][0]) * vxv + LC2VFU(rix[3][0]) * vnv) / (vxv + vnv); + vvv = vxv * vnv / (vxv + vnv); + // interpolated G-R(B) + muv = (xhv * vvv + xvv * vhv) / (vhv + vvv); + STC2VFU(rix[4][0], muv); + } + +#endif + + for (; cc < cc1 - 4; cc += 2) { + rix[0] = qix[0] + rr * cc1 + cc; + rix[1] = qix[1] + rr * cc1 + cc; + rix[2] = qix[2] + rr * cc1 + cc; + rix[3] = qix[3] + rr * cc1 + cc; + rix[4] = qix[4] + rr * cc1 + cc; + // horizontal + float p1 = rix[2][-4]; + float p2 = rix[2][-3]; + float p3 = rix[2][-2]; + float p4 = rix[2][-1]; + float p5 = rix[2][ 0]; + float p6 = rix[2][ 1]; + float p7 = rix[2][ 2]; + float p8 = rix[2][ 3]; + float p9 = rix[2][ 4]; + float mu = (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9) / 9.f; + float vx = 1e-7 + SQR(p1 - mu) + SQR(p2 - mu) + SQR(p3 - mu) + SQR(p4 - mu) + SQR(p5 - mu) + SQR(p6 - mu) + SQR(p7 - mu) + SQR(p8 - mu) + SQR(p9 - mu); + p1 -= rix[0][-4]; + p2 -= rix[0][-3]; + p3 -= rix[0][-2]; + p4 -= rix[0][-1]; + p5 -= rix[0][ 0]; + p6 -= rix[0][ 1]; + p7 -= rix[0][ 2]; + p8 -= rix[0][ 3]; + p9 -= rix[0][ 4]; + float vn = 1e-7 + SQR(p1) + SQR(p2) + SQR(p3) + SQR(p4) + SQR(p5) + SQR(p6) + SQR(p7) + SQR(p8) + SQR(p9); + float xh = (rix[0][0] * vx + rix[2][0] * vn) / (vx + vn); + float vh = vx * vn / (vx + vn); + + // vertical + p1 = rix[3][-w4]; + p2 = rix[3][-w3]; + p3 = rix[3][-w2]; + p4 = rix[3][-w1]; + p5 = rix[3][ 0]; + p6 = rix[3][ w1]; + p7 = rix[3][ w2]; + p8 = rix[3][ w3]; + p9 = rix[3][ w4]; + mu = (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9) / 9.f; + vx = 1e-7 + SQR(p1 - mu) + SQR(p2 - mu) + SQR(p3 - mu) + SQR(p4 - mu) + SQR(p5 - mu) + SQR(p6 - mu) + SQR(p7 - mu) + SQR(p8 - mu) + SQR(p9 - mu); + p1 -= rix[1][-w4]; + p2 -= rix[1][-w3]; + p3 -= rix[1][-w2]; + p4 -= rix[1][-w1]; + p5 -= rix[1][ 0]; + p6 -= rix[1][ w1]; + p7 -= rix[1][ w2]; + p8 -= rix[1][ w3]; + p9 -= rix[1][ w4]; + vn = 1e-7 + SQR(p1) + SQR(p2) + SQR(p3) + SQR(p4) + SQR(p5) + SQR(p6) + SQR(p7) + SQR(p8) + SQR(p9); + float xv = (rix[1][0] * vx + rix[3][0] * vn) / (vx + vn); + float vv = vx * vn / (vx + vn); + // interpolated G-R(B) + rix[4][0] = (xh * vv + xv * vh) / (vh + vv); + } + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.4); + } + } + + // copy CFA values +#ifdef _OPENMP + #pragma omp for +#endif + + for (int rr = 0; rr < rr1; rr++) + for (int cc = 0, row = rr - ba; cc < cc1; cc++) { + int col = cc - ba; + int c = FC(rr, cc); + rix[c] = qix[c] + rr * cc1 + cc; + + if ((row >= 0) & (row < height) & (col >= 0) & (col < width)) { + rix[c][0] = (*gamtab)[rawData[row][col]]; + } else { + rix[c][0] = 0.f; + } + + if (c != 1) { + rix[1] = qix[1] + rr * cc1 + cc; + rix[4] = qix[4] + rr * cc1 + cc; + rix[1][0] = rix[c][0] + rix[4][0]; + } + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.5); + } + } + + // bilinear interpolation for R/B + // interpolate R/B at G location +#ifdef _OPENMP + #pragma omp for +#endif + + for (int rr = 1; rr < rr1 - 1; rr++) + for (int cc = 1 + (FC(rr, 2) & 1), c = FC(rr, cc + 1); cc < cc1 - 1; cc += 2) { + rix[c] = qix[c] + rr * cc1 + cc; + rix[1] = qix[1] + rr * cc1 + cc; + rix[c][0] = rix[1][0] + xdiv2f(rix[c][ -1] - rix[1][ -1] + rix[c][ 1] - rix[1][ 1]); + c = 2 - c; + rix[c] = qix[c] + rr * cc1 + cc; + rix[c][0] = rix[1][0] + xdiv2f(rix[c][-w1] - rix[1][-w1] + rix[c][w1] - rix[1][w1]); + c = 2 - c; + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.6); + } + } + + // interpolate R/B at B/R location +#ifdef _OPENMP + #pragma omp for +#endif + + for (int rr = 1; rr < rr1 - 1; rr++) + for (int cc = 1 + (FC(rr, 1) & 1), c = 2 - FC(rr, cc); cc < cc1 - 1; cc += 2) { + rix[c] = qix[c] + rr * cc1 + cc; + rix[1] = qix[1] + rr * cc1 + cc; + rix[c][0] = rix[1][0] + x0250(rix[c][-w1] - rix[1][-w1] + rix[c][ -1] - rix[1][ -1] + rix[c][ 1] - rix[1][ 1] + rix[c][ w1] - rix[1][ w1]); + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.7); + } + } + + }// End of parallelization 1 + + // median filter/ + for (int pass = 0; pass < iter; pass++) { + // Apply 3x3 median filter + // Compute median(R-G) and median(B-G) + +#ifdef _OPENMP + #pragma omp parallel for private(rix) +#endif + + for (int rr = 1; rr < rr1 - 1; rr++) { + for (int c = 0; c < 3; c += 2) { + int d = c + 3 - (c == 0 ? 0 : 1); + int cc = 1; +#ifdef __SSE2__ + + for (; cc < cc1 - 4; cc += 4) { + rix[d] = qix[d] + rr * cc1 + cc; + rix[c] = qix[c] + rr * cc1 + cc; + rix[1] = qix[1] + rr * cc1 + cc; + // Assign 3x3 differential color values + const std::array p = { + LVFU(rix[c][-w1 - 1]) - LVFU(rix[1][-w1 - 1]), + LVFU(rix[c][-w1]) - LVFU(rix[1][-w1]), + LVFU(rix[c][-w1 + 1]) - LVFU(rix[1][-w1 + 1]), + LVFU(rix[c][ -1]) - LVFU(rix[1][ -1]), + LVFU(rix[c][ 0]) - LVFU(rix[1][ 0]), + LVFU(rix[c][ 1]) - LVFU(rix[1][ 1]), + LVFU(rix[c][ w1 - 1]) - LVFU(rix[1][ w1 - 1]), + LVFU(rix[c][ w1]) - LVFU(rix[1][ w1]), + LVFU(rix[c][ w1 + 1]) - LVFU(rix[1][ w1 + 1]) + }; + _mm_storeu_ps(&rix[d][0], median(p)); + } + +#endif + + for (; cc < cc1 - 1; cc++) { + rix[d] = qix[d] + rr * cc1 + cc; + rix[c] = qix[c] + rr * cc1 + cc; + rix[1] = qix[1] + rr * cc1 + cc; + // Assign 3x3 differential color values + const std::array p = { + rix[c][-w1 - 1] - rix[1][-w1 - 1], + rix[c][-w1] - rix[1][-w1], + rix[c][-w1 + 1] - rix[1][-w1 + 1], + rix[c][ -1] - rix[1][ -1], + rix[c][ 0] - rix[1][ 0], + rix[c][ 1] - rix[1][ 1], + rix[c][ w1 - 1] - rix[1][ w1 - 1], + rix[c][ w1] - rix[1][ w1], + rix[c][ w1 + 1] - rix[1][ w1 + 1] + }; + rix[d][0] = median(p); + } + } + } + + // red/blue at GREEN pixel locations & red/blue and green at BLUE/RED pixel locations +#ifdef _OPENMP + #pragma omp parallel for private (rix) +#endif + + for (int rr = 0; rr < rr1; rr++) { + rix[0] = qix[0] + rr * cc1; + rix[1] = qix[1] + rr * cc1; + rix[2] = qix[2] + rr * cc1; + rix[3] = qix[3] + rr * cc1; + rix[4] = qix[4] + rr * cc1; + int c0 = FC(rr, 0); + int c1 = FC(rr, 1); + + if(c0 == 1) { + c1 = 2 - c1; + int d = c1 + 3 - (c1 == 0 ? 0 : 1); + int cc; + + for (cc = 0; cc < cc1 - 1; cc += 2) { + rix[0][0] = rix[1][0] + rix[3][0]; + rix[2][0] = rix[1][0] + rix[4][0]; + rix[0]++; + rix[1]++; + rix[2]++; + rix[3]++; + rix[4]++; + rix[c1][0] = rix[1][0] + rix[d][0]; + rix[1][0] = 0.5f * (rix[0][0] - rix[3][0] + rix[2][0] - rix[4][0]); + rix[0]++; + rix[1]++; + rix[2]++; + rix[3]++; + rix[4]++; + } + + if(cc < cc1) { // remaining pixel, only if width is odd + rix[0][0] = rix[1][0] + rix[3][0]; + rix[2][0] = rix[1][0] + rix[4][0]; + } + } else { + c0 = 2 - c0; + int d = c0 + 3 - (c0 == 0 ? 0 : 1); + int cc; + + for (cc = 0; cc < cc1 - 1; cc += 2) { + rix[c0][0] = rix[1][0] + rix[d][0]; + rix[1][0] = 0.5f * (rix[0][0] - rix[3][0] + rix[2][0] - rix[4][0]); + rix[0]++; + rix[1]++; + rix[2]++; + rix[3]++; + rix[4]++; + rix[0][0] = rix[1][0] + rix[3][0]; + rix[2][0] = rix[1][0] + rix[4][0]; + rix[0]++; + rix[1]++; + rix[2]++; + rix[3]++; + rix[4]++; + } + + if(cc < cc1) { // remaining pixel, only if width is odd + rix[c0][0] = rix[1][0] + rix[d][0]; + rix[1][0] = 0.5f * (rix[0][0] - rix[3][0] + rix[2][0] - rix[4][0]); + } + } + } + } + + if (plistener) { + plistener->setProgress (0.8); + } + + if(applyGamma) { + gamtab = &(Color::igammatab_24_17); + } else { + gamtab->makeIdentity(); + } + + array2D* rgb[3]; + rgb[0] = &red; + rgb[1] = &green; + rgb[2] = &blue; + + // copy result back to image matrix +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int row = 0; row < height; row++) { + for (int col = 0, rr = row + ba; col < width; col++) { + int cc = col + ba; + int c = FC(row, col); + + for (int ii = 0; ii < 3; ii++) + if (ii != c) { + float *rix = qix[ii] + rr * cc1 + cc; + (*(rgb[ii]))[row][col] = (*gamtab)[65535.f * rix[0]]; + } else { + (*(rgb[ii]))[row][col] = CLIP(rawData[row][col]); + } + } + } + + if (plistener) { + plistener->setProgress (1.0); + } + + if(buffer) { + free(buffer); + } else + for(int i = 0; i < 5; i++) { + free(qix[i]); + } + + if(!applyGamma) { + delete gamtab; + } + + if(iterations > 4 && iterations <= 6) { + refinement(passref); + } else if(iterations > 6) { + refinement_lassus(passref); + } + +} + +/*** +* +* Bayer CFA Demosaicing using Integrated Gaussian Vector on Color Differences +* Revision 1.0 - 2013/02/28 +* +* Copyright (c) 2007-2013 Luis Sanz Rodriguez +* Using High Order Interpolation technique by Jim S, Jimmy Li, and Sharmil Randhawa +* +* Contact info: luis.sanz.rodriguez@gmail.com +* +* This code is distributed under a GNU General Public License, version 3. +* Visit for more information. +* +***/ +// Adapted to RawTherapee by Jacques Desmis 3/2013 +// SSE version by Ingo Weyrich 5/2013 +#ifdef __SSE2__ +#define CLIPV(a) LIMV(a,zerov,c65535v) +void RawImageSource::igv_interpolate(int winw, int winh) +{ + static const float eps = 1e-5f, epssq = 1e-5f; //mod epssq -10f =>-5f Jacques 3/2013 to prevent artifact (divide by zero) + + static const int h1 = 1, h2 = 2, h3 = 3, h5 = 5; + const int width = winw, height = winh; + const int v1 = 1 * width, v2 = 2 * width, v3 = 3 * width, v5 = 5 * width; + float* rgb[2]; + float* chr[4]; + float *rgbarray, *vdif, *hdif, *chrarray; + rgbarray = (float (*)) malloc((width * height) * sizeof( float ) ); + rgb[0] = rgbarray; + rgb[1] = rgbarray + (width * height) / 2; + + vdif = (float (*)) calloc( width * height / 2, sizeof * vdif ); + hdif = (float (*)) calloc( width * height / 2, sizeof * hdif ); + + chrarray = (float (*)) calloc( width * height, sizeof( float ) ); + chr[0] = chrarray; + chr[1] = chrarray + (width * height) / 2; + + // mapped chr[2] and chr[3] to hdif and hdif, because these are out of use, when chr[2] and chr[3] are used + chr[2] = hdif; + chr[3] = vdif; + + border_interpolate2(winw, winh, 7); + + if (plistener) { + plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::IGV))); + plistener->setProgress (0.0); + } + +#ifdef _OPENMP + #pragma omp parallel default(none) shared(rgb,vdif,hdif,chr) +#endif + { + __m128 ngv, egv, wgv, sgv, nvv, evv, wvv, svv, nwgv, negv, swgv, segv, nwvv, nevv, swvv, sevv, tempv, temp1v, temp2v, temp3v, temp4v, temp5v, temp6v, temp7v, temp8v; + __m128 epsv = _mm_set1_ps( eps ); + __m128 epssqv = _mm_set1_ps( epssq ); + __m128 c65535v = _mm_set1_ps( 65535.f ); + __m128 c23v = _mm_set1_ps( 23.f ); + __m128 c40v = _mm_set1_ps( 40.f ); + __m128 c51v = _mm_set1_ps( 51.f ); + __m128 c32v = _mm_set1_ps( 32.f ); + __m128 c8v = _mm_set1_ps( 8.f ); + __m128 c7v = _mm_set1_ps( 7.f ); + __m128 c6v = _mm_set1_ps( 6.f ); + __m128 c10v = _mm_set1_ps( 10.f ); + __m128 c21v = _mm_set1_ps( 21.f ); + __m128 c78v = _mm_set1_ps( 78.f ); + __m128 c69v = _mm_set1_ps( 69.f ); + __m128 c3145680v = _mm_set1_ps( 3145680.f ); + __m128 onev = _mm_set1_ps ( 1.f ); + __m128 zerov = _mm_set1_ps ( 0.f ); + __m128 d725v = _mm_set1_ps ( 0.725f ); + __m128 d1375v = _mm_set1_ps ( 0.1375f ); + + float *dest1, *dest2; + float ng, eg, wg, sg, nv, ev, wv, sv, nwg, neg, swg, seg, nwv, nev, swv, sev; +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 0; row < height - 0; row++) { + dest1 = rgb[FC(row, 0) & 1]; + dest2 = rgb[FC(row, 1) & 1]; + int col, indx; + + for (col = 0, indx = row * width + col; col < width - 7; col += 8, indx += 8) { + temp1v = LVFU( rawData[row][col] ); + temp1v = CLIPV( temp1v ); + temp2v = LVFU( rawData[row][col + 4] ); + temp2v = CLIPV( temp2v ); + tempv = _mm_shuffle_ps( temp1v, temp2v, _MM_SHUFFLE( 2, 0, 2, 0 ) ); + _mm_storeu_ps( &dest1[indx >> 1], tempv ); + tempv = _mm_shuffle_ps( temp1v, temp2v, _MM_SHUFFLE( 3, 1, 3, 1 ) ); + _mm_storeu_ps( &dest2[indx >> 1], tempv ); + } + + for (; col < width; col++, indx += 2) { + dest1[indx >> 1] = CLIP(rawData[row][col]); //rawData = RT datas + col++; + if(col < width) + dest2[indx >> 1] = CLIP(rawData[row][col]); //rawData = RT datas + } + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.13); + } + } + +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 5; row < height - 5; row++) { + int col, indx, indx1; + + for (col = 5 + (FC(row, 1) & 1), indx = row * width + col, indx1 = indx >> 1; col < width - 12; col += 8, indx += 8, indx1 += 4) { + //N,E,W,S Gradients + ngv = (epsv + (vabsf(LVFU(rgb[1][(indx - v1) >> 1]) - LVFU(rgb[1][(indx - v3) >> 1])) + vabsf(LVFU(rgb[0][indx1]) - LVFU(rgb[0][(indx1 - v1)]))) / c65535v); + egv = (epsv + (vabsf(LVFU(rgb[1][(indx + h1) >> 1]) - LVFU(rgb[1][(indx + h3) >> 1])) + vabsf(LVFU(rgb[0][indx1]) - LVFU(rgb[0][(indx1 + h1)]))) / c65535v); + wgv = (epsv + (vabsf(LVFU(rgb[1][(indx - h1) >> 1]) - LVFU(rgb[1][(indx - h3) >> 1])) + vabsf(LVFU(rgb[0][indx1]) - LVFU(rgb[0][(indx1 - h1)]))) / c65535v); + sgv = (epsv + (vabsf(LVFU(rgb[1][(indx + v1) >> 1]) - LVFU(rgb[1][(indx + v3) >> 1])) + vabsf(LVFU(rgb[0][indx1]) - LVFU(rgb[0][(indx1 + v1)]))) / c65535v); + //N,E,W,S High Order Interpolation (Li & Randhawa) + //N,E,W,S Hamilton Adams Interpolation + // (48.f * 65535.f) = 3145680.f + tempv = c40v * LVFU(rgb[0][indx1]); + nvv = LIMV(((c23v * LVFU(rgb[1][(indx - v1) >> 1]) + c23v * LVFU(rgb[1][(indx - v3) >> 1]) + LVFU(rgb[1][(indx - v5) >> 1]) + LVFU(rgb[1][(indx + v1) >> 1]) + tempv - c32v * LVFU(rgb[0][(indx1 - v1)]) - c8v * LVFU(rgb[0][(indx1 - v2)]))) / c3145680v, zerov, onev); + evv = LIMV(((c23v * LVFU(rgb[1][(indx + h1) >> 1]) + c23v * LVFU(rgb[1][(indx + h3) >> 1]) + LVFU(rgb[1][(indx + h5) >> 1]) + LVFU(rgb[1][(indx - h1) >> 1]) + tempv - c32v * LVFU(rgb[0][(indx1 + h1)]) - c8v * LVFU(rgb[0][(indx1 + h2)]))) / c3145680v, zerov, onev); + wvv = LIMV(((c23v * LVFU(rgb[1][(indx - h1) >> 1]) + c23v * LVFU(rgb[1][(indx - h3) >> 1]) + LVFU(rgb[1][(indx - h5) >> 1]) + LVFU(rgb[1][(indx + h1) >> 1]) + tempv - c32v * LVFU(rgb[0][(indx1 - h1)]) - c8v * LVFU(rgb[0][(indx1 - h2)]))) / c3145680v, zerov, onev); + svv = LIMV(((c23v * LVFU(rgb[1][(indx + v1) >> 1]) + c23v * LVFU(rgb[1][(indx + v3) >> 1]) + LVFU(rgb[1][(indx + v5) >> 1]) + LVFU(rgb[1][(indx - v1) >> 1]) + tempv - c32v * LVFU(rgb[0][(indx1 + v1)]) - c8v * LVFU(rgb[0][(indx1 + v2)]))) / c3145680v, zerov, onev); + //Horizontal and vertical color differences + tempv = LVFU( rgb[0][indx1] ) / c65535v; + _mm_storeu_ps( &vdif[indx1], (sgv * nvv + ngv * svv) / (ngv + sgv) - tempv ); + _mm_storeu_ps( &hdif[indx1], (wgv * evv + egv * wvv) / (egv + wgv) - tempv ); + } + + // borders without SSE + for (; col < width - 5; col += 2, indx += 2, indx1++) { + //N,E,W,S Gradients + ng = (eps + (fabsf(rgb[1][(indx - v1) >> 1] - rgb[1][(indx - v3) >> 1]) + fabsf(rgb[0][indx1] - rgb[0][(indx1 - v1)])) / 65535.f);; + eg = (eps + (fabsf(rgb[1][(indx + h1) >> 1] - rgb[1][(indx + h3) >> 1]) + fabsf(rgb[0][indx1] - rgb[0][(indx1 + h1)])) / 65535.f); + wg = (eps + (fabsf(rgb[1][(indx - h1) >> 1] - rgb[1][(indx - h3) >> 1]) + fabsf(rgb[0][indx1] - rgb[0][(indx1 - h1)])) / 65535.f); + sg = (eps + (fabsf(rgb[1][(indx + v1) >> 1] - rgb[1][(indx + v3) >> 1]) + fabsf(rgb[0][indx1] - rgb[0][(indx1 + v1)])) / 65535.f); + //N,E,W,S High Order Interpolation (Li & Randhawa) + //N,E,W,S Hamilton Adams Interpolation + // (48.f * 65535.f) = 3145680.f + nv = LIM(((23.0f * rgb[1][(indx - v1) >> 1] + 23.0f * rgb[1][(indx - v3) >> 1] + rgb[1][(indx - v5) >> 1] + rgb[1][(indx + v1) >> 1] + 40.0f * rgb[0][indx1] - 32.0f * rgb[0][(indx1 - v1)] - 8.0f * rgb[0][(indx1 - v2)])) / 3145680.f, 0.0f, 1.0f); + ev = LIM(((23.0f * rgb[1][(indx + h1) >> 1] + 23.0f * rgb[1][(indx + h3) >> 1] + rgb[1][(indx + h5) >> 1] + rgb[1][(indx - h1) >> 1] + 40.0f * rgb[0][indx1] - 32.0f * rgb[0][(indx1 + h1)] - 8.0f * rgb[0][(indx1 + h2)])) / 3145680.f, 0.0f, 1.0f); + wv = LIM(((23.0f * rgb[1][(indx - h1) >> 1] + 23.0f * rgb[1][(indx - h3) >> 1] + rgb[1][(indx - h5) >> 1] + rgb[1][(indx + h1) >> 1] + 40.0f * rgb[0][indx1] - 32.0f * rgb[0][(indx1 - h1)] - 8.0f * rgb[0][(indx1 - h2)])) / 3145680.f, 0.0f, 1.0f); + sv = LIM(((23.0f * rgb[1][(indx + v1) >> 1] + 23.0f * rgb[1][(indx + v3) >> 1] + rgb[1][(indx + v5) >> 1] + rgb[1][(indx - v1) >> 1] + 40.0f * rgb[0][indx1] - 32.0f * rgb[0][(indx1 + v1)] - 8.0f * rgb[0][(indx1 + v2)])) / 3145680.f, 0.0f, 1.0f); + //Horizontal and vertical color differences + vdif[indx1] = (sg * nv + ng * sv) / (ng + sg) - (rgb[0][indx1]) / 65535.f; + hdif[indx1] = (wg * ev + eg * wv) / (eg + wg) - (rgb[0][indx1]) / 65535.f; + } + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.26); + } + } +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 7; row < height - 7; row++) { + int col, d, indx1; + + for (col = 7 + (FC(row, 1) & 1), indx1 = (row * width + col) >> 1, d = FC(row, col) / 2; col < width - 14; col += 8, indx1 += 4) { + //H&V integrated gaussian vector over variance on color differences + //Mod Jacques 3/2013 + ngv = LIMV(epssqv + c78v * SQRV(LVFU(vdif[indx1])) + c69v * (SQRV(LVFU(vdif[indx1 - v1])) + SQRV(LVFU(vdif[indx1 + v1]))) + c51v * (SQRV(LVFU(vdif[indx1 - v2])) + SQRV(LVFU(vdif[indx1 + v2]))) + c21v * (SQRV(LVFU(vdif[indx1 - v3])) + SQRV(LVFU(vdif[indx1 + v3]))) - c6v * SQRV(LVFU(vdif[indx1 - v1]) + LVFU(vdif[indx1]) + LVFU(vdif[indx1 + v1])) + - c10v * (SQRV(LVFU(vdif[indx1 - v2]) + LVFU(vdif[indx1 - v1]) + LVFU(vdif[indx1])) + SQRV(LVFU(vdif[indx1]) + LVFU(vdif[indx1 + v1]) + LVFU(vdif[indx1 + v2]))) - c7v * (SQRV(LVFU(vdif[indx1 - v3]) + LVFU(vdif[indx1 - v2]) + LVFU(vdif[indx1 - v1])) + SQRV(LVFU(vdif[indx1 + v1]) + LVFU(vdif[indx1 + v2]) + LVFU(vdif[indx1 + v3]))), zerov, onev); + egv = LIMV(epssqv + c78v * SQRV(LVFU(hdif[indx1])) + c69v * (SQRV(LVFU(hdif[indx1 - h1])) + SQRV(LVFU(hdif[indx1 + h1]))) + c51v * (SQRV(LVFU(hdif[indx1 - h2])) + SQRV(LVFU(hdif[indx1 + h2]))) + c21v * (SQRV(LVFU(hdif[indx1 - h3])) + SQRV(LVFU(hdif[indx1 + h3]))) - c6v * SQRV(LVFU(hdif[indx1 - h1]) + LVFU(hdif[indx1]) + LVFU(hdif[indx1 + h1])) + - c10v * (SQRV(LVFU(hdif[indx1 - h2]) + LVFU(hdif[indx1 - h1]) + LVFU(hdif[indx1])) + SQRV(LVFU(hdif[indx1]) + LVFU(hdif[indx1 + h1]) + LVFU(hdif[indx1 + h2]))) - c7v * (SQRV(LVFU(hdif[indx1 - h3]) + LVFU(hdif[indx1 - h2]) + LVFU(hdif[indx1 - h1])) + SQRV(LVFU(hdif[indx1 + h1]) + LVFU(hdif[indx1 + h2]) + LVFU(hdif[indx1 + h3]))), zerov, onev); + //Limit chrominance using H/V neighbourhood + nvv = median(d725v * LVFU(vdif[indx1]) + d1375v * LVFU(vdif[indx1 - v1]) + d1375v * LVFU(vdif[indx1 + v1]), LVFU(vdif[indx1 - v1]), LVFU(vdif[indx1 + v1])); + evv = median(d725v * LVFU(hdif[indx1]) + d1375v * LVFU(hdif[indx1 - h1]) + d1375v * LVFU(hdif[indx1 + h1]), LVFU(hdif[indx1 - h1]), LVFU(hdif[indx1 + h1])); + //Chrominance estimation + tempv = (egv * nvv + ngv * evv) / (ngv + egv); + _mm_storeu_ps(&(chr[d][indx1]), tempv); + //Green channel population + temp1v = c65535v * tempv + LVFU(rgb[0][indx1]); + _mm_storeu_ps( &(rgb[0][indx1]), temp1v ); + } + + for (; col < width - 7; col += 2, indx1++) { + //H&V integrated gaussian vector over variance on color differences + //Mod Jacques 3/2013 + ng = LIM(epssq + 78.0f * SQR(vdif[indx1]) + 69.0f * (SQR(vdif[indx1 - v1]) + SQR(vdif[indx1 + v1])) + 51.0f * (SQR(vdif[indx1 - v2]) + SQR(vdif[indx1 + v2])) + 21.0f * (SQR(vdif[indx1 - v3]) + SQR(vdif[indx1 + v3])) - 6.0f * SQR(vdif[indx1 - v1] + vdif[indx1] + vdif[indx1 + v1]) + - 10.0f * (SQR(vdif[indx1 - v2] + vdif[indx1 - v1] + vdif[indx1]) + SQR(vdif[indx1] + vdif[indx1 + v1] + vdif[indx1 + v2])) - 7.0f * (SQR(vdif[indx1 - v3] + vdif[indx1 - v2] + vdif[indx1 - v1]) + SQR(vdif[indx1 + v1] + vdif[indx1 + v2] + vdif[indx1 + v3])), 0.f, 1.f); + eg = LIM(epssq + 78.0f * SQR(hdif[indx1]) + 69.0f * (SQR(hdif[indx1 - h1]) + SQR(hdif[indx1 + h1])) + 51.0f * (SQR(hdif[indx1 - h2]) + SQR(hdif[indx1 + h2])) + 21.0f * (SQR(hdif[indx1 - h3]) + SQR(hdif[indx1 + h3])) - 6.0f * SQR(hdif[indx1 - h1] + hdif[indx1] + hdif[indx1 + h1]) + - 10.0f * (SQR(hdif[indx1 - h2] + hdif[indx1 - h1] + hdif[indx1]) + SQR(hdif[indx1] + hdif[indx1 + h1] + hdif[indx1 + h2])) - 7.0f * (SQR(hdif[indx1 - h3] + hdif[indx1 - h2] + hdif[indx1 - h1]) + SQR(hdif[indx1 + h1] + hdif[indx1 + h2] + hdif[indx1 + h3])), 0.f, 1.f); + //Limit chrominance using H/V neighbourhood + nv = median(0.725f * vdif[indx1] + 0.1375f * vdif[indx1 - v1] + 0.1375f * vdif[indx1 + v1], vdif[indx1 - v1], vdif[indx1 + v1]); + ev = median(0.725f * hdif[indx1] + 0.1375f * hdif[indx1 - h1] + 0.1375f * hdif[indx1 + h1], hdif[indx1 - h1], hdif[indx1 + h1]); + //Chrominance estimation + chr[d][indx1] = (eg * nv + ng * ev) / (ng + eg); + //Green channel population + rgb[0][indx1] = rgb[0][indx1] + 65535.f * chr[d][indx1]; + } + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.39); + } + } +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 7; row < height - 7; row++) { + int col, indx, c; + + for (col = 7 + (FC(row, 1) & 1), indx = row * width + col, c = 1 - FC(row, col) / 2; col < width - 14; col += 8, indx += 8) { + //NW,NE,SW,SE Gradients + nwgv = onev / (epsv + vabsf(LVFU(chr[c][(indx - v1 - h1) >> 1]) - LVFU(chr[c][(indx - v3 - h3) >> 1])) + vabsf(LVFU(chr[c][(indx + v1 + h1) >> 1]) - LVFU(chr[c][(indx - v3 - h3) >> 1]))); + negv = onev / (epsv + vabsf(LVFU(chr[c][(indx - v1 + h1) >> 1]) - LVFU(chr[c][(indx - v3 + h3) >> 1])) + vabsf(LVFU(chr[c][(indx + v1 - h1) >> 1]) - LVFU(chr[c][(indx - v3 + h3) >> 1]))); + swgv = onev / (epsv + vabsf(LVFU(chr[c][(indx + v1 - h1) >> 1]) - LVFU(chr[c][(indx + v3 + h3) >> 1])) + vabsf(LVFU(chr[c][(indx - v1 + h1) >> 1]) - LVFU(chr[c][(indx + v3 - h3) >> 1]))); + segv = onev / (epsv + vabsf(LVFU(chr[c][(indx + v1 + h1) >> 1]) - LVFU(chr[c][(indx + v3 - h3) >> 1])) + vabsf(LVFU(chr[c][(indx - v1 - h1) >> 1]) - LVFU(chr[c][(indx + v3 + h3) >> 1]))); + //Limit NW,NE,SW,SE Color differences + nwvv = median(LVFU(chr[c][(indx - v1 - h1) >> 1]), LVFU(chr[c][(indx - v3 - h1) >> 1]), LVFU(chr[c][(indx - v1 - h3) >> 1])); + nevv = median(LVFU(chr[c][(indx - v1 + h1) >> 1]), LVFU(chr[c][(indx - v3 + h1) >> 1]), LVFU(chr[c][(indx - v1 + h3) >> 1])); + swvv = median(LVFU(chr[c][(indx + v1 - h1) >> 1]), LVFU(chr[c][(indx + v3 - h1) >> 1]), LVFU(chr[c][(indx + v1 - h3) >> 1])); + sevv = median(LVFU(chr[c][(indx + v1 + h1) >> 1]), LVFU(chr[c][(indx + v3 + h1) >> 1]), LVFU(chr[c][(indx + v1 + h3) >> 1])); + //Interpolate chrominance: R@B and B@R + tempv = (nwgv * nwvv + negv * nevv + swgv * swvv + segv * sevv) / (nwgv + negv + swgv + segv); + _mm_storeu_ps( &(chr[c][indx >> 1]), tempv); + } + + for (; col < width - 7; col += 2, indx += 2) { + //NW,NE,SW,SE Gradients + nwg = 1.0f / (eps + fabsf(chr[c][(indx - v1 - h1) >> 1] - chr[c][(indx - v3 - h3) >> 1]) + fabsf(chr[c][(indx + v1 + h1) >> 1] - chr[c][(indx - v3 - h3) >> 1])); + neg = 1.0f / (eps + fabsf(chr[c][(indx - v1 + h1) >> 1] - chr[c][(indx - v3 + h3) >> 1]) + fabsf(chr[c][(indx + v1 - h1) >> 1] - chr[c][(indx - v3 + h3) >> 1])); + swg = 1.0f / (eps + fabsf(chr[c][(indx + v1 - h1) >> 1] - chr[c][(indx + v3 + h3) >> 1]) + fabsf(chr[c][(indx - v1 + h1) >> 1] - chr[c][(indx + v3 - h3) >> 1])); + seg = 1.0f / (eps + fabsf(chr[c][(indx + v1 + h1) >> 1] - chr[c][(indx + v3 - h3) >> 1]) + fabsf(chr[c][(indx - v1 - h1) >> 1] - chr[c][(indx + v3 + h3) >> 1])); + //Limit NW,NE,SW,SE Color differences + nwv = median(chr[c][(indx - v1 - h1) >> 1], chr[c][(indx - v3 - h1) >> 1], chr[c][(indx - v1 - h3) >> 1]); + nev = median(chr[c][(indx - v1 + h1) >> 1], chr[c][(indx - v3 + h1) >> 1], chr[c][(indx - v1 + h3) >> 1]); + swv = median(chr[c][(indx + v1 - h1) >> 1], chr[c][(indx + v3 - h1) >> 1], chr[c][(indx + v1 - h3) >> 1]); + sev = median(chr[c][(indx + v1 + h1) >> 1], chr[c][(indx + v3 + h1) >> 1], chr[c][(indx + v1 + h3) >> 1]); + //Interpolate chrominance: R@B and B@R + chr[c][indx >> 1] = (nwg * nwv + neg * nev + swg * swv + seg * sev) / (nwg + neg + swg + seg); + } + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.65); + } + } +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 7; row < height - 7; row++) { + int col, indx; + + for (col = 7 + (FC(row, 0) & 1), indx = row * width + col; col < width - 14; col += 8, indx += 8) { + //N,E,W,S Gradients + ngv = onev / (epsv + vabsf(LVFU(chr[0][(indx - v1) >> 1]) - LVFU(chr[0][(indx - v3) >> 1])) + vabsf(LVFU(chr[0][(indx + v1) >> 1]) - LVFU(chr[0][(indx - v3) >> 1]))); + egv = onev / (epsv + vabsf(LVFU(chr[0][(indx + h1) >> 1]) - LVFU(chr[0][(indx + h3) >> 1])) + vabsf(LVFU(chr[0][(indx - h1) >> 1]) - LVFU(chr[0][(indx + h3) >> 1]))); + wgv = onev / (epsv + vabsf(LVFU(chr[0][(indx - h1) >> 1]) - LVFU(chr[0][(indx - h3) >> 1])) + vabsf(LVFU(chr[0][(indx + h1) >> 1]) - LVFU(chr[0][(indx - h3) >> 1]))); + sgv = onev / (epsv + vabsf(LVFU(chr[0][(indx + v1) >> 1]) - LVFU(chr[0][(indx + v3) >> 1])) + vabsf(LVFU(chr[0][(indx - v1) >> 1]) - LVFU(chr[0][(indx + v3) >> 1]))); + //Interpolate chrominance: R@G and B@G + tempv = ((ngv * LVFU(chr[0][(indx - v1) >> 1]) + egv * LVFU(chr[0][(indx + h1) >> 1]) + wgv * LVFU(chr[0][(indx - h1) >> 1]) + sgv * LVFU(chr[0][(indx + v1) >> 1])) / (ngv + egv + wgv + sgv)); + _mm_storeu_ps( &chr[0 + 2][indx >> 1], tempv); + } + + for (; col < width - 7; col += 2, indx += 2) { + //N,E,W,S Gradients + ng = 1.0f / (eps + fabsf(chr[0][(indx - v1) >> 1] - chr[0][(indx - v3) >> 1]) + fabsf(chr[0][(indx + v1) >> 1] - chr[0][(indx - v3) >> 1])); + eg = 1.0f / (eps + fabsf(chr[0][(indx + h1) >> 1] - chr[0][(indx + h3) >> 1]) + fabsf(chr[0][(indx - h1) >> 1] - chr[0][(indx + h3) >> 1])); + wg = 1.0f / (eps + fabsf(chr[0][(indx - h1) >> 1] - chr[0][(indx - h3) >> 1]) + fabsf(chr[0][(indx + h1) >> 1] - chr[0][(indx - h3) >> 1])); + sg = 1.0f / (eps + fabsf(chr[0][(indx + v1) >> 1] - chr[0][(indx + v3) >> 1]) + fabsf(chr[0][(indx - v1) >> 1] - chr[0][(indx + v3) >> 1])); + //Interpolate chrominance: R@G and B@G + chr[0 + 2][indx >> 1] = ((ng * chr[0][(indx - v1) >> 1] + eg * chr[0][(indx + h1) >> 1] + wg * chr[0][(indx - h1) >> 1] + sg * chr[0][(indx + v1) >> 1]) / (ng + eg + wg + sg)); + } + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.78); + } + } +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 7; row < height - 7; row++) { + int col, indx; + + for (col = 7 + (FC(row, 0) & 1), indx = row * width + col; col < width - 14; col += 8, indx += 8) { + //N,E,W,S Gradients + ngv = onev / (epsv + vabsf(LVFU(chr[1][(indx - v1) >> 1]) - LVFU(chr[1][(indx - v3) >> 1])) + vabsf(LVFU(chr[1][(indx + v1) >> 1]) - LVFU(chr[1][(indx - v3) >> 1]))); + egv = onev / (epsv + vabsf(LVFU(chr[1][(indx + h1) >> 1]) - LVFU(chr[1][(indx + h3) >> 1])) + vabsf(LVFU(chr[1][(indx - h1) >> 1]) - LVFU(chr[1][(indx + h3) >> 1]))); + wgv = onev / (epsv + vabsf(LVFU(chr[1][(indx - h1) >> 1]) - LVFU(chr[1][(indx - h3) >> 1])) + vabsf(LVFU(chr[1][(indx + h1) >> 1]) - LVFU(chr[1][(indx - h3) >> 1]))); + sgv = onev / (epsv + vabsf(LVFU(chr[1][(indx + v1) >> 1]) - LVFU(chr[1][(indx + v3) >> 1])) + vabsf(LVFU(chr[1][(indx - v1) >> 1]) - LVFU(chr[1][(indx + v3) >> 1]))); + //Interpolate chrominance: R@G and B@G + tempv = ((ngv * LVFU(chr[1][(indx - v1) >> 1]) + egv * LVFU(chr[1][(indx + h1) >> 1]) + wgv * LVFU(chr[1][(indx - h1) >> 1]) + sgv * LVFU(chr[1][(indx + v1) >> 1])) / (ngv + egv + wgv + sgv)); + _mm_storeu_ps( &chr[1 + 2][indx >> 1], tempv); + } + + for (; col < width - 7; col += 2, indx += 2) { + //N,E,W,S Gradients + ng = 1.0f / (eps + fabsf(chr[1][(indx - v1) >> 1] - chr[1][(indx - v3) >> 1]) + fabsf(chr[1][(indx + v1) >> 1] - chr[1][(indx - v3) >> 1])); + eg = 1.0f / (eps + fabsf(chr[1][(indx + h1) >> 1] - chr[1][(indx + h3) >> 1]) + fabsf(chr[1][(indx - h1) >> 1] - chr[1][(indx + h3) >> 1])); + wg = 1.0f / (eps + fabsf(chr[1][(indx - h1) >> 1] - chr[1][(indx - h3) >> 1]) + fabsf(chr[1][(indx + h1) >> 1] - chr[1][(indx - h3) >> 1])); + sg = 1.0f / (eps + fabsf(chr[1][(indx + v1) >> 1] - chr[1][(indx + v3) >> 1]) + fabsf(chr[1][(indx - v1) >> 1] - chr[1][(indx + v3) >> 1])); + //Interpolate chrominance: R@G and B@G + chr[1 + 2][indx >> 1] = ((ng * chr[1][(indx - v1) >> 1] + eg * chr[1][(indx + h1) >> 1] + wg * chr[1][(indx - h1) >> 1] + sg * chr[1][(indx + v1) >> 1]) / (ng + eg + wg + sg)); + } + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.91); + } + } + float *src1, *src2, *redsrc0, *redsrc1, *bluesrc0, *bluesrc1; +#ifdef _OPENMP + #pragma omp for +#endif + + for(int row = 7; row < height - 7; row++) { + int col, indx, fc; + fc = FC(row, 7) & 1; + src1 = rgb[fc]; + src2 = rgb[fc ^ 1]; + redsrc0 = chr[fc << 1]; + redsrc1 = chr[(fc ^ 1) << 1]; + bluesrc0 = chr[(fc << 1) + 1]; + bluesrc1 = chr[((fc ^ 1) << 1) + 1]; + + for(col = 7, indx = row * width + col; col < width - 14; col += 8, indx += 8) { + temp1v = LVFU( src1[indx >> 1] ); + temp2v = LVFU( src2[(indx + 1) >> 1] ); + tempv = _mm_shuffle_ps( temp1v, temp2v, _MM_SHUFFLE( 1, 0, 1, 0 ) ); + tempv = _mm_shuffle_ps( tempv, tempv, _MM_SHUFFLE( 3, 1, 2, 0 ) ); + _mm_storeu_ps( &green[row][col], CLIPV( tempv )); + temp5v = LVFU(redsrc0[indx >> 1]); + temp6v = LVFU(redsrc1[(indx + 1) >> 1]); + temp3v = _mm_shuffle_ps( temp5v, temp6v, _MM_SHUFFLE( 1, 0, 1, 0 ) ); + temp3v = _mm_shuffle_ps( temp3v, temp3v, _MM_SHUFFLE( 3, 1, 2, 0 ) ); + temp3v = CLIPV( tempv - c65535v * temp3v ); + _mm_storeu_ps( &red[row][col], temp3v); + temp7v = LVFU(bluesrc0[indx >> 1]); + temp8v = LVFU(bluesrc1[(indx + 1) >> 1]); + temp4v = _mm_shuffle_ps( temp7v, temp8v, _MM_SHUFFLE( 1, 0, 1, 0 ) ); + temp4v = _mm_shuffle_ps( temp4v, temp4v, _MM_SHUFFLE( 3, 1, 2, 0 ) ); + temp4v = CLIPV( tempv - c65535v * temp4v ); + _mm_storeu_ps( &blue[row][col], temp4v); + + tempv = _mm_shuffle_ps( temp1v, temp2v, _MM_SHUFFLE( 3, 2, 3, 2 ) ); + tempv = _mm_shuffle_ps( tempv, tempv, _MM_SHUFFLE( 3, 1, 2, 0 ) ); + _mm_storeu_ps( &green[row][col + 4], CLIPV( tempv )); + + temp3v = _mm_shuffle_ps( temp5v, temp6v, _MM_SHUFFLE( 3, 2, 3, 2 ) ); + temp3v = _mm_shuffle_ps( temp3v, temp3v, _MM_SHUFFLE( 3, 1, 2, 0 ) ); + temp3v = CLIPV( tempv - c65535v * temp3v ); + _mm_storeu_ps( &red[row][col + 4], temp3v); + temp4v = _mm_shuffle_ps( temp7v, temp8v, _MM_SHUFFLE( 3, 2, 3, 2 ) ); + temp4v = _mm_shuffle_ps( temp4v, temp4v, _MM_SHUFFLE( 3, 1, 2, 0 ) ); + temp4v = CLIPV( tempv - c65535v * temp4v ); + _mm_storeu_ps( &blue[row][col + 4], temp4v); + } + + for(; col < width - 7; col++, indx += 2) { + red [row][col] = CLIP(src1[indx >> 1] - 65535.f * redsrc0[indx >> 1]); + green[row][col] = CLIP(src1[indx >> 1]); + blue [row][col] = CLIP(src1[indx >> 1] - 65535.f * bluesrc0[indx >> 1]); + col++; + red [row][col] = CLIP(src2[(indx + 1) >> 1] - 65535.f * redsrc1[(indx + 1) >> 1]); + green[row][col] = CLIP(src2[(indx + 1) >> 1]); + blue [row][col] = CLIP(src2[(indx + 1) >> 1] - 65535.f * bluesrc1[(indx + 1) >> 1]); + } + } + }// End of parallelization + + if (plistener) { + plistener->setProgress (1.0); + } + + free(chrarray); + free(rgbarray); + free(vdif); + free(hdif); +} +#undef CLIPV +#else +void RawImageSource::igv_interpolate(int winw, int winh) +{ + static const float eps = 1e-5f, epssq = 1e-5f; //mod epssq -10f =>-5f Jacques 3/2013 to prevent artifact (divide by zero) + static const int h1 = 1, h2 = 2, h3 = 3, h4 = 4, h5 = 5, h6 = 6; + const int width = winw, height = winh; + const int v1 = 1 * width, v2 = 2 * width, v3 = 3 * width, v4 = 4 * width, v5 = 5 * width, v6 = 6 * width; + float* rgb[3]; + float* chr[2]; + float (*rgbarray), *vdif, *hdif, (*chrarray); + + rgbarray = (float (*)) calloc(width * height * 3, sizeof( float)); + rgb[0] = rgbarray; + rgb[1] = rgbarray + (width * height); + rgb[2] = rgbarray + 2 * (width * height); + + chrarray = (float (*)) calloc(width * height * 2, sizeof( float)); + chr[0] = chrarray; + chr[1] = chrarray + (width * height); + + vdif = (float (*)) calloc(width * height / 2, sizeof * vdif); + hdif = (float (*)) calloc(width * height / 2, sizeof * hdif); + + border_interpolate2(winw, winh, 7); + + if (plistener) { + plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::IGV))); + plistener->setProgress (0.0); + } + +#ifdef _OPENMP + #pragma omp parallel default(none) shared(rgb,vdif,hdif,chr) +#endif + { + + float ng, eg, wg, sg, nv, ev, wv, sv, nwg, neg, swg, seg, nwv, nev, swv, sev; + +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 0; row < height - 0; row++) + for (int col = 0, indx = row * width + col; col < width - 0; col++, indx++) { + int c = FC(row, col); + rgb[c][indx] = CLIP(rawData[row][col]); //rawData = RT datas + } + +// border_interpolate2(7, rgb); + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.13); + } + } + +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 5; row < height - 5; row++) + for (int col = 5 + (FC(row, 1) & 1), indx = row * width + col, c = FC(row, col); col < width - 5; col += 2, indx += 2) { + //N,E,W,S Gradients + ng = (eps + (fabsf(rgb[1][indx - v1] - rgb[1][indx - v3]) + fabsf(rgb[c][indx] - rgb[c][indx - v2])) / 65535.f);; + eg = (eps + (fabsf(rgb[1][indx + h1] - rgb[1][indx + h3]) + fabsf(rgb[c][indx] - rgb[c][indx + h2])) / 65535.f); + wg = (eps + (fabsf(rgb[1][indx - h1] - rgb[1][indx - h3]) + fabsf(rgb[c][indx] - rgb[c][indx - h2])) / 65535.f); + sg = (eps + (fabsf(rgb[1][indx + v1] - rgb[1][indx + v3]) + fabsf(rgb[c][indx] - rgb[c][indx + v2])) / 65535.f); + //N,E,W,S High Order Interpolation (Li & Randhawa) + //N,E,W,S Hamilton Adams Interpolation + // (48.f * 65535.f) = 3145680.f + nv = LIM(((23.0f * rgb[1][indx - v1] + 23.0f * rgb[1][indx - v3] + rgb[1][indx - v5] + rgb[1][indx + v1] + 40.0f * rgb[c][indx] - 32.0f * rgb[c][indx - v2] - 8.0f * rgb[c][indx - v4])) / 3145680.f, 0.0f, 1.0f); + ev = LIM(((23.0f * rgb[1][indx + h1] + 23.0f * rgb[1][indx + h3] + rgb[1][indx + h5] + rgb[1][indx - h1] + 40.0f * rgb[c][indx] - 32.0f * rgb[c][indx + h2] - 8.0f * rgb[c][indx + h4])) / 3145680.f, 0.0f, 1.0f); + wv = LIM(((23.0f * rgb[1][indx - h1] + 23.0f * rgb[1][indx - h3] + rgb[1][indx - h5] + rgb[1][indx + h1] + 40.0f * rgb[c][indx] - 32.0f * rgb[c][indx - h2] - 8.0f * rgb[c][indx - h4])) / 3145680.f, 0.0f, 1.0f); + sv = LIM(((23.0f * rgb[1][indx + v1] + 23.0f * rgb[1][indx + v3] + rgb[1][indx + v5] + rgb[1][indx - v1] + 40.0f * rgb[c][indx] - 32.0f * rgb[c][indx + v2] - 8.0f * rgb[c][indx + v4])) / 3145680.f, 0.0f, 1.0f); + //Horizontal and vertical color differences + vdif[indx >> 1] = (sg * nv + ng * sv) / (ng + sg) - (rgb[c][indx]) / 65535.f; + hdif[indx >> 1] = (wg * ev + eg * wv) / (eg + wg) - (rgb[c][indx]) / 65535.f; + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.26); + } + } + +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 7; row < height - 7; row++) + for (int col = 7 + (FC(row, 1) & 1), indx = row * width + col, c = FC(row, col), d = c / 2; col < width - 7; col += 2, indx += 2) { + //H&V integrated gaussian vector over variance on color differences + //Mod Jacques 3/2013 + ng = LIM(epssq + 78.0f * SQR(vdif[indx >> 1]) + 69.0f * (SQR(vdif[(indx - v2) >> 1]) + SQR(vdif[(indx + v2) >> 1])) + 51.0f * (SQR(vdif[(indx - v4) >> 1]) + SQR(vdif[(indx + v4) >> 1])) + 21.0f * (SQR(vdif[(indx - v6) >> 1]) + SQR(vdif[(indx + v6) >> 1])) - 6.0f * SQR(vdif[(indx - v2) >> 1] + vdif[indx >> 1] + vdif[(indx + v2) >> 1]) + - 10.0f * (SQR(vdif[(indx - v4) >> 1] + vdif[(indx - v2) >> 1] + vdif[indx >> 1]) + SQR(vdif[indx >> 1] + vdif[(indx + v2) >> 1] + vdif[(indx + v4) >> 1])) - 7.0f * (SQR(vdif[(indx - v6) >> 1] + vdif[(indx - v4) >> 1] + vdif[(indx - v2) >> 1]) + SQR(vdif[(indx + v2) >> 1] + vdif[(indx + v4) >> 1] + vdif[(indx + v6) >> 1])), 0.f, 1.f); + eg = LIM(epssq + 78.0f * SQR(hdif[indx >> 1]) + 69.0f * (SQR(hdif[(indx - h2) >> 1]) + SQR(hdif[(indx + h2) >> 1])) + 51.0f * (SQR(hdif[(indx - h4) >> 1]) + SQR(hdif[(indx + h4) >> 1])) + 21.0f * (SQR(hdif[(indx - h6) >> 1]) + SQR(hdif[(indx + h6) >> 1])) - 6.0f * SQR(hdif[(indx - h2) >> 1] + hdif[indx >> 1] + hdif[(indx + h2) >> 1]) + - 10.0f * (SQR(hdif[(indx - h4) >> 1] + hdif[(indx - h2) >> 1] + hdif[indx >> 1]) + SQR(hdif[indx >> 1] + hdif[(indx + h2) >> 1] + hdif[(indx + h4) >> 1])) - 7.0f * (SQR(hdif[(indx - h6) >> 1] + hdif[(indx - h4) >> 1] + hdif[(indx - h2) >> 1]) + SQR(hdif[(indx + h2) >> 1] + hdif[(indx + h4) >> 1] + hdif[(indx + h6) >> 1])), 0.f, 1.f); + //Limit chrominance using H/V neighbourhood + nv = median(0.725f * vdif[indx >> 1] + 0.1375f * vdif[(indx - v2) >> 1] + 0.1375f * vdif[(indx + v2) >> 1], vdif[(indx - v2) >> 1], vdif[(indx + v2) >> 1]); + ev = median(0.725f * hdif[indx >> 1] + 0.1375f * hdif[(indx - h2) >> 1] + 0.1375f * hdif[(indx + h2) >> 1], hdif[(indx - h2) >> 1], hdif[(indx + h2) >> 1]); + //Chrominance estimation + chr[d][indx] = (eg * nv + ng * ev) / (ng + eg); + //Green channel population + rgb[1][indx] = rgb[c][indx] + 65535.f * chr[d][indx]; + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.39); + } + } + +// free(vdif); free(hdif); +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 7; row < height - 7; row += 2) + for (int col = 7 + (FC(row, 1) & 1), indx = row * width + col, c = 1 - FC(row, col) / 2; col < width - 7; col += 2, indx += 2) { + //NW,NE,SW,SE Gradients + nwg = 1.0f / (eps + fabsf(chr[c][indx - v1 - h1] - chr[c][indx - v3 - h3]) + fabsf(chr[c][indx + v1 + h1] - chr[c][indx - v3 - h3])); + neg = 1.0f / (eps + fabsf(chr[c][indx - v1 + h1] - chr[c][indx - v3 + h3]) + fabsf(chr[c][indx + v1 - h1] - chr[c][indx - v3 + h3])); + swg = 1.0f / (eps + fabsf(chr[c][indx + v1 - h1] - chr[c][indx + v3 + h3]) + fabsf(chr[c][indx - v1 + h1] - chr[c][indx + v3 - h3])); + seg = 1.0f / (eps + fabsf(chr[c][indx + v1 + h1] - chr[c][indx + v3 - h3]) + fabsf(chr[c][indx - v1 - h1] - chr[c][indx + v3 + h3])); + //Limit NW,NE,SW,SE Color differences + nwv = median(chr[c][indx - v1 - h1], chr[c][indx - v3 - h1], chr[c][indx - v1 - h3]); + nev = median(chr[c][indx - v1 + h1], chr[c][indx - v3 + h1], chr[c][indx - v1 + h3]); + swv = median(chr[c][indx + v1 - h1], chr[c][indx + v3 - h1], chr[c][indx + v1 - h3]); + sev = median(chr[c][indx + v1 + h1], chr[c][indx + v3 + h1], chr[c][indx + v1 + h3]); + //Interpolate chrominance: R@B and B@R + chr[c][indx] = (nwg * nwv + neg * nev + swg * swv + seg * sev) / (nwg + neg + swg + seg); + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.52); + } + } +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 8; row < height - 7; row += 2) + for (int col = 7 + (FC(row, 1) & 1), indx = row * width + col, c = 1 - FC(row, col) / 2; col < width - 7; col += 2, indx += 2) { + //NW,NE,SW,SE Gradients + nwg = 1.0f / (eps + fabsf(chr[c][indx - v1 - h1] - chr[c][indx - v3 - h3]) + fabsf(chr[c][indx + v1 + h1] - chr[c][indx - v3 - h3])); + neg = 1.0f / (eps + fabsf(chr[c][indx - v1 + h1] - chr[c][indx - v3 + h3]) + fabsf(chr[c][indx + v1 - h1] - chr[c][indx - v3 + h3])); + swg = 1.0f / (eps + fabsf(chr[c][indx + v1 - h1] - chr[c][indx + v3 + h3]) + fabsf(chr[c][indx - v1 + h1] - chr[c][indx + v3 - h3])); + seg = 1.0f / (eps + fabsf(chr[c][indx + v1 + h1] - chr[c][indx + v3 - h3]) + fabsf(chr[c][indx - v1 - h1] - chr[c][indx + v3 + h3])); + //Limit NW,NE,SW,SE Color differences + nwv = median(chr[c][indx - v1 - h1], chr[c][indx - v3 - h1], chr[c][indx - v1 - h3]); + nev = median(chr[c][indx - v1 + h1], chr[c][indx - v3 + h1], chr[c][indx - v1 + h3]); + swv = median(chr[c][indx + v1 - h1], chr[c][indx + v3 - h1], chr[c][indx + v1 - h3]); + sev = median(chr[c][indx + v1 + h1], chr[c][indx + v3 + h1], chr[c][indx + v1 + h3]); + //Interpolate chrominance: R@B and B@R + chr[c][indx] = (nwg * nwv + neg * nev + swg * swv + seg * sev) / (nwg + neg + swg + seg); + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.65); + } + } +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 7; row < height - 7; row++) + for (int col = 7 + (FC(row, 0) & 1), indx = row * width + col; col < width - 7; col += 2, indx += 2) { + //N,E,W,S Gradients + ng = 1.0f / (eps + fabsf(chr[0][indx - v1] - chr[0][indx - v3]) + fabsf(chr[0][indx + v1] - chr[0][indx - v3])); + eg = 1.0f / (eps + fabsf(chr[0][indx + h1] - chr[0][indx + h3]) + fabsf(chr[0][indx - h1] - chr[0][indx + h3])); + wg = 1.0f / (eps + fabsf(chr[0][indx - h1] - chr[0][indx - h3]) + fabsf(chr[0][indx + h1] - chr[0][indx - h3])); + sg = 1.0f / (eps + fabsf(chr[0][indx + v1] - chr[0][indx + v3]) + fabsf(chr[0][indx - v1] - chr[0][indx + v3])); + //Interpolate chrominance: R@G and B@G + chr[0][indx] = ((ng * chr[0][indx - v1] + eg * chr[0][indx + h1] + wg * chr[0][indx - h1] + sg * chr[0][indx + v1]) / (ng + eg + wg + sg)); + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.78); + } + } +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 7; row < height - 7; row++) + for (int col = 7 + (FC(row, 0) & 1), indx = row * width + col; col < width - 7; col += 2, indx += 2) { + + //N,E,W,S Gradients + ng = 1.0f / (eps + fabsf(chr[1][indx - v1] - chr[1][indx - v3]) + fabsf(chr[1][indx + v1] - chr[1][indx - v3])); + eg = 1.0f / (eps + fabsf(chr[1][indx + h1] - chr[1][indx + h3]) + fabsf(chr[1][indx - h1] - chr[1][indx + h3])); + wg = 1.0f / (eps + fabsf(chr[1][indx - h1] - chr[1][indx - h3]) + fabsf(chr[1][indx + h1] - chr[1][indx - h3])); + sg = 1.0f / (eps + fabsf(chr[1][indx + v1] - chr[1][indx + v3]) + fabsf(chr[1][indx - v1] - chr[1][indx + v3])); + //Interpolate chrominance: R@G and B@G + chr[1][indx] = ((ng * chr[1][indx - v1] + eg * chr[1][indx + h1] + wg * chr[1][indx - h1] + sg * chr[1][indx + v1]) / (ng + eg + wg + sg)); + } + +#ifdef _OPENMP + #pragma omp single +#endif + { + if (plistener) { + plistener->setProgress (0.91); + } + + //Interpolate borders +// border_interpolate2(7, rgb); + } + /* +#ifdef _OPENMP + #pragma omp for +#endif + for (int row=0; row < height; row++) //borders + for (int col=0; col < width; col++) { + if (col==7 && row >= 7 && row < height-7) + col = width-7; + int indxc=row*width+col; + red [row][col] = rgb[indxc][0]; + green[row][col] = rgb[indxc][1]; + blue [row][col] = rgb[indxc][2]; + } + */ + +#ifdef _OPENMP + #pragma omp for +#endif + + for(int row = 7; row < height - 7; row++) + for(int col = 7, indx = row * width + col; col < width - 7; col++, indx++) { + red [row][col] = CLIP(rgb[1][indx] - 65535.f * chr[0][indx]); + green[row][col] = CLIP(rgb[1][indx]); + blue [row][col] = CLIP(rgb[1][indx] - 65535.f * chr[1][indx]); + } + }// End of parallelization + + if (plistener) { + plistener->setProgress (1.0); + } + + free(chrarray); + free(rgbarray); + free(vdif); + free(hdif); +} +#endif + + +/* + Adaptive Homogeneity-Directed interpolation is based on + the work of Keigo Hirakawa, Thomas Parks, and Paul Lee. + */ +#define TS 256 /* Tile Size */ +#define FORC(cnt) for (c=0; c < cnt; c++) +#define FORC3 FORC(3) + +void RawImageSource::ahd_demosaic() +{ + int i, j, k, top, left, row, col, tr, tc, c, d, val, hm[2]; + float (*pix)[4], (*rix)[3]; + static const int dir[4] = { -1, 1, -TS, TS }; + float ldiff[2][4], abdiff[2][4], leps, abeps; + float xyz[3], xyz_cam[3][4]; + float* cbrt; + float (*rgb)[TS][TS][3]; + float (*lab)[TS][TS][3]; + float (*lix)[3]; + char (*homo)[TS][TS], *buffer; + double r; + + int width = W, height = H; + float (*image)[4]; + unsigned int colors = 3; + + const double xyz_rgb[3][3] = { /* XYZ from RGB */ + { 0.412453, 0.357580, 0.180423 }, + { 0.212671, 0.715160, 0.072169 }, + { 0.019334, 0.119193, 0.950227 } + }; + + const float d65_white[3] = { 0.950456, 1, 1.088754 }; + + if (plistener) { + plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::AHD))); + plistener->setProgress (0.0); + } + + image = (float (*)[4]) calloc (H * W, sizeof * image); + + for (int ii = 0; ii < H; ii++) + for (int jj = 0; jj < W; jj++) { + image[ii * W + jj][fc(ii, jj)] = rawData[ii][jj]; + } + + cbrt = (float (*)) calloc (0x10000, sizeof * cbrt); + + for (i = 0; i < 0x10000; i++) { + r = (double)i / 65535.0; + cbrt[i] = r > 0.008856 ? std::cbrt(r) : 7.787 * r + 16 / 116.0; + } + + for (i = 0; i < 3; i++) + for (unsigned int j = 0; j < colors; j++) + for (xyz_cam[i][j] = k = 0; k < 3; k++) { + xyz_cam[i][j] += xyz_rgb[i][k] * imatrices.rgb_cam[k][j] / d65_white[i]; + } + + border_interpolate(5, image); + buffer = (char *) malloc (13 * TS * TS * sizeof(float)); /* 1664 kB */ + //merror (buffer, "ahd_interpolate()"); + rgb = (float(*)[TS][TS][3]) buffer; + lab = (float(*)[TS][TS][3])(buffer + 6 * TS * TS * sizeof(float)); + homo = (char (*)[TS][TS]) (buffer + 12 * TS * TS * sizeof(float)); + + // helper variables for progress indication + int n_tiles = ((height - 7 + (TS - 7)) / (TS - 6)) * ((width - 7 + (TS - 7)) / (TS - 6)); + int tile = 0; + + for (top = 2; top < height - 5; top += TS - 6) + for (left = 2; left < width - 5; left += TS - 6) { + /* Interpolate green horizontally and vertically: */ + for (row = top; row < top + TS && row < height - 2; row++) { + col = left + (FC(row, left) & 1); + + for (c = FC(row, col); col < left + TS && col < width - 2; col += 2) { + pix = image + (row * width + col); + val = 0.25 * ((pix[-1][1] + pix[0][c] + pix[1][1]) * 2 + - pix[-2][c] - pix[2][c]) ; + rgb[0][row - top][col - left][1] = median(static_cast(val), pix[-1][1], pix[1][1]); + val = 0.25 * ((pix[-width][1] + pix[0][c] + pix[width][1]) * 2 + - pix[-2 * width][c] - pix[2 * width][c]) ; + rgb[1][row - top][col - left][1] = median(static_cast(val), pix[-width][1], pix[width][1]); + } + } + + /* Interpolate red and blue, and convert to CIELab: */ + for (d = 0; d < 2; d++) + for (row = top + 1; row < top + TS - 1 && row < height - 3; row++) + for (col = left + 1; col < left + TS - 1 && col < width - 3; col++) { + pix = image + (row * width + col); + rix = &rgb[d][row - top][col - left]; + lix = &lab[d][row - top][col - left]; + + if ((c = 2 - FC(row, col)) == 1) { + c = FC(row + 1, col); + val = pix[0][1] + (0.5 * ( pix[-1][2 - c] + pix[1][2 - c] + - rix[-1][1] - rix[1][1] ) ); + rix[0][2 - c] = CLIP(val); + val = pix[0][1] + (0.5 * ( pix[-width][c] + pix[width][c] + - rix[-TS][1] - rix[TS][1] ) ); + } else + val = rix[0][1] + (0.25 * ( pix[-width - 1][c] + pix[-width + 1][c] + + pix[+width - 1][c] + pix[+width + 1][c] + - rix[-TS - 1][1] - rix[-TS + 1][1] + - rix[+TS - 1][1] - rix[+TS + 1][1]) ); + + rix[0][c] = CLIP(val); + c = FC(row, col); + rix[0][c] = pix[0][c]; + xyz[0] = xyz[1] = xyz[2] = 0.0; + FORCC { + xyz[0] += xyz_cam[0][c] * rix[0][c]; + xyz[1] += xyz_cam[1][c] * rix[0][c]; + xyz[2] += xyz_cam[2][c] * rix[0][c]; + } + + xyz[0] = CurveFactory::flinterp(cbrt, xyz[0]); + xyz[1] = CurveFactory::flinterp(cbrt, xyz[1]); + xyz[2] = CurveFactory::flinterp(cbrt, xyz[2]); + + //xyz[0] = xyz[0] > 0.008856 ? pow(xyz[0]/65535,1/3.0) : 7.787*xyz[0] + 16/116.0; + //xyz[1] = xyz[1] > 0.008856 ? pow(xyz[1]/65535,1/3.0) : 7.787*xyz[1] + 16/116.0; + //xyz[2] = xyz[2] > 0.008856 ? pow(xyz[2]/65535,1/3.0) : 7.787*xyz[2] + 16/116.0; + + lix[0][0] = (116 * xyz[1] - 16); + lix[0][1] = 500 * (xyz[0] - xyz[1]); + lix[0][2] = 200 * (xyz[1] - xyz[2]); + } + + /* Build homogeneity maps from the CIELab images: */ + memset (homo, 0, 2 * TS * TS); + + for (row = top + 2; row < top + TS - 2 && row < height - 4; row++) { + tr = row - top; + + for (col = left + 2; col < left + TS - 2 && col < width - 4; col++) { + tc = col - left; + + for (d = 0; d < 2; d++) { + lix = &lab[d][tr][tc]; + + for (i = 0; i < 4; i++) { + ldiff[d][i] = ABS(lix[0][0] - lix[dir[i]][0]); + abdiff[d][i] = SQR(lix[0][1] - lix[dir[i]][1]) + + SQR(lix[0][2] - lix[dir[i]][2]); + } + } + + leps = min(max(ldiff[0][0], ldiff[0][1]), + max(ldiff[1][2], ldiff[1][3])); + abeps = min(max(abdiff[0][0], abdiff[0][1]), + max(abdiff[1][2], abdiff[1][3])); + + for (d = 0; d < 2; d++) + for (i = 0; i < 4; i++) + if (ldiff[d][i] <= leps && abdiff[d][i] <= abeps) { + homo[d][tr][tc]++; + } + } + } + + /* Combine the most homogenous pixels for the final result: */ + for (row = top + 3; row < top + TS - 3 && row < height - 5; row++) { + tr = row - top; + + for (col = left + 3; col < left + TS - 3 && col < width - 5; col++) { + tc = col - left; + + for (d = 0; d < 2; d++) + for (hm[d] = 0, i = tr - 1; i <= tr + 1; i++) + for (j = tc - 1; j <= tc + 1; j++) { + hm[d] += homo[d][i][j]; + } + + if (hm[0] != hm[1]) { + FORC3 image[row * width + col][c] = rgb[hm[1] > hm[0]][tr][tc][c]; + } else + FORC3 image[row * width + col][c] = + 0.5 * (rgb[0][tr][tc][c] + rgb[1][tr][tc][c]) ; + } + } + + tile++; + + if(plistener) { + plistener->setProgress((double)tile / n_tiles); + } + } + + if(plistener) { + plistener->setProgress (1.0); + } + + free (buffer); + + for (int i = 0; i < H; i++) { + for (int j = 0; j < W; j++) { + red[i][j] = image[i * W + j][0]; + green[i][j] = image[i * W + j][1]; + blue[i][j] = image[i * W + j][2]; + } + } + + free (image); + free (cbrt); +} +#undef TS + +void RawImageSource::nodemosaic(bool bw) +{ + red(W, H); + green(W, H); + blue(W, H); + #pragma omp parallel for + + for (int i = 0; i < H; i++) { + for (int j = 0; j < W; j++) { + if (bw) { + red[i][j] = green[i][j] = blue[i][j] = rawData[i][j]; + } else if(ri->getSensorType() != ST_FUJI_XTRANS) { + switch( FC(i, j)) { + case 0: + red[i][j] = rawData[i][j]; + green[i][j] = blue[i][j] = 0; + break; + + case 1: + green[i][j] = rawData[i][j]; + red[i][j] = blue[i][j] = 0; + break; + + case 2: + blue[i][j] = rawData[i][j]; + red[i][j] = green[i][j] = 0; + break; + } + } else { + switch( ri->XTRANSFC(i, j)) { + case 0: + red[i][j] = rawData[i][j]; + green[i][j] = blue[i][j] = 0; + break; + + case 1: + green[i][j] = rawData[i][j]; + red[i][j] = blue[i][j] = 0; + break; + + case 2: + blue[i][j] = rawData[i][j]; + red[i][j] = green[i][j] = 0; + break; + } + } + } + } +} + +/* + Refinement based on EECI demosaicing algorithm by L. Chang and Y.P. Tan + Paul Lee + Adapted for RawTherapee - Jacques Desmis 04/2013 +*/ + +#ifdef __SSE2__ +#define CLIPV(a) LIMV(a,ZEROV,c65535v) +#endif +void RawImageSource::refinement(int PassCount) +{ + MyTime t1e, t2e; + t1e.set(); + + int width = W; + int height = H; + int w1 = width; + int w2 = 2 * w1; + + if (plistener) { + plistener->setProgressStr (M("TP_RAW_DMETHOD_PROGRESSBAR_REFINE")); + } + + array2D *rgb[3]; + rgb[0] = &red; + rgb[1] = &green; + rgb[2] = &blue; + + for (int b = 0; b < PassCount; b++) { + if (plistener) { + plistener->setProgress ((float)b / PassCount); + } + + +#ifdef _OPENMP + #pragma omp parallel +#endif + { + float *pix[3]; + + /* Reinforce interpolated green pixels on RED/BLUE pixel locations */ +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 2; row < height - 2; row++) { + int col = 2 + (FC(row, 2) & 1); + int c = FC(row, col); +#ifdef __SSE2__ + __m128 dLv, dRv, dUv, dDv, v0v; + __m128 onev = _mm_set1_ps(1.f); + __m128 zd5v = _mm_set1_ps(0.5f); + __m128 c65535v = _mm_set1_ps(65535.f); + + for (; col < width - 8; col += 8) { + int indx = row * width + col; + pix[c] = (float*)(*rgb[c]) + indx; + pix[1] = (float*)(*rgb[1]) + indx; + dLv = onev / (onev + vabsf(LC2VFU(pix[c][ -2]) - LC2VFU(pix[c][0])) + vabsf(LC2VFU(pix[1][ 1]) - LC2VFU(pix[1][ -1]))); + dRv = onev / (onev + vabsf(LC2VFU(pix[c][ 2]) - LC2VFU(pix[c][0])) + vabsf(LC2VFU(pix[1][ 1]) - LC2VFU(pix[1][ -1]))); + dUv = onev / (onev + vabsf(LC2VFU(pix[c][-w2]) - LC2VFU(pix[c][0])) + vabsf(LC2VFU(pix[1][w1]) - LC2VFU(pix[1][-w1]))); + dDv = onev / (onev + vabsf(LC2VFU(pix[c][ w2]) - LC2VFU(pix[c][0])) + vabsf(LC2VFU(pix[1][w1]) - LC2VFU(pix[1][-w1]))); + v0v = CLIPV(LC2VFU(pix[c][0]) + zd5v + ((LC2VFU(pix[1][-1]) - LC2VFU(pix[c][-1])) * dLv + (LC2VFU(pix[1][1]) - LC2VFU(pix[c][1])) * dRv + (LC2VFU(pix[1][-w1]) - LC2VFU(pix[c][-w1])) * dUv + (LC2VFU(pix[1][w1]) - LC2VFU(pix[c][w1])) * dDv ) / (dLv + dRv + dUv + dDv)); + STC2VFU(pix[1][0], v0v); + } + +#endif + + for (; col < width - 2; col += 2) { + int indx = row * width + col; + pix[c] = (float*)(*rgb[c]) + indx; + pix[1] = (float*)(*rgb[1]) + indx; + float dL = 1.f / (1.f + fabsf(pix[c][ -2] - pix[c][0]) + fabsf(pix[1][ 1] - pix[1][ -1])); + float dR = 1.f / (1.f + fabsf(pix[c][ 2] - pix[c][0]) + fabsf(pix[1][ 1] - pix[1][ -1])); + float dU = 1.f / (1.f + fabsf(pix[c][-w2] - pix[c][0]) + fabsf(pix[1][w1] - pix[1][-w1])); + float dD = 1.f / (1.f + fabsf(pix[c][ w2] - pix[c][0]) + fabsf(pix[1][w1] - pix[1][-w1])); + float v0 = (pix[c][0] + 0.5f + ((pix[1][ -1] - pix[c][ -1]) * dL + (pix[1][ 1] - pix[c][ 1]) * dR + (pix[1][-w1] - pix[c][-w1]) * dU + (pix[1][ w1] - pix[c][ w1]) * dD ) / (dL + dR + dU + dD)); + pix[1][0] = CLIP(v0); + } + } + + /* Reinforce interpolated red/blue pixels on GREEN pixel locations */ +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 2; row < height - 2; row++) { + int col = 2 + (FC(row, 3) & 1); + int c = FC(row, col + 1); +#ifdef __SSE2__ + __m128 dLv, dRv, dUv, dDv, v0v; + __m128 onev = _mm_set1_ps(1.f); + __m128 zd5v = _mm_set1_ps(0.5f); + __m128 c65535v = _mm_set1_ps(65535.f); + + for (; col < width - 8; col += 8) { + int indx = row * width + col; + pix[1] = (float*)(*rgb[1]) + indx; + + for (int i = 0; i < 2; c = 2 - c, i++) { + pix[c] = (float*)(*rgb[c]) + indx; + dLv = onev / (onev + vabsf(LC2VFU(pix[1][ -2]) - LC2VFU(pix[1][0])) + vabsf(LC2VFU(pix[c][ 1]) - LC2VFU(pix[c][ -1]))); + dRv = onev / (onev + vabsf(LC2VFU(pix[1][ 2]) - LC2VFU(pix[1][0])) + vabsf(LC2VFU(pix[c][ 1]) - LC2VFU(pix[c][ -1]))); + dUv = onev / (onev + vabsf(LC2VFU(pix[1][-w2]) - LC2VFU(pix[1][0])) + vabsf(LC2VFU(pix[c][w1]) - LC2VFU(pix[c][-w1]))); + dDv = onev / (onev + vabsf(LC2VFU(pix[1][ w2]) - LC2VFU(pix[1][0])) + vabsf(LC2VFU(pix[c][w1]) - LC2VFU(pix[c][-w1]))); + v0v = CLIPV(LC2VFU(pix[1][0]) + zd5v - ((LC2VFU(pix[1][-1]) - LC2VFU(pix[c][-1])) * dLv + (LC2VFU(pix[1][1]) - LC2VFU(pix[c][1])) * dRv + (LC2VFU(pix[1][-w1]) - LC2VFU(pix[c][-w1])) * dUv + (LC2VFU(pix[1][w1]) - LC2VFU(pix[c][w1])) * dDv ) / (dLv + dRv + dUv + dDv)); + STC2VFU(pix[c][0], v0v); + } + } + +#endif + + for (; col < width - 2; col += 2) { + int indx = row * width + col; + pix[1] = (float*)(*rgb[1]) + indx; + + for (int i = 0; i < 2; c = 2 - c, i++) { + pix[c] = (float*)(*rgb[c]) + indx; + float dL = 1.f / (1.f + fabsf(pix[1][ -2] - pix[1][0]) + fabsf(pix[c][ 1] - pix[c][ -1])); + float dR = 1.f / (1.f + fabsf(pix[1][ 2] - pix[1][0]) + fabsf(pix[c][ 1] - pix[c][ -1])); + float dU = 1.f / (1.f + fabsf(pix[1][-w2] - pix[1][0]) + fabsf(pix[c][w1] - pix[c][-w1])); + float dD = 1.f / (1.f + fabsf(pix[1][ w2] - pix[1][0]) + fabsf(pix[c][w1] - pix[c][-w1])); + float v0 = (pix[1][0] + 0.5f - ((pix[1][ -1] - pix[c][ -1]) * dL + (pix[1][ 1] - pix[c][ 1]) * dR + (pix[1][-w1] - pix[c][-w1]) * dU + (pix[1][ w1] - pix[c][ w1]) * dD ) / (dL + dR + dU + dD)); + pix[c][0] = CLIP(v0); + } + } + } + + /* Reinforce integrated red/blue pixels on BLUE/RED pixel locations */ +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 2; row < height - 2; row++) { + int col = 2 + (FC(row, 2) & 1); + int c = 2 - FC(row, col); +#ifdef __SSE2__ + __m128 dLv, dRv, dUv, dDv, v0v; + __m128 onev = _mm_set1_ps(1.f); + __m128 zd5v = _mm_set1_ps(0.5f); + __m128 c65535v = _mm_set1_ps(65535.f); + + for (; col < width - 8; col += 8) { + int indx = row * width + col; + pix[0] = (float*)(*rgb[0]) + indx; + pix[1] = (float*)(*rgb[1]) + indx; + pix[2] = (float*)(*rgb[2]) + indx; + int d = 2 - c; + dLv = onev / (onev + vabsf(LC2VFU(pix[d][ -2]) - LC2VFU(pix[d][0])) + vabsf(LC2VFU(pix[1][ 1]) - LC2VFU(pix[1][ -1]))); + dRv = onev / (onev + vabsf(LC2VFU(pix[d][ 2]) - LC2VFU(pix[d][0])) + vabsf(LC2VFU(pix[1][ 1]) - LC2VFU(pix[1][ -1]))); + dUv = onev / (onev + vabsf(LC2VFU(pix[d][-w2]) - LC2VFU(pix[d][0])) + vabsf(LC2VFU(pix[1][w1]) - LC2VFU(pix[1][-w1]))); + dDv = onev / (onev + vabsf(LC2VFU(pix[d][ w2]) - LC2VFU(pix[d][0])) + vabsf(LC2VFU(pix[1][w1]) - LC2VFU(pix[1][-w1]))); + v0v = CLIPV(LC2VFU(pix[1][0]) + zd5v - ((LC2VFU(pix[1][-1]) - LC2VFU(pix[c][-1])) * dLv + (LC2VFU(pix[1][1]) - LC2VFU(pix[c][1])) * dRv + (LC2VFU(pix[1][-w1]) - LC2VFU(pix[c][-w1])) * dUv + (LC2VFU(pix[1][w1]) - LC2VFU(pix[c][w1])) * dDv ) / (dLv + dRv + dUv + dDv)); + STC2VFU(pix[c][0], v0v); + } + +#endif + + for (; col < width - 2; col += 2) { + int indx = row * width + col; + pix[0] = (float*)(*rgb[0]) + indx; + pix[1] = (float*)(*rgb[1]) + indx; + pix[2] = (float*)(*rgb[2]) + indx; + int d = 2 - c; + float dL = 1.f / (1.f + fabsf(pix[d][ -2] - pix[d][0]) + fabsf(pix[1][ 1] - pix[1][ -1])); + float dR = 1.f / (1.f + fabsf(pix[d][ 2] - pix[d][0]) + fabsf(pix[1][ 1] - pix[1][ -1])); + float dU = 1.f / (1.f + fabsf(pix[d][-w2] - pix[d][0]) + fabsf(pix[1][w1] - pix[1][-w1])); + float dD = 1.f / (1.f + fabsf(pix[d][ w2] - pix[d][0]) + fabsf(pix[1][w1] - pix[1][-w1])); + float v0 = (pix[1][0] + 0.5f - ((pix[1][ -1] - pix[c][ -1]) * dL + (pix[1][ 1] - pix[c][ 1]) * dR + (pix[1][-w1] - pix[c][-w1]) * dU + (pix[1][ w1] - pix[c][ w1]) * dD ) / (dL + dR + dU + dD)); + pix[c][0] = CLIP(v0); + } + } + } // end parallel + } + + t2e.set(); + + if (settings->verbose) { + printf("Refinement Lee %d usec\n", t2e.etime(t1e)); + } +} +#ifdef __SSE2__ +#undef CLIPV +#endif + + +// Refinement based on EECI demozaicing algorithm by L. Chang and Y.P. Tan +// from "Lassus" : Luis Sanz Rodriguez, adapted by Jacques Desmis - JDC - and Oliver Duis for RawTherapee +// increases the signal to noise ratio (PSNR) # +1 to +2 dB : tested with Dcraw : eg: Lighthouse + AMaZE : whitout refinement:39.96dB, with refinement:41.86 dB +// reduce color artifacts, improves the interpolation +// but it's relatively slow +// +// Should be DISABLED if it decreases image quality by increases some image noise and generates blocky edges +void RawImageSource::refinement_lassus(int PassCount) +{ + // const int PassCount=1; + + // if (settings->verbose) printf("Refinement\n"); + + MyTime t1e, t2e; + t1e.set(); + int u = W, v = 2 * u, w = 3 * u, x = 4 * u, y = 5 * u; + float (*image)[3]; + image = (float(*)[3]) calloc(W * H, sizeof * image); +#ifdef _OPENMP + #pragma omp parallel shared(image) +#endif + { + // convert red, blue, green to image +#ifdef _OPENMP + #pragma omp for +#endif + + for (int i = 0; i < H; i++) { + for (int j = 0; j < W; j++) { + image[i * W + j][0] = red [i][j]; + image[i * W + j][1] = green[i][j]; + image[i * W + j][2] = blue [i][j]; + } + } + + for (int b = 0; b < PassCount; b++) { + if (plistener) { + plistener->setProgressStr (M("TP_RAW_DMETHOD_PROGRESSBAR_REFINE")); + plistener->setProgress ((float)b / PassCount); + } + + // Reinforce interpolated green pixels on RED/BLUE pixel locations +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 6; row < H - 6; row++) { + for (int col = 6 + (FC(row, 2) & 1), c = FC(row, col); col < W - 6; col += 2) { + float (*pix)[3] = image + row * W + col; + + // Cubic Spline Interpolation by Li and Randhawa, modified by Luis Sanz Rodriguez + + float f[4]; + f[0] = 1.0f / (1.0f + xmul2f(fabs(x1125(pix[-v][c]) - x0875(pix[0][c]) - x0250(pix[-x][c]))) + fabs(x0875(pix[u][1]) - x1125(pix[-u][1]) + x0250(pix[-w][1])) + fabs(x0875(pix[-w][1]) - x1125(pix[-u][1]) + x0250(pix[-y][1]))); + f[1] = 1.0f / (1.0f + xmul2f(fabs(x1125(pix[+2][c]) - x0875(pix[0][c]) - x0250(pix[+4][c]))) + fabs(x0875(pix[1][1]) - x1125(pix[-1][1]) + x0250(pix[+3][1])) + fabs(x0875(pix[+3][1]) - x1125(pix[+1][1]) + x0250(pix[+5][1]))); + f[2] = 1.0f / (1.0f + xmul2f(fabs(x1125(pix[-2][c]) - x0875(pix[0][c]) - x0250(pix[-4][c]))) + fabs(x0875(pix[1][1]) - x1125(pix[-1][1]) + x0250(pix[-3][1])) + fabs(x0875(pix[-3][1]) - x1125(pix[-1][1]) + x0250(pix[-5][1]))); + f[3] = 1.0f / (1.0f + xmul2f(fabs(x1125(pix[+v][c]) - x0875(pix[0][c]) - x0250(pix[+x][c]))) + fabs(x0875(pix[u][1]) - x1125(pix[-u][1]) + x0250(pix[+w][1])) + fabs(x0875(pix[+w][1]) - x1125(pix[+u][1]) + x0250(pix[+y][1]))); + + float g[4];//CLIREF avoid overflow + g[0] = pix[0][c] + (x0875(CLIREF(pix[-u][1] - pix[-u][c])) + x0125(CLIREF(pix[+u][1] - pix[+u][c]))); + g[1] = pix[0][c] + (x0875(CLIREF(pix[+1][1] - pix[+1][c])) + x0125(CLIREF(pix[-1][1] - pix[-1][c]))); + g[2] = pix[0][c] + (x0875(CLIREF(pix[-1][1] - pix[-1][c])) + x0125(CLIREF(pix[+1][1] - pix[+1][c]))); + g[3] = pix[0][c] + (x0875(CLIREF(pix[+u][1] - pix[+u][c])) + x0125(CLIREF(pix[-u][1] - pix[-u][c]))); + + pix[0][1] = (f[0] * g[0] + f[1] * g[1] + f[2] * g[2] + f[3] * g[3]) / (f[0] + f[1] + f[2] + f[3]); + + } + } + + // Reinforce interpolated red/blue pixels on GREEN pixel locations +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 6; row < H - 6; row++) { + for (int col = 6 + (FC(row, 3) & 1), c = FC(row, col + 1); col < W - 6; col += 2) { + float (*pix)[3] = image + row * W + col; + + for (int i = 0; i < 2; c = 2 - c, i++) { + float f[4]; + f[0] = 1.0f / (1.0f + xmul2f(fabs(x0875(pix[-v][1]) - x1125(pix[0][1]) + x0250(pix[-x][1]))) + fabs(pix[u] [c] - pix[-u][c]) + fabs(pix[-w][c] - pix[-u][c])); + f[1] = 1.0f / (1.0f + xmul2f(fabs(x0875(pix[+2][1]) - x1125(pix[0][1]) + x0250(pix[+4][1]))) + fabs(pix[+1][c] - pix[-1][c]) + fabs(pix[+3][c] - pix[+1][c])); + f[2] = 1.0f / (1.0f + xmul2f(fabs(x0875(pix[-2][1]) - x1125(pix[0][1]) + x0250(pix[-4][1]))) + fabs(pix[+1][c] - pix[-1][c]) + fabs(pix[-3][c] - pix[-1][c])); + f[3] = 1.0f / (1.0f + xmul2f(fabs(x0875(pix[+v][1]) - x1125(pix[0][1]) + x0250(pix[+x][1]))) + fabs(pix[u] [c] - pix[-u][c]) + fabs(pix[+w][c] - pix[+u][c])); + + float g[5];//CLIREF avoid overflow + g[0] = CLIREF(pix[-u][1] - pix[-u][c]); + g[1] = CLIREF(pix[+1][1] - pix[+1][c]); + g[2] = CLIREF(pix[-1][1] - pix[-1][c]); + g[3] = CLIREF(pix[+u][1] - pix[+u][c]); + g[4] = ((f[0] * g[0] + f[1] * g[1] + f[2] * g[2] + f[3] * g[3]) / (f[0] + f[1] + f[2] + f[3])); + pix[0][c] = pix[0][1] - (0.65f * g[4] + 0.35f * CLIREF(pix[0][1] - pix[0][c])); + } + } + } + + // Reinforce integrated red/blue pixels on BLUE/RED pixel locations +#ifdef _OPENMP + #pragma omp for +#endif + + for (int row = 6; row < H - 6; row++) { + for (int col = 6 + (FC(row, 2) & 1), c = 2 - FC(row, col), d = 2 - c; col < W - 6; col += 2) { + float (*pix)[3] = image + row * W + col; + + float f[4]; + f[0] = 1.0f / (1.0f + xmul2f(fabs(x1125(pix[-v][d]) - x0875(pix[0][d]) - x0250(pix[-x][d]))) + fabs(x0875(pix[u][1]) - x1125(pix[-u][1]) + x0250(pix[-w][1])) + fabs(x0875(pix[-w][1]) - x1125(pix[-u][1]) + x0250(pix[-y][1]))); + f[1] = 1.0f / (1.0f + xmul2f(fabs(x1125(pix[+2][d]) - x0875(pix[0][d]) - x0250(pix[+4][d]))) + fabs(x0875(pix[1][1]) - x1125(pix[-1][1]) + x0250(pix[+3][1])) + fabs(x0875(pix[+3][1]) - x1125(pix[+1][1]) + x0250(pix[+5][1]))); + f[2] = 1.0f / (1.0f + xmul2f(fabs(x1125(pix[-2][d]) - x0875(pix[0][d]) - x0250(pix[-4][d]))) + fabs(x0875(pix[1][1]) - x1125(pix[-1][1]) + x0250(pix[-3][1])) + fabs(x0875(pix[-3][1]) - x1125(pix[-1][1]) + x0250(pix[-5][1]))); + f[3] = 1.0f / (1.0f + xmul2f(fabs(x1125(pix[+v][d]) - x0875(pix[0][d]) - x0250(pix[+x][d]))) + fabs(x0875(pix[u][1]) - x1125(pix[-u][1]) + x0250(pix[+w][1])) + fabs(x0875(pix[+w][1]) - x1125(pix[+u][1]) + x0250(pix[+y][1]))); + + float g[5]; + g[0] = (x0875((pix[-u][1] - pix[-u][c])) + x0125((pix[-v][1] - pix[-v][c]))); + g[1] = (x0875((pix[+1][1] - pix[+1][c])) + x0125((pix[+2][1] - pix[+2][c]))); + g[2] = (x0875((pix[-1][1] - pix[-1][c])) + x0125((pix[-2][1] - pix[-2][c]))); + g[3] = (x0875((pix[+u][1] - pix[+u][c])) + x0125((pix[+v][1] - pix[+v][c]))); + + g[4] = (f[0] * g[0] + f[1] * g[1] + f[2] * g[2] + f[3] * g[3]) / (f[0] + f[1] + f[2] + f[3]); + + const std::array p = { + pix[-u - 1][1] - pix[-u - 1][c], + pix[-u + 0][1] - pix[-u + 0][c], + pix[-u + 1][1] - pix[-u + 1][c], + pix[+0 - 1][1] - pix[+0 - 1][c], + pix[+0 + 0][1] - pix[+0 + 0][c], + pix[+0 + 1][1] - pix[+0 + 1][c], + pix[+u - 1][1] - pix[+u - 1][c], + pix[+u + 0][1] - pix[+u + 0][c], + pix[+u + 1][1] - pix[+u + 1][c] + }; + + const float med = median(p); + + pix[0][c] = LIM(pix[0][1] - (1.30f * g[4] - 0.30f * (pix[0][1] - pix[0][c])), 0.99f * (pix[0][1] - med), 1.01f * (pix[0][1] - med)); + + } + } + + } + + // put modified values to red, green, blue +#ifdef _OPENMP + #pragma omp for +#endif + + for (int i = 0; i < H; i++) { + for (int j = 0; j < W; j++) { + red [i][j] = image[i * W + j][0]; + green[i][j] = image[i * W + j][1]; + blue [i][j] = image[i * W + j][2]; + } + } + } + + free(image); + + t2e.set(); + + if (settings->verbose) { + printf("Refinement Lassus %d usec\n", t2e.etime(t1e)); + } +} + + +/* + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * * Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above + * copyright notice, this list of conditions and the following disclaimer + * in the documentation and/or other materials provided with the + * distribution. + * * Neither the name of the author nor the names of its + * contributors may be used to endorse or promote products derived from + * this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +// If you want to use the code, you need to display name of the original authors in +// your software! + +/* DCB demosaicing by Jacek Gozdz (cuniek@kft.umcs.lublin.pl) + * the code is open source (BSD licence) +*/ + +#define TILESIZE 192 +#define TILEBORDER 10 +#define CACHESIZE (TILESIZE+2*TILEBORDER) + +inline void RawImageSource::dcb_initTileLimits(int &colMin, int &rowMin, int &colMax, int &rowMax, int x0, int y0, int border) +{ + rowMin = border; + colMin = border; + rowMax = CACHESIZE - border; + colMax = CACHESIZE - border; + + if(!y0 ) { + rowMin = TILEBORDER + border; + } + + if(!x0 ) { + colMin = TILEBORDER + border; + } + + if( y0 + TILESIZE + TILEBORDER >= H - border) { + rowMax = TILEBORDER + H - border - y0; + } + + if( x0 + TILESIZE + TILEBORDER >= W - border) { + colMax = TILEBORDER + W - border - x0; + } +} + +void RawImageSource::fill_raw( float (*cache )[3], int x0, int y0, float** rawData) +{ + int rowMin, colMin, rowMax, colMax; + dcb_initTileLimits(colMin, rowMin, colMax, rowMax, x0, y0, 0); + + for (int row = rowMin, y = y0 - TILEBORDER + rowMin; row < rowMax; row++, y++) + for (int col = colMin, x = x0 - TILEBORDER + colMin, indx = row * CACHESIZE + col; col < colMax; col++, x++, indx++) { + cache[indx][fc(y, x)] = rawData[y][x]; + } +} + +void RawImageSource::fill_border( float (*cache )[3], int border, int x0, int y0) +{ + unsigned f; + float sum[8]; + constexpr unsigned int colors = 3; // used in FORCC + + for (int row = y0; row < y0 + TILESIZE + TILEBORDER && row < H; row++) { + for (int col = x0; col < x0 + TILESIZE + TILEBORDER && col < W; col++) { + if (col >= border && col < W - border && row >= border && row < H - border) { + col = W - border; + + if(col >= x0 + TILESIZE + TILEBORDER ) { + break; + } + } + + memset(sum, 0, sizeof sum); + + for (int y = row - 1; y != row + 2; y++) + for (int x = col - 1; x != col + 2; x++) + if (y < H && y < y0 + TILESIZE + TILEBORDER && x < W && x < x0 + TILESIZE + TILEBORDER) { + f = fc(y, x); + sum[f] += cache[(y - y0 + TILEBORDER) * CACHESIZE + TILEBORDER + x - x0][f]; + sum[f + 4]++; + } + + f = fc(row, col); + FORCC + + if (c != f && sum[c + 4] > 0) { + cache[(row - y0 + TILEBORDER) * CACHESIZE + TILEBORDER + col - x0][c] = sum[c] / sum[c + 4]; + } + } + } +} + +// saves red and blue + +// change buffer[3] -> buffer[2], possibly to buffer[1] if split +// into two loops, one for R and another for B, could also be smaller because +// there is no need for green pixels pass +// this would decrease the amount of needed memory +// from megapixels*2 records to megapixels*0.5 +// also don't know if float is needed as data is 1-65536 integer (I believe!!) +// comment from Ingo: float is needed because rawdata in rt is float +void RawImageSource::copy_to_buffer( float (*buffer)[2], float (*image)[3]) +{ + for (int indx = 0; indx < CACHESIZE * CACHESIZE; indx++) { + buffer[indx][0] = image[indx][0]; //R + buffer[indx][1] = image[indx][2]; //B + } +} + +// restores red and blue + +// other comments like in copy_to_buffer +void RawImageSource::restore_from_buffer(float (*image)[3], float (*buffer)[2]) +{ + for (int indx = 0; indx < CACHESIZE * CACHESIZE; indx++) { + image[indx][0] = buffer[indx][0]; //R + image[indx][2] = buffer[indx][1]; //B + } +} + +// First pass green interpolation + +// remove entirely: bufferH and bufferV +void RawImageSource::dcb_hid(float (*image)[3], int x0, int y0) +{ + const int u = CACHESIZE; + int rowMin, colMin, rowMax, colMax; + dcb_initTileLimits(colMin, rowMin, colMax, rowMax, x0, y0, 2); + +// simple green bilinear in R and B pixels + for (int row = rowMin; row < rowMax; row++) + for (int col = colMin + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin) & 1), indx = row * CACHESIZE + col; col < colMax; col += 2, indx += 2) { + assert(indx - u - 1 >= 0 && indx + u + 1 < u * u); + + image[indx][1] = 0.25*(image[indx-1][1]+image[indx+1][1]+image[indx-u][1]+image[indx+u][1]); + } +} + +// missing colours are interpolated +void RawImageSource::dcb_color(float (*image)[3], int x0, int y0) +{ + const int u = CACHESIZE; + int rowMin, colMin, rowMax, colMax; + dcb_initTileLimits(colMin, rowMin, colMax, rowMax, x0, y0, 1); + + // red in blue pixel, blue in red pixel + for (int row = rowMin; row < rowMax; row++) + for (int col = colMin + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin) & 1), indx = row * CACHESIZE + col, c = 2 - FC(y0 - TILEBORDER + row, x0 - TILEBORDER + col); col < colMax; col += 2, indx += 2) { + assert(indx >= 0 && indx < u * u && c >= 0 && c < 4); + + +//Jacek comment: one multiplication less + image[indx][c] = image[indx][1] + + ( image[indx + u + 1][c] + image[indx + u - 1][c] + image[indx - u + 1][c] + image[indx - u - 1][c] + - (image[indx + u + 1][1] + image[indx + u - 1][1] + image[indx - u + 1][1] + image[indx - u - 1][1]) ) * 0.25f; + +/* original + image[indx][c] = ( 4.f * image[indx][1] + - image[indx + u + 1][1] - image[indx + u - 1][1] - image[indx - u + 1][1] - image[indx - u - 1][1] + + image[indx + u + 1][c] + image[indx + u - 1][c] + image[indx - u + 1][c] + image[indx - u - 1][c] ) * 0.25f; +*/ + } + + // red or blue in green pixels + for (int row = rowMin; row < rowMax; row++) + for (int col = colMin + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin + 1) & 1), indx = row * CACHESIZE + col, c = FC(y0 - TILEBORDER + row, x0 - TILEBORDER + col + 1), d = 2 - c; col < colMax; col += 2, indx += 2) { + assert(indx >= 0 && indx < u * u && c >= 0 && c < 4); + +//Jacek comment: two multiplications (in total) less + image[indx][c] = image[indx][1] + (image[indx + 1][c] + image[indx - 1][c] - (image[indx + 1][1] + image[indx - 1][1])) * 0.5f; + image[indx][d] = image[indx][1] + (image[indx + u][d] + image[indx - u][d] - (image[indx + u][1] + image[indx - u][1])) * 0.5f; + + +/* original + image[indx][c] = (2.f * image[indx][1] - image[indx + 1][1] - image[indx - 1][1] + image[indx + 1][c] + image[indx - 1][c]) * 0.5f; + image[indx][d] = (2.f * image[indx][1] - image[indx + u][1] - image[indx - u][1] + image[indx + u][d] + image[indx - u][d]) * 0.5f; +*/ + } +} + +// green correction +void RawImageSource::dcb_hid2(float (*image)[3], int x0, int y0) +{ + const int v = 2 * CACHESIZE; + int rowMin, colMin, rowMax, colMax; + dcb_initTileLimits(colMin, rowMin, colMax, rowMax, x0, y0, 2); + + for (int row = rowMin; row < rowMax; row++) { + for (int col = colMin + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin) & 1), indx = row * CACHESIZE + col, c = FC(y0 - TILEBORDER + row, x0 - TILEBORDER + col); col < colMax; col += 2, indx += 2) { + assert(indx - v >= 0 && indx + v < CACHESIZE * CACHESIZE); + +//Jacek comment: one multiplication less + image[indx][1] = image[indx][c] + + (image[indx + v][1] + image[indx - v][1] + image[indx - 2][1] + image[indx + 2][1] + - (image[indx + v][c] + image[indx - v][c] + image[indx - 2][c] + image[indx + 2][c])) * 0.25f; + +/* original + image[indx][1] = (image[indx + v][1] + image[indx - v][1] + image[indx - 2][1] + image[indx + 2][1]) * 0.25f + + image[indx][c] - ( image[indx + v][c] + image[indx - v][c] + image[indx - 2][c] + image[indx + 2][c]) * 0.25f; +*/ + } + } +} + +// green is used to create +// an interpolation direction map +// 1 = vertical +// 0 = horizontal +// saved in image[][3] + +// seems at least 2 persons implemented some code, as this one has different coding style, could be unified +// I don't know if *pix is faster than a loop working on image[] directly +void RawImageSource::dcb_map(float (*image)[3], uint8_t *map, int x0, int y0) +{ + const int u = 3 * CACHESIZE; + int rowMin, colMin, rowMax, colMax; + dcb_initTileLimits(colMin, rowMin, colMax, rowMax, x0, y0, 2); + + for (int row = rowMin; row < rowMax; row++) { + for (int col = colMin, indx = row * CACHESIZE + col; col < colMax; col++, indx++) { + float *pix = &(image[indx][1]); + + assert(indx >= 0 && indx < u * u); + + // comparing 4 * a to (b+c+d+e) instead of a to (b+c+d+e)/4 is faster because divisions are slow + if ( 4 * (*pix) > ( (pix[-3] + pix[+3]) + (pix[-u] + pix[+u])) ) { + map[indx] = ((min(pix[-3], pix[+3]) + (pix[-3] + pix[+3]) ) < (min(pix[-u], pix[+u]) + (pix[-u] + pix[+u]))); + } else { + map[indx] = ((max(pix[-3], pix[+3]) + (pix[-3] + pix[+3]) ) > (max(pix[-u], pix[+u]) + (pix[-u] + pix[+u]))); + } + } + } +} + +// interpolated green pixels are corrected using the map +void RawImageSource::dcb_correction(float (*image)[3], uint8_t *map, int x0, int y0) +{ + const int u = CACHESIZE, v = 2 * CACHESIZE; + int rowMin, colMin, rowMax, colMax; + dcb_initTileLimits(colMin, rowMin, colMax, rowMax, x0, y0, 2); + + for (int row = rowMin; row < rowMax; row++) { + for (int indx = row * CACHESIZE + colMin + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin) & 1); indx < row * CACHESIZE + colMax; indx += 2) { +// for (int col = colMin + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin) & 1), indx = row * CACHESIZE + col; col < colMax; col += 2, indx += 2) { + float current = 4 * map[indx] + + 2 * (map[indx + u] + map[indx - u] + map[indx + 1] + map[indx - 1]) + + map[indx + v] + map[indx - v] + map[indx + 2] + map[indx - 2]; + + assert(indx >= 0 && indx < u * u); + image[indx][1] = ((16.f - current) * (image[indx - 1][1] + image[indx + 1][1]) + current * (image[indx - u][1] + image[indx + u][1]) ) * 0.03125f; +// image[indx][1] = ((16.f - current) * (image[indx - 1][1] + image[indx + 1][1]) * 0.5f + current * (image[indx - u][1] + image[indx + u][1]) * 0.5f ) * 0.0625f; + } + } +} + +// R and B smoothing using green contrast, all pixels except 2 pixel wide border + +// again code with *pix, is this kind of calculating faster in C, than this what was commented? +void RawImageSource::dcb_pp(float (*image)[3], int x0, int y0) +{ + const int u = CACHESIZE; + int rowMin, colMin, rowMax, colMax; + dcb_initTileLimits(colMin, rowMin, colMax, rowMax, x0, y0, 2); + + for (int row = rowMin; row < rowMax; row++) + for (int col = colMin, indx = row * CACHESIZE + col; col < colMax; col++, indx++) { +// float r1 = image[indx-1][0] + image[indx+1][0] + image[indx-u][0] + image[indx+u][0] + image[indx-u-1][0] + image[indx+u+1][0] + image[indx-u+1][0] + image[indx+u-1][0]; +// float g1 = image[indx-1][1] + image[indx+1][1] + image[indx-u][1] + image[indx+u][1] + image[indx-u-1][1] + image[indx+u+1][1] + image[indx-u+1][1] + image[indx+u-1][1]; +// float b1 = image[indx-1][2] + image[indx+1][2] + image[indx-u][2] + image[indx+u][2] + image[indx-u-1][2] + image[indx+u+1][2] + image[indx-u+1][2] + image[indx+u-1][2]; + float (*pix)[3] = image + (indx - u - 1); + float r1 = (*pix)[0]; + float g1 = (*pix)[1]; + float b1 = (*pix)[2]; + pix++; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix++; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix += CACHESIZE - 2; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix += 2; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix += CACHESIZE - 2; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix++; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix++; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + r1 *= 0.125f; + g1 *= 0.125f; + b1 *= 0.125f; + r1 += ( image[indx][1] - g1 ); + b1 += ( image[indx][1] - g1 ); + + assert(indx >= 0 && indx < u * u); + image[indx][0] = r1; + image[indx][2] = b1; + } +} + +// interpolated green pixels are corrected using the map +// with correction +void RawImageSource::dcb_correction2(float (*image)[3], uint8_t *map, int x0, int y0) +{ + const int u = CACHESIZE, v = 2 * CACHESIZE; + int rowMin, colMin, rowMax, colMax; + dcb_initTileLimits(colMin, rowMin, colMax, rowMax, x0, y0, 4); + + for (int row = rowMin; row < rowMax; row++) { + for (int indx = row * CACHESIZE + colMin + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin) & 1), c = FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin) & 1)); indx < row * CACHESIZE + colMax; indx += 2) { + // map values are uint8_t either 0 or 1. Adding them using integer instructions is perfectly valid and fast. Final result is converted to float then + float current = 4 * map[indx] + + 2 * (map[indx + u] + map[indx - u] + map[indx + 1] + map[indx - 1]) + + map[indx + v] + map[indx - v] + map[indx + 2] + map[indx - 2]; + + assert(indx >= 0 && indx < u * u); + +// Jacek comment: works now, and has 3 float mults and 9 float adds + image[indx][1] = image[indx][c] + + ((16.f - current) * (image[indx - 1][1] + image[indx + 1][1] - (image[indx + 2][c] + image[indx - 2][c])) + + current * (image[indx - u][1] + image[indx + u][1] - (image[indx + v][c] + image[indx - v][c]))) * 0.03125f; + + + // 4 float mults and 9 float adds + // Jacek comment: not mathematically identical to original +/* image[indx][1] = 16.f * image[indx][c] + + ((16.f - current) * ((image[indx - 1][1] + image[indx + 1][1]) + - (image[indx + 2][c] + image[indx - 2][c])) + + current * ((image[indx - u][1] + image[indx + u][1]) - (image[indx + v][c] + image[indx - v][c]))) * 0.03125f; +*/ + // 7 float mults and 10 float adds + // original code +/* + image[indx][1] = ((16.f - current) * ((image[indx - 1][1] + image[indx + 1][1]) * 0.5f + + image[indx][c] - (image[indx + 2][c] + image[indx - 2][c]) * 0.5f) + + current * ((image[indx - u][1] + image[indx + u][1]) * 0.5f + image[indx][c] - (image[indx + v][c] + image[indx - v][c]) * 0.5f)) * 0.0625f; +*/ + } + } +} + +// image refinement +void RawImageSource::dcb_refinement(float (*image)[3], uint8_t *map, int x0, int y0) +{ + const int u = CACHESIZE, v = 2 * CACHESIZE; + int rowMin, colMin, rowMax, colMax; + dcb_initTileLimits(colMin, rowMin, colMax, rowMax, x0, y0, 4); + + float f0, f1, f2, g1, h0, h1, h2, g2; + + for (int row = rowMin; row < rowMax; row++) + for (int col = colMin + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin) & 1), indx = row * CACHESIZE + col, c = FC(y0 - TILEBORDER + row, x0 - TILEBORDER + col); col < colMax; col += 2, indx += 2) { + + float current = 4 * map[indx] + + 2 * (map[indx + u] + map[indx - u] + map[indx + 1] + map[indx - 1]) + + map[indx + v] + map[indx - v] + map[indx - 2] + map[indx + 2]; + + float currPix = image[indx][c]; + + f0 = (float)(image[indx - u][1] + image[indx + u][1]) / (1.f + 2.f * currPix); + f1 = 2.f * image[indx - u][1] / (1.f + image[indx - v][c] + currPix); + f2 = 2.f * image[indx + u][1] / (1.f + image[indx + v][c] + currPix); + + g1 = f0 + f1 + f2; + + h0 = (float)(image[indx - 1][1] + image[indx + 1][1]) / (1.f + 2.f * currPix); + h1 = 2.f * image[indx - 1][1] / (1.f + image[indx - 2][c] + currPix); + h2 = 2.f * image[indx + 1][1] / (1.f + image[indx + 2][c] + currPix); + + g2 = h0 + h1 + h2; + + // new green value + assert(indx >= 0 && indx < u * u); + currPix *= (current * g1 + (16.f - current) * g2) / 48.f; + + // get rid of the overshot pixels + float minVal = min(image[indx - 1][1], min(image[indx + 1][1], min(image[indx - u][1], image[indx + u][1]))); + float maxVal = max(image[indx - 1][1], max(image[indx + 1][1], max(image[indx - u][1], image[indx + u][1]))); + + image[indx][1] = LIM(currPix, minVal, maxVal); + + } +} + +// missing colours are interpolated using high quality algorithm by Luis Sanz Rodriguez +void RawImageSource::dcb_color_full(float (*image)[3], int x0, int y0, float (*chroma)[2]) +{ + const int u = CACHESIZE, w = 3 * CACHESIZE; + int rowMin, colMin, rowMax, colMax; + dcb_initTileLimits(colMin, rowMin, colMax, rowMax, x0, y0, 3); + + float f[4], g[4]; + + for (int row = 1; row < CACHESIZE - 1; row++) + for (int col = 1 + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + 1) & 1), indx = row * CACHESIZE + col, c = FC(y0 - TILEBORDER + row, x0 - TILEBORDER + col), d = c / 2; col < CACHESIZE - 1; col += 2, indx += 2) { + assert(indx >= 0 && indx < u * u && c >= 0 && c < 4); + chroma[indx][d] = image[indx][c] - image[indx][1]; + } + + for (int row = rowMin; row < rowMax; row++) + for (int col = colMin + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin) & 1), indx = row * CACHESIZE + col, c = 1 - FC(y0 - TILEBORDER + row, x0 - TILEBORDER + col) / 2; col < colMax; col += 2, indx += 2) { + f[0] = 1.f / (float)(1.f + fabs(chroma[indx - u - 1][c] - chroma[indx + u + 1][c]) + fabs(chroma[indx - u - 1][c] - chroma[indx - w - 3][c]) + fabs(chroma[indx + u + 1][c] - chroma[indx - w - 3][c])); + f[1] = 1.f / (float)(1.f + fabs(chroma[indx - u + 1][c] - chroma[indx + u - 1][c]) + fabs(chroma[indx - u + 1][c] - chroma[indx - w + 3][c]) + fabs(chroma[indx + u - 1][c] - chroma[indx - w + 3][c])); + f[2] = 1.f / (float)(1.f + fabs(chroma[indx + u - 1][c] - chroma[indx - u + 1][c]) + fabs(chroma[indx + u - 1][c] - chroma[indx + w + 3][c]) + fabs(chroma[indx - u + 1][c] - chroma[indx + w - 3][c])); + f[3] = 1.f / (float)(1.f + fabs(chroma[indx + u + 1][c] - chroma[indx - u - 1][c]) + fabs(chroma[indx + u + 1][c] - chroma[indx + w - 3][c]) + fabs(chroma[indx - u - 1][c] - chroma[indx + w + 3][c])); + g[0] = 1.325f * chroma[indx - u - 1][c] - 0.175f * chroma[indx - w - 3][c] - 0.075f * (chroma[indx - w - 1][c] + chroma[indx - u - 3][c]); + g[1] = 1.325f * chroma[indx - u + 1][c] - 0.175f * chroma[indx - w + 3][c] - 0.075f * (chroma[indx - w + 1][c] + chroma[indx - u + 3][c]); + g[2] = 1.325f * chroma[indx + u - 1][c] - 0.175f * chroma[indx + w - 3][c] - 0.075f * (chroma[indx + w - 1][c] + chroma[indx + u - 3][c]); + g[3] = 1.325f * chroma[indx + u + 1][c] - 0.175f * chroma[indx + w + 3][c] - 0.075f * (chroma[indx + w + 1][c] + chroma[indx + u + 3][c]); + +// g[0] = 1.325f * chroma[indx - u - 1][c] - 0.175f * chroma[indx - w - 3][c] - 0.075f * chroma[indx - w - 1][c] - 0.075f * chroma[indx - u - 3][c]; +// g[1] = 1.325f * chroma[indx - u + 1][c] - 0.175f * chroma[indx - w + 3][c] - 0.075f * chroma[indx - w + 1][c] - 0.075f * chroma[indx - u + 3][c]; +// g[2] = 1.325f * chroma[indx + u - 1][c] - 0.175f * chroma[indx + w - 3][c] - 0.075f * chroma[indx + w - 1][c] - 0.075f * chroma[indx + u - 3][c]; +// g[3] = 1.325f * chroma[indx + u + 1][c] - 0.175f * chroma[indx + w + 3][c] - 0.075f * chroma[indx + w + 1][c] - 0.075f * chroma[indx + u + 3][c]; + + assert(indx >= 0 && indx < u * u && c >= 0 && c < 2); + chroma[indx][c] = (f[0] * g[0] + f[1] * g[1] + f[2] * g[2] + f[3] * g[3]) / (f[0] + f[1] + f[2] + f[3]); + } + + for (int row = rowMin; row < rowMax; row++) + for (int col = colMin + (FC(y0 - TILEBORDER + row, x0 - TILEBORDER + colMin + 1) & 1), indx = row * CACHESIZE + col, c = FC(y0 - TILEBORDER + row, x0 - TILEBORDER + col + 1) / 2; col < colMax; col += 2, indx += 2) + for(int d = 0; d <= 1; c = 1 - c, d++) { + f[0] = 1.f / (1.f + fabs(chroma[indx - u][c] - chroma[indx + u][c]) + fabs(chroma[indx - u][c] - chroma[indx - w][c]) + fabs(chroma[indx + u][c] - chroma[indx - w][c])); + f[1] = 1.f / (1.f + fabs(chroma[indx + 1][c] - chroma[indx - 1][c]) + fabs(chroma[indx + 1][c] - chroma[indx + 3][c]) + fabs(chroma[indx - 1][c] - chroma[indx + 3][c])); + f[2] = 1.f / (1.f + fabs(chroma[indx - 1][c] - chroma[indx + 1][c]) + fabs(chroma[indx - 1][c] - chroma[indx - 3][c]) + fabs(chroma[indx + 1][c] - chroma[indx - 3][c])); + f[3] = 1.f / (1.f + fabs(chroma[indx + u][c] - chroma[indx - u][c]) + fabs(chroma[indx + u][c] - chroma[indx + w][c]) + fabs(chroma[indx - u][c] - chroma[indx + w][c])); + + g[0] = intp(0.875f, chroma[indx - u][c], chroma[indx - w][c]); + g[1] = intp(0.875f, chroma[indx + 1][c], chroma[indx + 3][c]); + g[2] = intp(0.875f, chroma[indx - 1][c], chroma[indx - 3][c]); + g[3] = intp(0.875f, chroma[indx + u][c], chroma[indx + w][c]); + +// g[0] = 0.875f * chroma[indx - u][c] + 0.125f * chroma[indx - w][c]; +// g[1] = 0.875f * chroma[indx + 1][c] + 0.125f * chroma[indx + 3][c]; +// g[2] = 0.875f * chroma[indx - 1][c] + 0.125f * chroma[indx - 3][c]; +// g[3] = 0.875f * chroma[indx + u][c] + 0.125f * chroma[indx + w][c]; + + assert(indx >= 0 && indx < u * u && c >= 0 && c < 2); + chroma[indx][c] = (f[0] * g[0] + f[1] * g[1] + f[2] * g[2] + f[3] * g[3]) / (f[0] + f[1] + f[2] + f[3]); + } + + for(int row = rowMin; row < rowMax; row++) + for(int col = colMin, indx = row * CACHESIZE + col; col < colMax; col++, indx++) { + assert(indx >= 0 && indx < u * u); + + image[indx][0] = chroma[indx][0] + image[indx][1]; + image[indx][2] = chroma[indx][1] + image[indx][1]; + } +} + +// DCB demosaicing main routine +void RawImageSource::dcb_demosaic(int iterations, bool dcb_enhance) +{ +BENCHFUN + double currentProgress = 0.0; + + if(plistener) { + plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::DCB))); + plistener->setProgress (currentProgress); + } + + int wTiles = W / TILESIZE + (W % TILESIZE ? 1 : 0); + int hTiles = H / TILESIZE + (H % TILESIZE ? 1 : 0); + int numTiles = wTiles * hTiles; + int tilesDone = 0; + constexpr int cldf = 2; // factor to multiply cache line distance. 1 = 64 bytes, 2 = 128 bytes ... + +#ifdef _OPENMP + #pragma omp parallel +#endif +{ + // assign working space + char *buffer0 = (char *) malloc(5 * sizeof(float) * CACHESIZE * CACHESIZE + sizeof(uint8_t) * CACHESIZE * CACHESIZE + 3 * cldf * 64 + 63); + // aligned to 64 byte boundary + char *data = (char*)( ( uintptr_t(buffer0) + uintptr_t(63)) / 64 * 64); + + float (*tile)[3] = (float(*)[3]) data; + float (*buffer)[2] = (float(*)[2]) ((char*)tile + sizeof(float) * CACHESIZE * CACHESIZE * 3 + cldf * 64); + float (*chrm)[2] = (float(*)[2]) (buffer); // No overlap in usage of buffer and chrm means we can reuse buffer + uint8_t *map = (uint8_t*) ((char*)buffer + sizeof(float) * CACHESIZE * CACHESIZE * 2 + cldf * 64); + +#ifdef _OPENMP + #pragma omp for schedule(dynamic) nowait +#endif + + for( int iTile = 0; iTile < numTiles; iTile++) { + int xTile = iTile % wTiles; + int yTile = iTile / wTiles; + int x0 = xTile * TILESIZE; + int y0 = yTile * TILESIZE; + + memset(tile, 0, CACHESIZE * CACHESIZE * sizeof * tile); + memset(map, 0, CACHESIZE * CACHESIZE * sizeof * map); + + fill_raw( tile, x0, y0, rawData ); + + if( !xTile || !yTile || xTile == wTiles - 1 || yTile == hTiles - 1) { + fill_border(tile, 6, x0, y0); + } + + copy_to_buffer(buffer, tile); + dcb_hid(tile, x0, y0); + + for (int i = iterations; i > 0; i--) { + dcb_hid2(tile, x0, y0); + dcb_hid2(tile, x0, y0); + dcb_hid2(tile, x0, y0); + dcb_map(tile, map, x0, y0); + dcb_correction(tile, map, x0, y0); + } + + dcb_color(tile, x0, y0); + dcb_pp(tile, x0, y0); + dcb_map(tile, map, x0, y0); + dcb_correction2(tile, map, x0, y0); + dcb_map(tile, map, x0, y0); + dcb_correction(tile, map, x0, y0); + dcb_color(tile, x0, y0); + dcb_map(tile, map, x0, y0); + dcb_correction(tile, map, x0, y0); + dcb_map(tile, map, x0, y0); + dcb_correction(tile, map, x0, y0); + dcb_map(tile, map, x0, y0); + restore_from_buffer(tile, buffer); + + if (!dcb_enhance) + dcb_color(tile, x0, y0); + else + { + memset(chrm, 0, CACHESIZE * CACHESIZE * sizeof * chrm); + dcb_refinement(tile, map, x0, y0); + dcb_color_full(tile, x0, y0, chrm); + } + + /* + dcb_hid(tile, buffer, buffer2, x0, y0); + dcb_color(tile, x0, y0); + + copy_to_buffer(buffer, tile); + + for (int i = iterations; i > 0; i--) { + dcb_hid2(tile, x0, y0); + dcb_hid2(tile, x0, y0); + dcb_hid2(tile, x0, y0); + dcb_map(tile, x0, y0); + dcb_correction(tile, x0, y0); + } + + dcb_color(tile, x0, y0); + dcb_pp(tile, x0, y0); + dcb_map(tile, x0, y0); + dcb_correction2(tile, x0, y0); + dcb_map(tile, x0, y0); + dcb_correction(tile, x0, y0); + dcb_color(tile, x0, y0); + dcb_map(tile, x0, y0); + dcb_correction(tile, x0, y0); + dcb_map(tile, x0, y0); + dcb_correction(tile, x0, y0); + dcb_map(tile, x0, y0); + restore_from_buffer(tile, buffer); + dcb_color_full(tile, x0, y0, chrm); + + if (dcb_enhance) { + dcb_refinement(tile, x0, y0); + dcb_color_full(tile, x0, y0, chrm); + } +*/ + for(int y = 0; y < TILESIZE && y0 + y < H; y++) { + for (int j = 0; j < TILESIZE && x0 + j < W; j++) { + red[y0 + y][x0 + j] = tile[(y + TILEBORDER) * CACHESIZE + TILEBORDER + j][0]; + green[y0 + y][x0 + j] = tile[(y + TILEBORDER) * CACHESIZE + TILEBORDER + j][1]; + blue[y0 + y][x0 + j] = tile[(y + TILEBORDER) * CACHESIZE + TILEBORDER + j][2]; + } + } + +#ifdef _OPENMP + + if(omp_get_thread_num() == 0) +#endif + { + if( plistener && double(tilesDone) / numTiles > currentProgress) { + currentProgress += 0.1; // Show progress each 10% + plistener->setProgress (currentProgress); + } + } + +#ifdef _OPENMP + #pragma omp atomic +#endif + tilesDone++; + } + free(buffer0); +} + + if(plistener) { + plistener->setProgress (1.0); + } +} + +#undef TILEBORDER +#undef TILESIZE +#undef CACHESIZE +} /* namespace */ diff --git a/rtengine/pixelshift.cc b/rtengine/pixelshift.cc index 488538e8a..248336407 100644 --- a/rtengine/pixelshift.cc +++ b/rtengine/pixelshift.cc @@ -323,19 +323,22 @@ BENCHFUN if(motionDetection) { if(!showOnlyMask) { if(bayerParams.pixelShiftMedian) { // We need the demosaiced frames for motion correction - if(bayerParams.pixelShiftLmmse) { + if (bayerParams.pixelShiftDemosaicMethod == bayerParams.getPSDemosaicMethodString(bayerParams.PSDemosaicMethod::LMMSE)) { lmmse_interpolate_omp(winw, winh, *(rawDataFrames[0]), red, green, blue, bayerParams.lmmse_iterations); + } else if (bayerParams.pixelShiftDemosaicMethod == bayerParams.getPSDemosaicMethodString(bayerParams.PSDemosaicMethod::AMAZEVNG4)) { + amaze_vng4_demosaic_RT (winw, winh, *(rawDataFrames[0]), red, green, blue, bayerParams.dualDemosaicContrast); } else { amaze_demosaic_RT(winx, winy, winw, winh, *(rawDataFrames[0]), red, green, blue); } - multi_array2D redTmp(winw, winh); multi_array2D greenTmp(winw, winh); multi_array2D blueTmp(winw, winh); for(int i = 0; i < 3; i++) { - if(bayerParams.pixelShiftLmmse) { + if (bayerParams.pixelShiftDemosaicMethod == bayerParams.getPSDemosaicMethodString(bayerParams.PSDemosaicMethod::LMMSE)) { lmmse_interpolate_omp(winw, winh, *(rawDataFrames[i + 1]), redTmp[i], greenTmp[i], blueTmp[i], bayerParams.lmmse_iterations); + } else if (bayerParams.pixelShiftDemosaicMethod == bayerParams.getPSDemosaicMethodString(bayerParams.PSDemosaicMethod::AMAZEVNG4)) { + amaze_vng4_demosaic_RT (winw, winh, *(rawDataFrames[i + 1]), redTmp[i], greenTmp[i], blueTmp[i], bayerParams.dualDemosaicContrast); } else { amaze_demosaic_RT(winx, winy, winw, winh, *(rawDataFrames[i + 1]), redTmp[i], greenTmp[i], blueTmp[i]); } @@ -359,8 +362,10 @@ BENCHFUN } } } else { - if(bayerParams.pixelShiftLmmse) { + if (bayerParams.pixelShiftDemosaicMethod == bayerParams.getPSDemosaicMethodString(bayerParams.PSDemosaicMethod::LMMSE)) { lmmse_interpolate_omp(winw, winh, rawData, red, green, blue, bayerParams.lmmse_iterations); + } else if (bayerParams.pixelShiftDemosaicMethod == bayerParams.getPSDemosaicMethodString(bayerParams.PSDemosaicMethod::AMAZEVNG4)) { + amaze_vng4_demosaic_RT (winw, winh, rawData, red, green, blue, bayerParams.dualDemosaicContrast); } else { amaze_demosaic_RT(winx, winy, winw, winh, rawData, red, green, blue); } diff --git a/rtengine/procparams.cc b/rtengine/procparams.cc index 3a244f5e8..7fc004aa0 100644 --- a/rtengine/procparams.cc +++ b/rtengine/procparams.cc @@ -2367,10 +2367,10 @@ RAWParams::BayerSensor::BayerSensor() : pixelShiftGreen(true), pixelShiftBlur(true), pixelShiftSmoothFactor(0.7), - pixelShiftLmmse(false), pixelShiftEqualBright(false), pixelShiftEqualBrightChannel(false), pixelShiftNonGreenCross(true), + pixelShiftDemosaicMethod(getPSDemosaicMethodString(PSDemosaicMethod::AMAZE)), dcb_enhance(true), pdafLinesFilter(false) { @@ -2403,10 +2403,10 @@ bool RAWParams::BayerSensor::operator ==(const BayerSensor& other) const && pixelShiftGreen == other.pixelShiftGreen && pixelShiftBlur == other.pixelShiftBlur && pixelShiftSmoothFactor == other.pixelShiftSmoothFactor - && pixelShiftLmmse == other.pixelShiftLmmse && pixelShiftEqualBright == other.pixelShiftEqualBright && pixelShiftEqualBrightChannel == other.pixelShiftEqualBrightChannel && pixelShiftNonGreenCross == other.pixelShiftNonGreenCross + && pixelShiftDemosaicMethod == other.pixelShiftDemosaicMethod && dcb_enhance == other.dcb_enhance && pdafLinesFilter == other.pdafLinesFilter; } @@ -2426,10 +2426,10 @@ void RAWParams::BayerSensor::setPixelShiftDefaults() pixelShiftGreen = true; pixelShiftBlur = true; pixelShiftSmoothFactor = 0.7; - pixelShiftLmmse = false; pixelShiftEqualBright = false; pixelShiftEqualBrightChannel = false; pixelShiftNonGreenCross = true; + pixelShiftDemosaicMethod = getPSDemosaicMethodString(PSDemosaicMethod::AMAZE); } const std::vector& RAWParams::BayerSensor::getMethodStrings() @@ -2458,6 +2458,23 @@ Glib::ustring RAWParams::BayerSensor::getMethodString(Method method) return getMethodStrings()[toUnderlying(method)]; } +const std::vector& RAWParams::BayerSensor::getPSDemosaicMethodStrings() +{ + static const std::vector method_strings { + "amaze", + "amazevng4", + "lmmse" + }; + return method_strings; +} + +Glib::ustring RAWParams::BayerSensor::getPSDemosaicMethodString(PSDemosaicMethod method) +{ + return getPSDemosaicMethodStrings()[toUnderlying(method)]; +} + + + RAWParams::XTransSensor::XTransSensor() : method(getMethodString(Method::THREE_PASS)), ccSteps(0), @@ -3345,10 +3362,10 @@ int ProcParams::save(const Glib::ustring& fname, const Glib::ustring& fname2, bo saveToKeyfile(!pedited || pedited->raw.bayersensor.pixelShiftGreen, "RAW Bayer", "pixelShiftGreen", raw.bayersensor.pixelShiftGreen, keyFile); saveToKeyfile(!pedited || pedited->raw.bayersensor.pixelShiftBlur, "RAW Bayer", "pixelShiftBlur", raw.bayersensor.pixelShiftBlur, keyFile); saveToKeyfile(!pedited || pedited->raw.bayersensor.pixelShiftSmooth, "RAW Bayer", "pixelShiftSmoothFactor", raw.bayersensor.pixelShiftSmoothFactor, keyFile); - saveToKeyfile(!pedited || pedited->raw.bayersensor.pixelShiftLmmse, "RAW Bayer", "pixelShiftLmmse", raw.bayersensor.pixelShiftLmmse, keyFile); saveToKeyfile(!pedited || pedited->raw.bayersensor.pixelShiftEqualBright, "RAW Bayer", "pixelShiftEqualBright", raw.bayersensor.pixelShiftEqualBright, keyFile); saveToKeyfile(!pedited || pedited->raw.bayersensor.pixelShiftEqualBrightChannel, "RAW Bayer", "pixelShiftEqualBrightChannel", raw.bayersensor.pixelShiftEqualBrightChannel, keyFile); saveToKeyfile(!pedited || pedited->raw.bayersensor.pixelShiftNonGreenCross, "RAW Bayer", "pixelShiftNonGreenCross", raw.bayersensor.pixelShiftNonGreenCross, keyFile); + saveToKeyfile(!pedited || pedited->raw.bayersensor.pixelShiftDemosaicMethod, "RAW Bayer", "pixelShiftDemosaicMethod", raw.bayersensor.pixelShiftDemosaicMethod, keyFile); saveToKeyfile(!pedited || pedited->raw.bayersensor.pdafLinesFilter, "RAW Bayer", "PDAFLinesFilter", raw.bayersensor.pdafLinesFilter, keyFile); saveToKeyfile(!pedited || pedited->raw.xtranssensor.method, "RAW X-Trans", "Method", raw.xtranssensor.method, keyFile); saveToKeyfile(!pedited || pedited->raw.xtranssensor.ccSteps, "RAW X-Trans", "CcSteps", raw.xtranssensor.ccSteps, keyFile); @@ -4682,10 +4699,26 @@ int ProcParams::load(const Glib::ustring& fname, ParamsEdited* pedited) assignFromKeyfile(keyFile, "RAW Bayer", "pixelShiftGreen", pedited, raw.bayersensor.pixelShiftGreen, pedited->raw.bayersensor.pixelShiftGreen); assignFromKeyfile(keyFile, "RAW Bayer", "pixelShiftBlur", pedited, raw.bayersensor.pixelShiftBlur, pedited->raw.bayersensor.pixelShiftBlur); assignFromKeyfile(keyFile, "RAW Bayer", "pixelShiftSmoothFactor", pedited, raw.bayersensor.pixelShiftSmoothFactor, pedited->raw.bayersensor.pixelShiftSmooth); - assignFromKeyfile(keyFile, "RAW Bayer", "pixelShiftLmmse", pedited, raw.bayersensor.pixelShiftLmmse, pedited->raw.bayersensor.pixelShiftLmmse); assignFromKeyfile(keyFile, "RAW Bayer", "pixelShiftEqualBright", pedited, raw.bayersensor.pixelShiftEqualBright, pedited->raw.bayersensor.pixelShiftEqualBright); assignFromKeyfile(keyFile, "RAW Bayer", "pixelShiftEqualBrightChannel", pedited, raw.bayersensor.pixelShiftEqualBrightChannel, pedited->raw.bayersensor.pixelShiftEqualBrightChannel); assignFromKeyfile(keyFile, "RAW Bayer", "pixelShiftNonGreenCross", pedited, raw.bayersensor.pixelShiftNonGreenCross, pedited->raw.bayersensor.pixelShiftNonGreenCross); + + if (ppVersion < 336) { + if (keyFile.has_key("RAW Bayer", "pixelShiftLmmse")) { + bool useLmmse = keyFile.get_boolean ("RAW Bayer", "pixelShiftLmmse"); + if (useLmmse) { + raw.bayersensor.pixelShiftDemosaicMethod = raw.bayersensor.getPSDemosaicMethodString(raw.bayersensor.PSDemosaicMethod::LMMSE); + } else { + raw.bayersensor.pixelShiftDemosaicMethod = raw.bayersensor.getPSDemosaicMethodString(raw.bayersensor.PSDemosaicMethod::AMAZE); + } + if (pedited) { + pedited->raw.bayersensor.pixelShiftDemosaicMethod = true; + } + } + } else { + assignFromKeyfile(keyFile, "RAW Bayer", "pixelShiftDemosaicMethod", pedited, raw.bayersensor.pixelShiftDemosaicMethod, pedited->raw.bayersensor.pixelShiftDemosaicMethod); + } + assignFromKeyfile(keyFile, "RAW Bayer", "PDAFLinesFilter", pedited, raw.bayersensor.pdafLinesFilter, pedited->raw.bayersensor.pdafLinesFilter); } diff --git a/rtengine/procparams.h b/rtengine/procparams.h index be7784e1b..9b52cc828 100644 --- a/rtengine/procparams.h +++ b/rtengine/procparams.h @@ -1245,6 +1245,12 @@ struct RAWParams { CUSTOM }; + enum class PSDemosaicMethod { + AMAZE, + AMAZEVNG4, + LMMSE + }; + Glib::ustring method; int imageNum; int ccSteps; @@ -1275,10 +1281,10 @@ struct RAWParams { bool pixelShiftGreen; bool pixelShiftBlur; double pixelShiftSmoothFactor; - bool pixelShiftLmmse; bool pixelShiftEqualBright; bool pixelShiftEqualBrightChannel; bool pixelShiftNonGreenCross; + Glib::ustring pixelShiftDemosaicMethod; bool dcb_enhance; bool pdafLinesFilter; @@ -1291,6 +1297,9 @@ struct RAWParams { static const std::vector& getMethodStrings(); static Glib::ustring getMethodString(Method method); + + static const std::vector& getPSDemosaicMethodStrings(); + static Glib::ustring getPSDemosaicMethodString(PSDemosaicMethod method); }; /** diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index ce300bcd7..c5af124e5 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -2068,7 +2068,7 @@ void RawImageSource::demosaic(const RAWParams &raw) if ( raw.bayersensor.method == RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::HPHD) ) { hphd_demosaic (); } else if (raw.bayersensor.method == RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::VNG4) ) { - vng4_demosaic (); + vng4_demosaic (rawData, red, green, blue); } else if (raw.bayersensor.method == RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::AHD) ) { ahd_demosaic (); } else if (raw.bayersensor.method == RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::AMAZE) ) { diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h index e2467dba6..878886fed 100644 --- a/rtengine/rawimagesource.h +++ b/rtengine/rawimagesource.h @@ -262,13 +262,13 @@ protected: void nodemosaic(bool bw); void eahd_demosaic(); void hphd_demosaic(); - void vng4_demosaic(); + void vng4_demosaic(const array2D &rawData, array2D &red, array2D &green, array2D &blue); void ppg_demosaic(); void jdl_interpolate_omp(); void igv_interpolate(int winw, int winh); void lmmse_interpolate_omp(int winw, int winh, array2D &rawData, array2D &red, array2D &green, array2D &blue, int iterations); - void amaze_demosaic_RT(int winx, int winy, int winw, int winh, array2D &rawData, array2D &red, array2D &green, array2D &blue);//Emil's code for AMaZE - void amaze_vng4_demosaic_RT(int winw, int winh, array2D &rawData, array2D &red, array2D &green, array2D &blue, double contrast = 0.0); + void amaze_demosaic_RT(int winx, int winy, int winw, int winh, const array2D &rawData, array2D &red, array2D &green, array2D &blue);//Emil's code for AMaZE + void amaze_vng4_demosaic_RT(int winw, int winh, const array2D &rawData, array2D &red, array2D &green, array2D &blue, double contrast = 0.0); void fast_demosaic();//Emil's code for fast demosaicing void dcb_demosaic(int iterations, bool dcb_enhance); void ahd_demosaic(); diff --git a/rtengine/simpleprocess.cc b/rtengine/simpleprocess.cc index 7818871a7..ad622c6bd 100644 --- a/rtengine/simpleprocess.cc +++ b/rtengine/simpleprocess.cc @@ -1465,8 +1465,6 @@ private: } if (params.raw.bayersensor.method == procparams::RAWParams::BayerSensor::getMethodString(procparams::RAWParams::BayerSensor::Method::PIXELSHIFT)) { - params.raw.bayersensor.method = procparams::RAWParams::BayerSensor::getMethodString(params.raw.bayersensor.pixelShiftLmmse ? procparams::RAWParams::BayerSensor::Method::LMMSE : procparams::RAWParams::BayerSensor::Method::AMAZE); - } else if (params.raw.bayersensor.method == procparams::RAWParams::BayerSensor::getMethodString(procparams::RAWParams::BayerSensor::Method::AMAZE)) { params.raw.bayersensor.method = procparams::RAWParams::BayerSensor::getMethodString(procparams::RAWParams::BayerSensor::Method::RCD); } diff --git a/rtgui/bayerprocess.cc b/rtgui/bayerprocess.cc index 339666aa0..23e848fc6 100644 --- a/rtgui/bayerprocess.cc +++ b/rtgui/bayerprocess.cc @@ -29,6 +29,7 @@ BayerProcess::BayerProcess () : FoldableToolPanel(this, "bayerprocess", M("TP_RA auto m = ProcEventMapper::getInstance(); EvDemosaicContrast = m->newEvent(DEMOSAIC, "HISTORY_MSG_DUALDEMOSAIC_CONTRAST"); + EvDemosaicPixelshiftDemosaicMethod = m->newEvent(DEMOSAIC, "HISTORY_MSG_PIXELSHIFT_DEMOSAIC"); Gtk::HBox* hb1 = Gtk::manage (new Gtk::HBox ()); hb1->pack_start (*Gtk::manage (new Gtk::Label ( M("TP_RAW_DMETHOD") + ": ")), Gtk::PACK_SHRINK, 4); @@ -149,6 +150,17 @@ BayerProcess::BayerProcess () : FoldableToolPanel(this, "bayerprocess", M("TP_RA pixelShiftShowMotionMaskOnly->set_tooltip_text (M("TP_RAW_PIXELSHIFTSHOWMOTIONMASKONLY_TOOLTIP")); pixelShiftFrame->pack_start(*pixelShiftShowMotionMaskOnly); + + Gtk::HBox* hb4 = Gtk::manage (new Gtk::HBox ()); + hb4->pack_start (*Gtk::manage (new Gtk::Label ( M("TP_RAW_PIXELSHIFTDMETHOD") + ": ")), Gtk::PACK_SHRINK, 4); + pixelShiftDemosaicMethod = Gtk::manage (new MyComboBoxText ()); + for(const auto method_string : procparams::RAWParams::BayerSensor::getPSDemosaicMethodStrings()) { + pixelShiftDemosaicMethod->append(M("TP_RAW_" + Glib::ustring(method_string).uppercase())); + } + pixelShiftDemosaicMethod->set_active(0); + hb4->pack_start(*pixelShiftDemosaicMethod); + pixelShiftOptions->pack_start(*hb4); + pixelShiftGreen = Gtk::manage (new CheckBox(M("TP_RAW_PIXELSHIFTGREEN"), multiImage)); pixelShiftGreen->setCheckBoxListener (this); pixelShiftOptions->pack_start(*pixelShiftGreen); @@ -207,13 +219,6 @@ BayerProcess::BayerProcess () : FoldableToolPanel(this, "bayerprocess", M("TP_RA pixelShiftMedian->set_tooltip_text (M("TP_RAW_PIXELSHIFTMEDIAN_TOOLTIP")); pixelShiftOptions->pack_start(*pixelShiftMedian); - - pixelShiftLmmse = Gtk::manage (new CheckBox(M("TP_RAW_PIXELSHIFTLMMSE"), multiImage)); - pixelShiftLmmse->setCheckBoxListener (this); - pixelShiftLmmse->set_tooltip_text (M("TP_RAW_PIXELSHIFTLMMSE_TOOLTIP")); - pixelShiftOptions->pack_start(*pixelShiftLmmse); - - pixelShiftFrame->pack_start(*pixelShiftOptions); pixelShiftOptions->hide(); @@ -222,6 +227,7 @@ BayerProcess::BayerProcess () : FoldableToolPanel(this, "bayerprocess", M("TP_RA method->connect(method->signal_changed().connect( sigc::mem_fun(*this, &BayerProcess::methodChanged) )); imageNumber->connect(imageNumber->signal_changed().connect( sigc::mem_fun(*this, &BayerProcess::imageNumberChanged) )); pixelShiftMotionMethod->connect(pixelShiftMotionMethod->signal_changed().connect( sigc::mem_fun(*this, &BayerProcess::pixelShiftMotionMethodChanged) )); + pixelShiftDemosaicMethod->connect(pixelShiftDemosaicMethod->signal_changed().connect( sigc::mem_fun(*this, &BayerProcess::pixelShiftDemosaicMethodChanged) )); } @@ -230,6 +236,7 @@ void BayerProcess::read(const rtengine::procparams::ProcParams* pp, const Params disableListener (); method->block (true); imageNumber->block (true); + pixelShiftDemosaicMethod->block(true); //allEnhconn.block (true); method->set_active(std::numeric_limits::max()); @@ -242,8 +249,14 @@ void BayerProcess::read(const rtengine::procparams::ProcParams* pp, const Params break; } } + pixelShiftDemosaicMethod->set_active(std::numeric_limits::max()); + for (size_t i = 0; i < procparams::RAWParams::BayerSensor::getPSDemosaicMethodStrings().size(); ++i) { + if (pp->raw.bayersensor.pixelShiftDemosaicMethod == procparams::RAWParams::BayerSensor::getPSDemosaicMethodStrings()[i]) { + pixelShiftDemosaicMethod->set_active(i); + break; + } + } - //allEnhance->set_active(pp->raw.bayersensor.all_enhance); dcbIterations->setValue (pp->raw.bayersensor.dcb_iterations); dcbEnhance->setValue (pp->raw.bayersensor.dcb_enhance); @@ -260,7 +273,6 @@ void BayerProcess::read(const rtengine::procparams::ProcParams* pp, const Params pixelShiftSmooth->set_sensitive (pp->raw.bayersensor.pixelShiftBlur); } pixelShiftSmooth->setValue (pp->raw.bayersensor.pixelShiftSmoothFactor); - pixelShiftLmmse->setValue (pp->raw.bayersensor.pixelShiftLmmse); pixelShiftEqualBright->setValue (pp->raw.bayersensor.pixelShiftEqualBright); pixelShiftEqualBrightChannel->set_sensitive (pp->raw.bayersensor.pixelShiftEqualBright); pixelShiftEqualBrightChannel->setValue (pp->raw.bayersensor.pixelShiftEqualBrightChannel); @@ -286,7 +298,6 @@ void BayerProcess::read(const rtengine::procparams::ProcParams* pp, const Params pixelShiftGreen->setEdited (pedited->raw.bayersensor.pixelShiftGreen); pixelShiftBlur->setEdited (pedited->raw.bayersensor.pixelShiftBlur); pixelShiftSmooth->setEditedState ( pedited->raw.bayersensor.pixelShiftSmooth ? Edited : UnEdited); - pixelShiftLmmse->setEdited (pedited->raw.bayersensor.pixelShiftLmmse); pixelShiftEqualBright->setEdited (pedited->raw.bayersensor.pixelShiftEqualBright); pixelShiftEqualBrightChannel->setEdited (pedited->raw.bayersensor.pixelShiftEqualBrightChannel); pixelShiftNonGreenCross->setEdited (pedited->raw.bayersensor.pixelShiftNonGreenCross); @@ -298,12 +309,19 @@ void BayerProcess::read(const rtengine::procparams::ProcParams* pp, const Params if(!pedited->raw.bayersensor.method) { method->set_active(std::numeric_limits::max()); // No name } + if(!pedited->raw.bayersensor.imageNum) { imageNumber->set_active_text(M("GENERAL_UNCHANGED")); } + if(!pedited->raw.bayersensor.pixelShiftMotionCorrectionMethod) { pixelShiftMotionMethod->set_active_text(M("GENERAL_UNCHANGED")); } + + if(!pedited->raw.bayersensor.pixelShiftDemosaicMethod) { + pixelShiftDemosaicMethod->set_active(std::numeric_limits::max()); // No name + } + } if (!batchMode) { @@ -330,6 +348,7 @@ void BayerProcess::read(const rtengine::procparams::ProcParams* pp, const Params method->block (false); imageNumber->block (false); + pixelShiftDemosaicMethod->block(false); //allEnhconn.block (false); enableListener (); @@ -353,7 +372,6 @@ void BayerProcess::write( rtengine::procparams::ProcParams* pp, ParamsEdited* pe pp->raw.bayersensor.pixelShiftGreen = pixelShiftGreen->getLastActive (); pp->raw.bayersensor.pixelShiftBlur = pixelShiftBlur->getLastActive (); pp->raw.bayersensor.pixelShiftSmoothFactor = pixelShiftSmooth->getValue(); - pp->raw.bayersensor.pixelShiftLmmse = pixelShiftLmmse->getLastActive (); pp->raw.bayersensor.pixelShiftEqualBright = pixelShiftEqualBright->getLastActive (); pp->raw.bayersensor.pixelShiftEqualBrightChannel = pixelShiftEqualBrightChannel->getLastActive (); pp->raw.bayersensor.pixelShiftNonGreenCross = pixelShiftNonGreenCross->getLastActive (); @@ -368,6 +386,11 @@ void BayerProcess::write( rtengine::procparams::ProcParams* pp, ParamsEdited* pe pp->raw.bayersensor.imageNum = currentRow; } + currentRow = pixelShiftDemosaicMethod->get_active_row_number(); + if( currentRow >= 0 && currentRow < std::numeric_limits::max()) { + pp->raw.bayersensor.pixelShiftDemosaicMethod = procparams::RAWParams::BayerSensor::getPSDemosaicMethodString(RAWParams::BayerSensor::PSDemosaicMethod(currentRow)); + } + if (pedited) { pedited->raw.bayersensor.ccSteps = ccSteps->getEditedState (); @@ -379,6 +402,7 @@ void BayerProcess::write( rtengine::procparams::ProcParams* pp, ParamsEdited* pe pedited->raw.bayersensor.lmmseIterations = lmmseIterations->getEditedState (); pedited->raw.bayersensor.dualDemosaicContrast = dualDemosaicContrast->getEditedState (); pedited->raw.bayersensor.pixelShiftMotionCorrectionMethod = pixelShiftMotionMethod->get_active_text() != M("GENERAL_UNCHANGED"); + pedited->raw.bayersensor.pixelShiftDemosaicMethod = pixelShiftDemosaicMethod->get_active_row_number() != std::numeric_limits::max(); pedited->raw.bayersensor.pixelShiftEperIso = pixelShiftEperIso->getEditedState (); pedited->raw.bayersensor.pixelShiftSigma = pixelShiftSigma->getEditedState (); pedited->raw.bayersensor.pixelShiftShowMotion = !pixelShiftShowMotion->get_inconsistent(); @@ -388,7 +412,6 @@ void BayerProcess::write( rtengine::procparams::ProcParams* pp, ParamsEdited* pe pedited->raw.bayersensor.pixelShiftGreen = !pixelShiftGreen->get_inconsistent(); pedited->raw.bayersensor.pixelShiftBlur = !pixelShiftBlur->get_inconsistent(); pedited->raw.bayersensor.pixelShiftSmooth = pixelShiftSmooth->getEditedState(); - pedited->raw.bayersensor.pixelShiftLmmse = !pixelShiftLmmse->get_inconsistent(); pedited->raw.bayersensor.pixelShiftEqualBright = !pixelShiftEqualBright->get_inconsistent(); pedited->raw.bayersensor.pixelShiftEqualBrightChannel = !pixelShiftEqualBrightChannel->get_inconsistent(); pedited->raw.bayersensor.pixelShiftNonGreenCross = !pixelShiftNonGreenCross->get_inconsistent(); @@ -423,6 +446,8 @@ void BayerProcess::setBatchMode(bool batchMode) method->set_active(std::numeric_limits::max()); // No name pixelShiftMotionMethod->append (M("GENERAL_UNCHANGED")); pixelShiftMotionMethod->set_active_text (M("GENERAL_UNCHANGED")); + pixelShiftDemosaicMethod->append (M("GENERAL_UNCHANGED")); + pixelShiftDemosaicMethod->set_active(std::numeric_limits::max()); // No name imageNumber->append (M("GENERAL_UNCHANGED")); imageNumber->set_active_text (M("GENERAL_UNCHANGED")); ToolPanel::setBatchMode (batchMode); @@ -535,6 +560,19 @@ void BayerProcess::methodChanged () } } +void BayerProcess::pixelShiftDemosaicMethodChanged () +{ + + if (listener) { + const int curSelection = pixelShiftDemosaicMethod->get_active_row_number(); + if (curSelection >= 0 && curSelection < std::numeric_limits::max()) { + const RAWParams::BayerSensor::PSDemosaicMethod method = RAWParams::BayerSensor::PSDemosaicMethod(curSelection); + Glib::ustring methodName = procparams::RAWParams::BayerSensor::getPSDemosaicMethodString(method);; + listener->panelChanged (EvDemosaicPixelshiftDemosaicMethod, methodName); + } + } +} + void BayerProcess::imageNumberChanged () { if (listener) { @@ -579,10 +617,6 @@ void BayerProcess::checkBoxToggled (CheckBox* c, CheckValue newval) if (listener) { listener->panelChanged (EvPixelShiftBlur, pixelShiftBlur->getValueAsStr ()); } - } else if (c == pixelShiftLmmse) { - if (listener) { - listener->panelChanged (EvPixelShiftLmmse, pixelShiftLmmse->getValueAsStr ()); - } } else if (c == pixelShiftEqualBright) { if (!batchMode) { pixelShiftEqualBrightChannel->set_sensitive(newval != CheckValue::off); diff --git a/rtgui/bayerprocess.h b/rtgui/bayerprocess.h index ee532f27c..bfcb99878 100644 --- a/rtgui/bayerprocess.h +++ b/rtgui/bayerprocess.h @@ -42,6 +42,7 @@ protected: Gtk::VBox *pixelShiftFrame; Gtk::VBox *pixelShiftOptions; MyComboBoxText* pixelShiftMotionMethod; + MyComboBoxText* pixelShiftDemosaicMethod; CheckBox* pixelShiftShowMotion; CheckBox* pixelShiftShowMotionMaskOnly; CheckBox* pixelShiftNonGreenCross; @@ -49,7 +50,6 @@ protected: CheckBox* pixelShiftBlur; CheckBox* pixelShiftHoleFill; CheckBox* pixelShiftMedian; - CheckBox* pixelShiftLmmse; CheckBox* pixelShiftEqualBright; CheckBox* pixelShiftEqualBrightChannel; Adjuster* pixelShiftSmooth; @@ -61,6 +61,7 @@ protected: IdleRegister idle_register; rtengine::ProcEvent EvDemosaicContrast; + rtengine::ProcEvent EvDemosaicPixelshiftDemosaicMethod; public: BayerProcess (); @@ -77,6 +78,7 @@ public: void adjusterChanged(Adjuster* a, double newval); void checkBoxToggled(CheckBox* c, CheckValue newval); void pixelShiftMotionMethodChanged(); + void pixelShiftDemosaicMethodChanged(); void FrameCountChanged(int n, int frameNum); }; diff --git a/rtgui/paramsedited.cc b/rtgui/paramsedited.cc index bf5a63fa8..9e6d3b6e7 100644 --- a/rtgui/paramsedited.cc +++ b/rtgui/paramsedited.cc @@ -414,10 +414,10 @@ void ParamsEdited::set (bool v) raw.bayersensor.pixelShiftGreen = v; raw.bayersensor.pixelShiftBlur = v; raw.bayersensor.pixelShiftSmooth = v; - raw.bayersensor.pixelShiftLmmse = v; raw.bayersensor.pixelShiftEqualBright = v; raw.bayersensor.pixelShiftEqualBrightChannel = v; raw.bayersensor.pixelShiftNonGreenCross = v; + raw.bayersensor.pixelShiftDemosaicMethod = v; raw.bayersensor.greenEq = v; raw.bayersensor.linenoise = v; raw.bayersensor.linenoiseDirection = v; @@ -964,10 +964,10 @@ void ParamsEdited::initFrom (const std::vector raw.bayersensor.pixelShiftGreen = raw.bayersensor.pixelShiftGreen && p.raw.bayersensor.pixelShiftGreen == other.raw.bayersensor.pixelShiftGreen; raw.bayersensor.pixelShiftBlur = raw.bayersensor.pixelShiftBlur && p.raw.bayersensor.pixelShiftBlur == other.raw.bayersensor.pixelShiftBlur; raw.bayersensor.pixelShiftSmooth = raw.bayersensor.pixelShiftSmooth && p.raw.bayersensor.pixelShiftSmoothFactor == other.raw.bayersensor.pixelShiftSmoothFactor; - raw.bayersensor.pixelShiftLmmse = raw.bayersensor.pixelShiftLmmse && p.raw.bayersensor.pixelShiftLmmse == other.raw.bayersensor.pixelShiftLmmse; raw.bayersensor.pixelShiftEqualBright = raw.bayersensor.pixelShiftEqualBright && p.raw.bayersensor.pixelShiftEqualBright == other.raw.bayersensor.pixelShiftEqualBright; raw.bayersensor.pixelShiftEqualBrightChannel = raw.bayersensor.pixelShiftEqualBrightChannel && p.raw.bayersensor.pixelShiftEqualBrightChannel == other.raw.bayersensor.pixelShiftEqualBrightChannel; raw.bayersensor.pixelShiftNonGreenCross = raw.bayersensor.pixelShiftNonGreenCross && p.raw.bayersensor.pixelShiftNonGreenCross == other.raw.bayersensor.pixelShiftNonGreenCross; + raw.bayersensor.pixelShiftDemosaicMethod = raw.bayersensor.pixelShiftDemosaicMethod && p.raw.bayersensor.pixelShiftDemosaicMethod == other.raw.bayersensor.pixelShiftDemosaicMethod; raw.bayersensor.greenEq = raw.bayersensor.greenEq && p.raw.bayersensor.greenthresh == other.raw.bayersensor.greenthresh; raw.bayersensor.linenoise = raw.bayersensor.linenoise && p.raw.bayersensor.linenoise == other.raw.bayersensor.linenoise; raw.bayersensor.linenoiseDirection = raw.bayersensor.linenoiseDirection && p.raw.bayersensor.linenoiseDirection == other.raw.bayersensor.linenoiseDirection; @@ -2542,10 +2542,6 @@ void ParamsEdited::combine (rtengine::procparams::ProcParams& toEdit, const rten toEdit.raw.bayersensor.pixelShiftSmoothFactor = mods.raw.bayersensor.pixelShiftSmoothFactor; } - if (raw.bayersensor.pixelShiftLmmse) { - toEdit.raw.bayersensor.pixelShiftLmmse = mods.raw.bayersensor.pixelShiftLmmse; - } - if (raw.bayersensor.pixelShiftEqualBright) { toEdit.raw.bayersensor.pixelShiftEqualBright = mods.raw.bayersensor.pixelShiftEqualBright; } @@ -2558,6 +2554,10 @@ void ParamsEdited::combine (rtengine::procparams::ProcParams& toEdit, const rten toEdit.raw.bayersensor.pixelShiftNonGreenCross = mods.raw.bayersensor.pixelShiftNonGreenCross; } + if (raw.bayersensor.pixelShiftDemosaicMethod) { + toEdit.raw.bayersensor.pixelShiftDemosaicMethod = mods.raw.bayersensor.pixelShiftDemosaicMethod; + } + if (raw.bayersensor.greenEq) { toEdit.raw.bayersensor.greenthresh = dontforceSet && options.baBehav[ADDSET_PREPROCESS_GREENEQUIL] ? toEdit.raw.bayersensor.greenthresh + mods.raw.bayersensor.greenthresh : mods.raw.bayersensor.greenthresh; } @@ -3080,7 +3080,7 @@ bool RAWParamsEdited::BayerSensor::isUnchanged() const { return method && imageNum && dcbIterations && dcbEnhance && lmmseIterations && dualDemosaicContrast /*&& allEnhance*/ && greenEq && pixelShiftMotionCorrectionMethod && pixelShiftEperIso && pixelShiftSigma && pixelShiftShowMotion && pixelShiftShowMotionMaskOnly - && pixelShiftHoleFill && pixelShiftMedian && pixelShiftNonGreenCross && pixelShiftGreen && pixelShiftBlur && pixelShiftSmooth && pixelShiftLmmse && pixelShiftEqualBright && pixelShiftEqualBrightChannel + && pixelShiftHoleFill && pixelShiftMedian && pixelShiftNonGreenCross && pixelShiftDemosaicMethod && pixelShiftGreen && pixelShiftBlur && pixelShiftSmooth && pixelShiftEqualBright && pixelShiftEqualBrightChannel && linenoise && linenoiseDirection && pdafLinesFilter && exBlack0 && exBlack1 && exBlack2 && exBlack3 && exTwoGreen; } diff --git a/rtgui/paramsedited.h b/rtgui/paramsedited.h index 5b3370b9a..1db298dfd 100644 --- a/rtgui/paramsedited.h +++ b/rtgui/paramsedited.h @@ -742,10 +742,10 @@ public: bool pixelShiftGreen; bool pixelShiftBlur; bool pixelShiftSmooth; - bool pixelShiftLmmse; bool pixelShiftEqualBright; bool pixelShiftEqualBrightChannel; bool pixelShiftNonGreenCross; + bool pixelShiftDemosaicMethod; //bool allEnhance; bool greenEq; diff --git a/rtgui/partialpastedlg.cc b/rtgui/partialpastedlg.cc index c9ed1cd14..99411e46c 100644 --- a/rtgui/partialpastedlg.cc +++ b/rtgui/partialpastedlg.cc @@ -852,7 +852,7 @@ void PartialPasteDlg::applyPaste (rtengine::procparams::ProcParams* dstPP, Param filterPE.raw.bayersensor.pixelShiftEqualBrightChannel = falsePE.raw.bayersensor.pixelShiftEqualBrightChannel; filterPE.raw.bayersensor.pixelShiftGreen = falsePE.raw.bayersensor.pixelShiftGreen; filterPE.raw.bayersensor.pixelShiftHoleFill = falsePE.raw.bayersensor.pixelShiftHoleFill; - filterPE.raw.bayersensor.pixelShiftLmmse = falsePE.raw.bayersensor.pixelShiftLmmse; + filterPE.raw.bayersensor.pixelShiftDemosaicMethod = falsePE.raw.bayersensor.pixelShiftDemosaicMethod; filterPE.raw.bayersensor.pixelShiftMedian = falsePE.raw.bayersensor.pixelShiftMedian; filterPE.raw.bayersensor.pixelShiftMotionCorrectionMethod = falsePE.raw.bayersensor.pixelShiftMotionCorrectionMethod; filterPE.raw.bayersensor.pixelShiftNonGreenCross = falsePE.raw.bayersensor.pixelShiftNonGreenCross; diff --git a/rtgui/ppversion.h b/rtgui/ppversion.h index 2b6bc49de..1fb34ca0a 100644 --- a/rtgui/ppversion.h +++ b/rtgui/ppversion.h @@ -1,11 +1,13 @@ #pragma once // This number has to be incremented whenever the PP3 file format is modified or the behaviour of a tool changes -#define PPVERSION 335 +#define PPVERSION 336 #define PPVERSION_AEXP 301 //value of PPVERSION when auto exposure algorithm was modified /* Log of version changes + 336 2018-06-01 + new demosaic method combobox for pixelshift 335 2018-05-30 new contrast adjuster in Bayer process tool 334 2018-05-13