/*
 *  This file is part of RawTherapee.
 *
 *  Copyright (c) 2017-2018 Luis Sanz Rodriguez (luis.sanz.rodriguez(at)gmail(dot)com) and Ingo Weyrich (heckflosse67@gmx.de)
 *
 *  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 <http://www.gnu.org/licenses/>.
 */
#include <cmath>

#include "rawimagesource.h"
#include "rt_math.h"
#include "../rtgui/multilangmgr.h"
#include "opthelper.h"
#define BENCHMARK
#include "StopWatch.h"

using namespace std;

namespace rtengine
{

/**
* RATIO CORRECTED DEMOSAICING
* Luis Sanz Rodriguez (luis.sanz.rodriguez(at)gmail(dot)com)
*
* Release 2.3 @ 171125
*
* Original code from https://github.com/LuisSR/RCD-Demosaicing
* Licensed under the GNU GPL version 3
*/
// Tiled version by Ingo Weyrich (heckflosse67@gmx.de)
void RawImageSource::rcd_demosaic()
{
    BENCHFUN

    if (plistener) {
        plistener->setProgressStr(Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), "rcd"));
        plistener->setProgress(0);
    }
    
    const int width = W, height = H;
    constexpr int tileBorder = 8;
    constexpr int tileSize = 228;
    constexpr int tileSizeN = tileSize - 2 * tileBorder;
    const int numTh = H / (tileSizeN) + ((H % (tileSizeN)) ? 1 : 0);
    const int numTw = W / (tileSizeN) + ((W % (tileSizeN)) ? 1 : 0);
    constexpr int w1 = tileSize, w2 = 2 * tileSize, w3 = 3 * tileSize, w4 = 4 * tileSize;
    //Tolerance to avoid dividing by zero
    static constexpr float eps = 1e-5f;
    static constexpr float epssq = 1e-10f;

