diff --git a/rtdata/languages/default b/rtdata/languages/default
index 8320a351b..ba9a11466 100644
--- a/rtdata/languages/default
+++ b/rtdata/languages/default
@@ -1788,6 +1788,7 @@ TP_RAW_PIXELSHIFTSMOOTH_TOOLTIP;Smooth transitions between areas with motion and
 TP_RAW_PIXELSHIFTSTDDEVFACTORBLUE;StdDev factor Blue
 TP_RAW_PIXELSHIFTSTDDEVFACTORGREEN;StdDev factor Green
 TP_RAW_PIXELSHIFTSTDDEVFACTORRED;StdDev factor Red
+TP_RAW_RCD;RCD
 TP_RAW_SENSOR_BAYER_LABEL;Sensor with Bayer Matrix
 TP_RAW_SENSOR_XTRANS_DMETHOD_TOOLTIP;3-pass gives best results (recommended for low ISO images).\n1-pass is almost undistinguishable from 3-pass for high ISO images and is faster.
 TP_RAW_SENSOR_XTRANS_LABEL;Sensor with X-Trans Matrix
diff --git a/rtengine/demosaic_algos.cc b/rtengine/demosaic_algos.cc
index 88d369969..61258fca2 100644
--- a/rtengine/demosaic_algos.cc
+++ b/rtengine/demosaic_algos.cc
@@ -3927,6 +3927,286 @@ 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.2 @ 171117
+*/
+void RawImageSource::rcd_demosaic()
+{
+    // RT ---------------------------------------------------------------------
+    if (plistener) {
+        plistener->setProgressStr(Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), "rcd"));
+        plistener->setProgress(0);
+    }
+    
+    int width = W, height = H;
+
+    std::vector<float> cfa(width * height);
+    std::vector<std::array<float, 3>> rgb(width * height);
+
+#ifdef _OPENMP
+    #pragma omp parallel for
+#endif
+    for (int row = 0; row < height; row++) {
+        for (int col = 0, indx = row * width + col; col < width; col++, indx++) {
+            int c = FC(row, col);
+            cfa[indx] = rgb[indx][c] = CLIP(rawData[row][col]) / 65535.f; 
+        }
+    }
+
+    if (plistener) {
+        plistener->setProgress(0.05);
+    }
+    // ------------------------------------------------------------------------
+/* RT    
+    int row, col, indx, c;
+*/    
+    int w1 = width, w2 = 2 * width, w3 = 3 * width, w4 = 4 * width;
+
+    //Tolerance to avoid dividing by zero
+    static const float eps = 1e-5, epssq = 1e-10;
+
+/* RT    
+    //Gradients
+    float N_Grad, E_Grad, W_Grad, S_Grad, NW_Grad, NE_Grad, SW_Grad, SE_Grad;
+
+    //Pixel estimation
+    float N_Est, E_Est, W_Est, S_Est, NW_Est, NE_Est, SW_Est, SE_Est, V_Est, H_Est, P_Est, Q_Est;
+
+    //Directional discrimination
+    //float V_Stat, H_Stat, P_Stat, Q_Stat;
+    float VH_Central_Value, VH_Neighbour_Value, PQ_Central_Value, PQ_Neighbour_Value;
+*/    
+    float ( *VH_Dir ), ( *VH_Disc ), ( *PQ_Dir ), ( *PQ_Disc );
+
+    //Low pass filter
+    float ( *lpf );
+
+
+    /**
+    * STEP 1: Find cardinal and diagonal interpolation directions
+    */
+
+    VH_Dir = ( float ( * ) ) calloc( width * height, sizeof *VH_Dir ); //merror ( VH_Dir, "rcd_demosaicing_171117()" );
+    PQ_Dir = ( float ( * ) ) calloc( width * height, sizeof *PQ_Dir ); //merror ( PQ_Dir, "rcd_demosaicing_171117()" );
+
+#ifdef _OPENMP
+    #pragma omp parallel for
+#endif
+    for (int row = 4; row < height - 4; row++ ) {
+        for (int col = 4, indx = row * width + col; col < width - 4; col++, indx++ ) {
+
+            //Calculate h/v local discrimination
+            float V_Stat = 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 = 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);
+
+            //Calculate m/p local discrimination
+            float P_Stat = epssq - 18.0f * cfa[indx] * cfa[indx - w1 - 1] - 18.0f * cfa[indx] * cfa[indx + w1 + 1] - 36.0f * cfa[indx] * cfa[indx - w2 - 2] - 36.0f * cfa[indx] * cfa[indx + w2 + 2] + 18.0f * cfa[indx] * cfa[indx - w3 - 3] + 18.0f * cfa[indx] * cfa[indx + w3 + 3] - 2.0f * cfa[indx] * cfa[indx - w4 - 4] - 2.0f * cfa[indx] * cfa[indx + w4 + 4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx - w1 - 1] * cfa[indx + w1 + 1] - 12.0f * cfa[indx - w1 - 1] * cfa[indx - w2 - 2] + 24.0f * cfa[indx - w1 - 1] * cfa[indx + w2 + 2] - 38.0f * cfa[indx - w1 - 1] * cfa[indx - w3 - 3] + 16.0f * cfa[indx - w1 - 1] * cfa[indx + w3 + 3] + 12.0f * cfa[indx - w1 - 1] * cfa[indx - w4 - 4] - 6.0f * cfa[indx - w1 - 1] * cfa[indx + w4 + 4] + 46.0f * cfa[indx - w1 - 1] * cfa[indx - w1 - 1] + 24.0f * cfa[indx + w1 + 1] * cfa[indx - w2 - 2] - 12.0f * cfa[indx + w1 + 1] * cfa[indx + w2 + 2] + 16.0f * cfa[indx + w1 + 1] * cfa[indx - w3 - 3] - 38.0f * cfa[indx + w1 + 1] * cfa[indx + w3 + 3] - 6.0f * cfa[indx + w1 + 1] * cfa[indx - w4 - 4] + 12.0f * cfa[indx + w1 + 1] * cfa[indx + w4 + 4] + 46.0f * cfa[indx + w1 + 1] * cfa[indx + w1 + 1] + 14.0f * cfa[indx - w2 - 2] * cfa[indx + w2 + 2] - 12.0f * cfa[indx - w2 - 2] * cfa[indx + w3 + 3] - 2.0f * cfa[indx - w2 - 2] * cfa[indx - w4 - 4] + 2.0f * cfa[indx - w2 - 2] * cfa[indx + w4 + 4] + 11.0f * cfa[indx - w2 - 2] * cfa[indx - w2 - 2] - 12.0f * cfa[indx + w2 + 2] * cfa[indx - w3 - 3] + 2 * cfa[indx + w2 + 2] * cfa[indx - w4 - 4] - 2.0f * cfa[indx + w2 + 2] * cfa[indx + w4 + 4] + 11.0f * cfa[indx + w2 + 2] * cfa[indx + w2 + 2] + 2.0f * cfa[indx - w3 - 3] * cfa[indx + w3 + 3] - 6.0f * cfa[indx - w3 - 3] * cfa[indx - w4 - 4] + 10.0f * cfa[indx - w3 - 3] * cfa[indx - w3 - 3] - 6.0f * cfa[indx + w3 + 3] * cfa[indx + w4 + 4] + 10.0f * cfa[indx + w3 + 3] * cfa[indx + w3 + 3] + 1.0f * cfa[indx - w4 - 4] * cfa[indx - w4 - 4] + 1.0f * cfa[indx + w4 + 4] * cfa[indx + w4 + 4];
+            float Q_Stat = epssq - 18.0f * cfa[indx] * cfa[indx + w1 - 1] - 18.0f * cfa[indx] * cfa[indx - w1 + 1] - 36.0f * cfa[indx] * cfa[indx + w2 - 2] - 36.0f * cfa[indx] * cfa[indx - w2 + 2] + 18.0f * cfa[indx] * cfa[indx + w3 - 3] + 18.0f * cfa[indx] * cfa[indx - w3 + 3] - 2.0f * cfa[indx] * cfa[indx + w4 - 4] - 2.0f * cfa[indx] * cfa[indx - w4 + 4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx + w1 - 1] * cfa[indx - w1 + 1] - 12.0f * cfa[indx + w1 - 1] * cfa[indx + w2 - 2] + 24.0f * cfa[indx + w1 - 1] * cfa[indx - w2 + 2] - 38.0f * cfa[indx + w1 - 1] * cfa[indx + w3 - 3] + 16.0f * cfa[indx + w1 - 1] * cfa[indx - w3 + 3] + 12.0f * cfa[indx + w1 - 1] * cfa[indx + w4 - 4] - 6.0f * cfa[indx + w1 - 1] * cfa[indx - w4 + 4] + 46.0f * cfa[indx + w1 - 1] * cfa[indx + w1 - 1] + 24.0f * cfa[indx - w1 + 1] * cfa[indx + w2 - 2] - 12.0f * cfa[indx - w1 + 1] * cfa[indx - w2 + 2] + 16.0f * cfa[indx - w1 + 1] * cfa[indx + w3 - 3] - 38.0f * cfa[indx - w1 + 1] * cfa[indx - w3 + 3] - 6.0f * cfa[indx - w1 + 1] * cfa[indx + w4 - 4] + 12.0f * cfa[indx - w1 + 1] * cfa[indx - w4 + 4] + 46.0f * cfa[indx - w1 + 1] * cfa[indx - w1 + 1] + 14.0f * cfa[indx + w2 - 2] * cfa[indx - w2 + 2] - 12.0f * cfa[indx + w2 - 2] * cfa[indx - w3 + 3] - 2.0f * cfa[indx + w2 - 2] * cfa[indx + w4 - 4] + 2.0f * cfa[indx + w2 - 2] * cfa[indx - w4 + 4] + 11.0f * cfa[indx + w2 - 2] * cfa[indx + w2 - 2] - 12.0f * cfa[indx - w2 + 2] * cfa[indx + w3 - 3] + 2 * cfa[indx - w2 + 2] * cfa[indx + w4 - 4] - 2.0f * cfa[indx - w2 + 2] * cfa[indx - w4 + 4] + 11.0f * cfa[indx - w2 + 2] * cfa[indx - w2 + 2] + 2.0f * cfa[indx + w3 - 3] * cfa[indx - w3 + 3] - 6.0f * cfa[indx + w3 - 3] * cfa[indx + w4 - 4] + 10.0f * cfa[indx + w3 - 3] * cfa[indx + w3 - 3] - 6.0f * cfa[indx - w3 + 3] * cfa[indx - w4 + 4] + 10.0f * cfa[indx - w3 + 3] * cfa[indx - w3 + 3] + 1.0f * cfa[indx + w4 - 4] * cfa[indx + w4 - 4] + 1.0f * cfa[indx - w4 + 4] * cfa[indx - w4 + 4];
+
+            PQ_Dir[indx] = P_Stat / ( P_Stat + Q_Stat );
+
+        }
+    }
+
+    // RT ---------------------------------------------------------------------
+    if (plistener) {
+        plistener->setProgress(0.2);
+    }
+    // -------------------------------------------------------------------------
+    
+
+    VH_Disc = ( float ( * ) ) calloc( width * height, sizeof *VH_Disc ); //merror ( VH_Disc, "rcd_demosaicing_171117()" );
+    PQ_Disc = ( float ( * ) ) calloc( width * height, sizeof *PQ_Disc ); //merror ( PQ_Disc, "rcd_demosaicing_171117()" );
+
+#ifdef _OPENMP
+    #pragma omp parallel for
+#endif
+    for ( int row = 4; row < height - 4; row++ ) {
+        for ( int col = 4, indx = row * width + col; col < width - 4; col++, indx++ ) {
+
+            //Refined h/v local discrimination
+            float VH_Central_Value   = VH_Dir[indx];
+            float VH_Neighbour_Value = 0.25f * (VH_Dir[indx - w1 - 1] + VH_Dir[indx - w1 + 1] + VH_Dir[indx + w1 - 1] + VH_Dir[indx + w1 + 1]);
+
+            VH_Disc[indx] = ( fabs( 0.5f - VH_Central_Value ) < fabs( 0.5f - VH_Neighbour_Value ) ) ? VH_Neighbour_Value : VH_Central_Value;
+
+            //Refined m/p local discrimination
+            float PQ_Central_Value   = PQ_Dir[indx];
+            float PQ_Neighbour_Value = 0.25f * (PQ_Dir[indx - w1 - 1] + PQ_Dir[indx - w1 + 1] + PQ_Dir[indx + w1 - 1] + PQ_Dir[indx + w1 + 1]);
+
+            PQ_Disc[indx] = ( fabs( 0.5f - PQ_Central_Value ) < fabs( 0.5f - PQ_Neighbour_Value ) ) ? PQ_Neighbour_Value : PQ_Central_Value;
+
+        }
+    }
+
+    free( VH_Dir );
+    free( PQ_Dir );
+
+    // RT ---------------------------------------------------------------------
+    if (plistener) {
+        plistener->setProgress(0.4);
+    }
+    // ------------------------------------------------------------------------
+
+    /**
+    * STEP 2: Calculate the low pass filter
+    */
+
+    lpf = ( float ( * ) ) calloc( width * height, sizeof *lpf ); //merror ( lpf, "rcd_demosaicing_171117()" );
+
+#ifdef _OPENMP
+    #pragma omp parallel for
+#endif    
+    for ( int row = 1; row < height - 1; row++ ) {
+        for ( int col = 1, indx = row * width + col; col < width - 1; col++, indx++ ) {
+
+            //Low pass filter incorporating red and blue local samples
+            lpf[indx] = 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] );
+
+        }
+    }
+
+    // RT ---------------------------------------------------------------------
+    if (plistener) {
+        plistener->setProgress(0.5);
+    }
+    // ------------------------------------------------------------------------
+
+    /**
+    * STEP 3: Populate the green channel
+    */
+#ifdef _OPENMP
+    #pragma omp parallel for
+#endif
+    for ( int row = 4; row < height - 4; row++ ) {
+        for ( int col = 4 + ( FC( row, 0 )&1 ), indx = row * width + col; col < width - 4; col += 2, indx += 2 ) {
+
+            //Cardinal gradients
+            float N_Grad = eps + fabs( cfa[indx - w1] - cfa[indx + w1] ) + fabs( cfa[indx] - cfa[indx - w2] ) + fabs( cfa[indx - w1] - cfa[indx - w3] ) + fabs( cfa[indx - w2] - cfa[indx - w4] );
+            float S_Grad = eps + fabs( cfa[indx + w1] - cfa[indx - w1] ) + fabs( cfa[indx] - cfa[indx + w2] ) + fabs( cfa[indx + w1] - cfa[indx + w3] ) + fabs( cfa[indx + w2] - cfa[indx + w4] );
+            float W_Grad = eps + fabs( cfa[indx -  1] - cfa[indx +  1] ) + fabs( cfa[indx] - cfa[indx -  2] ) + fabs( cfa[indx -  1] - cfa[indx -  3] ) + fabs( cfa[indx -  2] - cfa[indx -  4] );
+            float E_Grad = eps + fabs( cfa[indx +  1] - cfa[indx -  1] ) + fabs( cfa[indx] - cfa[indx +  2] ) + fabs( cfa[indx +  1] - cfa[indx +  3] ) + fabs( cfa[indx +  2] - cfa[indx +  4] );
+
+            //Cardinal pixel estimations
+            float N_Est = cfa[indx - w1] * ( 1.f + ( lpf[indx] - lpf[indx - w2] ) / ( eps + lpf[indx] + lpf[indx - w2] ) );
+            float S_Est = cfa[indx + w1] * ( 1.f + ( lpf[indx] - lpf[indx + w2] ) / ( eps + lpf[indx] + lpf[indx + w2] ) );
+            float W_Est = cfa[indx -  1] * ( 1.f + ( lpf[indx] - lpf[indx -  2] ) / ( eps + lpf[indx] + lpf[indx -  2] ) );
+            float E_Est = cfa[indx +  1] * ( 1.f + ( lpf[indx] - lpf[indx +  2] ) / ( eps + lpf[indx] + lpf[indx +  2] ) );
+
+            //Interpolate G@R & G@B
+            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 );
+
+            rgb[indx][1] = LIM( VH_Disc[indx] * H_Est + ( 1.0f - VH_Disc[indx] ) * V_Est, 0.f, 1.f );
+
+        }
+    }
+
+    free( lpf );
+
+    // RT ---------------------------------------------------------------------
+    if (plistener) {
+        plistener->setProgress(0.7);
+    }
+    // -------------------------------------------------------------------------
+
+    /**
+    * STEP 4: Populate the red and blue channel
+    */
+#ifdef _OPENMP
+    #pragma omp parallel for
+#endif
+    for ( int row = 4; row < height - 4; row++ ) {
+        for ( int col = 4 + ( FC( row, 0 )&1 ), indx = row * width + col, c = 2 - FC( row, col ); col < width - 4; col += 2, indx += 2 ) {
+
+            //Diagonal gradients
+            float NW_Grad = eps + fabs( rgb[indx - w1 - 1][c] - rgb[indx + w1 + 1][c] ) + fabs( rgb[indx - w1 - 1][c] - rgb[indx - w3 - 3][c] ) + fabs( rgb[indx][1] - rgb[indx - w2 - 2][1] );
+            float NE_Grad = eps + fabs( rgb[indx - w1 + 1][c] - rgb[indx + w1 - 1][c] ) + fabs( rgb[indx - w1 + 1][c] - rgb[indx - w3 + 3][c] ) + fabs( rgb[indx][1] - rgb[indx - w2 + 2][1] );
+            float SW_Grad = eps + fabs( rgb[indx + w1 - 1][c] - rgb[indx - w1 + 1][c] ) + fabs( rgb[indx + w1 - 1][c] - rgb[indx + w3 - 3][c] ) + fabs( rgb[indx][1] - rgb[indx + w2 - 2][1] );
+            float SE_Grad = eps + fabs( rgb[indx + w1 + 1][c] - rgb[indx - w1 - 1][c] ) + fabs( rgb[indx + w1 + 1][c] - rgb[indx + w3 + 3][c] ) + fabs( rgb[indx][1] - rgb[indx + w2 + 2][1] );
+
+            //Diagonal colour differences
+            float NW_Est = rgb[indx - w1 - 1][c] - rgb[indx - w1 - 1][1];
+            float NE_Est = rgb[indx - w1 + 1][c] - rgb[indx - w1 + 1][1];
+            float SW_Est = rgb[indx + w1 - 1][c] - rgb[indx + w1 - 1][1];
+            float SE_Est = rgb[indx + w1 + 1][c] - rgb[indx + w1 + 1][1];
+
+            //Interpolate R@B and B@R
+            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 );
+
+            rgb[indx][c] = LIM( rgb[indx][1] + ( 1.0f - PQ_Disc[indx] ) * P_Est + PQ_Disc[indx] * Q_Est, 0.f, 1.f );
+
+        }
+    }
+
+    // RT ---------------------------------------------------------------------
+    if (plistener) {
+        plistener->setProgress(0.825);
+    }
+    // -------------------------------------------------------------------------
+    
+#ifdef _OPENMP
+    #pragma omp parallel for
+#endif
+    for ( int row = 4; row < height - 4; row++ ) {
+        for ( int col = 4 + ( FC( row, 1 )&1 ), indx = row * width + col; col < width - 4; col += 2, indx += 2 ) {
+
+            for ( int c = 0; c <= 2; c += 2 ) {
+
+                //Cardinal gradients
+                float N_Grad = eps + fabs( rgb[indx][1] - rgb[indx - w2][1] ) + fabs( rgb[indx - w1][c] - rgb[indx + w1][c] ) + fabs( rgb[indx - w1][c] - rgb[indx - w3][c] );
+                float S_Grad = eps + fabs( rgb[indx][1] - rgb[indx + w2][1] ) + fabs( rgb[indx + w1][c] - rgb[indx - w1][c] ) + fabs( rgb[indx + w1][c] - rgb[indx + w3][c] );
+                float W_Grad = eps + fabs( rgb[indx][1] - rgb[indx -  2][1] ) + fabs( rgb[indx -  1][c] - rgb[indx +  1][c] ) + fabs( rgb[indx -  1][c] - rgb[indx -  3][c] );
+                float E_Grad = eps + fabs( rgb[indx][1] - rgb[indx +  2][1] ) + fabs( rgb[indx +  1][c] - rgb[indx -  1][c] ) + fabs( rgb[indx +  1][c] - rgb[indx +  3][c] );
+
+                //Cardinal colour differences
+                float N_Est = rgb[indx - w1][c] - rgb[indx - w1][1];
+                float S_Est = rgb[indx + w1][c] - rgb[indx + w1][1];
+                float W_Est = rgb[indx -  1][c] - rgb[indx -  1][1];
+                float E_Est = rgb[indx +  1][c] - rgb[indx +  1][1];
+
+                //Interpolate R@G and B@G
+                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 );
+
+                rgb[indx][c] = LIM( rgb[indx][1] + ( 1.0f - VH_Disc[indx] ) * V_Est + VH_Disc[indx] * H_Est, 0.f, 1.f );
+
+            }
+        }
+    }
+
+    free( PQ_Disc );
+    free( VH_Disc );
+
+    // RT ---------------------------------------------------------------------
+    if (plistener) {
+        plistener->setProgress(0.95);
+    }
+
+#ifdef _OPENMP
+    #pragma omp parallel for
+#endif
+    for (int row = 0; row < height; ++row) {
+        for (int col = 0, idx = row * width + col ; col < width; ++col, ++idx) {
+            red[row][col] = CLIP(rgb[idx][0] * 65535.f);
+            green[row][col] = CLIP(rgb[idx][1] * 65535.f);
+            blue[row][col] = CLIP(rgb[idx][2] * 65535.f);
+        }
+    }
+
+    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/procparams.cc b/rtengine/procparams.cc
index 0ada68a62..0447da364 100644
--- a/rtengine/procparams.cc
+++ b/rtengine/procparams.cc
@@ -53,7 +53,7 @@ const int tl = (int) options.rtSettings.top_left;
 const int bl = (int) options.rtSettings.bot_left;
 
 const char *LensProfParams::methodstring[static_cast<size_t>(LensProfParams::LcMode::LCP) + 1u] = {"none", "lfauto", "lfmanual", "lcp"};
