From 7c10f92ace3ba69ccccc314ceb49f6d9a5b16f7d Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Fri, 12 Oct 2018 16:01:48 +0200 Subject: [PATCH] dehaze: improved use of the guided filter for less halos --- rtengine/ipdehaze.cc | 152 ++++++++++++++++++++++++++----------------- 1 file changed, 94 insertions(+), 58 deletions(-) diff --git a/rtengine/ipdehaze.cc b/rtengine/ipdehaze.cc index 7fa1988ce..355d7f843 100644 --- a/rtengine/ipdehaze.cc +++ b/rtengine/ipdehaze.cc @@ -58,28 +58,28 @@ namespace { #endif -int get_dark_channel(const Imagefloat &src, array2D &dst, +int get_dark_channel(const array2D &R, const array2D &G, const array2D &B, array2D &dst, int patchsize, float *ambient, bool multithread) { - const int w = src.getWidth(); - const int h = src.getHeight(); + const int W = R.width(); + const int H = R.height(); int npatches = 0; #ifdef _OPENMP #pragma omp parallel for if (multithread) #endif - for (int y = 0; y < src.getHeight(); y += patchsize) { - int pH = std::min(y+patchsize, h); - for (int x = 0; x < src.getWidth(); x += patchsize, ++npatches) { + for (int y = 0; y < H; y += patchsize) { + int pH = std::min(y+patchsize, H); + for (int x = 0; x < W; x += patchsize, ++npatches) { float val = RT_INFINITY_F; - int pW = std::min(x+patchsize, w); + int pW = std::min(x+patchsize, W); for (int yy = y; yy < pH; ++yy) { float yval = RT_INFINITY_F; for (int xx = x; xx < pW; ++xx) { - float r = src.r(yy, xx); - float g = src.g(yy, xx); - float b = src.b(yy, xx); + float r = R[yy][xx]; + float g = G[yy][xx]; + float b = B[yy][xx]; if (ambient) { r /= ambient[0]; g /= ambient[1]; @@ -89,14 +89,16 @@ int get_dark_channel(const Imagefloat &src, array2D &dst, } val = min(val, yval); } + val = LIM01(val); for (int yy = y; yy < pH; ++yy) { std::fill(dst[yy]+x, dst[yy]+pW, val); } + float val2 = RT_INFINITY_F; for (int yy = y; yy < pH; ++yy) { for (int xx = x; xx < pW; ++xx) { - float r = src.r(yy, xx); - float g = src.g(yy, xx); - float b = src.b(yy, xx); + float r = R[yy][xx]; + float g = G[yy][xx]; + float b = B[yy][xx]; if (ambient) { r /= ambient[0]; g /= ambient[1]; @@ -104,7 +106,18 @@ int get_dark_channel(const Imagefloat &src, array2D &dst, } float l = min(r, g, b); if (l >= 2.f * val) { - dst[yy][xx] = l; + val2 = min(val2, l); + dst[yy][xx] = -1; + } + } + } + if (val2 < RT_INFINITY_F) { + val2 = LIM01(val2); + for (int yy = y; yy < pH; ++yy) { + for (int xx = x; xx < pW; ++xx) { + if (dst[yy][xx] < 0.f) { + dst[yy][xx] = val2; + } } } } @@ -115,10 +128,10 @@ int get_dark_channel(const Imagefloat &src, array2D &dst, } -int estimate_ambient_light(const Imagefloat *img, const array2D &dark, const array2D &Y, int patchsize, int npatches, float ambient[3]) +int estimate_ambient_light(const array2D &R, const array2D &G, const array2D &B, const array2D &dark, const array2D &Y, int patchsize, int npatches, float ambient[3]) { - const int W = img->getWidth(); - const int H = img->getHeight(); + const int W = R.width(); + const int H = R.height(); const auto get_percentile = [](std::priority_queue &q, float prcnt) -> float @@ -183,9 +196,9 @@ int estimate_ambient_light(const Imagefloat *img, const array2D &dark, co for (int y = p.second; y < pH; ++y) { for (int x = p.first; x < pW; ++x) { if (Y[y][x] >= lim) { - float r = img->r(y, x); - float g = img->g(y, x); - float b = img->b(y, x); + float r = R[y][x]; + float g = G[y][x]; + float b = B[y][x]; rr += r; gg += g; bb += b; @@ -218,41 +231,25 @@ void get_luminance(Imagefloat *img, array2D &Y, TMatrix ws, bool multithr } -void apply_contrast(array2D &dark, int contrast, double scale, bool multithread) +void apply_contrast(array2D &dark, float ambient, int contrast, double scale, bool multithread) { if (contrast) { const int W = dark.width(); const int H = dark.height(); - double tot = 0.0; -#ifdef _OPENMP - #pragma omp parallel for if (multithread) -#endif - for (int y = 0; y < H; ++y) { - double ytot = 0.0; - for (int x = 0; x < W; ++x) { - ytot += dark[y][x]; - } -#ifdef _OPENMP - #pragma omp critical -#endif - { - tot += ytot; - } - } - - float avg = tot / (W * H); + float avg = ambient * 0.25f; + float c = contrast * 0.3f; std::vector pts = { DCT_NURBS, 0, //black point. Value in [0 ; 1] range 0, //black point. Value in [0 ; 1] range - avg - avg * (0.6 - contrast / 250.0), //toe point - avg - avg * (0.6 + contrast / 250.0), //value at toe point + avg - avg * (0.6 - c / 250.0), //toe point + avg - avg * (0.6 + c / 250.0), //value at toe point - avg + (1 - avg) * (0.6 - contrast / 250.0), //shoulder point - avg + (1 - avg) * (0.6 + contrast / 250.0), //value at shoulder point + avg + (1 - avg) * (0.6 - c / 250.0), //shoulder point + avg + (1 - avg) * (0.6 + c / 250.0), //value at shoulder point 1., // white point 1. // value at white point @@ -271,6 +268,28 @@ void apply_contrast(array2D &dark, int contrast, double scale, bool multi } } + +void extract_channels(Imagefloat *img, const array2D &Y, array2D &r, array2D &g, array2D &b, int radius, float epsilon, bool multithread) +{ + const int W = img->getWidth(); + const int H = img->getHeight(); + +#ifdef _OPENMP + #pragma omp parallel for if (multithread) +#endif + for (int y = 0; y < H; ++y) { + for (int x = 0; x < W; ++x) { + r[y][x] = img->r(y, x); + g[y][x] = img->g(y, x); + b[y][x] = img->b(y, x); + } + } + + guidedFilter(Y, r, r, radius, epsilon, multithread); + guidedFilter(Y, g, g, radius, epsilon, multithread); + guidedFilter(Y, b, b, radius, epsilon, multithread); +} + } // namespace @@ -289,18 +308,24 @@ void ImProcFunctions::dehaze(Imagefloat *img) if (options.rtSettings.verbose) { std::cout << "dehaze: strength = " << strength << std::endl; } - - array2D dark(W, H); - const int patchsize = std::max(W / (200 + max(params->dehaze.detail, 0)), 2); - int npatches = get_dark_channel(*img, dark, patchsize, nullptr, multiThread); - DEBUG_DUMP(dark); TMatrix ws = ICCStore::getInstance()->workingSpaceMatrix(params->icm.workingProfile); array2D Y(W, H); get_luminance(img, Y, ws, multiThread); + + array2D R(W, H); + array2D G(W, H); + array2D B(W, H); + int patchsize = max(int(20 / scale), 2); + extract_channels(img, Y, R, G, B, patchsize, 1e-1, multiThread); + array2D dark(W, H); + patchsize = std::max(W / (200 + params->dehaze.detail * (SGN(params->dehaze.detail) > 0 ? 4 : 1)), 2); + int npatches = get_dark_channel(R, G, B, dark, patchsize, nullptr, multiThread); + DEBUG_DUMP(dark); + float ambient[3]; - int n = estimate_ambient_light(img, dark, Y, patchsize, npatches, ambient); + int n = estimate_ambient_light(R, G, B, dark, Y, patchsize, npatches, ambient); float ambient_Y = Color::rgbLuminance(ambient[0], ambient[1], ambient[2], ws); if (options.rtSettings.verbose) { @@ -320,8 +345,8 @@ void ImProcFunctions::dehaze(Imagefloat *img) } array2D &t_tilde = dark; - get_dark_channel(*img, dark, patchsize, ambient, multiThread); - apply_contrast(dark, params->dehaze.depth, scale, multiThread); + get_dark_channel(R, G, B, dark, patchsize, ambient, multiThread); + apply_contrast(dark, ambient_Y, params->dehaze.depth, scale, multiThread); DEBUG_DUMP(t_tilde); if (!params->dehaze.showDepthMap) { @@ -344,7 +369,8 @@ void ImProcFunctions::dehaze(Imagefloat *img) const int radius = max(int(patchsize * mult), 1); const float epsilon = 2.5e-4; array2D &t = t_tilde; - + + if (!params->dehaze.showDepthMap) guidedFilter(Y, t_tilde, t, radius, epsilon, multiThread); DEBUG_DUMP(t); @@ -362,16 +388,26 @@ void ImProcFunctions::dehaze(Imagefloat *img) return; } - const float t0 = 0.01; + const float t0 = 0.1; + const float teps = 1e-3; #ifdef _OPENMP #pragma omp parallel for if (multiThread) #endif for (int y = 0; y < H; ++y) { for (int x = 0; x < W; ++x) { - float mt = std::max(t[y][x], t0); - float r = (img->r(y, x) - ambient[0]) / mt + ambient[0]; - float g = (img->g(y, x) - ambient[1]) / mt + ambient[1]; - float b = (img->b(y, x) - ambient[2]) / mt + ambient[2]; + float rgb[3] = { img->r(y, x), img->g(y, x), img->b(y, x) }; + float tl = 1.f - min(rgb[0]/ambient[0], rgb[1]/ambient[1], rgb[2]/ambient[2]); + float tu = t0 - teps; + for (int c = 0; c < 3; ++c) { + if (ambient[c] < 1) { + tu = max(tu, (rgb[c] - ambient[c])/(1.f - ambient[c])); + } + } + float mt = max(t[y][x], t0, tl + teps, tu + teps); + float r = (rgb[0] - ambient[0]) / mt + ambient[0]; + float g = (rgb[1] - ambient[1]) / mt + ambient[1]; + float b = (rgb[2] - ambient[2]) / mt + ambient[2]; + img->r(y, x) = r; img->g(y, x) = g; img->b(y, x) = b; @@ -388,7 +424,7 @@ void ImProcFunctions::dehaze(Imagefloat *img) if (newmed > 1e-5f) { const float f1 = oldmed / newmed; - const float f = /* f1 * */ 65535.f; + const float f = f1 * 65535.f; #ifdef _OPENMP #pragma omp parallel for if (multiThread) #endif