diff --git a/rtengine/demosaic_algos.cc b/rtengine/demosaic_algos.cc index cc0b99beb..7d5c74264 100644 --- a/rtengine/demosaic_algos.cc +++ b/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]; - threshold = (int)(0.008856 * MAXVALD); + cache = new float[maxindex]; + 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), - homh (W, 3), homv (W, 3); + lLv (W, 3), lav (W, 3), lbv (W, 3); + JaggedArray 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]; - int dLmapv[9]; - int dCamaph[9]; - int dCamapv[9]; - int dCbmaph[9]; - int dCbmapv[9]; + float dLmaph[9]; + float dLmapv[9]; + float dCamaph[9]; + float dCamapv[9]; + float dCbmaph[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 <= 1; x++) { - idx = (i + x) % 3; + for (int x = -1; x <= 0; x++) { + int idx = mod[x + 1]; for (int y = -1; y <= 1; y++) { // compute distance in a, b, and L if (dmi < 4) { - sh = homh[idx][j + y]; - sv = homv[idx][j + y]; + int sh = homh[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])); - int 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 eL = min(max(dLmaph[3], dLmaph[5]), max(dLmapv[1], dLmapv[7])); + float eCa = min(max(dCamaph[3], dCamaph[5]), max(dCamapv[1], dCamapv[7])); + float eCb = min(max(dCbmaph[3], dCbmaph[5]), max(dCbmapv[1], dCbmapv[7])); int wh = 0; - for (int dmi = 0; dmi < 9; dmi++) - if (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++; - } + for (int dmi = 0; dmi < 9; dmi++) { + wh += (dLmaph[dmi] <= eL) * (dCamaph[dmi] <= eCa) * (dCbmaph[dmi] <= eCb); + } 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]; - vc = homv[imx][j]; + int hc = homh[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]; - vc = homv[(i - 1) % 3][j]; + int hc = homh[(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); diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h index 53d3e0dc1..0c335b4ab 100644 --- a/rtengine/rawimagesource.h +++ b/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; - double* cache; - int threshold; + float lc00, lc01, lc02, lc10, lc11, lc12, lc20, lc21, lc22; + float* cache; + float threshold; array2D rawData; // holds preprocessed pixel values, rowData[i][j] corresponds to the ith row and jth column array2D *rawDataFrames[4] = {nullptr}; diff --git a/rtengine/rawimagesource_i.h b/rtengine/rawimagesource_i.h index 91d62ecab..fb85e4933 100644 --- a/rtengine/rawimagesource_i.h +++ b/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]; - double g = ag[j]; - double b = ab[j]; + float r = ar[j]; + float g = ag[j]; + float b = ab[j]; - double x = lc00 * r + lc01 * g + lc02 * b; - double y = lc10 * r + lc11 * g + lc12 * b; - double z = lc20 * r + lc21 * g + lc22 * b; + float x = lc00 * r + lc01 * g + lc02 * b; + float y = lc10 * r + lc11 * g + lc12 * 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))); - 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))); + 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] = 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))) { - // RGRGR or GRGRGR line - 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; + float *nonGreen1 = ar; + float *nonGreen2 = ab; - if (i > 0 && j > 0) { - b += rawData[i - 1][j - 1] - getPg(j - 1); - n++; - } + if ((ri->ISBLUE(i, 0) || ri->ISBLUE(i, 1))) { + std::swap(nonGreen1, nonGreen2); + } + int j = FC(i, 0) & 1; + if (j) { + // linear R-G interp. horizontally + float val1; - if (i > 0 && j < W - 1) { - b += rawData[i - 1][j + 1] - getPg(j + 1); - n++; - } + val1 = cg[0] + rawData[i][1] - cg[1]; - if (i < H - 1 && j > 0) { - b += rawData[i + 1][j - 1] - getNg(j - 1); - n++; - } + nonGreen1[0] = CLIP(val1); + // linear B-G interp. vertically + float val2; - if (i < H - 1 && j < W - 1) { - b += rawData[i + 1][j + 1] - getNg(j + 1); - n++; - } - - b = cg[j] + b / std::max(1, n); - 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; - } + if (i == 0) { + val2 = getNg(0) + rawData[1][0] - cg[0]; + } else if (i == H - 1) { + val2 = getPg(0) + rawData[H - 2][0] - cg[0]; + } else { + val2 = cg[0] + (rawData[i - 1][0] - getPg(0) + rawData[i + 1][0] - getNg(0)) / 2; } - } 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) { - r += rawData[i - 1][j - 1] - getPg(j - 1); - n++; - } + nonGreen2[0] = val2; + } + // 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) { - r += 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; + if (i > 0) { + if (j > 0) { + nonGreen += rawData[i - 1][j - 1] - getPg(j - 1); + n++; } + 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; } }