-const char *RAWParams::BayerSensor::methodstring[RAWParams::BayerSensor::numMethods] = {"amaze", "igv", "lmmse", "eahd", "hphd", "vng4", "dcb", "ahd", "fast", "mono", "none", "pixelshift" };
+const char *RAWParams::BayerSensor::methodstring[RAWParams::BayerSensor::numMethods] = {"amaze", "igv", "lmmse", "eahd", "hphd", "vng4", "dcb", "ahd", "rcd", "fast", "mono", "none", "pixelshift" };
 const char *RAWParams::XTransSensor::methodstring[RAWParams::XTransSensor::numMethods] = {"3-pass (best)", "1-pass (medium)", "fast", "mono", "none" };
 
 const char *RAWParams::ff_BlurTypestring[RAWParams::numFlatFileBlurTypes] = {/*"Parametric",*/ "Area Flatfield", "Vertical Flatfield", "Horizontal Flatfield", "V+H Flatfield"};
diff --git a/rtengine/procparams.h b/rtengine/procparams.h
index b3139cdcb..d7ebd1921 100644
--- a/rtengine/procparams.h
+++ b/rtengine/procparams.h
@@ -1262,7 +1262,7 @@ public:
     public:
         //enum eMethod{ eahd,hphd,vng4,dcb,amaze,ahd,IGV_noise,fast,
         //numMethods }; // This MUST be the last enum
