EAHD demosaic: further speedup, #4727

rtengine/demosaic_algos.cc

@@ -85,8 +85,8 @@ void RawImageSource::eahd_demosaic ()
     lc22 = (0.019334 * imatrices.rgb_cam[2][0] + 0.119193 * imatrices.rgb_cam[2][1] + 0.950227 * imatrices.rgb_cam[2][2]) ;// / 1.088754;
 
     int maxindex = 3 * 65536; //2*65536 3 = avoid crash 3/2013 J.Desmis
-    cache = new double[maxindex];
+    cache = new float[maxindex];
-    threshold = (int)(0.008856 * MAXVALD);
+    threshold = 0.008856 * MAXVALD;
 
     for (int i = 0; i < maxindex; i++) {
         cache[i] = std::cbrt(double(i) / MAXVALD);
@@ -98,8 +98,8 @@ void RawImageSource::eahd_demosaic ()
     rh (W, 3), gh (W, 4), bh (W, 3),
     rv (W, 3), gv (W, 4), bv (W, 3),
     lLh (W, 3), lah (W, 3), lbh (W, 3),
-    lLv (W, 3), lav (W, 3), lbv (W, 3),
+    lLv (W, 3), lav (W, 3), lbv (W, 3);
-    homh (W, 3), homv (W, 3);
+    JaggedArray<uint16_t> homh (W, 3), homv (W, 3);
 
     // interpolate first two lines
     interpolate_row_g (gh[0], gv[0], 0);
@@ -122,14 +122,15 @@ void RawImageSource::eahd_demosaic ()
         homv[1][j] = 0;
     }
 
-    int dLmaph[9];
+    float dLmaph[9];
-    int dLmapv[9];
+    float dLmapv[9];
-    int dCamaph[9];
+    float dCamaph[9];
-    int dCamapv[9];
+    float dCamapv[9];
-    int dCbmaph[9];
+    float dCbmaph[9];
-    int dCbmapv[9];
+    float dCbmapv[9];
 
     for (int i = 1; i < H - 1; i++) {
+        int mod[3] = {(i-1) % 3, i % 3, (i+1) %3};
         int ix = i % 3;
         int imx = (i - 1) % 3;
         int ipx = (i + 1) % 3;
@@ -151,19 +152,16 @@ void RawImageSource::eahd_demosaic ()
             homv[ipx][j] = 0;
         }
 
