"Fossies" - the Fresh Open Source Software Archive

Member "rawtherapee-5.7/rtengine/rcd_demosaic.cc" (10 Sep 2019, 19877 Bytes) of package /linux/misc/rawtherapee-5.7.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "rcd_demosaic.cc" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 5.6_vs_5.7.

    1 /*
    2  *  This file is part of RawTherapee.
    3  *
    4  *  Copyright (c) 2017-2018 Luis Sanz Rodriguez (luis.sanz.rodriguez(at)gmail(dot)com) and Ingo Weyrich (heckflosse67@gmx.de)
    5  *
    6  *  RawTherapee is free software: you can redistribute it and/or modify
    7  *  it under the terms of the GNU General Public License as published by
    8  *  the Free Software Foundation, either version 3 of the License, or
    9  *  (at your option) any later version.
   10  *
   11  *  RawTherapee is distributed in the hope that it will be useful,
   12  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   13  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   14  *  GNU General Public License for more details.
   15  *
   16  *  You should have received a copy of the GNU General Public License
   17  *  along with RawTherapee.  If not, see <https://www.gnu.org/licenses/>.
   18  */
   19 #include <cmath>
   20 
   21 #include "rawimagesource.h"
   22 #include "rt_math.h"
   23 #include "procparams.h"
   24 #include "../rtgui/multilangmgr.h"
   25 #include "opthelper.h"
   26 #include "StopWatch.h"
   27 
   28 using namespace std;
   29 
   30 namespace rtengine
   31 {
   32 
   33 /**
   34 * RATIO CORRECTED DEMOSAICING
   35 * Luis Sanz Rodriguez (luis.sanz.rodriguez(at)gmail(dot)com)
   36 *
   37 * Release 2.3 @ 171125
   38 *
   39 * Original code from https://github.com/LuisSR/RCD-Demosaicing
   40 * Licensed under the GNU GPL version 3
   41 */
   42 // Tiled version by Ingo Weyrich (heckflosse67@gmx.de)
   43 void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
   44 {
   45     std::unique_ptr<StopWatch> stop;
   46 
   47     if (measure) {
   48         std::cout << "Demosaicing " << W << "x" << H << " image using rcd with " << chunkSize << " tiles per thread" << std::endl;
   49         stop.reset(new StopWatch("rcd demosaic"));
   50     }
   51 
   52     double progress = 0.0;
   53 
   54     if (plistener) {
   55         plistener->setProgressStr(Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), M("TP_RAW_RCD")));
   56         plistener->setProgress(progress);
   57     }
   58     
   59     constexpr int rcdBorder = 9;
   60     constexpr int tileSize = 214;
   61     constexpr int tileSizeN = tileSize - 2 * rcdBorder;
   62     const int numTh = H / (tileSizeN) + ((H % (tileSizeN)) ? 1 : 0);
   63     const int numTw = W / (tileSizeN) + ((W % (tileSizeN)) ? 1 : 0);
   64     constexpr int w1 = tileSize, w2 = 2 * tileSize, w3 = 3 * tileSize, w4 = 4 * tileSize;
   65     //Tolerance to avoid dividing by zero
   66     static constexpr float eps = 1e-5f;
   67     static constexpr float epssq = 1e-10f;
   68 
   69 #ifdef _OPENMP
   70 #pragma omp parallel
   71 #endif
   72 {
   73     int progresscounter = 0;
   74     float *cfa = (float*) calloc(tileSize * tileSize, sizeof *cfa);
   75     float (*rgb)[tileSize * tileSize] = (float (*)[tileSize * tileSize])malloc(3 * sizeof *rgb);
   76     float *VH_Dir = (float*) calloc(tileSize * tileSize, sizeof *VH_Dir);
   77     float *PQ_Dir = (float*) calloc(tileSize * tileSize, sizeof *PQ_Dir);
   78     float *lpf = PQ_Dir; // reuse buffer, they don't overlap in usage
   79 
   80 #ifdef _OPENMP
   81     #pragma omp for schedule(dynamic, chunkSize) collapse(2) nowait
   82 #endif
   83     for(int tr = 0; tr < numTh; ++tr) {
   84         for(int tc = 0; tc < numTw; ++tc) {
   85             const int rowStart = tr * tileSizeN;
   86             const int rowEnd = std::min(rowStart + tileSize, H);
   87             if(rowStart + rcdBorder == rowEnd - rcdBorder) {
   88                 continue;
   89             }
   90             const int colStart = tc * tileSizeN;
   91             const int colEnd = std::min(colStart + tileSize, W);
   92             if(colStart + rcdBorder == colEnd - rcdBorder) {
   93                 continue;
   94             }
   95 
   96             const int tileRows = std::min(rowEnd - rowStart, tileSize);
   97             const int tilecols = std::min(colEnd - colStart, tileSize);
   98 
   99             for (int row = rowStart; row < rowEnd; row++) {
  100                 int indx = (row - rowStart) * tileSize;
  101                 int c0 = FC(row, colStart);
  102                 int c1 = FC(row, colStart + 1);
  103                 int col = colStart;
  104 
  105                 for (; col < colEnd - 1; col+=2, indx+=2) {
  106                     cfa[indx] = rgb[c0][indx] = LIM01(rawData[row][col] / 65535.f);
  107                     cfa[indx + 1] = rgb[c1][indx + 1] = LIM01(rawData[row][col + 1] / 65535.f);
  108                 }
  109                 if(col < colEnd) {
  110                     cfa[indx] = rgb[c0][indx] = LIM01(rawData[row][col] / 65535.f);
  111                 }
  112             }
  113 
  114             /**
  115             * STEP 1: Find cardinal and diagonal interpolation directions
  116             */
  117 
  118             for (int row = 4; row < tileRows - 4; row++) {
  119                 for (int col = 4, indx = row * tileSize + col; col < tilecols - 4; col++, indx++) {
  120                     const float cfai = cfa[indx];
  121                     //Calculate h/v local discrimination
  122                     float V_Stat = max(epssq, -18.f * cfai * (cfa[indx - w1] + cfa[indx + w1] + 2.f * (cfa[indx - w2] + cfa[indx + w2]) - cfa[indx - w3] - cfa[indx + w3]) - 2.f * cfai * (cfa[indx - w4] + cfa[indx + w4] - 19.f * cfai) - cfa[indx - w1] * (70.f * cfa[indx + w1] + 12.f * cfa[indx - w2] - 24.f * cfa[indx + w2] + 38.f * cfa[indx - w3] - 16.f * cfa[indx + w3] - 12.f * cfa[indx - w4] + 6.f * cfa[indx + w4] - 46.f * cfa[indx - w1]) + cfa[indx + w1] * (24.f * cfa[indx - w2] - 12.f * cfa[indx + w2] + 16.f * cfa[indx - w3] - 38.f * cfa[indx + w3] - 6.f * cfa[indx - w4] + 12.f * cfa[indx + w4] + 46.f * cfa[indx + w1]) + cfa[indx - w2] * (14.f * cfa[indx + w2] - 12.f * cfa[indx + w3] - 2.f * cfa[indx - w4] + 2.f * cfa[indx + w4] + 11.f * cfa[indx - w2]) + cfa[indx + w2] * (-12.f * cfa[indx - w3] + 2.f * (cfa[indx - w4] - cfa[indx + w4]) + 11.f * cfa[indx + w2]) + cfa[indx - w3] * (2.f * cfa[indx + w3] - 6.f * cfa[indx - w4] + 10.f * cfa[indx - w3]) + cfa[indx + w3] * (-6.f * cfa[indx + w4] + 10.f * cfa[indx + w3]) + cfa[indx - w4] * cfa[indx - w4] + cfa[indx + w4] * cfa[indx + w4]);
  123                     float H_Stat = max(epssq, -18.f * cfai * (cfa[indx -  1] + cfa[indx +  1] + 2.f * (cfa[indx -  2] + cfa[indx +  2]) - cfa[indx -  3] - cfa[indx +  3]) - 2.f * cfai * (cfa[indx -  4] + cfa[indx +  4] - 19.f * cfai) - cfa[indx -  1] * (70.f * cfa[indx +  1] + 12.f * cfa[indx -  2] - 24.f * cfa[indx +  2] + 38.f * cfa[indx -  3] - 16.f * cfa[indx +  3] - 12.f * cfa[indx -  4] + 6.f * cfa[indx +  4] - 46.f * cfa[indx -  1]) + cfa[indx +  1] * (24.f * cfa[indx -  2] - 12.f * cfa[indx +  2] + 16.f * cfa[indx -  3] - 38.f * cfa[indx +  3] - 6.f * cfa[indx -  4] + 12.f * cfa[indx +  4] + 46.f * cfa[indx +  1]) + cfa[indx -  2] * (14.f * cfa[indx +  2] - 12.f * cfa[indx +  3] - 2.f * cfa[indx -  4] + 2.f * cfa[indx +  4] + 11.f * cfa[indx -  2]) + cfa[indx +  2] * (-12.f * cfa[indx -  3] + 2.f * (cfa[indx -  4] - cfa[indx +  4]) + 11.f * cfa[indx +  2]) + cfa[indx -  3] * (2.f * cfa[indx +  3] - 6.f * cfa[indx -  4] + 10.f * cfa[indx -  3]) + cfa[indx +  3] * (-6.f * cfa[indx +  4] + 10.f * cfa[indx +  3]) + cfa[indx -  4] * cfa[indx -  4] + cfa[indx +  4] * cfa[indx +  4]);
  124 
  125                     VH_Dir[indx] = V_Stat / (V_Stat + H_Stat);
  126                 }
  127             }
  128 
  129             /**
  130             * STEP 2: Calculate the low pass filter
  131             */
  132             // Step 2.1: Low pass filter incorporating green, red and blue local samples from the raw data
  133 
  134             for (int row = 2; row < tileRows - 2; row++) {
  135                 for (int col = 2 + (FC(row, 0) & 1), indx = row * tileSize + col; col < tilecols - 2; col += 2, indx += 2) {
  136                     lpf[indx>>1] = 0.25f * cfa[indx] + 0.125f * (cfa[indx - w1] + cfa[indx + w1] + cfa[indx - 1] + cfa[indx + 1]) + 0.0625f * (cfa[indx - w1 - 1] + cfa[indx - w1 + 1] + cfa[indx + w1 - 1] + cfa[indx + w1 + 1]);
  137                 }
  138             }
  139 
  140             /**
  141             * STEP 3: Populate the green channel
  142             */
  143             // Step 3.1: Populate the green channel at blue and red CFA positions
  144            for (int row = 4; row < tileRows - 4; row++) {
  145                for (int col = 4 + (FC(row, 0) & 1), indx = row * tileSize + col; col < tilecols - 4; col += 2, indx += 2) {
  146 
  147                     // Refined vertical and horizontal local discrimination
  148                     float VH_Central_Value = VH_Dir[indx];
  149                     float VH_Neighbourhood_Value = 0.25f * ((VH_Dir[indx - w1 - 1] + VH_Dir[indx - w1 + 1]) + (VH_Dir[indx + w1 - 1] + VH_Dir[indx + w1 + 1]));
  150 
  151                     float VH_Disc = std::fabs(0.5f - VH_Central_Value) < std::fabs(0.5f - VH_Neighbourhood_Value) ? VH_Neighbourhood_Value : VH_Central_Value;
  152 
  153                     // Cardinal gradients
  154                     float N_Grad = eps + std::fabs(cfa[indx - w1] - cfa[indx + w1]) + std::fabs(cfa[indx] - cfa[indx - w2]) + std::fabs(cfa[indx - w1] - cfa[indx - w3]) + std::fabs(cfa[indx - w2] - cfa[indx - w4]);
  155                     float S_Grad = eps + std::fabs(cfa[indx - w1] - cfa[indx + w1]) + std::fabs(cfa[indx] - cfa[indx + w2]) + std::fabs(cfa[indx + w1] - cfa[indx + w3]) + std::fabs(cfa[indx + w2] - cfa[indx + w4]);
  156                     float W_Grad = eps + std::fabs(cfa[indx -  1] - cfa[indx +  1]) + std::fabs(cfa[indx] - cfa[indx -  2]) + std::fabs(cfa[indx -  1] - cfa[indx -  3]) + std::fabs(cfa[indx -  2] - cfa[indx -  4]);
  157                     float E_Grad = eps + std::fabs(cfa[indx -  1] - cfa[indx +  1]) + std::fabs(cfa[indx] - cfa[indx +  2]) + std::fabs(cfa[indx +  1] - cfa[indx +  3]) + std::fabs(cfa[indx +  2] - cfa[indx +  4]);
  158 
  159                     // Cardinal pixel estimations
  160                     float N_Est = cfa[indx - w1] * (1.f + (lpf[indx>>1] - lpf[(indx - w2)>>1]) / (eps + lpf[indx>>1] + lpf[(indx - w2)>>1]));
  161                     float S_Est = cfa[indx + w1] * (1.f + (lpf[indx>>1] - lpf[(indx + w2)>>1]) / (eps + lpf[indx>>1] + lpf[(indx + w2)>>1]));
  162                     float W_Est = cfa[indx -  1] * (1.f + (lpf[indx>>1] - lpf[(indx -  2)>>1]) / (eps + lpf[indx>>1] + lpf[(indx -  2)>>1]));
  163                     float E_Est = cfa[indx +  1] * (1.f + (lpf[indx>>1] - lpf[(indx +  2)>>1]) / (eps + lpf[indx>>1] + lpf[(indx +  2)>>1]));
  164 
  165                     // Vertical and horizontal estimations
  166                     float V_Est = (S_Grad * N_Est + N_Grad * S_Est) / (N_Grad + S_Grad);
  167                     float H_Est = (W_Grad * E_Est + E_Grad * W_Est) / (E_Grad + W_Grad);
  168 
  169                     // G@B and G@R interpolation
  170                     rgb[1][indx] = VH_Disc * H_Est + (1.f - VH_Disc) * V_Est;
  171 
  172                 }
  173             }
  174             /**
  175             * STEP 4: Populate the red and blue channels
  176             */
  177 
  178             // Step 4.1: Calculate P/Q diagonal local discrimination
  179             for (int row = rcdBorder - 4; row < tileRows - rcdBorder + 4; row++) {
  180                 for (int col = rcdBorder - 4 + (FC(row, rcdBorder) & 1), indx = row * tileSize + col; col < tilecols - rcdBorder + 4; col += 2, indx += 2) {
  181                     const float cfai = cfa[indx];
  182 
  183                     float P_Stat = max(epssq, - 18.f * cfai * (cfa[indx - w1 - 1] + cfa[indx + w1 + 1] + 2.f * (cfa[indx - w2 - 2] + cfa[indx + w2 + 2]) - cfa[indx - w3 - 3] - cfa[indx + w3 + 3]) - 2.f * cfai * (cfa[indx - w4 - 4] + cfa[indx + w4 + 4] - 19.f * cfai) - cfa[indx - w1 - 1] * (70.f * cfa[indx + w1 + 1] - 12.f * cfa[indx - w2 - 2] + 24.f * cfa[indx + w2 + 2] - 38.f * cfa[indx - w3 - 3] + 16.f * cfa[indx + w3 + 3] + 12.f * cfa[indx - w4 - 4] - 6.f * cfa[indx + w4 + 4] + 46.f * cfa[indx - w1 - 1]) + cfa[indx + w1 + 1] * (24.f * cfa[indx - w2 - 2] - 12.f * cfa[indx + w2 + 2] + 16.f * cfa[indx - w3 - 3] - 38.f * cfa[indx + w3 + 3] - 6.f * cfa[indx - w4 - 4] + 12.f * cfa[indx + w4 + 4] + 46.f * cfa[indx + w1 + 1]) + cfa[indx - w2 - 2] * (14.f * cfa[indx + w2 + 2] - 12.f * cfa[indx + w3 + 3] - 2.f * (cfa[indx - w4 - 4] - cfa[indx + w4 + 4]) + 11.f * cfa[indx - w2 - 2]) - cfa[indx + w2 + 2] * (12.f * cfa[indx - w3 - 3] + 2.f * (cfa[indx - w4 - 4] - cfa[indx + w4 + 4]) + 11.f * cfa[indx + w2 + 2]) + cfa[indx - w3 - 3] * (2.f * cfa[indx + w3 + 3] - 6.f * cfa[indx - w4 - 4] + 10.f * cfa[indx - w3 - 3]) - cfa[indx + w3 + 3] * (6.f * cfa[indx + w4 + 4] + 10.f * cfa[indx + w3 + 3]) + cfa[indx - w4 - 4] * cfa[indx - w4 - 4] + cfa[indx + w4 + 4] * cfa[indx + w4 + 4]);
  184                     float Q_Stat = max(epssq, - 18.f * cfai * (cfa[indx + w1 - 1] + cfa[indx - w1 + 1] + 2.f * (cfa[indx + w2 - 2] + cfa[indx - w2 + 2]) - cfa[indx + w3 - 3] - cfa[indx - w3 + 3]) - 2.f * cfai * (cfa[indx + w4 - 4] + cfa[indx - w4 + 4] - 19.f * cfai) - cfa[indx + w1 - 1] * (70.f * cfa[indx - w1 + 1] - 12.f * cfa[indx + w2 - 2] + 24.f * cfa[indx - w2 + 2] - 38.f * cfa[indx + w3 - 3] + 16.f * cfa[indx - w3 + 3] + 12.f * cfa[indx + w4 - 4] - 6.f * cfa[indx - w4 + 4] + 46.f * cfa[indx + w1 - 1]) + cfa[indx - w1 + 1] * (24.f * cfa[indx + w2 - 2] - 12.f * cfa[indx - w2 + 2] + 16.f * cfa[indx + w3 - 3] - 38.f * cfa[indx - w3 + 3] - 6.f * cfa[indx + w4 - 4] + 12.f * cfa[indx - w4 + 4] + 46.f * cfa[indx - w1 + 1]) + cfa[indx + w2 - 2] * (14.f * cfa[indx - w2 + 2] - 12.f * cfa[indx - w3 + 3] - 2.f * (cfa[indx + w4 - 4] - cfa[indx - w4 + 4]) + 11.f * cfa[indx + w2 - 2]) - cfa[indx - w2 + 2] * (12.f * cfa[indx + w3 - 3] + 2.f * (cfa[indx + w4 - 4] - cfa[indx - w4 + 4]) + 11.f * cfa[indx - w2 + 2]) + cfa[indx + w3 - 3] * (2.f * cfa[indx - w3 + 3] - 6.f * cfa[indx + w4 - 4] + 10.f * cfa[indx + w3 - 3]) - cfa[indx - w3 + 3] * (6.f * cfa[indx - w4 + 4] + 10.f * cfa[indx - w3 + 3]) + cfa[indx + w4 - 4] * cfa[indx + w4 - 4] + cfa[indx - w4 + 4] * cfa[indx - w4 + 4]);
  185 
  186                     PQ_Dir[indx] = P_Stat / (P_Stat + Q_Stat);
  187 
  188                 }
  189             }
  190 
  191             // Step 4.2: Populate the red and blue channels at blue and red CFA positions
  192             for (int row = rcdBorder - 3; row < tileRows - rcdBorder + 3; row++) {
  193                 for (int col = rcdBorder - 3 + (FC(row, rcdBorder - 1) & 1), indx = row * tileSize + col, c = 2 - FC(row, col); col < tilecols - rcdBorder + 3; col += 2, indx += 2) {
  194 
  195                     // Refined P/Q diagonal local discrimination
  196                     float PQ_Central_Value   = PQ_Dir[indx];
  197                     float PQ_Neighbourhood_Value = 0.25f * (PQ_Dir[indx - w1 - 1] + PQ_Dir[indx - w1 + 1] + PQ_Dir[indx + w1 - 1] + PQ_Dir[indx + w1 + 1]);
  198 
  199                     float PQ_Disc = (std::fabs(0.5f - PQ_Central_Value) < std::fabs(0.5f - PQ_Neighbourhood_Value)) ? PQ_Neighbourhood_Value : PQ_Central_Value;
  200 
  201                     // Diagonal gradients
  202                     float NW_Grad = eps + std::fabs(rgb[c][indx - w1 - 1] - rgb[c][indx + w1 + 1]) + std::fabs(rgb[c][indx - w1 - 1] - rgb[c][indx - w3 - 3]) + std::fabs(rgb[1][indx] - rgb[1][indx - w2 - 2]);
  203                     float NE_Grad = eps + std::fabs(rgb[c][indx - w1 + 1] - rgb[c][indx + w1 - 1]) + std::fabs(rgb[c][indx - w1 + 1] - rgb[c][indx - w3 + 3]) + std::fabs(rgb[1][indx] - rgb[1][indx - w2 + 2]);
  204                     float SW_Grad = eps + std::fabs(rgb[c][indx - w1 + 1] - rgb[c][indx + w1 - 1]) + std::fabs(rgb[c][indx + w1 - 1] - rgb[c][indx + w3 - 3]) + std::fabs(rgb[1][indx] - rgb[1][indx + w2 - 2]);
  205                     float SE_Grad = eps + std::fabs(rgb[c][indx - w1 - 1] - rgb[c][indx + w1 + 1]) + std::fabs(rgb[c][indx + w1 + 1] - rgb[c][indx + w3 + 3]) + std::fabs(rgb[1][indx] - rgb[1][indx + w2 + 2]);
  206 
  207                     // Diagonal colour differences
  208                     float NW_Est = rgb[c][indx - w1 - 1] - rgb[1][indx - w1 - 1];
  209                     float NE_Est = rgb[c][indx - w1 + 1] - rgb[1][indx - w1 + 1];
  210                     float SW_Est = rgb[c][indx + w1 - 1] - rgb[1][indx + w1 - 1];
  211                     float SE_Est = rgb[c][indx + w1 + 1] - rgb[1][indx + w1 + 1];
  212 
  213                     // P/Q estimations
  214                     float P_Est = (NW_Grad * SE_Est + SE_Grad * NW_Est) / (NW_Grad + SE_Grad);
  215                     float Q_Est = (NE_Grad * SW_Est + SW_Grad * NE_Est) / (NE_Grad + SW_Grad);
  216 
  217                     // R@B and B@R interpolation
  218                     rgb[c][indx] = rgb[1][indx] + (1.f - PQ_Disc) * P_Est + PQ_Disc * Q_Est;
  219 
  220                 }
  221             }
  222 
  223             // Step 4.3: Populate the red and blue channels at green CFA positions
  224             for (int row = rcdBorder; row < tileRows - rcdBorder; row++) {
  225                 for (int col = rcdBorder + (FC(row, rcdBorder - 1) & 1), indx = row * tileSize + col; col < tilecols - rcdBorder; col += 2, indx += 2) {
  226 
  227                     // Refined vertical and horizontal local discrimination
  228                     float VH_Central_Value   = VH_Dir[indx];
  229                     float VH_Neighbourhood_Value = 0.25f * ((VH_Dir[indx - w1 - 1] + VH_Dir[indx - w1 + 1]) + (VH_Dir[indx + w1 - 1] + VH_Dir[indx + w1 + 1]));
  230 
  231                     float VH_Disc = (std::fabs(0.5f - VH_Central_Value) < std::fabs(0.5f - VH_Neighbourhood_Value)) ? VH_Neighbourhood_Value : VH_Central_Value;
  232                     float rgb1 = rgb[1][indx];
  233                     float N1 = eps + std::fabs(rgb1 - rgb[1][indx - w2]);
  234                     float S1 = eps + std::fabs(rgb1 - rgb[1][indx + w2]);
  235                     float W1 = eps + std::fabs(rgb1 - rgb[1][indx -  2]);
  236                     float E1 = eps + std::fabs(rgb1 - rgb[1][indx +  2]);
  237 
  238                     float rgb1mw1 = rgb[1][indx - w1];
  239                     float rgb1pw1 = rgb[1][indx + w1];
  240                     float rgb1m1 = rgb[1][indx - 1];
  241                     float rgb1p1 = rgb[1][indx + 1];
  242                     for (int c = 0; c <= 2; c += 2) {
  243                         // Cardinal gradients
  244                         float SNabs = std::fabs(rgb[c][indx - w1] - rgb[c][indx + w1]);
  245                         float EWabs = std::fabs(rgb[c][indx -  1] - rgb[c][indx +  1]);
  246                         float N_Grad = N1 + SNabs + std::fabs(rgb[c][indx - w1] - rgb[c][indx - w3]);
  247                         float S_Grad = S1 + SNabs + std::fabs(rgb[c][indx + w1] - rgb[c][indx + w3]);
  248                         float W_Grad = W1 + EWabs + std::fabs(rgb[c][indx -  1] - rgb[c][indx -  3]);
  249                         float E_Grad = E1 + EWabs + std::fabs(rgb[c][indx +  1] - rgb[c][indx +  3]);
  250 
  251                         // Cardinal colour differences
  252                         float N_Est = rgb[c][indx - w1] - rgb1mw1;
  253                         float S_Est = rgb[c][indx + w1] - rgb1pw1;
  254                         float W_Est = rgb[c][indx -  1] - rgb1m1;
  255                         float E_Est = rgb[c][indx +  1] - rgb1p1;
  256 
  257                         // Vertical and horizontal estimations
  258                         float V_Est = (N_Grad * S_Est + S_Grad * N_Est) / (N_Grad + S_Grad);
  259                         float H_Est = (E_Grad * W_Est + W_Grad * E_Est) / (E_Grad + W_Grad);
  260 
  261                         // R@G and B@G interpolation
  262                         rgb[c][indx] = rgb1 + (1.f - VH_Disc) * V_Est + VH_Disc * H_Est;
  263 
  264                     }
  265                 }
  266             }
  267 
  268             for (int row = rowStart + rcdBorder; row < rowEnd - rcdBorder; ++row) {
  269                 for (int col = colStart + rcdBorder; col < colEnd - rcdBorder; ++col) {
  270                     int idx = (row - rowStart) * tileSize + col - colStart ;
  271                     red[row][col] = CLIP(rgb[0][idx] * 65535.f);
  272                     green[row][col] = CLIP(rgb[1][idx] * 65535.f);
  273                     blue[row][col] = CLIP(rgb[2][idx] * 65535.f);
  274                 }
  275             }
  276 
  277             if(plistener) {
  278                 progresscounter++;
  279                 if(progresscounter % 32 == 0) {
  280 #ifdef _OPENMP
  281                     #pragma omp critical (rcdprogress)
  282 #endif
  283                     {
  284                         progress += (double)32 * ((tileSizeN) * (tileSizeN)) / (H * W);
  285                         progress = progress > 1.0 ? 1.0 : progress;
  286                         plistener->setProgress(progress);
  287                     }
  288                 }
  289             }
  290 
  291         }
  292     }
  293 
  294     free(cfa);
  295     free(rgb);
  296     free(VH_Dir);
  297     free(PQ_Dir);
  298 }
  299 
  300     border_interpolate2(W, H, rcdBorder, rawData, red, green, blue);
  301 
  302     if (plistener) {
  303         plistener->setProgress(1);
  304     }
  305     // -------------------------------------------------------------------------
  306 }
  307 
  308 } /* namespace */