#ifdef _OPENMP
#pragma omp parallel
#endif
{
    float *cfa = (float*) calloc( tileSize * tileSize, sizeof *cfa );
    float (*rgb)[tileSize * tileSize] = (float (*)[tileSize * tileSize])malloc(3 * sizeof *rgb);
    float *VH_Dir = (float*) calloc( tileSize * tileSize, sizeof *VH_Dir );
    float *PQ_Dir = (float*) calloc( tileSize * tileSize, sizeof *PQ_Dir );
    float *lpf = PQ_Dir; // reuse buffer, they don't overlap in usage

#ifdef _OPENMP
    #pragma omp for schedule(dynamic) collapse(2) nowait
#endif
    for(int tr = 0; tr < numTh; ++tr) {
        for(int tc = 0; tc < numTw; ++tc) {
            int rowStart = tr * tileSizeN;
            int rowEnd = std::min(tr * tileSizeN + tileSize, H);
            int colStart = tc * tileSizeN;
            int colEnd = std::min(tc * tileSizeN + tileSize, W);

            for (int row = rowStart; row < rowEnd; row++) {
                int indx = (row - rowStart) * tileSize;
                int c0 = FC(row, colStart);
                int c1 = FC(row, colStart + 1);
                int col = colStart;

                for (; col < colEnd - 1; col+=2, indx+=2) {
                    cfa[indx] = rgb[c0][indx] = LIM01(rawData[row][col] / 65535.f);
                    cfa[indx + 1] = rgb[c1][indx + 1] = LIM01(rawData[row][col + 1] / 65535.f);
                }
                if(col < colEnd) {
                    cfa[indx] = rgb[c0][indx] = LIM01(rawData[row][col] / 65535.f);
                }
            }

            /**
            * STEP 1: Find cardinal and diagonal interpolation directions
            */

            for (int row = 4; row < tileSize - 4; row++ ) {
                for (int col = 4, indx = row * tileSize + col; col < tileSize - 4; col++, indx++ ) {
                    //Calculate h/v local discrimination
                    float V_Stat = max(epssq, - 18.0f  *  cfa[indx] * cfa[indx - w1] - 18.0f * cfa[indx] * cfa[indx + w1] - 36.0f * cfa[indx] * cfa[indx - w2] - 36.0f * cfa[indx] * cfa[indx + w2] + 18.0f * cfa[indx] * cfa[indx - w3] + 18.0f * cfa[indx] * cfa[indx + w3] - 2.0f * cfa[indx] * cfa[indx - w4] - 2.0f * cfa[indx] * cfa[indx + w4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx - w1] * cfa[indx + w1] - 12.0f * cfa[indx - w1] * cfa[indx - w2] + 24.0f * cfa[indx - w1] * cfa[indx + w2] - 38.0f * cfa[indx - w1] * cfa[indx - w3] + 16.0f * cfa[indx - w1] * cfa[indx + w3] + 12.0f * cfa[indx - w1] * cfa[indx - w4] - 6.0f * cfa[indx - w1] * cfa[indx + w4] + 46.0f * cfa[indx - w1] * cfa[indx - w1] + 24.0f * cfa[indx + w1] * cfa[indx - w2] - 12.0f * cfa[indx + w1] * cfa[indx + w2] + 16.0f * cfa[indx + w1] * cfa[indx - w3] - 38.0f * cfa[indx + w1] * cfa[indx + w3] - 6.0f * cfa[indx + w1] * cfa[indx - w4] + 12.0f * cfa[indx + w1] * cfa[indx + w4] + 46.0f * cfa[indx + w1] * cfa[indx + w1] + 14.0f * cfa[indx - w2] * cfa[indx + w2] - 12.0f * cfa[indx - w2] * cfa[indx + w3] - 2.0f * cfa[indx - w2] * cfa[indx - w4] + 2.0f * cfa[indx - w2] * cfa[indx + w4] + 11.0f * cfa[indx - w2] * cfa[indx - w2] - 12.0f * cfa[indx + w2] * cfa[indx - w3] + 2.0f * cfa[indx + w2] * cfa[indx - w4] - 2.0f * cfa[indx + w2] * cfa[indx + w4] + 11.0f * cfa[indx + w2] * cfa[indx + w2] + 2.0f * cfa[indx - w3] * cfa[indx + w3] - 6.0f * cfa[indx - w3] * cfa[indx - w4] + 10.0f * cfa[indx - w3] * cfa[indx - w3] - 6.0f * cfa[indx + w3] * cfa[indx + w4] + 10.0f * cfa[indx + w3] * cfa[indx + w3] + 1.0f * cfa[indx - w4] * cfa[indx - w4] + 1.0f * cfa[indx + w4] * cfa[indx + w4]);

                    float H_Stat = max(epssq, - 18.0f  *  cfa[indx] * cfa[indx -  1] - 18.0f * cfa[indx] * cfa[indx +  1] - 36.0f * cfa[indx] * cfa[indx -  2] - 36.0f * cfa[indx] * cfa[indx +  2] + 18.0f * cfa[indx] * cfa[indx -  3] + 18.0f * cfa[indx] * cfa[indx +  3] - 2.0f * cfa[indx] * cfa[indx -  4] - 2.0f * cfa[indx] * cfa[indx +  4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx -  1] * cfa[indx +  1] - 12.0f * cfa[indx -  1] * cfa[indx -  2] + 24.0f * cfa[indx -  1] * cfa[indx +  2] - 38.0f * cfa[indx -  1] * cfa[indx -  3] + 16.0f * cfa[indx -  1] * cfa[indx +  3] + 12.0f * cfa[indx -  1] * cfa[indx -  4] - 6.0f * cfa[indx -  1] * cfa[indx +  4] + 46.0f * cfa[indx -  1] * cfa[indx -  1] + 24.0f * cfa[indx +  1] * cfa[indx -  2] - 12.0f * cfa[indx +  1] * cfa[indx +  2] + 16.0f * cfa[indx +  1] * cfa[indx -  3] - 38.0f * cfa[indx +  1] * cfa[indx +  3] - 6.0f * cfa[indx +  1] * cfa[indx -  4] + 12.0f * cfa[indx +  1] * cfa[indx +  4] + 46.0f * cfa[indx +  1] * cfa[indx +  1] + 14.0f * cfa[indx -  2] * cfa[indx +  2] - 12.0f * cfa[indx -  2] * cfa[indx +  3] - 2.0f * cfa[indx -  2] * cfa[indx -  4] + 2.0f * cfa[indx -  2] * cfa[indx +  4] + 11.0f * cfa[indx -  2] * cfa[indx -  2] - 12.0f * cfa[indx +  2] * cfa[indx -  3] + 2.0f * cfa[indx +  2] * cfa[indx -  4] - 2.0f * cfa[indx +  2] * cfa[indx +  4] + 11.0f * cfa[indx +  2] * cfa[indx +  2] + 2.0f * cfa[indx -  3] * cfa[indx +  3] - 6.0f * cfa[indx -  3] * cfa[indx -  4] + 10.0f * cfa[indx -  3] * cfa[indx -  3] - 6.0f * cfa[indx +  3] * cfa[indx +  4] + 10.0f * cfa[indx +  3] * cfa[indx +  3] + 1.0f * cfa[indx -  4] * cfa[indx -  4] + 1.0f * cfa[indx +  4] * cfa[indx +  4]);

                    VH_Dir[indx] = V_Stat / (V_Stat + H_Stat);
                }
            }

            /**
            * STEP 2: Calculate the low pass filter
            */
            // Step 2.1: Low pass filter incorporating green, red and blue local samples from the raw data

            for ( int row = 2; row < tileSize - 2; row++ ) {
                for ( int col = 2 + (FC( row, 0 ) & 1), indx = row * tileSize + col; col < tileSize - 2; col += 2, indx += 2 ) {
                    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] );
                }
            }

            /**
            * STEP 3: Populate the green channel
            */
            // Step 3.1: Populate the green channel at blue and red CFA positions
           for ( int row = 4; row < tileSize - 4; row++ ) {
               for ( int col = 4 + (FC( row, 0 ) & 1), indx = row * tileSize + col; col < tileSize - 4; col += 2, indx += 2 ) {

                    // Refined vertical and horizontal local discrimination
                    float VH_Central_Value = VH_Dir[indx];
                    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]));

                    float VH_Disc = std::fabs(0.5f - VH_Central_Value) < std::fabs(0.5f - VH_Neighbourhood_Value) ? VH_Neighbourhood_Value : VH_Central_Value;

                    // Cardinal gradients
                    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] );
                    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] );
                    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] );
                    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] );

                    // Cardinal pixel estimations
                    float N_Est = cfa[indx - w1] * ( 1.f + ( lpf[indx>>1] - lpf[(indx - w2)>>1] ) / ( eps + lpf[indx>>1] + lpf[(indx - w2)>>1] ) );
                    float S_Est = cfa[indx + w1] * ( 1.f + ( lpf[indx>>1] - lpf[(indx + w2)>>1] ) / ( eps + lpf[indx>>1] + lpf[(indx + w2)>>1] ) );
                    float W_Est = cfa[indx -  1] * ( 1.f + ( lpf[indx>>1] - lpf[(indx -  2)>>1] ) / ( eps + lpf[indx>>1] + lpf[(indx -  2)>>1] ) );
                    float E_Est = cfa[indx +  1] * ( 1.f + ( lpf[indx>>1] - lpf[(indx +  2)>>1] ) / ( eps + lpf[indx>>1] + lpf[(indx +  2)>>1] ) );

                    // Vertical and horizontal estimations
                    float V_Est = ( S_Grad * N_Est + N_Grad * S_Est ) / (N_Grad + S_Grad );
                    float H_Est = ( W_Grad * E_Est + E_Grad * W_Est ) / (E_Grad + W_Grad );

                    // G@B and G@R interpolation
                    rgb[1][indx] = VH_Disc * H_Est + ( 1.f - VH_Disc ) * V_Est;

                }
            }
            /**
            * STEP 4: Populate the red and blue channels
            */

            // Step 4.1: Calculate P/Q diagonal local discrimination
            for ( int row = 4; row < tileSize - 4; row++ ) {
                for ( int col = 4 + (FC( row, 0 ) & 1), indx = row * tileSize + col; col < tileSize - 4; col += 2, indx += 2 ) {

                    float P_Stat = max( - 18.f * cfa[indx] * cfa[indx - w1 - 1] - 18.f * cfa[indx] * cfa[indx + w1 + 1] - 36.f * cfa[indx] * cfa[indx - w2 - 2] - 36.f * cfa[indx] * cfa[indx + w2 + 2] + 18.f * cfa[indx] * cfa[indx - w3 - 3] + 18.f * cfa[indx] * cfa[indx + w3 + 3] - 2.f * cfa[indx] * cfa[indx - w4 - 4] - 2.f * cfa[indx] * cfa[indx + w4 + 4] + 38.f * cfa[indx] * cfa[indx] - 70.f * cfa[indx - w1 - 1] * cfa[indx + w1 + 1] - 12.f * cfa[indx - w1 - 1] * cfa[indx - w2 - 2] + 24.f * cfa[indx - w1 - 1] * cfa[indx + w2 + 2] - 38.f * cfa[indx - w1 - 1] * cfa[indx - w3 - 3] + 16.f * cfa[indx - w1 - 1] * cfa[indx + w3 + 3] + 12.f * cfa[indx - w1 - 1] * cfa[indx - w4 - 4] - 6.f * cfa[indx - w1 - 1] * cfa[indx + w4 + 4] + 46.f * cfa[indx - w1 - 1] * cfa[indx - w1 - 1] + 24.f * cfa[indx + w1 + 1] * cfa[indx - w2 - 2] - 12.f * cfa[indx + w1 + 1] * cfa[indx + w2 + 2] + 16.f * cfa[indx + w1 + 1] * cfa[indx - w3 - 3] - 38.f * cfa[indx + w1 + 1] * cfa[indx + w3 + 3] - 6.f * cfa[indx + w1 + 1] * cfa[indx - w4 - 4] + 12.f * cfa[indx + w1 + 1] * cfa[indx + w4 + 4] + 46.f * cfa[indx + w1 + 1] * cfa[indx + w1 + 1] + 14.f * cfa[indx - w2 - 2] * cfa[indx + w2 + 2] - 12.f * cfa[indx - w2 - 2] * cfa[indx + w3 + 3] - 2.f * cfa[indx - w2 - 2] * cfa[indx - w4 - 4] + 2.f * cfa[indx - w2 - 2] * cfa[indx + w4 + 4] + 11.f * cfa[indx - w2 - 2] * cfa[indx - w2 - 2] - 12.f * cfa[indx + w2 + 2] * cfa[indx - w3 - 3] + 2 * cfa[indx + w2 + 2] * cfa[indx - w4 - 4] - 2.f * cfa[indx + w2 + 2] * cfa[indx + w4 + 4] + 11.f * cfa[indx + w2 + 2] * cfa[indx + w2 + 2] + 2.f * cfa[indx - w3 - 3] * cfa[indx + w3 + 3] - 6.f * cfa[indx - w3 - 3] * cfa[indx - w4 - 4] + 10.f * cfa[indx - w3 - 3] * cfa[indx - w3 - 3] - 6.f * cfa[indx + w3 + 3] * cfa[indx + w4 + 4] + 10.f * cfa[indx + w3 + 3] * cfa[indx + w3 + 3] + 1.f * cfa[indx - w4 - 4] * cfa[indx - w4 - 4] + 1.f * cfa[indx + w4 + 4] * cfa[indx + w4 + 4], epssq );
                    float Q_Stat = max( - 18.f * cfa[indx] * cfa[indx + w1 - 1] - 18.f * cfa[indx] * cfa[indx - w1 + 1] - 36.f * cfa[indx] * cfa[indx + w2 - 2] - 36.f * cfa[indx] * cfa[indx - w2 + 2] + 18.f * cfa[indx] * cfa[indx + w3 - 3] + 18.f * cfa[indx] * cfa[indx - w3 + 3] - 2.f * cfa[indx] * cfa[indx + w4 - 4] - 2.f * cfa[indx] * cfa[indx - w4 + 4] + 38.f * cfa[indx] * cfa[indx] - 70.f * cfa[indx + w1 - 1] * cfa[indx - w1 + 1] - 12.f * cfa[indx + w1 - 1] * cfa[indx + w2 - 2] + 24.f * cfa[indx + w1 - 1] * cfa[indx - w2 + 2] - 38.f * cfa[indx + w1 - 1] * cfa[indx + w3 - 3] + 16.f * cfa[indx + w1 - 1] * cfa[indx - w3 + 3] + 12.f * cfa[indx + w1 - 1] * cfa[indx + w4 - 4] - 6.f * cfa[indx + w1 - 1] * cfa[indx - w4 + 4] + 46.f * cfa[indx + w1 - 1] * cfa[indx + w1 - 1] + 24.f * cfa[indx - w1 + 1] * cfa[indx + w2 - 2] - 12.f * cfa[indx - w1 + 1] * cfa[indx - w2 + 2] + 16.f * cfa[indx - w1 + 1] * cfa[indx + w3 - 3] - 38.f * cfa[indx - w1 + 1] * cfa[indx - w3 + 3] - 6.f * cfa[indx - w1 + 1] * cfa[indx + w4 - 4] + 12.f * cfa[indx - w1 + 1] * cfa[indx - w4 + 4] + 46.f * cfa[indx - w1 + 1] * cfa[indx - w1 + 1] + 14.f * cfa[indx + w2 - 2] * cfa[indx - w2 + 2] - 12.f * cfa[indx + w2 - 2] * cfa[indx - w3 + 3] - 2.f * cfa[indx + w2 - 2] * cfa[indx + w4 - 4] + 2.f * cfa[indx + w2 - 2] * cfa[indx - w4 + 4] + 11.f * cfa[indx + w2 - 2] * cfa[indx + w2 - 2] - 12.f * cfa[indx - w2 + 2] * cfa[indx + w3 - 3] + 2 * cfa[indx - w2 + 2] * cfa[indx + w4 - 4] - 2.f * cfa[indx - w2 + 2] * cfa[indx - w4 + 4] + 11.f * cfa[indx - w2 + 2] * cfa[indx - w2 + 2] + 2.f * cfa[indx + w3 - 3] * cfa[indx - w3 + 3] - 6.f * cfa[indx + w3 - 3] * cfa[indx + w4 - 4] + 10.f * cfa[indx + w3 - 3] * cfa[indx + w3 - 3] - 6.f * cfa[indx - w3 + 3] * cfa[indx - w4 + 4] + 10.f * cfa[indx - w3 + 3] * cfa[indx - w3 + 3] + 1.f * cfa[indx + w4 - 4] * cfa[indx + w4 - 4] + 1.f * cfa[indx - w4 + 4] * cfa[indx - w4 + 4], epssq );

                    PQ_Dir[indx] = P_Stat / ( P_Stat + Q_Stat );

                }
            }

            // Step 4.2: Populate the red and blue channels at blue and red CFA positions
            for ( int row = 4; row < tileSize - 4; row++ ) {
                for ( int col = 4 + (FC( row, 0 ) & 1), indx = row * tileSize + col, c = 2 - FC( row, col ); col < tileSize - 4; col += 2, indx += 2 ) {

                    // Refined P/Q diagonal local discrimination
                    float PQ_Central_Value   = PQ_Dir[indx];
                    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] );

                    float PQ_Disc = ( std::fabs( 0.5f - PQ_Central_Value ) < std::fabs( 0.5f - PQ_Neighbourhood_Value ) ) ? PQ_Neighbourhood_Value : PQ_Central_Value;

                    // Diagonal gradients
                    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] );
                    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] );
                    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] );
                    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] );

                    // Diagonal colour differences
                    float NW_Est = rgb[c][indx - w1 - 1] - rgb[1][indx - w1 - 1];
                    float NE_Est = rgb[c][indx - w1 + 1] - rgb[1][indx - w1 + 1];
                    float SW_Est = rgb[c][indx + w1 - 1] - rgb[1][indx + w1 - 1];
                    float SE_Est = rgb[c][indx + w1 + 1] - rgb[1][indx + w1 + 1];

                    // P/Q estimations
                    float P_Est = ( NW_Grad * SE_Est + SE_Grad * NW_Est ) / (NW_Grad + SE_Grad );
                    float Q_Est = ( NE_Grad * SW_Est + SW_Grad * NE_Est ) / (NE_Grad + SW_Grad );

                    // R@B and B@R interpolation
                    rgb[c][indx] = rgb[1][indx] + ( 1.f - PQ_Disc ) * P_Est + PQ_Disc * Q_Est;

                }
            }

            // Step 4.3: Populate the red and blue channels at green CFA positions
            for ( int row = 4; row < tileSize - 4; row++ ) {
                for ( int col = 4 + (FC( row, 1 ) & 1), indx = row * tileSize + col; col < tileSize - 4; col += 2, indx += 2 ) {

                    // Refined vertical and horizontal local discrimination
                    float VH_Central_Value   = VH_Dir[indx];
                    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]) );

                    float VH_Disc = ( std::fabs( 0.5f - VH_Central_Value ) < std::fabs( 0.5f - VH_Neighbourhood_Value ) ) ? VH_Neighbourhood_Value : VH_Central_Value;
                    float N1 = eps + std::fabs( rgb[1][indx] - rgb[1][indx - w2] );
                    float S1 = eps + std::fabs( rgb[1][indx] - rgb[1][indx + w2] );
                    float W1 = eps + std::fabs( rgb[1][indx] - rgb[1][indx -  2] );
                    float E1 = eps + std::fabs( rgb[1][indx] - rgb[1][indx +  2] );
                    for ( int c = 0; c <= 2; c += 2 ) {


                        // Cardinal gradients
                        float N_Grad = N1 + std::fabs( rgb[c][indx - w1] - rgb[c][indx + w1] ) + std::fabs( rgb[c][indx - w1] - rgb[c][indx - w3] );
                        float S_Grad = S1 + std::fabs( rgb[c][indx - w1] - rgb[c][indx + w1] ) + std::fabs( rgb[c][indx + w1] - rgb[c][indx + w3] );
                        float W_Grad = W1 + std::fabs( rgb[c][indx -  1] - rgb[c][indx +  1] ) + std::fabs( rgb[c][indx -  1] - rgb[c][indx -  3] );
                        float E_Grad = E1 + std::fabs( rgb[c][indx -  1] - rgb[c][indx +  1] ) + std::fabs( rgb[c][indx +  1] - rgb[c][indx +  3] );

                        // Cardinal colour differences
                        float N_Est = rgb[c][indx - w1] - rgb[1][indx - w1];
                        float S_Est = rgb[c][indx + w1] - rgb[1][indx + w1];
                        float W_Est = rgb[c][indx -  1] - rgb[1][indx -  1];
                        float E_Est = rgb[c][indx +  1] - rgb[1][indx +  1];

                        // Vertical and horizontal estimations
                        float V_Est = ( N_Grad * S_Est + S_Grad * N_Est ) / (N_Grad + S_Grad );
                        float H_Est = ( E_Grad * W_Est + W_Grad * E_Est ) / (E_Grad + W_Grad );

                        // R@G and B@G interpolation
                        rgb[c][indx] = rgb[1][indx] + ( 1.f - VH_Disc ) * V_Est + VH_Disc * H_Est;

                    }
                }
            }

            for (int row = rowStart + tileBorder; row < rowEnd - tileBorder; ++row) {
                for (int col = colStart + tileBorder; col < colEnd - tileBorder; ++col) {
                    int idx = (row - rowStart) * tileSize + col - colStart ;
                    red[row][col] = CLIP(rgb[0][idx] * 65535.f);
                    green[row][col] = CLIP(rgb[1][idx] * 65535.f);
                    blue[row][col] = CLIP(rgb[2][idx] * 65535.f);
                }
            }
        }
    }

    free(cfa);
    free(rgb);
    free(VH_Dir);
    free(PQ_Dir);
}

    border_interpolate2(width, height, 8);

    if (plistener) {
        plistener->setProgress(1);
    }
    // -------------------------------------------------------------------------
}

} /* namespace */