-        int sh, sv, idx;
-
         for (int j = 1; j < W - 1; j++) {
             int dmi = 0;
-
+            for (int x = -1; x <= 0; x++) {
-            for (int x = -1; x <= 1; x++) {
+                int idx = mod[x + 1];
-                idx = (i + x) % 3;
 
                 for (int y = -1; y <= 1; y++) {
                     // compute distance in a, b, and L
                     if (dmi < 4) {
-                        sh = homh[idx][j + y];
+                        int sh = homh[idx][j + y];
-                        sv = homv[idx][j + y];
+                        int sv = homv[idx][j + y];
 
                         if (sh > sv) { // fixate horizontal pixel
                             dLmaph[dmi]  = DIST(lLh[ix][j], lLh[idx][j + y]);
@@ -199,25 +197,29 @@ void RawImageSource::eahd_demosaic ()
                     dmi++;
                 }
             }
+            int idx = mod[2];
+
+            for (int y = -1; y <= 1; y++) {
+                // compute distance in a, b, and L
+                dLmaph[dmi]  = DIST(lLh[ix][j], lLh[idx][j + y]);
+                dCamaph[dmi] = DIST(lah[ix][j], lah[idx][j + y]);
+                dCbmaph[dmi] = DIST(lbh[ix][j], lbh[idx][j + y]);
+                dLmapv[dmi]  = DIST(lLv[ix][j], lLv[idx][j + y]);
+                dCamapv[dmi] = DIST(lav[ix][j], lav[idx][j + y]);
+                dCbmapv[dmi] = DIST(lbv[ix][j], lbv[idx][j + y]);
+                dmi++;
+            }
 
             // compute eL & eC
-            int eL = min(max(dLmaph[3], dLmaph[5]), max(dLmapv[1], dLmapv[7]));
+            float eL = min(max(dLmaph[3], dLmaph[5]), max(dLmapv[1], dLmapv[7]));
-            int eCa = min(max(dCamaph[3], dCamaph[5]), max(dCamapv[1], dCamapv[7]));
+            float eCa = min(max(dCamaph[3], dCamaph[5]), max(dCamapv[1], dCamapv[7]));
-            int eCb = min(max(dCbmaph[3], dCbmaph[5]), max(dCbmapv[1], dCbmapv[7]));
+            float eCb = min(max(dCbmaph[3], dCbmaph[5]), max(dCbmapv[1], dCbmapv[7]));
 
             int wh = 0;
 
-            for (int dmi = 0; dmi < 9; dmi++)
+            for (int dmi = 0; dmi < 9; dmi++) {
-                if (dLmaph[dmi] <= eL && dCamaph[dmi] <= eCa && dCbmaph[dmi] <= eCb) {
+                wh += (dLmaph[dmi] <= eL) * (dCamaph[dmi] <= eCa) * (dCbmaph[dmi] <= eCb);
-                    wh++;
+            }
-                }
-
-            int wv = 0;
-
-            for (int dmi = 0; dmi < 9; dmi++)
-                if (dLmapv[dmi] <= eL && dCamapv[dmi] <= eCa && dCbmapv[dmi] <= eCb) {
-                    wv++;
-                }
 
             homh[imx][j - 1] += wh;
             homh[imx][j]  += wh;
@@ -229,6 +231,12 @@ void RawImageSource::eahd_demosaic ()
             homh[ipx][j]  += wh;
             homh[ipx][j + 1] += wh;
 
+            int wv = 0;
+
+            for (int dmi = 0; dmi < 9; dmi++) {
+                wv += (dLmapv[dmi] <= eL) * (dCamapv[dmi] <= eCa) * (dCbmapv[dmi] <= eCb);
+            }
+
             homv[imx][j - 1] += wv;
             homv[imx][j]  += wv;
             homv[imx][j + 1] += wv;
@@ -240,16 +248,13 @@ void RawImageSource::eahd_demosaic ()
             homv[ipx][j + 1] += wv;
         }
 
-//}
         // finalize image
-        int hc, vc;
-
         for (int j = 0; j < W; j++) {
             if (ri->ISGREEN(i - 1, j)) {
                 green[i - 1][j] = rawData[i - 1][j];
             } else {
-                hc = homh[imx][j];
+                int hc = homh[imx][j];
-                vc = homv[imx][j];
+                int vc = homv[imx][j];
 
                 if (hc > vc) {
                     green[i - 1][j] = gh[(i - 1) % 4][j];
@@ -267,12 +272,10 @@ void RawImageSource::eahd_demosaic ()
     }
 
     // finish H-2th and H-1th row, homogenity value is still valailable
-    int hc, vc;
-
     for (int i = H - 1; i < H + 1; i++)
         for (int j = 0; j < W; j++) {
-            hc = homh[(i - 1) % 3][j];
+            int hc = homh[(i - 1) % 3][j];
-            vc = homv[(i - 1) % 3][j];
+            int vc = homv[(i - 1) % 3][j];
 
             if (hc > vc) {
                 green[i - 1][j] = gh[(i - 1) % 4][j];
@@ -284,6 +287,7 @@ void RawImageSource::eahd_demosaic ()
         }
 
     // Interpolate R and B
+    #pragma omp parallel for
     for (int i = 0; i < H; i++) {
         if (i == 0) {
             interpolate_row_rb_mul_pp (rawData, red[i], blue[i], nullptr, green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1);

rtengine/rawimagesource.h

@@ -76,9 +76,9 @@ protected:
     unsigned int numFrames = 0;
 
     // to accelerate CIELAB conversion:
-    double lc00, lc01, lc02, lc10, lc11, lc12, lc20, lc21, lc22;
+    float lc00, lc01, lc02, lc10, lc11, lc12, lc20, lc21, lc22;
-    double* cache;
+    float* cache;
-    int threshold;
+    float threshold;
 
     array2D<float> rawData;  // holds preprocessed pixel values, rowData[i][j] corresponds to the ith row and jth column
     array2D<float> *rawDataFrames[4] = {nullptr};

rtengine/rawimagesource_i.h

@@ -62,22 +62,22 @@ inline void RawImageSource::convert_to_cielab_row (float* ar, float* ag, float*
 {
 
     for (int j = 0; j < W; j++) {
-        double r = ar[j];
+        float r = ar[j];
-        double g = ag[j];
+        float g = ag[j];
-        double b = ab[j];
+        float b = ab[j];
 
-        double x = lc00 * r + lc01 * g + lc02 * b;
+        float x = lc00 * r + lc01 * g + lc02 * b;
-        double y = lc10 * r + lc11 * g + lc12 * b;
+        float y = lc10 * r + lc11 * g + lc12 * b;
-        double z = lc20 * r + lc21 * g + lc22 * b;
+        float z = lc20 * r + lc21 * g + lc22 * b;
 
         if (y > threshold) {
             oL[j] = cache[(int)y];
         } else {
-            oL[j] = float(903.3 * y / MAXVALD);
+            oL[j] = float(903.3f * y / MAXVALF);
         }
 
-        oa[j] = float(500.0 * ((x > threshold ? cache[(int)x] : 7.787 * x / MAXVALD + 16.0 / 116.0) - (y > threshold ? cache[(int)y] : 7.787 * y / MAXVALD + 16.0 / 116.0)));
+        oa[j] = 500.f * ((x > threshold ? cache[(int)x] : 7.787f * x / MAXVALF + 16.f / 116.f) - (y > threshold ? cache[(int)y] : 7.787f * y / MAXVALF + 16.f / 116.f));
-        ob[j] = float(200.0 * ((y > threshold ? cache[(int)y] : 7.787 * y / MAXVALD + 16.0 / 116.0) - (z > threshold ? cache[(int)z] : 7.787 * z / MAXVALD + 16.0 / 116.0)));
+        ob[j] = 200.f * ((y > threshold ? cache[(int)y] : 7.787f * y / MAXVALF + 16.f / 116.f) - (z > threshold ? cache[(int)z] : 7.787f * z / MAXVALF + 16.f / 116.f));
     }
 }
 
@@ -154,125 +154,103 @@ inline void RawImageSource::interpolate_row_rb (float* ar, float* ab, float* pg,
                 : 0.f;
     };
 
-    if ((ri->ISRED(i, 0) || ri->ISRED(i, 1))) {
+    float *nonGreen1 = ar;
-        // RGRGR or GRGRGR line
+    float *nonGreen2 = ab;
-        for (int j = 0; j < W; j++) {
-            if (ri->ISRED(i, j)) {
-                // red is simple
-                ar[j] = rawData[i][j];
-                // blue: cross interpolation
-                int b = 0;
-                int n = 0;
 
-                if (i > 0 && j > 0) {
+    if ((ri->ISBLUE(i, 0) || ri->ISBLUE(i, 1))) {
-                    b += rawData[i - 1][j - 1] - getPg(j - 1);
+        std::swap(nonGreen1, nonGreen2);
-                    n++;
+    }
-                }
+    int j = FC(i, 0) & 1;
+    if (j) {
+        // linear R-G interp. horizontally
+        float val1;
 
-                if (i > 0 && j < W - 1) {
+        val1 = cg[0] + rawData[i][1] - cg[1];
-                    b += rawData[i - 1][j + 1] - getPg(j + 1);
-                    n++;
-                }
 
-                if (i < H - 1 && j > 0) {
+        nonGreen1[0] = CLIP(val1);
-                    b += rawData[i + 1][j - 1] - getNg(j - 1);
+        // linear B-G interp. vertically
-                    n++;
+        float val2;
-                }
 
-                if (i < H - 1 && j < W - 1) {
+        if (i == 0) {
-                    b += rawData[i + 1][j + 1] - getNg(j + 1);
+            val2 = getNg(0) + rawData[1][0] - cg[0];
-                    n++;
+        } else if (i == H - 1) {
-                }
+            val2 = getPg(0) + rawData[H - 2][0] - cg[0];
-
+        } else {
-                b = cg[j] + b / std::max(1, n);
+            val2 = cg[0] + (rawData[i - 1][0] - getPg(0) + rawData[i + 1][0] - getNg(0)) / 2;
-                ab[j] = b;
-            } else {
-                // linear R-G interp. horizontally
-                int r;
-
-                if (j == 0) {
-                    r = cg[0] + rawData[i][1] - cg[1];
-                } else if (j == W - 1) {
-                    r = cg[W - 1] + rawData[i][W - 2] - cg[W - 2];
-                } else {
-                    r = cg[j] + (rawData[i][j - 1] - cg[j - 1] + rawData[i][j + 1] - cg[j + 1]) / 2;
-                }
-
-                ar[j] = CLIP(r);
-                // linear B-G interp. vertically
-                int b;
-
-                if (i == 0) {
-                    b = getNg(j) + rawData[1][j] - cg[j];
-                } else if (i == H - 1) {
-                    b = getPg(j) + rawData[H - 2][j] - cg[j];
-                } else {
-                    b = cg[j] + (rawData[i - 1][j] - getPg(j) + rawData[i + 1][j] - getNg(j)) / 2;
-                }
-
-                ab[j] = b;
-            }
         }
-    } else {
-        // BGBGB or GBGBGB line
-        for (int j = 0; j < W; j++) {
-            if (ri->ISBLUE(i, j)) {
-                // red is simple
-                ab[j] = rawData[i][j];
-                // blue: cross interpolation
-                int r = 0;
-                int n = 0;
 
-                if (i > 0 && j > 0) {
+        nonGreen2[0] = val2;
-                    r += rawData[i - 1][j - 1] - getPg(j - 1);
+    }
-                    n++;
+    // RGRGR or GRGRGR line
-                }
+    for (; j < W - 1; j += 2) {
+        // nonGreen of row is simple
+        nonGreen1[j] = rawData[i][j];
+        // non green of next row: cross interpolation
+        float nonGreen = 0.f;
+        int n = 0;
 
-                if (i > 0 && j < W - 1) {
+        if (i > 0) {
-                    r += rawData[i - 1][j + 1] - getPg(j + 1);
+            if (j > 0) {
-                    n++;
+                nonGreen += rawData[i - 1][j - 1] - getPg(j - 1);
-                }
+                n++;
-
-                if (i < H - 1 && j > 0) {
-                    r += rawData[i + 1][j - 1] - getNg(j - 1);
-                    n++;
-                }
-
-                if (i < H - 1 && j < W - 1) {
-                    r += rawData[i + 1][j + 1] - getNg(j + 1);
-                    n++;
-                }
-
-                r = cg[j] + r / std::max(n, 1);
-
-                ar[j] = r;
-            } else {
-                // linear B-G interp. horizontally
-                int b;
-
-                if (j == 0) {
-                    b = cg[0] + rawData[i][1] - cg[1];
-                } else if (j == W - 1) {
-                    b = cg[W - 1] + rawData[i][W - 2] - cg[W - 2];
-                } else {
-                    b = cg[j] + (rawData[i][j - 1] - cg[j - 1] + rawData[i][j + 1] - cg[j + 1]) / 2;
-                }
-
-                ab[j] = CLIP(b);
-                // linear R-G interp. vertically
-                int r;
-
-                if (i == 0) {
-                    r = getNg(j) + rawData[1][j] - cg[j];
-                } else if (i == H - 1) {
-                    r = getPg(j) + rawData[H - 2][j] - cg[j];
-                } else {
-                    r = cg[j] + (rawData[i - 1][j] - getPg(j) + rawData[i + 1][j] - getNg(j)) / 2;
-                }
-
-                ar[j] = r;
             }
+            nonGreen += rawData[i - 1][j + 1] - getPg(j + 1);
+            n++;
         }
+
+        if (i < H - 1) {
+            if (j > 0) {
+                nonGreen += rawData[i + 1][j - 1] - getNg(j - 1);
+                n++;
+            }
+            nonGreen += rawData[i + 1][j + 1] - getNg(j + 1);
+            n++;
+        }
+
+        nonGreen2[j] = cg[j] + nonGreen / n;
+
+        // linear R-G interp. horizontally
+        float val1;
+
+        if (j == W - 2) {
+            val1 = cg[W - 1] + rawData[i][W - 2] - cg[W - 2];
+        } else {
+            val1 = cg[j + 1] + (rawData[i][j] - cg[j] + rawData[i][j + 2] - cg[j + 2]) / 2;
+        }
+
+        nonGreen1[j + 1] = CLIP(val1);
+        // linear B-G interp. vertically
+        float val2;
+
+        if (i == 0) {
+            val2 = getNg(j + 1) + rawData[1][j + 1] - cg[j + 1];
+        } else if (i == H - 1) {
+            val2 = getPg(j + 1) + rawData[H - 2][j + 1] - cg[j + 1];
+        } else {
+            val2 = cg[j + 1] + (rawData[i - 1][j + 1] - getPg(j + 1) + rawData[i + 1][j + 1] - getNg(j + 1)) / 2;
+        }
+
+        nonGreen2[j + 1] = val2;
+    }
+
+    if(j == W - 1) {
+        // nonGreen of row is simple
+        nonGreen1[j] = rawData[i][j];
+        // non green of next row: cross interpolation
+        float nonGreen = 0.f;
+        int n = 0;
+
+        if (i > 0) {
+            nonGreen += rawData[i - 1][j - 1] - getPg(j - 1);
+            n++;
+        }
+
+        if (i < H - 1) {
+            nonGreen += rawData[i + 1][j - 1] - getNg(j - 1);
+            n++;
+        }
+
+        nonGreen2[j] = cg[j] + nonGreen / n;
     }
 }