-        enum eMethod { amaze, igv, lmmse, eahd, hphd, vng4, dcb, ahd, fast, mono, none, pixelshift,
+        enum eMethod { amaze, igv, lmmse, eahd, hphd, vng4, dcb, ahd, rcd, fast, mono, none, pixelshift,
                        numMethods
                      }; // This MUST be the last enum
         enum ePSMotionCorrection {
diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc
index 5eb7aa7b9..a6b6d4b53 100644
--- a/rtengine/rawimagesource.cc
+++ b/rtengine/rawimagesource.cc
@@ -2049,6 +2049,8 @@ void RawImageSource::demosaic(const RAWParams &raw)
             fast_demosaic();
         } else if (raw.bayersensor.method == RAWParams::BayerSensor::methodstring[RAWParams::BayerSensor::mono] ) {
             nodemosaic(true);
+        } else if (raw.bayersensor.method == RAWParams::BayerSensor::methodstring[RAWParams::BayerSensor::rcd] ) {
+            rcd_demosaic ();
         } else {
             nodemosaic(false);
         }
diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h
index d2ce77fed..4924b955a 100644
--- a/rtengine/rawimagesource.h
+++ b/rtengine/rawimagesource.h
@@ -248,6 +248,7 @@ protected:
     void fast_demosaic();//Emil's code for fast demosaicing
     void dcb_demosaic(int iterations, bool dcb_enhance);
     void ahd_demosaic();
+    void rcd_demosaic();
     void border_interpolate(unsigned int border, float (*image)[4], unsigned int start = 0, unsigned int end = 0);
     void border_interpolate2(int winw, int winh, int lborders);
     void dcb_initTileLimits(int &colMin, int &rowMin, int &colMax, int &rowMax, int x0, int y0, int border);