From bad28bb0add0f6c04bdea910c23be6f076e903d2 Mon Sep 17 00:00:00 2001 From: heckflosse Date: Tue, 27 Feb 2018 19:31:05 +0100 Subject: [PATCH] Made own compilation unit for rcd_demosaic() --- rtengine/CMakeLists.txt | 1 + rtengine/demosaic_algos.cc | 239 -------------------------------- rtengine/rcd_demosaic.cc | 270 +++++++++++++++++++++++++++++++++++++ 3 files changed, 271 insertions(+), 239 deletions(-) create mode 100644 rtengine/rcd_demosaic.cc diff --git a/rtengine/CMakeLists.txt b/rtengine/CMakeLists.txt index 080a76410..230c0957a 100644 --- a/rtengine/CMakeLists.txt +++ b/rtengine/CMakeLists.txt @@ -105,6 +105,7 @@ set(RTENGINESOURCEFILES profilestore.cc rawimage.cc rawimagesource.cc + rcd_demosaic.cc refreshmap.cc rt_algo.cc rtthumbnail.cc diff --git a/rtengine/demosaic_algos.cc b/rtengine/demosaic_algos.cc index 8d0a39401..99d59e9d1 100644 --- a/rtengine/demosaic_algos.cc +++ b/rtengine/demosaic_algos.cc @@ -37,7 +37,6 @@ #include "sleef.c" #include "opthelper.h" #include "median.h" -#define BENCHMARK #include "StopWatch.h" #ifdef _OPENMP #include @@ -3927,244 +3926,6 @@ void RawImageSource::cielab (const float (*rgb)[3], float* l, float* a, float *b } } - -/** -* 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); - } - // ------------------------------------------------------------------------- -} - #define fcol(row,col) xtrans[(row)%6][(col)%6] #define isgreen(row,col) (xtrans[(row)%3][(col)%3]&1) diff --git a/rtengine/rcd_demosaic.cc b/rtengine/rcd_demosaic.cc new file mode 100644 index 000000000..decbe498b --- /dev/null +++ b/rtengine/rcd_demosaic.cc @@ -0,0 +1,270 @@ +/* + * 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 . + */ +#include + +#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 */