diff --git a/rtengine/CMakeLists.txt b/rtengine/CMakeLists.txt index 385ba921e..9ac09fbcc 100644 --- a/rtengine/CMakeLists.txt +++ b/rtengine/CMakeLists.txt @@ -30,15 +30,19 @@ link_directories("${PROJECT_SOURCE_DIR}/rtexif" set(CAMCONSTSFILE "camconst.json") set(RTENGINESOURCEFILES + badpixels.cc + CA_correct_RT.cc + EdgePreservingDecomposition.cc + FTblockDN.cc + PF_correct_RT.cc ahd_demosaic_RT.cc amaze_demosaic_RT.cc - CA_correct_RT.cc + cJSON.c calc_distort.cc camconst.cc cfa_linedn_RT.cc ciecam02.cc cieimage.cc - cJSON.c clutstore.cc color.cc colortemp.cc @@ -55,18 +59,13 @@ set(RTENGINESOURCEFILES dual_demosaic_RT.cc dynamicprofile.cc eahd_demosaic.cc - EdgePreservingDecomposition.cc fast_demo.cc ffmanager.cc filmnegativeproc.cc flatcurves.cc - FTblockDN.cc - gamutwarning.cc gauss.cc green_equil_RT.cc - guidedfilter.cc hilite_recon.cc - histmatching.cc hphd_demosaic_RT.cc iccjpeg.cc iccstore.cc @@ -81,15 +80,10 @@ set(RTENGINESOURCEFILES improcfun.cc impulse_denoise.cc init.cc - ipdehaze.cc iplab2rgb.cc - iplabregions.cc - iplocalcontrast.cc ipresize.cc ipretinex.cc - ipshadowshighlights.cc ipsharpen.cc - ipsoftlight.cc iptransform.cc ipvibrance.cc ipwavelet.cc @@ -97,8 +91,8 @@ set(RTENGINESOURCEFILES jpeg_ijg/jpeg_memsrc.cc klt/convolve.cc klt/error.cc - klt/klt_util.cc klt/klt.cc + klt/klt_util.cc klt/pnmio.cc klt/pyramid.cc klt/selectGoodFeatures.cc @@ -107,11 +101,8 @@ set(RTENGINESOURCEFILES klt/writeFeatures.cc labimage.cc lcp.cc - lj92.c loadinitial.cc myfile.cc - pdaflinesfilter.cc - PF_correct_RT.cc pipettebuffer.cc pixelshift.cc previewimage.cc @@ -123,17 +114,27 @@ set(RTENGINESOURCEFILES rcd_demosaic.cc refreshmap.cc rt_algo.cc - rtlensfun.cc rtthumbnail.cc shmap.cc simpleprocess.cc slicer.cc stdimagesource.cc - tmo_fattal02.cc utils.cc - vng4_demosaic_RT.cc + rtlensfun.cc + tmo_fattal02.cc + iplocalcontrast.cc + histmatching.cc + pdaflinesfilter.cc + gamutwarning.cc + ipshadowshighlights.cc xtrans_demosaic.cc -) + vng4_demosaic_RT.cc + ipsoftlight.cc + guidedfilter.cc + ipdehaze.cc + iplabregions.cc + lj92.c + ) if(LENSFUN_HAS_LOAD_DIRECTORY) set_source_files_properties(rtlensfun.cc PROPERTIES COMPILE_DEFINITIONS RT_LENSFUN_HAS_LOAD_DIRECTORY) diff --git a/rtengine/badpixels.cc b/rtengine/badpixels.cc new file mode 100644 index 000000000..895294cae --- /dev/null +++ b/rtengine/badpixels.cc @@ -0,0 +1,576 @@ +/* + * This file is part of RawTherapee. + * + * Copyright (c) 2004-2019 Gabor Horvath + * + * 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 "array2D.h" +#include "median.h" +#include "pixelsmap.h" +#include "rawimagesource.h" + +namespace rtengine +{ + +/* interpolateBadPixelsBayer: correct raw pixels looking at the bitmap + * takes into consideration if there are multiple bad pixels in the neighborhood + */ +int RawImageSource::interpolateBadPixelsBayer(const PixelsMap &bitmapBads, array2D &rawData) +{ + constexpr float eps = 1.f; + int counter = 0; + +#ifdef _OPENMP + #pragma omp parallel for reduction(+:counter) schedule(dynamic,16) +#endif + + for (int row = 2; row < H - 2; ++row) { + for (int col = 2; col < W - 2; ++col) { + const int sk = bitmapBads.skipIfZero(col, row); //optimization for a stripe all zero + + if (sk) { + col += sk - 1; //-1 is because of col++ in cycle + continue; + } + + if (!bitmapBads.get(col, row)) { + continue; + } + + float wtdsum = 0.f, norm = 0.f; + + // diagonal interpolation + if (FC(row, col) == 1) { + // green channel. We can use closer pixels than for red or blue channel. Distance to center pixel is sqrt(2) => weighting is 0.70710678 + // For green channel following pixels will be used for interpolation. Pixel to be interpolated is in center. + // 1 means that pixel is used in this step, if itself and his counterpart are not marked bad + // 0 0 0 0 0 + // 0 1 0 1 0 + // 0 0 0 0 0 + // 0 1 0 1 0 + // 0 0 0 0 0 + for (int dx = -1; dx <= 1; dx += 2) { + if (bitmapBads.get(col + dx, row - 1) || bitmapBads.get(col - dx, row + 1)) { + continue; + } + + const float dirwt = 0.70710678f / (fabsf(rawData[row - 1][col + dx] - rawData[row + 1][col - dx]) + eps); + wtdsum += dirwt * (rawData[row - 1][col + dx] + rawData[row + 1][col - dx]); + norm += dirwt; + } + } else { + // red and blue channel. Distance to center pixel is sqrt(8) => weighting is 0.35355339 + // For red and blue channel following pixels will be used for interpolation. Pixel to be interpolated is in center. + // 1 means that pixel is used in this step, if itself and his counterpart are not marked bad + // 1 0 0 0 1 + // 0 0 0 0 0 + // 0 0 0 0 0 + // 0 0 0 0 0 + // 1 0 0 0 1 + for (int dx = -2; dx <= 2; dx += 4) { + if (bitmapBads.get(col + dx, row - 2) || bitmapBads.get(col - dx, row + 2)) { + continue; + } + + const float dirwt = 0.35355339f / (fabsf(rawData[row - 2][col + dx] - rawData[row + 2][col - dx]) + eps); + wtdsum += dirwt * (rawData[row - 2][col + dx] + rawData[row + 2][col - dx]); + norm += dirwt; + } + } + + // channel independent. Distance to center pixel is 2 => weighting is 0.5 + // Additionally for all channel following pixels will be used for interpolation. Pixel to be interpolated is in center. + // 1 means that pixel is used in this step, if itself and his counterpart are not marked bad + // 0 0 1 0 0 + // 0 0 0 0 0 + // 1 0 0 0 1 + // 0 0 0 0 0 + // 0 0 1 0 0 + + // horizontal interpolation + if (!(bitmapBads.get(col - 2, row) || bitmapBads.get(col + 2, row))) { + const float dirwt = 0.5f / (fabsf(rawData[row][col - 2] - rawData[row][col + 2]) + eps); + wtdsum += dirwt * (rawData[row][col - 2] + rawData[row][col + 2]); + norm += dirwt; + } + + // vertical interpolation + if (!(bitmapBads.get(col, row - 2) || bitmapBads.get(col, row + 2))) { + const float dirwt = 0.5f / (fabsf(rawData[row - 2][col] - rawData[row + 2][col]) + eps); + wtdsum += dirwt * (rawData[row - 2][col] + rawData[row + 2][col]); + norm += dirwt; + } + + if (LIKELY(norm > 0.f)) { // This means, we found at least one pair of valid pixels in the steps above, likelihood of this case is about 99.999% + rawData[row][col] = wtdsum / (2.f * norm); //gradient weighted average, Factor of 2.f is an optimization to avoid multiplications in former steps + counter++; + } else { //backup plan -- simple average. Same method for all channels. We could improve this, but it's really unlikely that this case happens + int tot = 0; + float sum = 0.f; + + for (int dy = -2; dy <= 2; dy += 2) { + for (int dx = -2; dx <= 2; dx += 2) { + if (bitmapBads.get(col + dx, row + dy)) { + continue; + } + + sum += rawData[row + dy][col + dx]; + tot++; + } + } + + if (tot > 0) { + rawData[row][col] = sum / tot; + counter ++; + } + } + } + } + + return counter; // Number of interpolated pixels. +} + +/* interpolateBadPixelsNcolors: correct raw pixels looking at the bitmap + * takes into consideration if there are multiple bad pixels in the neighborhood + */ +int RawImageSource::interpolateBadPixelsNColours(const PixelsMap &bitmapBads, const int colors) +{ + constexpr float eps = 1.f; + int counter = 0; + +#ifdef _OPENMP + #pragma omp parallel for reduction(+:counter) schedule(dynamic,16) +#endif + + for (int row = 2; row < H - 2; ++row) { + for (int col = 2; col < W - 2; ++col) { + const int sk = bitmapBads.skipIfZero(col, row); //optimization for a stripe all zero + + if (sk) { + col += sk - 1; //-1 is because of col++ in cycle + continue; + } + + if (!bitmapBads.get(col, row)) { + continue; + } + + float wtdsum[colors] = {}; + float norm[colors] = {}; + + // diagonal interpolation + for (int dx = -1; dx <= 1; dx += 2) { + if (bitmapBads.get(col + dx, row - 1) || bitmapBads.get(col - dx, row + 1)) { + continue; + } + + for (int c = 0; c < colors; ++c) { + const float dirwt = 0.70710678f / (fabsf(rawData[row - 1][(col + dx) * colors + c] - rawData[row + 1][(col - dx) * colors + c]) + eps); + wtdsum[c] += dirwt * (rawData[row - 1][(col + dx) * colors + c] + rawData[row + 1][(col - dx) * colors + c]); + norm[c] += dirwt; + } + } + + // horizontal interpolation + if (!(bitmapBads.get(col - 1, row) || bitmapBads.get(col + 1, row))) { + for (int c = 0; c < colors; ++c) { + const float dirwt = 1.f / (fabsf(rawData[row][(col - 1) * colors + c] - rawData[row][(col + 1) * colors + c]) + eps); + wtdsum[c] += dirwt * (rawData[row][(col - 1) * colors + c] + rawData[row][(col + 1) * colors + c]); + norm[c] += dirwt; + } + } + + // vertical interpolation + if (!(bitmapBads.get(col, row - 1) || bitmapBads.get(col, row + 1))) { + for (int c = 0; c < colors; ++c) { + const float dirwt = 1.f / (fabsf(rawData[row - 1][col * colors + c] - rawData[row + 1][col * colors + c]) + eps); + wtdsum[c] += dirwt * (rawData[row - 1][col * colors + c] + rawData[row + 1][col * colors + c]); + norm[c] += dirwt; + } + } + + if (LIKELY(norm[0] > 0.f)) { // This means, we found at least one pair of valid pixels in the steps above, likelihood of this case is about 99.999% + for (int c = 0; c < colors; ++c) { + rawData[row][col * colors + c] = wtdsum[c] / (2.f * norm[c]); //gradient weighted average, Factor of 2.f is an optimization to avoid multiplications in former steps + } + + counter++; + } else { //backup plan -- simple average. Same method for all channels. We could improve this, but it's really unlikely that this case happens + int tot = 0; + float sum[colors] = {}; + + for (int dy = -2; dy <= 2; dy += 2) { + for (int dx = -2; dx <= 2; dx += 2) { + if (bitmapBads.get(col + dx, row + dy)) { + continue; + } + + for (int c = 0; c < colors; ++c) { + sum[c] += rawData[row + dy][(col + dx) * colors + c]; + } + + tot++; + } + } + + if (tot > 0) { + for (int c = 0; c < colors; ++c) { + rawData[row][col * colors + c] = sum[c] / tot; + } + + counter ++; + } + } + } + } + + return counter; // Number of interpolated pixels. +} + +/* interpolateBadPixelsXtrans: correct raw pixels looking at the bitmap + * takes into consideration if there are multiple bad pixels in the neighborhood + */ +int RawImageSource::interpolateBadPixelsXtrans(const PixelsMap &bitmapBads) +{ + constexpr float eps = 1.f; + int counter = 0; + +#ifdef _OPENMP + #pragma omp parallel for reduction(+:counter) schedule(dynamic,16) +#endif + + for (int row = 2; row < H - 2; ++row) { + for (int col = 2; col < W - 2; ++col) { + const int skip = bitmapBads.skipIfZero(col, row); //optimization for a stripe all zero + + if (skip) { + col += skip - 1; //-1 is because of col++ in cycle + continue; + } + + if (!bitmapBads.get(col, row)) { + continue; + } + + float wtdsum = 0.f, norm = 0.f; + const unsigned int pixelColor = ri->XTRANSFC(row, col); + + if (pixelColor == 1) { + // green channel. A green pixel can either be a solitary green pixel or a member of a 2x2 square of green pixels + if (ri->XTRANSFC(row, col - 1) == ri->XTRANSFC(row, col + 1)) { + // If left and right neighbor have same color, then this is a solitary green pixel + // For these the following pixels will be used for interpolation. Pixel to be interpolated is in center and marked with a P. + // Pairs of pixels used in this step are numbered. A pair will be used if none of the pixels of the pair is marked bad + // 0 means, the pixel has a different color and will not be used + // 0 1 0 2 0 + // 3 5 0 6 4 + // 0 0 P 0 0 + // 4 6 0 5 3 + // 0 2 0 1 0 + for (int dx = -1; dx <= 1; dx += 2) { // pixels marked 5 or 6 in above example. Distance to P is sqrt(2) => weighting is 0.70710678f + if (bitmapBads.get(col + dx, row - 1) || bitmapBads.get(col - dx, row + 1)) { + continue; + } + + const float dirwt = 0.70710678f / (fabsf(rawData[row - 1][col + dx] - rawData[row + 1][col - dx]) + eps); + wtdsum += dirwt * (rawData[row - 1][col + dx] + rawData[row + 1][col - dx]); + norm += dirwt; + } + + for (int dx = -1; dx <= 1; dx += 2) { // pixels marked 1 or 2 on above example. Distance to P is sqrt(5) => weighting is 0.44721359f + if (bitmapBads.get(col + dx, row - 2) || bitmapBads.get(col - dx, row + 2)) { + continue; + } + + const float dirwt = 0.44721359f / (fabsf(rawData[row - 2][col + dx] - rawData[row + 2][col - dx]) + eps); + wtdsum += dirwt * (rawData[row - 2][col + dx] + rawData[row + 2][col - dx]); + norm += dirwt; + } + + for (int dx = -2; dx <= 2; dx += 4) { // pixels marked 3 or 4 on above example. Distance to P is sqrt(5) => weighting is 0.44721359f + if (bitmapBads.get(col + dx, row - 1) || bitmapBads.get(col - dx, row + 1)) { + continue; + } + + const float dirwt = 0.44721359f / (fabsf(rawData[row - 1][col + dx] - rawData[row + 1][col - dx]) + eps); + wtdsum += dirwt * (rawData[row - 1][col + dx] + rawData[row + 1][col - dx]); + norm += dirwt; + } + } else { + // this is a member of a 2x2 square of green pixels + // For these the following pixels will be used for interpolation. Pixel to be interpolated is at position P in the example. + // Pairs of pixels used in this step are numbered. A pair will be used if none of the pixels of the pair is marked bad + // 0 means, the pixel has a different color and will not be used + // 1 0 0 3 + // 0 P 2 0 + // 0 2 1 0 + // 3 0 0 0 + + // pixels marked 1 in above example. Distance to P is sqrt(2) => weighting is 0.70710678f + const int offset1 = ri->XTRANSFC(row - 1, col - 1) == ri->XTRANSFC(row + 1, col + 1) ? 1 : -1; + + if (!(bitmapBads.get(col - offset1, row - 1) || bitmapBads.get(col + offset1, row + 1))) { + const float dirwt = 0.70710678f / (fabsf(rawData[row - 1][col - offset1] - rawData[row + 1][col + offset1]) + eps); + wtdsum += dirwt * (rawData[row - 1][col - offset1] + rawData[row + 1][col + offset1]); + norm += dirwt; + } + + // pixels marked 2 in above example. Distance to P is 1 => weighting is 1.f + int offsety = ri->XTRANSFC(row - 1, col) != 1 ? 1 : -1; + int offsetx = offset1 * offsety; + + if (!(bitmapBads.get(col + offsetx, row) || bitmapBads.get(col, row + offsety))) { + const float dirwt = 1.f / (fabsf(rawData[row][col + offsetx] - rawData[row + offsety][col]) + eps); + wtdsum += dirwt * (rawData[row][col + offsetx] + rawData[row + offsety][col]); + norm += dirwt; + } + + const int offsety2 = -offsety; + const int offsetx2 = -offsetx; + offsetx *= 2; + offsety *= 2; + + // pixels marked 3 in above example. Distance to P is sqrt(5) => weighting is 0.44721359f + if (!(bitmapBads.get(col + offsetx, row + offsety2) || bitmapBads.get(col + offsetx2, row + offsety))) { + const float dirwt = 0.44721359f / (fabsf(rawData[row + offsety2][col + offsetx] - rawData[row + offsety][col + offsetx2]) + eps); + wtdsum += dirwt * (rawData[row + offsety2][col + offsetx] + rawData[row + offsety][col + offsetx2]); + norm += dirwt; + } + } + } else { + // red and blue channel. + // Each red or blue pixel has exactly one neighbor of same color in distance 2 and four neighbors of same color which can be reached by a move of a knight in chess. + // For the distance 2 pixel (marked with an X) we generate a virtual counterpart (marked with a V) + // For red and blue channel following pixels will be used for interpolation. Pixel to be interpolated is in center and marked with a P. + // Pairs of pixels used in this step are numbered except for distance 2 pixels which are marked X and V. A pair will be used if none of the pixels of the pair is marked bad + // 0 1 0 0 0 0 0 X 0 0 remaining cases are symmetric + // 0 0 0 0 2 1 0 0 0 2 + // X 0 P 0 V 0 0 P 0 0 + // 0 0 0 0 1 0 0 0 0 0 + // 0 2 0 0 0 0 2 V 1 0 + + // Find two knight moves landing on a pixel of same color as the pixel to be interpolated. + // If we look at first and last row of 5x5 square, we will find exactly two knight pixels. + // Additionally we know that the column of this pixel has 1 or -1 horizontal distance to the center pixel + // When we find a knight pixel, we get its counterpart, which has distance (+-3,+-3), where the signs of distance depend on the corner of the found knight pixel. + // These pixels are marked 1 or 2 in above examples. Distance to P is sqrt(5) => weighting is 0.44721359f + // The following loop simply scans the four possible places. To keep things simple, it does not stop after finding two knight pixels, because it will not find more than two + for (int d1 = -2, offsety = 3; d1 <= 2; d1 += 4, offsety -= 6) { + for (int d2 = -1, offsetx = 3; d2 < 1; d2 += 2, offsetx -= 6) { + if (ri->XTRANSFC(row + d1, col + d2) == pixelColor) { + if (!(bitmapBads.get(col + d2, row + d1) || bitmapBads.get(col + d2 + offsetx, row + d1 + offsety))) { + const float dirwt = 0.44721359f / (fabsf(rawData[row + d1][col + d2] - rawData[row + d1 + offsety][col + d2 + offsetx]) + eps); + wtdsum += dirwt * (rawData[row + d1][col + d2] + rawData[row + d1 + offsety][col + d2 + offsetx]); + norm += dirwt; + } + } + } + } + + // now scan for the pixel of same color in distance 2 in each direction (marked with an X in above examples). + bool distance2PixelFound = false; + int dx, dy; + + // check horizontal + for (dx = -2, dy = 0; dx <= 2; dx += 4) { + if (ri->XTRANSFC(row, col + dx) == pixelColor) { + distance2PixelFound = true; + break; + } + } + + if (!distance2PixelFound) { + // no distance 2 pixel on horizontal, check vertical + for (dx = 0, dy = -2; dy <= 2; dy += 4) { + if (ri->XTRANSFC(row + dy, col) == pixelColor) { + distance2PixelFound = true; + break; + } + } + } + + // calculate the value of its virtual counterpart (marked with a V in above examples) + float virtualPixel; + + if (dy == 0) { + virtualPixel = 0.5f * (rawData[row - 1][col - dx] + rawData[row + 1][col - dx]); + } else { + virtualPixel = 0.5f * (rawData[row - dy][col - 1] + rawData[row - dy][col + 1]); + } + + // and weight as usual. Distance to P is 2 => weighting is 0.5f + const float dirwt = 0.5f / (fabsf(virtualPixel - rawData[row + dy][col + dx]) + eps); + wtdsum += dirwt * (virtualPixel + rawData[row + dy][col + dx]); + norm += dirwt; + } + + if (LIKELY(norm > 0.f)) { // This means, we found at least one pair of valid pixels in the steps above, likelihood of this case is about 99.999% + rawData[row][col] = wtdsum / (2.f * norm); //gradient weighted average, Factor of 2.f is an optimization to avoid multiplications in former steps + counter++; + } + } + } + + return counter; // Number of interpolated pixels. +} + +/* Search for hot or dead pixels in the image and update the map + * For each pixel compare its value to the average of similar color surrounding + * (Taken from Emil Martinec idea) + * (Optimized by Ingo Weyrich 2013 and 2015) + */ +int RawImageSource::findHotDeadPixels(PixelsMap &bpMap, const float thresh, const bool findHotPixels, const bool findDeadPixels) const +{ + const float varthresh = (20.0 * (thresh / 100.0) + 1.0) / 24.f; + + // allocate temporary buffer + float* cfablur = new float[H * W]; + + // counter for dead or hot pixels + int counter = 0; + +#ifdef _OPENMP + #pragma omp parallel +#endif + { +#ifdef _OPENMP + #pragma omp for schedule(dynamic,16) nowait +#endif + + for (int i = 2; i < H - 2; i++) { + for (int j = 2; j < W - 2; j++) { + const float temp = median(rawData[i - 2][j - 2], rawData[i - 2][j], rawData[i - 2][j + 2], + rawData[i][j - 2], rawData[i][j], rawData[i][j + 2], + rawData[i + 2][j - 2], rawData[i + 2][j], rawData[i + 2][j + 2]); + cfablur[i * W + j] = rawData[i][j] - temp; + } + } + + // process borders. Former version calculated the median using mirrored border which does not make sense because the original pixel loses weight + // Setting the difference between pixel and median for border pixels to zero should do the job not worse then former version +#ifdef _OPENMP + #pragma omp single +#endif + { + for (int i = 0; i < 2; ++i) { + for (int j = 0; j < W; ++j) { + cfablur[i * W + j] = 0.f; + } + } + + for (int i = 2; i < H - 2; ++i) { + for (int j = 0; j < 2; ++j) { + cfablur[i * W + j] = 0.f; + } + + for (int j = W - 2; j < W; ++j) { + cfablur[i * W + j] = 0.f; + } + } + + for (int i = H - 2; i < H; ++i) { + for (int j = 0; j < W; ++j) { + cfablur[i * W + j] = 0.f; + } + } + } + +#ifdef _OPENMP + #pragma omp barrier // barrier because of nowait clause above + + #pragma omp for reduction(+:counter) schedule(dynamic,16) +#endif + + //cfa pixel heat/death evaluation + for (int rr = 2; rr < H - 2; ++rr) { + for (int cc = 2, rrmWpcc = rr * W + 2; cc < W - 2; ++cc, ++rrmWpcc) { + //evaluate pixel for heat/death + float pixdev = cfablur[rrmWpcc]; + + if (pixdev == 0.f) { + continue; + } + + if ((!findDeadPixels) && pixdev < 0) { + continue; + } + + if ((!findHotPixels) && pixdev > 0) { + continue; + } + + pixdev = fabsf(pixdev); + float hfnbrave = -pixdev; + +#ifdef __SSE2__ + // sum up 5*4 = 20 values using SSE + // 10 fabs function calls and 10 float additions with SSE + vfloat sum = vabsf(LVFU(cfablur[(rr - 2) * W + cc - 2])) + vabsf(LVFU(cfablur[(rr - 1) * W + cc - 2])); + sum += vabsf(LVFU(cfablur[(rr) * W + cc - 2])); + sum += vabsf(LVFU(cfablur[(rr + 1) * W + cc - 2])); + sum += vabsf(LVFU(cfablur[(rr + 2) * W + cc - 2])); + // horizontally add the values and add the result to hfnbrave + hfnbrave += vhadd(sum); + + // add remaining 5 values of last column + for (int mm = rr - 2; mm <= rr + 2; ++mm) { + hfnbrave += fabsf(cfablur[mm * W + cc + 2]); + } + +#else + + // 25 fabs function calls and 25 float additions without SSE + for (int mm = rr - 2; mm <= rr + 2; ++mm) { + for (int nn = cc - 2; nn <= cc + 2; ++nn) { + hfnbrave += fabsf(cfablur[mm * W + nn]); + } + } + +#endif + + if (pixdev > varthresh * hfnbrave) { + // mark the pixel as "bad" + bpMap.set(cc, rr); + counter++; + } + }//end of pixel evaluation + } + }//end of parallel processing + delete [] cfablur; + return counter; +} + +int RawImageSource::findZeroPixels(PixelsMap &bpMap) const +{ + int counter = 0; + +#ifdef _OPENMP + #pragma omp parallel for reduction(+:counter) schedule(dynamic,16) +#endif + + for (int i = 0; i < H; ++i) { + for (int j = 0; j < W; ++j) { + if (ri->data[i][j] == 0.f) { + bpMap.set(j, i); + counter++; + } + } + } + return counter; +} + + +} diff --git a/rtengine/dfmanager.h b/rtengine/dfmanager.h index 1cc22723b..62379187b 100644 --- a/rtengine/dfmanager.h +++ b/rtengine/dfmanager.h @@ -20,6 +20,7 @@ #include #include #include +#include "pixelsmap.h" #include "rawimage.h" namespace rtengine diff --git a/rtengine/ffmanager.h b/rtengine/ffmanager.h index 4a65c2ed7..7022d1641 100644 --- a/rtengine/ffmanager.h +++ b/rtengine/ffmanager.h @@ -86,7 +86,6 @@ public: protected: typedef std::multimap ffList_t; - typedef std::map > bpList_t; ffList_t ffList; bool initialized; Glib::ustring currentPath; diff --git a/rtengine/pixelsmap.h b/rtengine/pixelsmap.h new file mode 100644 index 000000000..9089ac91c --- /dev/null +++ b/rtengine/pixelsmap.h @@ -0,0 +1,98 @@ +#pragma once + +/* + * This file is part of RawTherapee. + * + * Copyright (c) 2004-2019 Gabor Horvath + * + * 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 + +#include "noncopyable.h" + +namespace rtengine +{ + +struct badPix { + uint16_t x; + uint16_t y; + badPix(uint16_t xc, uint16_t yc): x(xc), y(yc) {} +}; + +class PixelsMap : + public NonCopyable +{ + int w; // line width in base_t units + int h; // height + typedef unsigned long base_t; + static constexpr size_t base_t_size = sizeof(base_t); + base_t *pm; + +public: + PixelsMap(int width, int height) + : w((width / (base_t_size * 8)) + 1), h(height), pm(new base_t[h * w]) + { + clear(); + } + + ~PixelsMap() + { + delete [] pm; + } + int width() const + { + return w; + } + int height() const + { + return h; + } + + // if a pixel is set returns true + bool get(int x, int y) const + { + return (pm[y * w + x / (base_t_size * 8)] & (base_t)1 << (x % (base_t_size * 8))) != 0; + } + + // set a pixel + void set(int x, int y) + { + pm[y * w + x / (base_t_size * 8)] |= (base_t)1 << (x % (base_t_size * 8)) ; + } + + // set pixels from a list + int set(const std::vector &bp) + { + for (std::vector::const_iterator iter = bp.begin(); iter != bp.end(); ++iter) { + set(iter->x, iter->y); + } + + return bp.size(); + } + + void clear() + { + memset(pm, 0, h * w * base_t_size); + } + // return 0 if at least one pixel in the word(base_t) is set, otherwise return the number of pixels to skip to the next word base_t + int skipIfZero(int x, int y) const + { + return pm[y * w + x / (base_t_size * 8)] == 0 ? base_t_size * 8 - x % (base_t_size * 8) : 0; + } +}; + +} diff --git a/rtengine/procparams.cc b/rtengine/procparams.cc index ade3b0581..4ed97f767 100644 --- a/rtengine/procparams.cc +++ b/rtengine/procparams.cc @@ -2763,23 +2763,23 @@ ProcParams::ProcParams() void ProcParams::setDefaults() { - toneCurve = ToneCurveParams(); + toneCurve = {}; - labCurve = LCurveParams(); + labCurve = {}; - rgbCurves = RGBCurvesParams(); + rgbCurves = {}; - localContrast = LocalContrastParams(); + localContrast = {}; - colorToning = ColorToningParams(); + colorToning = {}; - sharpenEdge = SharpenEdgeParams(); + sharpenEdge = {}; - sharpenMicro = SharpenMicroParams(); + sharpenMicro = {}; - sharpening = SharpeningParams(); + sharpening = {}; - prsharpening = SharpeningParams(); + prsharpening = {}; prsharpening.contrast = 15.0; prsharpening.method = "rld"; prsharpening.deconvamount = 100; @@ -2787,69 +2787,69 @@ void ProcParams::setDefaults() prsharpening.deconviter = 100; prsharpening.deconvdamping = 0; - vibrance = VibranceParams(); + vibrance = {}; - wb = WBParams(); + wb = {}; - colorappearance = ColorAppearanceParams(); + colorappearance = {}; - defringe = DefringeParams(); + defringe = {}; - impulseDenoise = ImpulseDenoiseParams(); + impulseDenoise = {}; - dirpyrDenoise = DirPyrDenoiseParams(); + dirpyrDenoise = {}; - epd = EPDParams(); + epd = {}; - fattal = FattalToneMappingParams(); + fattal = {}; - sh = SHParams(); + sh = {}; - crop = CropParams(); + crop = {}; - coarse = CoarseTransformParams(); + coarse = {}; - commonTrans = CommonTransformParams(); + commonTrans = {}; - rotate = RotateParams(); + rotate = {}; - distortion = DistortionParams(); + distortion = {}; - lensProf = LensProfParams(); + lensProf = {}; - perspective = PerspectiveParams(); + perspective = {}; - gradient = GradientParams(); + gradient = {}; - pcvignette = PCVignetteParams(); + pcvignette = {}; - vignetting = VignettingParams(); + vignetting = {}; - chmixer = ChannelMixerParams(); + chmixer = {}; - blackwhite = BlackWhiteParams(); + blackwhite = {}; - cacorrection = CACorrParams(); + cacorrection = {}; - resize = ResizeParams(); + resize = {}; - icm = ColorManagementParams(); + icm = {}; - wavelet = WaveletParams(); + wavelet = {}; - dirpyrequalizer = DirPyrEqualizerParams(); + dirpyrequalizer = {}; - hsvequalizer = HSVEqualizerParams(); + hsvequalizer = {}; - filmSimulation = FilmSimulationParams(); + filmSimulation = {}; - softlight = SoftLightParams(); + softlight = {}; - dehaze = DehazeParams(); + dehaze = {}; - raw = RAWParams(); + raw = {}; - metadata = MetaDataParams(); + metadata = {}; exif.clear(); iptc.clear(); diff --git a/rtengine/rawimage.h b/rtengine/rawimage.h index 4ff6d79c9..ae0bcea84 100644 --- a/rtengine/rawimage.h +++ b/rtengine/rawimage.h @@ -25,82 +25,10 @@ #include "dcraw.h" #include "imageformat.h" -#include "noncopyable.h" namespace rtengine { -struct badPix { - uint16_t x; - uint16_t y; - badPix( uint16_t xc, uint16_t yc ): x(xc), y(yc) {} -}; - -class PixelsMap : - public NonCopyable -{ - int w; // line width in base_t units - int h; // height - typedef unsigned long base_t; - static const size_t base_t_size = sizeof(base_t); - base_t *pm; - -public: - PixelsMap(int width, int height ) - : h(height) - { - w = (width / (base_t_size * 8)) + 1; - pm = new base_t [h * w ]; - memset(pm, 0, h * w * base_t_size ); - } - - ~PixelsMap() - { - delete [] pm; - } - int width() const - { - return w; - } - int height() const - { - return h; - } - - // if a pixel is set returns true - bool get(int x, int y) - { - return (pm[y * w + x / (base_t_size * 8) ] & (base_t)1 << (x % (base_t_size * 8)) ) != 0; - } - - // set a pixel - void set(int x, int y) - { - pm[y * w + x / (base_t_size * 8) ] |= (base_t)1 << (x % (base_t_size * 8)) ; - } - - // set pixels from a list - int set( std::vector &bp) - { - for(std::vector::iterator iter = bp.begin(); iter != bp.end(); ++iter) { - set( iter->x, iter->y); - } - - return bp.size(); - } - - void clear() - { - memset(pm, 0, h * w * base_t_size ); - } - // return 0 if at least one pixel in the word(base_t) is set, otherwise return the number of pixels to skip to the next word base_t - int skipIfZero(int x, int y) - { - return pm[y * w + x / (base_t_size * 8) ] == 0 ? base_t_size * 8 - x % (base_t_size * 8) : 0; - } -}; - - class RawImage: public DCraw { public: diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index b8c422545..3b2738405 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -930,542 +930,6 @@ void RawImageSource::convertColorSpace(Imagefloat* image, const ColorManagementP colorSpaceConversion (image, cmp, wb, pre_mul, embProfile, camProfile, imatrices.xyz_cam, (static_cast(getMetaData()))->getCamera()); } -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -/* interpolateBadPixelsBayer: correct raw pixels looking at the bitmap - * takes into consideration if there are multiple bad pixels in the neighbourhood - */ -int RawImageSource::interpolateBadPixelsBayer( PixelsMap &bitmapBads, array2D &rawData ) -{ - static const float eps = 1.f; - int counter = 0; -#ifdef _OPENMP - #pragma omp parallel for reduction(+:counter) schedule(dynamic,16) -#endif - - for( int row = 2; row < H - 2; row++ ) { - for(int col = 2; col < W - 2; col++ ) { - int sk = bitmapBads.skipIfZero(col, row); //optimization for a stripe all zero - - if( sk ) { - col += sk - 1; //-1 is because of col++ in cycle - continue; - } - - if(!bitmapBads.get(col, row)) { - continue; - } - - float wtdsum = 0.f, norm = 0.f; - - // diagonal interpolation - if(FC(row, col) == 1) { - // green channel. We can use closer pixels than for red or blue channel. Distance to centre pixel is sqrt(2) => weighting is 0.70710678 - // For green channel following pixels will be used for interpolation. Pixel to be interpolated is in centre. - // 1 means that pixel is used in this step, if itself and his counterpart are not marked bad - // 0 0 0 0 0 - // 0 1 0 1 0 - // 0 0 0 0 0 - // 0 1 0 1 0 - // 0 0 0 0 0 - for( int dx = -1; dx <= 1; dx += 2) { - if( bitmapBads.get(col + dx, row - 1) || bitmapBads.get(col - dx, row + 1)) { - continue; - } - - float dirwt = 0.70710678f / ( fabsf( rawData[row - 1][col + dx] - rawData[row + 1][col - dx]) + eps); - wtdsum += dirwt * (rawData[row - 1][col + dx] + rawData[row + 1][col - dx]); - norm += dirwt; - } - } else { - // red and blue channel. Distance to centre pixel is sqrt(8) => weighting is 0.35355339 - // For red and blue channel following pixels will be used for interpolation. Pixel to be interpolated is in centre. - // 1 means that pixel is used in this step, if itself and his counterpart are not marked bad - // 1 0 0 0 1 - // 0 0 0 0 0 - // 0 0 0 0 0 - // 0 0 0 0 0 - // 1 0 0 0 1 - for( int dx = -2; dx <= 2; dx += 4) { - if( bitmapBads.get(col + dx, row - 2) || bitmapBads.get(col - dx, row + 2)) { - continue; - } - - float dirwt = 0.35355339f / ( fabsf( rawData[row - 2][col + dx] - rawData[row + 2][col - dx]) + eps); - wtdsum += dirwt * (rawData[row - 2][col + dx] + rawData[row + 2][col - dx]); - norm += dirwt; - } - } - - // channel independent. Distance to centre pixel is 2 => weighting is 0.5 - // Additionally for all channel following pixels will be used for interpolation. Pixel to be interpolated is in centre. - // 1 means that pixel is used in this step, if itself and his counterpart are not marked bad - // 0 0 1 0 0 - // 0 0 0 0 0 - // 1 0 0 0 1 - // 0 0 0 0 0 - // 0 0 1 0 0 - - // horizontal interpolation - if(!(bitmapBads.get(col - 2, row) || bitmapBads.get(col + 2, row))) { - float dirwt = 0.5f / ( fabsf( rawData[row][col - 2] - rawData[row][col + 2]) + eps); - wtdsum += dirwt * (rawData[row][col - 2] + rawData[row][col + 2]); - norm += dirwt; - } - - // vertical interpolation - if(!(bitmapBads.get(col, row - 2) || bitmapBads.get(col, row + 2))) { - float dirwt = 0.5f / ( fabsf( rawData[row - 2][col] - rawData[row + 2][col]) + eps); - wtdsum += dirwt * (rawData[row - 2][col] + rawData[row + 2][col]); - norm += dirwt; - } - - if (LIKELY(norm > 0.f)) { // This means, we found at least one pair of valid pixels in the steps above, likelihood of this case is about 99.999% - rawData[row][col] = wtdsum / (2.f * norm); //gradient weighted average, Factor of 2.f is an optimization to avoid multiplications in former steps - counter++; - } else { //backup plan -- simple average. Same method for all channels. We could improve this, but it's really unlikely that this case happens - int tot = 0; - float sum = 0; - - for( int dy = -2; dy <= 2; dy += 2) { - for( int dx = -2; dx <= 2; dx += 2) { - if(bitmapBads.get(col + dx, row + dy)) { - continue; - } - - sum += rawData[row + dy][col + dx]; - tot++; - } - } - - if (tot > 0) { - rawData[row][col] = sum / tot; - counter ++; - } - } - } - } - - return counter; // Number of interpolated pixels. -} - -/* interpolateBadPixels3Colours: correct raw pixels looking at the bitmap - * takes into consideration if there are multiple bad pixels in the neighbourhood - */ -int RawImageSource::interpolateBadPixelsNColours( PixelsMap &bitmapBads, const int colours ) -{ - static const float eps = 1.f; - int counter = 0; -#ifdef _OPENMP - #pragma omp parallel for reduction(+:counter) schedule(dynamic,16) -#endif - - for( int row = 2; row < H - 2; row++ ) { - for(int col = 2; col < W - 2; col++ ) { - int sk = bitmapBads.skipIfZero(col, row); //optimization for a stripe all zero - - if( sk ) { - col += sk - 1; //-1 is because of col++ in cycle - continue; - } - - if(!bitmapBads.get(col, row)) { - continue; - } - - float wtdsum[colours], norm[colours]; - - for (int i = 0; i < colours; ++i) { - wtdsum[i] = norm[i] = 0.f; - } - - // diagonal interpolation - for( int dx = -1; dx <= 1; dx += 2) { - if( bitmapBads.get(col + dx, row - 1) || bitmapBads.get(col - dx, row + 1)) { - continue; - } - - for(int c = 0; c < colours; c++) { - float dirwt = 0.70710678f / ( fabsf( rawData[row - 1][(col + dx) * colours + c] - rawData[row + 1][(col - dx) * colours + c]) + eps); - wtdsum[c] += dirwt * (rawData[row - 1][(col + dx) * colours + c] + rawData[row + 1][(col - dx) * colours + c]); - norm[c] += dirwt; - } - } - - // horizontal interpolation - if(!(bitmapBads.get(col - 1, row) || bitmapBads.get(col + 1, row))) { - for(int c = 0; c < colours; c++) { - float dirwt = 1.f / ( fabsf( rawData[row][(col - 1) * colours + c] - rawData[row][(col + 1) * colours + c]) + eps); - wtdsum[c] += dirwt * (rawData[row][(col - 1) * colours + c] + rawData[row][(col + 1) * colours + c]); - norm[c] += dirwt; - } - } - - // vertical interpolation - if(!(bitmapBads.get(col, row - 1) || bitmapBads.get(col, row + 1))) { - for(int c = 0; c < colours; c++) { - float dirwt = 1.f / ( fabsf( rawData[row - 1][col * colours + c] - rawData[row + 1][col * colours + c]) + eps); - wtdsum[c] += dirwt * (rawData[row - 1][col * colours + c] + rawData[row + 1][col * colours + c]); - norm[c] += dirwt; - } - } - - if (LIKELY(norm[0] > 0.f)) { // This means, we found at least one pair of valid pixels in the steps above, likelihood of this case is about 99.999% - for(int c = 0; c < colours; c++) { - rawData[row][col * colours + c] = wtdsum[c] / (2.f * norm[c]); //gradient weighted average, Factor of 2.f is an optimization to avoid multiplications in former steps - } - - counter++; - } else { //backup plan -- simple average. Same method for all channels. We could improve this, but it's really unlikely that this case happens - int tot = 0; - float sum[colours]; - - for (int i = 0; i < colours; ++i) { - sum[i] = 0.f; - } - - for( int dy = -2; dy <= 2; dy += 2) { - for( int dx = -2; dx <= 2; dx += 2) { - if(bitmapBads.get(col + dx, row + dy)) { - continue; - } - - for(int c = 0; c < colours; c++) { - sum[c] += rawData[row + dy][(col + dx) * colours + c]; - } - - tot++; - } - } - - if (tot > 0) { - for(int c = 0; c < colours; c++) { - rawData[row][col * colours + c] = sum[c] / tot; - } - - counter ++; - } - } - } - } - - return counter; // Number of interpolated pixels. -} -/* interpolateBadPixelsXtrans: correct raw pixels looking at the bitmap - * takes into consideration if there are multiple bad pixels in the neighbourhood - */ -int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) -{ - static const float eps = 1.f; - int counter = 0; -#ifdef _OPENMP - #pragma omp parallel for reduction(+:counter) schedule(dynamic,16) -#endif - - for( int row = 2; row < H - 2; row++ ) { - for(int col = 2; col < W - 2; col++ ) { - int skip = bitmapBads.skipIfZero(col, row); //optimization for a stripe all zero - - if( skip ) { - col += skip - 1; //-1 is because of col++ in cycle - continue; - } - - if(!bitmapBads.get(col, row)) { - continue; - } - - float wtdsum = 0.f, norm = 0.f; - unsigned int pixelColor = ri->XTRANSFC(row, col); - - if(pixelColor == 1) { - // green channel. A green pixel can either be a solitary green pixel or a member of a 2x2 square of green pixels - if(ri->XTRANSFC(row, col - 1) == ri->XTRANSFC(row, col + 1)) { - // If left and right neighbour have same colour, then this is a solitary green pixel - // For these the following pixels will be used for interpolation. Pixel to be interpolated is in centre and marked with a P. - // Pairs of pixels used in this step are numbered. A pair will be used if none of the pixels of the pair is marked bad - // 0 means, the pixel has a different colour and will not be used - // 0 1 0 2 0 - // 3 5 0 6 4 - // 0 0 P 0 0 - // 4 6 0 5 3 - // 0 2 0 1 0 - for( int dx = -1; dx <= 1; dx += 2) { // pixels marked 5 or 6 in above example. Distance to P is sqrt(2) => weighting is 0.70710678f - if( bitmapBads.get(col + dx, row - 1) || bitmapBads.get(col - dx, row + 1)) { - continue; - } - - float dirwt = 0.70710678f / ( fabsf( rawData[row - 1][col + dx] - rawData[row + 1][col - dx]) + eps); - wtdsum += dirwt * (rawData[row - 1][col + dx] + rawData[row + 1][col - dx]); - norm += dirwt; - } - - for( int dx = -1; dx <= 1; dx += 2) { // pixels marked 1 or 2 on above example. Distance to P is sqrt(5) => weighting is 0.44721359f - if( bitmapBads.get(col + dx, row - 2) || bitmapBads.get(col - dx, row + 2)) { - continue; - } - - float dirwt = 0.44721359f / ( fabsf( rawData[row - 2][col + dx] - rawData[row + 2][col - dx]) + eps); - wtdsum += dirwt * (rawData[row - 2][col + dx] + rawData[row + 2][col - dx]); - norm += dirwt; - } - - for( int dx = -2; dx <= 2; dx += 4) { // pixels marked 3 or 4 on above example. Distance to P is sqrt(5) => weighting is 0.44721359f - if( bitmapBads.get(col + dx, row - 1) || bitmapBads.get(col - dx, row + 1)) { - continue; - } - - float dirwt = 0.44721359f / ( fabsf( rawData[row - 1][col + dx] - rawData[row + 1][col - dx]) + eps); - wtdsum += dirwt * (rawData[row - 1][col + dx] + rawData[row + 1][col - dx]); - norm += dirwt; - } - } else { - // this is a member of a 2x2 square of green pixels - // For these the following pixels will be used for interpolation. Pixel to be interpolated is at position P in the example. - // Pairs of pixels used in this step are numbered. A pair will be used if none of the pixels of the pair is marked bad - // 0 means, the pixel has a different colour and will not be used - // 1 0 0 3 - // 0 P 2 0 - // 0 2 1 0 - // 3 0 0 0 - - // pixels marked 1 in above example. Distance to P is sqrt(2) => weighting is 0.70710678f - int offset1 = ri->XTRANSFC(row - 1, col - 1) == ri->XTRANSFC(row + 1, col + 1) ? 1 : -1; - - if( !(bitmapBads.get(col - offset1, row - 1) || bitmapBads.get(col + offset1, row + 1))) { - float dirwt = 0.70710678f / ( fabsf( rawData[row - 1][col - offset1] - rawData[row + 1][col + offset1]) + eps); - wtdsum += dirwt * (rawData[row - 1][col - offset1] + rawData[row + 1][col + offset1]); - norm += dirwt; - } - - // pixels marked 2 in above example. Distance to P is 1 => weighting is 1.f - int offsety = (ri->XTRANSFC(row - 1, col) != 1 ? 1 : -1); - int offsetx = offset1 * offsety; - - if( !(bitmapBads.get(col + offsetx, row) || bitmapBads.get(col, row + offsety))) { - float dirwt = 1.f / ( fabsf( rawData[row][col + offsetx] - rawData[row + offsety][col]) + eps); - wtdsum += dirwt * (rawData[row][col + offsetx] + rawData[row + offsety][col]); - norm += dirwt; - } - - int offsety2 = -offsety; - int offsetx2 = -offsetx; - offsetx *= 2; - offsety *= 2; - - // pixels marked 3 in above example. Distance to P is sqrt(5) => weighting is 0.44721359f - if( !(bitmapBads.get(col + offsetx, row + offsety2) || bitmapBads.get(col + offsetx2, row + offsety))) { - float dirwt = 0.44721359f / ( fabsf( rawData[row + offsety2][col + offsetx] - rawData[row + offsety][col + offsetx2]) + eps); - wtdsum += dirwt * (rawData[row + offsety2][col + offsetx] + rawData[row + offsety][col + offsetx2]); - norm += dirwt; - } - } - } else { - // red and blue channel. - // Each red or blue pixel has exactly one neighbour of same colour in distance 2 and four neighbours of same colour which can be reached by a move of a knight in chess. - // For the distance 2 pixel (marked with an X) we generate a virtual counterpart (marked with a V) - // For red and blue channel following pixels will be used for interpolation. Pixel to be interpolated is in centre and marked with a P. - // Pairs of pixels used in this step are numbered except for distance 2 pixels which are marked X and V. A pair will be used if none of the pixels of the pair is marked bad - // 0 1 0 0 0 0 0 X 0 0 remaining cases are symmetric - // 0 0 0 0 2 1 0 0 0 2 - // X 0 P 0 V 0 0 P 0 0 - // 0 0 0 0 1 0 0 0 0 0 - // 0 2 0 0 0 0 2 V 1 0 - - // Find two knight moves landing on a pixel of same colour as the pixel to be interpolated. - // If we look at first and last row of 5x5 square, we will find exactly two knight pixels. - // Additionally we know that the column of this pixel has 1 or -1 horizontal distance to the centre pixel - // When we find a knight pixel, we get its counterpart, which has distance (+-3,+-3), where the signs of distance depend on the corner of the found knight pixel. - // These pixels are marked 1 or 2 in above examples. Distance to P is sqrt(5) => weighting is 0.44721359f - // The following loop simply scans the four possible places. To keep things simple, it does not stop after finding two knight pixels, because it will not find more than two - for(int d1 = -2, offsety = 3; d1 <= 2; d1 += 4, offsety -= 6) { - for(int d2 = -1, offsetx = 3; d2 < 1; d2 += 2, offsetx -= 6) { - if(ri->XTRANSFC(row + d1, col + d2) == pixelColor) { - if( !(bitmapBads.get(col + d2, row + d1) || bitmapBads.get(col + d2 + offsetx, row + d1 + offsety))) { - float dirwt = 0.44721359f / ( fabsf( rawData[row + d1][col + d2] - rawData[row + d1 + offsety][col + d2 + offsetx]) + eps); - wtdsum += dirwt * (rawData[row + d1][col + d2] + rawData[row + d1 + offsety][col + d2 + offsetx]); - norm += dirwt; - } - } - } - } - - // now scan for the pixel of same colour in distance 2 in each direction (marked with an X in above examples). - bool distance2PixelFound = false; - int dx, dy; - - // check horizontal - for(dx = -2, dy = 0; dx <= 2 && !distance2PixelFound; dx += 4) - if(ri->XTRANSFC(row, col + dx) == pixelColor) { - distance2PixelFound = true; - break; - } - - if(!distance2PixelFound) - - // no distance 2 pixel on horizontal, check vertical - for(dx = 0, dy = -2; dy <= 2 && !distance2PixelFound; dy += 4) - if(ri->XTRANSFC(row + dy, col) == pixelColor) { - distance2PixelFound = true; - break; - } - - // calculate the value of its virtual counterpart (marked with a V in above examples) - float virtualPixel; - - if(dy == 0) { - virtualPixel = 0.5f * (rawData[row - 1][col - dx] + rawData[row + 1][col - dx]); - } else { - virtualPixel = 0.5f * (rawData[row - dy][col - 1] + rawData[row - dy][col + 1]); - } - - // and weight as usual. Distance to P is 2 => weighting is 0.5f - float dirwt = 0.5f / ( fabsf( virtualPixel - rawData[row + dy][col + dx]) + eps); - wtdsum += dirwt * (virtualPixel + rawData[row + dy][col + dx]); - norm += dirwt; - } - - if (LIKELY(norm > 0.f)) { // This means, we found at least one pair of valid pixels in the steps above, likelihood of this case is about 99.999% - rawData[row][col] = wtdsum / (2.f * norm); //gradient weighted average, Factor of 2.f is an optimization to avoid multiplications in former steps - counter++; - } - } - } - - return counter; // Number of interpolated pixels. -} -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -/* Search for hot or dead pixels in the image and update the map - * For each pixel compare its value to the average of similar colour surrounding - * (Taken from Emil Martinec idea) - * (Optimized by Ingo Weyrich 2013 and 2015) - */ -int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thresh, bool findHotPixels, bool findDeadPixels ) -{ - float varthresh = (20.0 * (thresh / 100.0) + 1.0 ) / 24.f; - - // allocate temporary buffer - float* cfablur; - cfablur = (float (*)) malloc (H * W * sizeof * cfablur); - - // counter for dead or hot pixels - int counter = 0; - -#ifdef _OPENMP - #pragma omp parallel -#endif - { -#ifdef _OPENMP - #pragma omp for schedule(dynamic,16) nowait -#endif - - for (int i = 2; i < H - 2; i++) { - for (int j = 2; j < W - 2; j++) { - const float& temp = median(rawData[i - 2][j - 2], rawData[i - 2][j], rawData[i - 2][j + 2], - rawData[i][j - 2], rawData[i][j], rawData[i][j + 2], - rawData[i + 2][j - 2], rawData[i + 2][j], rawData[i + 2][j + 2]); - cfablur[i * W + j] = rawData[i][j] - temp; - } - } - - // process borders. Former version calculated the median using mirrored border which does not make sense because the original pixel loses weight - // Setting the difference between pixel and median for border pixels to zero should do the job not worse then former version -#ifdef _OPENMP - #pragma omp single -#endif - { - for(int i = 0; i < 2; i++) { - for(int j = 0; j < W; j++) { - cfablur[i * W + j] = 0.f; - } - } - - for(int i = 2; i < H - 2; i++) { - for(int j = 0; j < 2; j++) { - cfablur[i * W + j] = 0.f; - } - - for(int j = W - 2; j < W; j++) { - cfablur[i * W + j] = 0.f; - } - } - - for(int i = H - 2; i < H; i++) { - for(int j = 0; j < W; j++) { - cfablur[i * W + j] = 0.f; - } - } - } -#ifdef _OPENMP - #pragma omp barrier // barrier because of nowait clause above - - #pragma omp for reduction(+:counter) schedule(dynamic,16) -#endif - - //cfa pixel heat/death evaluation - for (int rr = 2; rr < H - 2; rr++) { - int rrmWpcc = rr * W + 2; - - for (int cc = 2; cc < W - 2; cc++, rrmWpcc++) { - //evaluate pixel for heat/death - float pixdev = cfablur[rrmWpcc]; - - if(pixdev == 0.f) { - continue; - } - - if((!findDeadPixels) && pixdev < 0) { - continue; - } - - if((!findHotPixels) && pixdev > 0) { - continue; - } - - pixdev = fabsf(pixdev); - float hfnbrave = -pixdev; - -#ifdef __SSE2__ - // sum up 5*4 = 20 values using SSE - // 10 fabs function calls and float 10 additions with SSE - vfloat sum = vabsf(LVFU(cfablur[(rr - 2) * W + cc - 2])) + vabsf(LVFU(cfablur[(rr - 1) * W + cc - 2])); - sum += vabsf(LVFU(cfablur[(rr) * W + cc - 2])); - sum += vabsf(LVFU(cfablur[(rr + 1) * W + cc - 2])); - sum += vabsf(LVFU(cfablur[(rr + 2) * W + cc - 2])); - // horizontally add the values and add the result to hfnbrave - hfnbrave += vhadd(sum); - - // add remaining 5 values of last column - for (int mm = rr - 2; mm <= rr + 2; mm++) { - hfnbrave += fabsf(cfablur[mm * W + cc + 2]); - } - -#else - - // 25 fabs function calls and 25 float additions without SSE - for (int mm = rr - 2; mm <= rr + 2; mm++) { - for (int nn = cc - 2; nn <= cc + 2; nn++) { - hfnbrave += fabsf(cfablur[mm * W + nn]); - } - } - -#endif - - if (pixdev > varthresh * hfnbrave) { - // mark the pixel as "bad" - bpMap.set(cc, rr); - counter++; - } - }//end of pixel evaluation - } - }//end of parallel processing - free (cfablur); - return counter; -} - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - void RawImageSource::getFullSize (int& w, int& h, int tr) { @@ -1741,23 +1205,13 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le printf( "Subtracting Darkframe:%s\n", rid->get_filename().c_str()); } - PixelsMap *bitmapBads = nullptr; + std::unique_ptr bitmapBads; int totBP = 0; // Hold count of bad pixels to correct if(ri->zeroIsBad()) { // mark all pixels with value zero as bad, has to be called before FF and DF. dcraw sets this flag only for some cameras (mainly Panasonic and Leica) - bitmapBads = new PixelsMap(W, H); -#ifdef _OPENMP - #pragma omp parallel for reduction(+:totBP) schedule(dynamic,16) -#endif - - for(int i = 0; i < H; i++) - for(int j = 0; j < W; j++) { - if(ri->data[i][j] == 0.f) { - bitmapBads->set(j, i); - totBP++; - } - } + bitmapBads.reset(new PixelsMap(W, H)); + totBP = findZeroPixels(*(bitmapBads.get())); if( settings->verbose) { printf( "%d pixels with value zero marked as bad pixels\n", totBP); @@ -1821,10 +1275,10 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le if( bp ) { if(!bitmapBads) { - bitmapBads = new PixelsMap(W, H); + bitmapBads.reset(new PixelsMap(W, H)); } - totBP += bitmapBads->set( *bp ); + totBP += bitmapBads->set(*bp); if( settings->verbose ) { std::cout << "Correcting " << bp->size() << " pixels from .badpixels" << std::endl; @@ -1842,10 +1296,10 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le if(bp) { if(!bitmapBads) { - bitmapBads = new PixelsMap(W, H); + bitmapBads.reset(new PixelsMap(W, H)); } - totBP += bitmapBads->set( *bp ); + totBP += bitmapBads->set(*bp); if( settings->verbose && !bp->empty()) { std::cout << "Correcting " << bp->size() << " hotpixels from darkframe" << std::endl; @@ -1917,10 +1371,10 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le } if(!bitmapBads) { - bitmapBads = new PixelsMap(W, H); + bitmapBads.reset(new PixelsMap(W, H)); } - int nFound = findHotDeadPixels( *bitmapBads, raw.hotdeadpix_thresh, raw.hotPixelFilter, raw.deadPixelFilter ); + int nFound = findHotDeadPixels(*(bitmapBads.get()), raw.hotdeadpix_thresh, raw.hotPixelFilter, raw.deadPixelFilter ); totBP += nFound; if( settings->verbose && nFound > 0) { @@ -1932,10 +1386,10 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le PDAFLinesFilter f(ri); if (!bitmapBads) { - bitmapBads = new PixelsMap(W, H); + bitmapBads.reset(new PixelsMap(W, H)); } - int n = f.mark(rawData, *bitmapBads); + int n = f.mark(rawData, *(bitmapBads.get())); totBP += n; if (n > 0) { @@ -1999,15 +1453,15 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le if ( ri->getSensorType() == ST_BAYER ) { if(numFrames == 4) { for(int i = 0; i < 4; ++i) { - interpolateBadPixelsBayer( *bitmapBads, *rawDataFrames[i] ); + interpolateBadPixelsBayer(*(bitmapBads.get()), *rawDataFrames[i]); } } else { - interpolateBadPixelsBayer( *bitmapBads, rawData ); + interpolateBadPixelsBayer(*(bitmapBads.get()), rawData); } } else if ( ri->getSensorType() == ST_FUJI_XTRANS ) { - interpolateBadPixelsXtrans( *bitmapBads ); + interpolateBadPixelsXtrans(*(bitmapBads.get())); } else { - interpolateBadPixelsNColours( *bitmapBads, ri->get_colors() ); + interpolateBadPixelsNColours(*(bitmapBads.get()), ri->get_colors()); } } @@ -2060,10 +1514,6 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le printf("Preprocessing: %d usec\n", t2.etime(t1)); } - if(bitmapBads) { - delete bitmapBads; - } - rawDirty = true; return; } diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h index c04da5b66..b33b19c08 100644 --- a/rtengine/rawimagesource.h +++ b/rtengine/rawimagesource.h @@ -28,6 +28,7 @@ #include "dcp.h" #include "iimage.h" #include "imagesource.h" +#include "pixelsmap.h" #define HR_SCALE 2 @@ -100,7 +101,7 @@ protected: void transformRect (const PreviewProps &pp, int tran, int &sx1, int &sy1, int &width, int &height, int &fw); void transformPosition (int x, int y, int tran, int& tx, int& ty); - unsigned FC(int row, int col) + unsigned FC(int row, int col) const { return ri->FC(row, col); } @@ -258,11 +259,11 @@ protected: ); void ddct8x8s(int isgn, float a[8][8]); - int interpolateBadPixelsBayer( PixelsMap &bitmapBads, array2D &rawData ); - int interpolateBadPixelsNColours( PixelsMap &bitmapBads, const int colours ); - int interpolateBadPixelsXtrans( PixelsMap &bitmapBads ); - int findHotDeadPixels( PixelsMap &bpMap, float thresh, bool findHotPixels, bool findDeadPixels ); - + int interpolateBadPixelsBayer(const PixelsMap &bitmapBads, array2D &rawData); + int interpolateBadPixelsNColours(const PixelsMap &bitmapBads, int colours); + int interpolateBadPixelsXtrans(const PixelsMap &bitmapBads); + int findHotDeadPixels(PixelsMap &bpMap, float thresh, bool findHotPixels, bool findDeadPixels) const; + int findZeroPixels(PixelsMap &bpMap) const; void cfa_linedn (float linenoiselevel, bool horizontal, bool vertical, const CFALineDenoiseRowBlender &rowblender);//Emil's line denoise void green_equilibrate_global (array2D &rawData);