lmmse.c

Go to the documentation of this file.

/*
    This file is part of darktable,
    Copyright (C) 2021 Hanno Schwalm.
    Copyright (C) 2021 luzpaz.
    Copyright (C) 2021 Pascal Obry.
    Copyright (C) 2022 Martin Bařinka.
    Copyright (C) 2023, 2026 Aurélien PIERRE.
    
    darktable 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.
    
    darktable 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 darktable.  If not, see <http://www.gnu.org/licenses/>.
*/
 
/*
    The lmmse code base used for the darktable port has been taken from rawtherapee derived librtprocess.
    Adapt for dt and tiling - hanno schwalm 06/2021
 
    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
*/
 
/*
   Refinement based on EECI demosaicing algorithm by L. Chang and Y.P. Tan
   Paul Lee
   Adapted for RawTherapee - Jacques Desmis 04/2013
*/
 
/* Why tiling?
   The internal tiling vastly reduces memory footprint and allows data processing to be done mostly
   with in-cache data thus increasing performance.
 
   The performance has been tested on a E-2288G for 45mpix images, tiling improves performance > 2-fold.
   times in sec: basic (0.5->0.15), median (0.6->0.18), 3xmedian (0.8->0.22), 3xmedian + 2x refine (1.2->0.30)
   The default is now 2 times slower than RCD and 2 times faster than AMaZE
*/
 
#ifndef LMMSE_GRP
  #define LMMSE_GRP 136
#endif
 
#ifdef __GNUC__
  #pragma GCC push_options
  #pragma GCC optimize ("fast-math", "fp-contract=fast", "finite-math-only", "no-math-errno")
#endif
 
#define LMMSE_OVERLAP 8
#define BORDER_AROUND 4
#define LMMSE_TILESIZE (LMMSE_GRP - 2 * BORDER_AROUND)
#define LMMSE_TILEVALID (LMMSE_TILESIZE - 2 * LMMSE_OVERLAP)
#define w1 (LMMSE_GRP)
#define w2 (LMMSE_GRP * 2)
#define w3 (LMMSE_GRP * 3)
#define w4 (LMMSE_GRP * 4)
 
static INLINE float limf(float x, float min, float max)
{
  return fmaxf(min, fminf(x, max));
}
 
static INLINE float median3f(float x0, float x1, float x2)
{
  return fmaxf(fminf(x0,x1), fminf(x2, fmaxf(x0,x1)));
}
 
static INLINE float median9f(float a0, float a1, float a2, float a3, float a4, float a5, float a6, float a7, float a8)
{
  float tmp;
  tmp = fminf(a1, a2);
  a2  = fmaxf(a1, a2);
  a1  = tmp;
  tmp = fminf(a4, a5);
  a5  = fmaxf(a4, a5);
  a4  = tmp;
  tmp = fminf(a7, a8);
  a8  = fmaxf(a7, a8);
  a7  = tmp;
  tmp = fminf(a0, a1);
  a1  = fmaxf(a0, a1);
  a0  = tmp;
  tmp = fminf(a3, a4);
  a4  = fmaxf(a3, a4);
  a3  = tmp;
  tmp = fminf(a6, a7);
  a7  = fmaxf(a6, a7);
  a6  = tmp;
  tmp = fminf(a1, a2);
  a2  = fmaxf(a1, a2);
  a1  = tmp;
  tmp = fminf(a4, a5);
  a5  = fminf(a4, a5);
  a4  = tmp;
  tmp = fminf(a7, a8);
  a8  = fmaxf(a7, a8);
  a3  = fmaxf(a0, a3);
  a5  = fminf(a5, a8);
  a7  = fmaxf(a4, tmp);
  tmp = fminf(a4, tmp);
  a6  = fmaxf(a3, a6);
  a4  = fmaxf(a1, tmp);
  a2  = fminf(a2, a5);
  a4  = fminf(a4, a7);
  tmp = fminf(a4, a2);
  a2  = fmaxf(a4, a2);
  a4  = fmaxf(a6, tmp);
  return fminf(a4, a2);
}
 
static INLINE float calc_gamma(float val, float *table)
{
  const float index = val * 65535.0f;
  if(index < 0.0f)      return 0.0f;
  if(index > 65534.99f) return 1.0f;
  const int idx = (int)index;
 
  const float diff = index - (float)idx;
  const float p1 = table[idx];
  const float p2 = table[idx+1] - p1;
  return (p1 + p2 * diff);
}
 
