diff --git a/data/kernels/dwt.cl b/data/kernels/dwt.cl index d25bfce8ccb..64238301640 100644 --- a/data/kernels/dwt.cl +++ b/data/kernels/dwt.cl @@ -58,9 +58,11 @@ dwt_hat_transform_row(global float4 *lpass, global float4 *hpass, int width, int if(x < sc) { const int h_idx = mad24(y, width, x); + // the right neighbour also needs reflection once 2*sc > size + const int rightpos = (x + sc < size) ? (x + sc) : (2 * size - 2 - (x + sc)); lpass[h_idx] = hat_mult * hpass[mad24(y, width, x)] + hpass[mad24(y, width, (sc - x))] - + hpass[mad24(y, width, (x + sc))]; + + hpass[mad24(y, width, rightpos)]; } else if(x + sc < size) { @@ -92,8 +94,11 @@ dwt_hat_transform_col(global float4 *lpass, int width, int height, const int sc, if(y < sc) { + // the lower neighbour also needs reflection once 2*sc > size + const int belowpos = (y + sc < size) ? (y + sc) : (2 * size - 2 - (y + sc)); + temp_buffer[mad24(y, width, x)] = (hat_mult * lpass[mad24(y, width, x)] + lpass[mad24((sc - y), width, x)] - + lpass[mad24((y + sc), width, x)]) + + lpass[mad24(belowpos, width, x)]) * lpass_mult; } else if(y + sc < size) diff --git a/src/common/dwt.c b/src/common/dwt.c index eaf312f7267..b13d249ef32 100644 --- a/src/common/dwt.c +++ b/src/common/dwt.c @@ -151,7 +151,11 @@ static void dwt_decompose_horiz( const size_t width, const size_t lev) { - const int hscale = MIN(1 << lev, width); //(int because we need a signed difference below) + // cap at width-1 (not width): the right-edge reflection below computes + // 2*width-2-(col+hscale), which underflows for the last column when hscale==width. + // this mirrors the vertical pass, which caps vscale at height-1 for the same reason. + const int hscale = + MIN(1 << lev, (int)width - 1); //(int because we need a signed difference below) DT_OMP_FOR() for(int row = 0; row < height ; row++) { @@ -382,7 +386,10 @@ static void dwt_denoise_vert_1ch( const size_t width, const size_t lev) { - const int vscale = MIN(1 << lev, height); + // cap at height-1 (not height): the top-edge reflection below reads + // abs(row-vscale), which lands one row past the buffer at row 0 once + // vscale == height. (matches dwt_decompose_vert.) + const int vscale = MIN(1 << lev, (int)height - 1); DT_OMP_FOR() for(int rowid = 0; rowid < height ; rowid++) { @@ -418,7 +425,9 @@ static void dwt_denoise_horiz_1ch( const float thold, const int last) { - const int hscale = MIN(1 << lev, width); + // cap at width-1 (not width): the reflected neighbour indices below read out + // of bounds once hscale == width. (matches dwt_decompose_horiz.) + const int hscale = MIN(1 << lev, (int)width - 1); DT_OMP_FOR() for(int row = 0; row < height ; row++) { @@ -431,23 +440,31 @@ static void dwt_denoise_horiz_1ch( float *const restrict details = in + rowindex; float *const restrict coarse = out + rowindex; float *const restrict accum_row = accum + rowindex; - // handle reflection at left edge - DT_OMP_SIMD() - for(int col = 0; col < hscale; col++) + // the row splits into three regions: a left edge where the left neighbour + // (and, once 2*hscale > width, the right one too) needs reflection, a + // branch-free middle where both neighbours are in bounds, and a right edge + // where the right neighbour needs reflection. clamp the boundaries with + // MIN/MAX so the regions stay a non-overlapping partition for small images + // / large scales (2*hscale > width); otherwise they would overlap and + // double-process columns, and the middle loop would index out of bounds. + const int left_end = MIN(hscale, (int)width); + const int right_start = MAX((int)width - hscale, left_end); + // left edge: left neighbour reflects about column 0 (abs(col-hscale) == hscale-col + // here); guard the right neighbour explicitly since it can also fall off the end. + for(int col = 0; col < left_end; col++) { + const int rightpos = + (col + hscale < (int)width) ? (col + hscale) : 2 * ((int)width - 1) - (col + hscale); // add up left/center/right, and renormalize by dividing by the total weight of all numbers added together - const float hat = (2.f * coarse[col] + coarse[hscale-col] + coarse[col+hscale]) / 16.f; - // the normalized value is our 'coarse' result; 'diff' is the difference between original input and 'coarse' - // (which would ordinarily be stored as the details scale, but we don't need it any further) + const float hat = (2.f * coarse[col] + coarse[hscale - col] + coarse[rightpos]) / 16.f; + // the normalized value is our 'coarse' result; 'diff' is the difference between original + // input and 'coarse' const float diff = details[col] - hat; - details[col] = hat; // done with original input, so we can overwrite it with 'coarse' - // GCC8 won't vectorize if we use the following line, but it turns out that just adding the two conditional - // alternatives produces exactly the same result, and *that* does get vectorized - //const float excess = diff < 0.0 ? MIN(diff + thold, 0.0f) : MAX(diff - thold, 0.0f); + details[col] = hat; // done with original input, so we can overwrite it with 'coarse' accum_row[col] += MAX(diff - thold,0.0f) + MIN(diff + thold, 0.0f); } DT_OMP_SIMD() - for(int col = hscale; col < width - hscale; col++) + for(int col = left_end; col < right_start; col++) { // add up left/center/right, and renormalize by dividing by the total weight of all numbers added together const float hat = (2.f * coarse[col] + coarse[col-hscale] + coarse[col+hscale]) / 16.f; @@ -460,15 +477,15 @@ static void dwt_denoise_horiz_1ch( //const float excess = diff < 0.0 ? MIN(diff + thold, 0.0f) : MAX(diff - thold, 0.0f); accum_row[col] += MAX(diff - thold,0.0f) + MIN(diff + thold, 0.0f); } - // handle reflection at right edge - DT_OMP_SIMD() - for(int col = width - hscale; col < width; col++) + // right edge: right neighbour reflects about the last column. col >= right_start + // >= hscale here, so the left neighbour (col-hscale) is always in bounds. + for(int col = right_start; col < width; col++) { - const float right = coarse[2*width - 2 - (col+hscale)]; + const float right = coarse[2 * ((int)width - 1) - (col + hscale)]; // add up left/center/right, and renormalize by dividing by the total weight of all numbers added together const float hat = (2.f * coarse[col] + coarse[col-hscale] + right) / 16.f; - // the normalized value is our 'coarse' result; 'diff' is the difference between original input and 'coarse' - // (which would ordinarily be stored as the details scale, but we don't need it any further) + // the normalized value is our 'coarse' result; 'diff' is the difference between original + // input and 'coarse' const float diff = details[col] - hat; details[col] = hat; // done with original input, so we can overwrite it with 'coarse' accum_row[col] += MAX(diff - thold,0.0f) + MIN(diff + thold, 0.0f); @@ -686,7 +703,10 @@ static cl_int dwt_wavelet_decompose_cl(cl_mem img, dwt_params_cl_t *const p, _dw { int sc = 1 << lev; sc = (int)(sc * p->preview_scale); - if(sc > p->width) sc = p->width; + // cap at width-1 (not width): the kernel's reflected neighbour indices + // read out of bounds once sc == width. (matches the CPU path.) + if(sc > p->width - 1) + sc = p->width - 1; err = dt_opencl_enqueue_kernel_2d_args(devid, p->global->kernel_dwt_hat_transform_row, p->width, p->height, CLARG(temp), CLARG((buffer[hpass])), CLARG(p->width), CLARG(p->height), CLARG(sc)); @@ -697,7 +717,9 @@ static cl_int dwt_wavelet_decompose_cl(cl_mem img, dwt_params_cl_t *const p, _dw { int sc = 1 << lev; sc = (int)(sc * p->preview_scale); - if(sc > p->height) sc = p->height; + // cap at height-1 (not height): see the row pass above. + if(sc > p->height - 1) + sc = p->height - 1; const float lpass_mult = (1.f / 16.f); err = dt_opencl_enqueue_kernel_2d_args(devid, p->global->kernel_dwt_hat_transform_col, p->width, p->height,