#ifdef _OPENMP
  #pragma omp declare simd aligned(in, out, gamma_in, gamma_out)
#endif
static void lmmse_demosaic(const dt_dev_pixelpipe_iop_t *piece, float *const restrict out, const float *const restrict in, dt_iop_roi_t *const roi_out,
                                   const dt_iop_roi_t *const roi_in, const uint32_t filters, const uint32_t mode, float *const restrict gamma_in, float *const restrict gamma_out)
{
  const int width = roi_in->width;
  const int height = roi_in->height;
 
  if((width < 16) || (height < 16))
  {
    dt_control_log(_("[lmmse_demosaic] too small area"));
    return;
  }
 
  float h0 = 1.0f;
  float h1 = expf( -1.0f / 8.0f);
  float h2 = expf( -4.0f / 8.0f);
  float h3 = expf( -9.0f / 8.0f);
  float h4 = expf(-16.0f / 8.0f);
  float hs = h0 + 2.0f * (h1 + h2 + h3 + h4);
  h0 /= hs;
  h1 /= hs;
  h2 /= hs;
  h3 /= hs;
  h4 /= hs;
 
  // median filter iterations
  const int medians = (mode < 2) ? mode : 3;
  // refinement steps
  const int refine = (mode > 2) ? mode - 2 : 0;
 
  const float scaler = fmaxf(piece->dsc_in.processed_maximum[0], fmaxf(piece->dsc_in.processed_maximum[1], piece->dsc_in.processed_maximum[2]));
  const float revscaler = 1.0f / scaler;
 
  const int num_vertical =   1 + (height - 2 * LMMSE_OVERLAP -1) / LMMSE_TILEVALID;
  const int num_horizontal = 1 + (width  - 2 * LMMSE_OVERLAP -1) / LMMSE_TILEVALID;
#ifdef _OPENMP
  #pragma omp parallel \
  dt_omp_firstprivate(width, height, out, in, scaler, revscaler, filters)
#endif
  {
    float *qix[6];
    float *buffer = dt_pixelpipe_cache_alloc_align_float_cache(LMMSE_GRP * LMMSE_GRP * 6, 0);
 
    qix[0] = buffer;
    for(int i = 1; i < 6; i++)
    {
      qix[i] = qix[i - 1] + LMMSE_GRP * LMMSE_GRP;
    }
    memset_zero(buffer, sizeof(float) * LMMSE_GRP * LMMSE_GRP * 6);
 
#ifdef _OPENMP
  #pragma omp for schedule(simd:dynamic, 6) collapse(2)
#endif
    for(int tile_vertical = 0; tile_vertical < num_vertical; tile_vertical++)
    {
      for(int tile_horizontal = 0; tile_horizontal < num_horizontal; tile_horizontal++)
      {
        const int rowStart = tile_vertical * LMMSE_TILEVALID;
        const int rowEnd = MIN(rowStart + LMMSE_TILESIZE, height);
 
        const int colStart = tile_horizontal * LMMSE_TILEVALID;
        const int colEnd = MIN(colStart + LMMSE_TILESIZE, width);
 
        const int tileRows = MIN(rowEnd - rowStart, LMMSE_TILESIZE);
        const int tileCols = MIN(colEnd - colStart, LMMSE_TILESIZE);
 
        // index limit; normally is LMMSE_GRP but maybe missing bottom lines or right columns for outermost tile
        const int last_rr = tileRows + 2 * BORDER_AROUND;
        const int last_cc = tileCols + 2 * BORDER_AROUND;
 
        for(int rrr = BORDER_AROUND, row = rowStart; rrr < tileRows + BORDER_AROUND; rrr++, row++)
        {
          float *cfa = qix[5] + rrr * LMMSE_GRP + BORDER_AROUND;
          int idx = row * width + colStart;
          for(int ccc = BORDER_AROUND, col = colStart; ccc < tileCols + BORDER_AROUND; ccc++, col++, cfa++, idx++)
          {
            cfa[0] = calc_gamma(revscaler * in[idx], gamma_in);
          }
        }
 
        // G-R(B)
        for(int rr = 2; rr < last_rr - 2; rr++)
        {
          // G-R(B) at R(B) location
          for(int cc = 2 + (FC(rr, 2, filters) & 1); cc < last_cc - 2; cc += 2)
          {
            float *cfa = qix[5] + rr * LMMSE_GRP + cc;
            const float v0 = 0.0625f * (cfa[-w1 - 1] + cfa[-w1 + 1] + cfa[w1 - 1] + cfa[w1 + 1]) + 0.25f * cfa[0];
            // horizontal
            float *hdiff = qix[0] + rr * LMMSE_GRP + cc;
            hdiff[0] = -0.25f * (cfa[ -2] + cfa[ 2]) + 0.5f * (cfa[ -1] + cfa[0] + cfa[ 1]);
            const float Y0 = v0 + 0.5f * hdiff[0];
            hdiff[0] = (cfa[0] > 1.75f * Y0) ? median3f(hdiff[0], cfa[ -1], cfa[ 1]) : limf(hdiff[0], 0.0f, 1.0f);
            hdiff[0] -= cfa[0];
 
            // vertical
            float *vdiff = qix[1] + rr * LMMSE_GRP + cc;
            vdiff[0] = -0.25f * (cfa[-w2] + cfa[w2]) + 0.5f * (cfa[-w1] + cfa[0] + cfa[w1]);
            const float Y1 = v0 + 0.5f * vdiff[0];
            vdiff[0] = (cfa[0] > 1.75f * Y1) ? median3f(vdiff[0], cfa[-w1], cfa[w1]) : limf(vdiff[0], 0.0f, 1.0f);
            vdiff[0] -= cfa[0];
          }
 
          // G-R(B) at G location
          for(int ccc = 2 + (FC(rr, 3, filters) & 1); ccc < last_cc - 2; ccc += 2)
          {
            float *cfa = qix[5] + rr * LMMSE_GRP + ccc;
            float *hdiff = qix[0] + rr * LMMSE_GRP + ccc;
            float *vdiff = qix[1] + rr * LMMSE_GRP + ccc;
            hdiff[0] = 0.25f * (cfa[ -2] + cfa[ 2]) - 0.5f * (cfa[ -1] + cfa[0] + cfa[ 1]);
            vdiff[0] = 0.25f * (cfa[-w2] + cfa[w2]) - 0.5f * (cfa[-w1] + cfa[0] + cfa[w1]);
            hdiff[0] = limf(hdiff[0], -1.0f, 0.0f) + cfa[0];
            vdiff[0] = limf(vdiff[0], -1.0f, 0.0f) + cfa[0];
          }
        }
 
        // apply low pass filter on differential colors
        for (int rr = 4; rr < last_rr - 4; rr++)
        {
          for(int cc = 4; cc < last_cc - 4; cc++)
          {
            float *hdiff = qix[0] + rr * LMMSE_GRP + cc;
            float *vdiff = qix[1] + rr * LMMSE_GRP + cc;
            float *hlp   = qix[2] + rr * LMMSE_GRP + cc;
            float *vlp   = qix[3] + rr * LMMSE_GRP + cc;
            hlp[0] = h0 * hdiff[0] + h1 * (hdiff[ -1] + hdiff[ 1]) + h2 * (hdiff[ -2] + hdiff[ 2]) + h3 * (hdiff[ -3] + hdiff[ 3]) + h4 * (hdiff[ -4] + hdiff[ 4]);
            vlp[0] = h0 * vdiff[0] + h1 * (vdiff[-w1] + vdiff[w1]) + h2 * (vdiff[-w2] + vdiff[w2]) + h3 * (vdiff[-w3] + vdiff[w3]) + h4 * (vdiff[-w4] + vdiff[w4]);
          }
        }
 
        for(int rr = 4; rr < last_rr - 4; rr++)
        {
          for(int cc = 4 + (FC(rr, 4, filters) & 1); cc < last_cc - 4; cc += 2)
          {
            float *hdiff = qix[0] + rr * LMMSE_GRP + cc;
            float *vdiff = qix[1] + rr * LMMSE_GRP + cc;
            float *hlp   = qix[2] + rr * LMMSE_GRP + cc;
            float *vlp   = qix[3] + rr * LMMSE_GRP + cc;
            float *interp = qix[4] + rr * LMMSE_GRP + cc;
            // horizontal
            float p1 = hlp[-4];
            float p2 = hlp[-3];
            float p3 = hlp[-2];
            float p4 = hlp[-1];
            float p5 = hlp[ 0];
            float p6 = hlp[ 1];
            float p7 = hlp[ 2];
            float p8 = hlp[ 3];
            float p9 = hlp[ 4];
            float mu = (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9) / 9.0f;
            float vx = 1e-7f + sqf(p1 - mu) + sqf(p2 - mu) + sqf(p3 - mu) + sqf(p4 - mu) + sqf(p5 - mu) + sqf(p6 - mu) + sqf(p7 - mu) + sqf(p8 - mu) + sqf(p9 - mu);
            p1 -= hdiff[-4];
            p2 -= hdiff[-3];
            p3 -= hdiff[-2];
            p4 -= hdiff[-1];
            p5 -= hdiff[ 0];
            p6 -= hdiff[ 1];
            p7 -= hdiff[ 2];
            p8 -= hdiff[ 3];
            p9 -= hdiff[ 4];
            float vn = 1e-7f + sqf(p1) + sqf(p2) + sqf(p3) + sqf(p4) + sqf(p5) + sqf(p6) + sqf(p7) + sqf(p8) + sqf(p9);
            float xh = (hdiff[0] * vx + hlp[0] * vn) / (vx + vn);
            float vh = vx * vn / (vx + vn);
 
            // vertical
            p1 = vlp[-w4];
            p2 = vlp[-w3];
            p3 = vlp[-w2];
            p4 = vlp[-w1];
            p5 = vlp[  0];
            p6 = vlp[ w1];
            p7 = vlp[ w2];
            p8 = vlp[ w3];
            p9 = vlp[ w4];
            mu = (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9) / 9.0f;
            vx = 1e-7f + sqf(p1 - mu) + sqf(p2 - mu) + sqf(p3 - mu) + sqf(p4 - mu) + sqf(p5 - mu) + sqf(p6 - mu) + sqf(p7 - mu) + sqf(p8 - mu) + sqf(p9 - mu);
            p1 -= vdiff[-w4];
            p2 -= vdiff[-w3];
            p3 -= vdiff[-w2];
            p4 -= vdiff[-w1];
            p5 -= vdiff[  0];
            p6 -= vdiff[ w1];
            p7 -= vdiff[ w2];
            p8 -= vdiff[ w3];
            p9 -= vdiff[ w4];
            vn = 1e-7f + sqf(p1) + sqf(p2) + sqf(p3) + sqf(p4) + sqf(p5) + sqf(p6) + sqf(p7) + sqf(p8) + sqf(p9);
            float xv = (vdiff[0] * vx + vlp[0] * vn) / (vx + vn);
            float vv = vx * vn / (vx + vn);
            // interpolated G-R(B)
            interp[0] = (xh * vv + xv * vh) / (vh + vv);
          }
        }
 
        // copy CFA values
        for(int rr = 0, row_in = rowStart - BORDER_AROUND; rr < last_rr; rr++, row_in++)
        {
          for(int cc = 0, col_in = colStart - BORDER_AROUND; cc < last_cc; cc++, col_in++)
          {
            const int c = FC(rr, cc, filters);
            const gboolean inside = ((row_in >= 0) && (row_in < height) && (col_in >= 0) && (col_in < width));
            float *colc = qix[c] + rr * LMMSE_GRP + cc;
            colc[0] = (inside) ? qix[5][rr * LMMSE_GRP + cc] : 0.0f;
            if(c != 1)
            {
              float *col1   = qix[1] + rr * LMMSE_GRP + cc;
              float *interp = qix[4] + rr * LMMSE_GRP + cc;
              col1[0] = (inside) ? colc[0] + interp[0] : 0.0f;
            }
          }
        }
 
        // bilinear interpolation for R/B
        // interpolate R/B at G location
        for(int rr = 1; rr < last_rr - 1; rr++)
        {
          for(int cc = 1 + (FC(rr, 2, filters) & 1), c = FC(rr, cc + 1, filters); cc < last_cc - 1; cc += 2)
          {
            float *colc = qix[c] + rr * LMMSE_GRP + cc;
            float *col1 = qix[1] + rr * LMMSE_GRP + cc;
            colc[0] = col1[0] + 0.5f * (colc[ -1] - col1[ -1] + colc[ 1] - col1[ 1]);
            c = 2 - c;
            colc = qix[c] + rr * LMMSE_GRP + cc;
            colc[0] = col1[0] + 0.5f * (colc[-w1] - col1[-w1] + colc[w1] - col1[w1]);
            c = 2 - c;
          }
        }
 
        // interpolate R/B at B/R location
        for(int rr = 1; rr < last_rr - 1; rr++)
        {
          for(int cc = 1 + (FC(rr, 1, filters) & 1), c = 2 - FC(rr, cc, filters); cc < last_cc - 1; cc += 2)
          {
            float *colc = qix[c] + rr * LMMSE_GRP + cc;
            float *col1 = qix[1] + rr * LMMSE_GRP + cc;
            colc[0] = col1[0] + 0.25f * (colc[-w1] - col1[-w1] + colc[ -1] - col1[ -1] + colc[  1] - col1[  1] + colc[ w1] - col1[ w1]);
          }
        }
 
        // for the median and refine corrections we need to specify other loop bounds
        // for inner vs outer tiles
        const int ccmin = (tile_horizontal == 0) ? 6 : 0 ;
        const int ccmax = last_cc - ((tile_horizontal == num_horizontal - 1) ? 6 : 0);
        const int rrmin = (tile_vertical == 0) ? 6 : 0 ;
        const int rrmax = last_rr - ((tile_vertical == num_vertical - 1) ? 6 : 0);
 
        // median filter/
        for(int pass = 0; pass < medians; pass++)
        {
          // Apply 3x3 median filter
          // Compute median(R-G) and median(B-G)
          for(int rr = 1; rr < last_rr - 1; rr++)
          {
            for(int c = 0; c < 3; c += 2)
            {
              const int d = c + 3 - (c == 0 ? 0 : 1);
              for(int cc = 1; cc < last_cc - 1; cc++)
              {
                float *corr = qix[d] + rr * LMMSE_GRP + cc;
                float *colc = qix[c] + rr * LMMSE_GRP + cc;
                float *col1 = qix[1] + rr * LMMSE_GRP + cc;
                // Assign 3x3 differential color values
                corr[0] = median9f(colc[-w1-1] - col1[-w1-1],
                                   colc[-w1  ] - col1[-w1  ],
                                   colc[-w1+1] - col1[-w1+1],
                                   colc[   -1] - col1[   -1],
                                   colc[    0] - col1[    0],
                                   colc[    1] - col1[    1],
                                   colc[ w1-1] - col1[ w1-1],
                                   colc[ w1  ] - col1[ w1  ],
                                   colc[ w1+1] - col1[ w1+1]);
              }
            }
          }
 
          // red/blue at GREEN pixel locations & red/blue and green at BLUE/RED pixel locations
          for(int rr = rrmin; rr < rrmax - 1; rr++)
          {
            float *col0 = qix[0] + rr * LMMSE_GRP + ccmin;
            float *col1 = qix[1] + rr * LMMSE_GRP + ccmin;
            float *col2 = qix[2] + rr * LMMSE_GRP + ccmin;
            float *corr3 = qix[3] + rr * LMMSE_GRP + ccmin;
            float *corr4 = qix[4] + rr * LMMSE_GRP + ccmin;
            int c0 = FC(rr, 0, filters);
            int c1 = FC(rr, 1, filters);
 
            if(c0 == 1)
            {
              c1 = 2 - c1;
              const int d = c1 + 3 - (c1 == 0 ? 0 : 1);
              int cc;
              float *col_c1 = qix[c1] + rr * LMMSE_GRP + ccmin;
              float *corr_d = qix[d] + rr * LMMSE_GRP + ccmin;
              for(cc = ccmin; cc < ccmax - 1; cc += 2)
              {
                col0[0] = col1[0] + corr3[0];
                col2[0] = col1[0] + corr4[0];
                col0++;
                col1++;
                col2++;
                corr3++;
                corr4++;
                col_c1++;
                corr_d++;
                col_c1[0] = col1[0] + corr_d[0];
                col1[0] = 0.5f * (col0[0] - corr3[0] + col2[0] - corr4[0]);
                col0++;
                col1++;
                col2++;
                corr3++;
                corr4++;
                col_c1++;
                corr_d++;
              }
 
              if(cc < ccmax)
              { // remaining pixel, only if width is odd
                col0[0] = col1[0] + corr3[0];
                col2[0] = col1[0] + corr4[0];
              }
            }
            else
            {
              c0 = 2 - c0;
              const int d = c0 + 3 - (c0 == 0 ? 0 : 1);
              float *col_c0 = qix[c0] + rr * LMMSE_GRP + ccmin;
              float *corr_d = qix[d] + rr * LMMSE_GRP + ccmin;
              int cc;
              for(cc = ccmin; cc < ccmax - 1; cc += 2)
              {
                col_c0[0] = col1[0] + corr_d[0];
                col1[0] = 0.5f * (col0[0] - corr3[0] + col2[0] - corr4[0]);
                col0++;
                col1++;
                col2++;
                corr3++;
                corr4++;
                col_c0++;
                corr_d++;
                col0[0] = col1[0] + corr3[0];
                col2[0] = col1[0] + corr4[0];
                col0++;
                col1++;
                col2++;
                corr3++;
                corr4++;
                col_c0++;
                corr_d++;
             }
 
              if(cc < ccmax)
              { // remaining pixel, only if width is odd
                col_c0[0] = col1[0] + corr_d[0];
                col1[0] = 0.5f * (col0[0] - corr3[0] + col2[0] - corr4[0]);
              }
            }
          }
        }
 
        // we fill the non-approximated color channels from gamma corrected cfa data
        for(int rrr = 4; rrr < last_rr - 4; rrr++)
        {
          for(int ccc = 4; ccc < last_cc - 4; ccc++)
          {
            const int idx = rrr * LMMSE_GRP + ccc;
            const int c = FC(rrr, ccc, filters);
            qix[c][idx] = qix[5][idx];
          }
        }
 
        // As we have the color channels fully available we can do the refinements here in tiled code
        for(int step = 0; step < refine; step++)
        {
          // Reinforce interpolated green pixels on RED/BLUE pixel locations
          for(int rr = rrmin + 2; rr < rrmax - 2; rr++)
          {
            for(int cc = ccmin + 2 + (FC(rr, 2, filters) & 1), c = FC(rr, cc, filters); cc < ccmax - 2; cc += 2)
            {
              float *rgb1 = qix[1] + rr * LMMSE_GRP + cc;
              float *rgbc = qix[c] + rr * LMMSE_GRP + cc;
 
              const float dL = 1.0f / (1.0f + fabsf(rgbc[ -2] - rgbc[0]) + fabsf(rgb1[ 1] - rgb1[ -1]));
              const float dR = 1.0f / (1.0f + fabsf(rgbc[  2] - rgbc[0]) + fabsf(rgb1[ 1] - rgb1[ -1]));
              const float dU = 1.0f / (1.0f + fabsf(rgbc[-w2] - rgbc[0]) + fabsf(rgb1[w1] - rgb1[-w1]));
              const float dD = 1.0f / (1.0f + fabsf(rgbc[ w2] - rgbc[0]) + fabsf(rgb1[w1] - rgb1[-w1]));
              rgb1[0] = (rgbc[0] + ((rgb1[-1] - rgbc[-1]) * dL + (rgb1[1] - rgbc[1]) * dR + (rgb1[-w1] - rgbc[-w1]) * dU + (rgb1[w1] - rgbc[w1]) * dD ) / (dL + dR + dU + dD));
            }
          }
          // Reinforce interpolated red/blue pixels on GREEN pixel locations
          for(int rr = rrmin + 2; rr < rrmax - 2; rr++)
          {
            for(int cc = ccmin + 2 + (FC(rr, 3, filters) & 1), c = FC(rr, cc + 1, filters); cc < ccmax - 2; cc += 2)
            {
              for(int i = 0; i < 2; c = 2 - c, i++)
              {
                float *rgb1 = qix[1] + rr * LMMSE_GRP + cc;
                float *rgbc = qix[c] + rr * LMMSE_GRP + cc;
 
                const float dL = 1.0f / (1.0f + fabsf(rgb1[ -2] - rgb1[0]) + fabsf(rgbc[ 1] - rgbc[ -1]));
                const float dR = 1.0f / (1.0f + fabsf(rgb1[  2] - rgb1[0]) + fabsf(rgbc[ 1] - rgbc[ -1]));
                const float dU = 1.0f / (1.0f + fabsf(rgb1[-w2] - rgb1[0]) + fabsf(rgbc[w1] - rgbc[-w1]));
                const float dD = 1.0f / (1.0f + fabsf(rgb1[ w2] - rgb1[0]) + fabsf(rgbc[w1] - rgbc[-w1]));
                rgbc[0] = (rgb1[0] - ((rgb1[-1] - rgbc[-1]) * dL + (rgb1[1] - rgbc[1]) * dR + (rgb1[-w1] - rgbc[-w1]) * dU + (rgb1[w1] - rgbc[w1]) * dD ) / (dL + dR + dU + dD));
              }
            }
          }
          // Reinforce integrated red/blue pixels on BLUE/RED pixel locations
          for(int rr = rrmin + 2; rr < rrmax - 2; rr++)
          {
            for(int cc = ccmin + 2 + (FC(rr, 2, filters) & 1), c = 2 - FC(rr, cc, filters); cc < ccmax - 2; cc += 2)
            {
              const int d = 2 - c;
              float *rgb1 = qix[1] + rr * LMMSE_GRP + cc;
              float *rgbc = qix[c] + rr * LMMSE_GRP + cc;
              float *rgbd = qix[d] + rr * LMMSE_GRP + cc;
 
              const float dL = 1.0f / (1.0f + fabsf(rgbd[ -2] - rgbd[0]) + fabsf(rgb1[ 1] - rgb1[ -1]));
              const float dR = 1.0f / (1.0f + fabsf(rgbd[  2] - rgbd[0]) + fabsf(rgb1[ 1] - rgb1[ -1]));
              const float dU = 1.0f / (1.0f + fabsf(rgbd[-w2] - rgbd[0]) + fabsf(rgb1[w1] - rgb1[-w1]));
              const float dD = 1.0f / (1.0f + fabsf(rgbd[ w2] - rgbd[0]) + fabsf(rgb1[w1] - rgb1[-w1]));
              rgbc[0] = (rgb1[0] - ((rgb1[-1] - rgbc[-1]) * dL + (rgb1[1] - rgbc[1]) * dR + (rgb1[-w1] - rgbc[-w1]) * dU + (rgb1[w1] - rgbc[w1]) * dD ) / (dL + dR + dU + dD));
            }
          }
        }
 
        // write result to out
        // For the outermost tiles in all directions we also write the otherwise overlapped area
        const int first_vertical =   rowStart + ((tile_vertical == 0)                    ? 0 : LMMSE_OVERLAP);
        const int last_vertical =    rowEnd   - ((tile_vertical == num_vertical - 1)     ? 0 : LMMSE_OVERLAP);
        const int first_horizontal = colStart + ((tile_horizontal == 0)                  ? 0 : LMMSE_OVERLAP);
        const int last_horizontal =  colEnd   - ((tile_horizontal == num_horizontal - 1) ? 0 : LMMSE_OVERLAP);
        for(int row = first_vertical, rr = row - rowStart + BORDER_AROUND; row < last_vertical; row++, rr++)
        {
          float *dest = out + 4 * (row * width + first_horizontal);
          const int idx = rr * LMMSE_GRP + first_horizontal - colStart + BORDER_AROUND;
          float *col0 = qix[0] + idx;
          float *col1 = qix[1] + idx;
          float *col2 = qix[2] + idx;
          for(int col = first_horizontal; col < last_horizontal; col++, dest +=4, col0++, col1++, col2++)
          {
            dest[0] = scaler * calc_gamma(col0[0], gamma_out);
            dest[1] = scaler * calc_gamma(col1[0], gamma_out);
            dest[2] = scaler * calc_gamma(col2[0], gamma_out);
            dest[3] = 0.0f;
          }
        }
      }
    }
    dt_pixelpipe_cache_free_align(buffer);
  }
}
 
// revert specific aggressive optimizing
#ifdef __GNUC__
  #pragma GCC pop_options
#endif
 
#undef LMMSE_TILESIZE
#undef LMMSE_OVERLAP
#undef BORDER_AROUND
#undef LMMSE_TILEVALID
#undef w1
#undef w2
#undef w3
#undef w4
 
// clang-format off
// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py
// vim: shiftwidth=2 expandtab tabstop=2 cindent
// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified;
// clang-format on