diff options
Diffstat (limited to 'src/imageplugins/noisereduction/noisereduction.cpp')
| -rw-r--r-- | src/imageplugins/noisereduction/noisereduction.cpp | 809 |
1 files changed, 809 insertions, 0 deletions
diff --git a/src/imageplugins/noisereduction/noisereduction.cpp b/src/imageplugins/noisereduction/noisereduction.cpp new file mode 100644 index 00000000..b9be4729 --- /dev/null +++ b/src/imageplugins/noisereduction/noisereduction.cpp @@ -0,0 +1,809 @@ +/* ============================================================ + * + * This file is a part of digiKam project + * http://www.digikam.org + * + * Date : 2005-05-25 + * Description : Noise Reduction threaded image filter. + * + * Copyright (C) 2005-2007 by Gilles Caulier <caulier dot gilles at gmail dot com> + * + * Original Noise Filter algorithm copyright (C) 2005 + * Peter Heckert <peter dot heckert at arcor dot de> + * from dcamnoise2 gimp plugin available at this url : + * http://home.arcor.de/peter.heckert/dcamnoise2-0.63.c + * + * This program 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 2, or (at your option) + * any later version. + * + * This program 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. + * + * ============================================================ */ + +#define IIR1(dest,src) (dest) = (d3 = ((((src) * b + d3) * b3 + d2) * b2 + d1) * b1) +#define IIR2(dest,src) (dest) = (d2 = ((((src) * b + d2) * b3 + d1) * b2 + d3) * b1) +#define IIR3(dest,src) (dest) = (d1 = ((((src) * b + d1) * b3 + d3) * b2 + d2) * b1) + +#define IIR1A(dest,src) (dest) = fabs(d3 = ((((src) * b + d3) * b3 + d2) * b2 + d1) * b1) +#define IIR2A(dest,src) (dest) = fabs(d2 = ((((src) * b + d2) * b3 + d1) * b2 + d3) * b1) +#define IIR3A(dest,src) (dest) = fabs(d1 = ((((src) * b + d1) * b3 + d3) * b2 + d2) * b1) + +#define FR 0.212671 +#define FG 0.715160 +#define FB 0.072169 + +// Local includes. + +#include "ddebug.h" +#include "dimg.h" +#include "dimgimagefilters.h" +#include "noisereduction.h" + +namespace DigikamNoiseReductionImagesPlugin +{ + +NoiseReduction::NoiseReduction(Digikam::DImg *orgImage, TQObject *parent, + double radius, double lsmooth, double effect, double texture, double sharp, + double csmooth, double lookahead, double gamma, double damping, double phase) + : Digikam::DImgThreadedFilter(orgImage, parent, "NoiseReduction") +{ + m_radius = radius; /* default radius default = 1.0 */ + m_sharp = sharp; /* Sharpness factor default = 0.25 */ + m_lsmooth = lsmooth; /* Luminance Tolerance default = 1.0 */ + m_effect = effect; /* Adaptive filter-effect threshold default = 0.08 */ + m_texture = texture; /* Texture Detail default = 0.0 */ + + m_csmooth = csmooth; /* RGB Tolerance default = 1.0 */ + m_lookahead = lookahead; /* Lookahead default = 2.0 */ + m_gamma = gamma; /* Filter gamma default = 1.0 */ + m_damping = damping; /* Phase jitter Damping default = 5.0 */ + m_phase = phase; /* Area Noise Clip default = 1.0 */ + + m_iir.B = 0.0; + m_iir.b1 = 0.0; + m_iir.b2 = 0.0; + m_iir.b3 = 0.0; + m_iir.b0 = 0.0; + m_iir.r = 0.0; + m_iir.q = 0.0; + m_iir.p = 0; + + m_clampMax = m_orgImage.sixteenBit() ? 65535 : 255; + + initFilter(); +} + +// Remove noise on the region, given a source region, dest. +// region, width and height of the regions, and corner coordinates of +// a subregion to act upon. Everything outside the subregion is unaffected. + +void NoiseReduction::filterImage(void) +{ + int bytes = m_orgImage.bytesDepth(); // Bytes per pixel sample + uchar *srcPR = m_orgImage.bits(); + uchar *destPR = m_destImage.bits(); + int width = m_orgImage.width(); + int height = m_orgImage.height(); + + int row, col, i, progress; + float prob = 0.0; + + int w = (int)((m_radius + m_lookahead + m_damping + m_phase) * 4.0 + 40.0); + + // NOTE: commented from original implementation + // if (radius < m_lookahead) w = m_lookahead * 4.0 + 40.0; + + float csmooth = m_csmooth; + + // Raw Filter preview + + if (csmooth >= 0.99) csmooth = 1.0; + + // Allocate and init buffers + + uchar *src = new uchar[ TQMAX (width, height) * bytes ]; + uchar *dest = new uchar[ TQMAX (width, height) * bytes ]; + float *data = new float[ TQMAX (width, height) + 2*w ]; + float *data2 = new float[ TQMAX (width, height) + 2*w ]; + float *buffer = new float[ TQMAX (width, height) + 2*w ]; + float *rbuf = new float[ TQMAX (width, height) + 2*w ]; + float *tbuf = new float[ TQMAX (width, height) + 2*w ]; + + memset (src, 0, TQMAX (width, height) * bytes); + memset (dest, 0, TQMAX (width, height) * bytes); + + for (i=0 ; i < TQMAX(width,height)+2*w-1 ; i++) + data[i] = data2[i] = buffer[i] = rbuf[i] = tbuf[i] = 0.0; + + // Initialize the damping filter coefficients + + iir_init(m_radius); + + // blur the rows + + for (row = 0 ; !m_cancel && (row < height) ; row++) + { + memcpy(src, srcPR + row*width*bytes, width*bytes); + memcpy(dest, src, width*bytes); + + blur_line (data+w, data2+w, buffer+w, rbuf+w, tbuf+w, src, dest, width); + + memcpy(destPR + row*width*bytes, dest, width*bytes); + + progress = (int)(((double)row * 20.0) / height); + if ( progress%2 == 0 ) + postProgress( progress ); + } + + // blur the cols + + for (col = 0 ; !m_cancel && (col < width) ; col++) + { + for (int n = 0 ; n < height ; n++) + memcpy(src + n*bytes, destPR + (col + width*n)*bytes, bytes); + + for (int n = 0 ; n < height ; n++) + memcpy(dest + n*bytes, srcPR + (col + width*n)*bytes, bytes); + + blur_line (data+w, data2+w, buffer+w, rbuf+w, tbuf+w, src, dest, height); + + for (int n = 0 ; n < height ; n++) + memcpy(destPR + (col + width*n)*bytes, dest + n*bytes, bytes); + + progress = (int)(20.0 + ((double)col * 20.0) / width); + if ( progress%2 == 0 ) + postProgress( progress ); + } + + // merge the source and destination (which currently contains + // the blurred version) images + + for (row = 0 ; !m_cancel && (row < height) ; row++) + { + uchar *s = src; + uchar *d = dest; + unsigned short *s16 = (unsigned short *)src; + unsigned short *d16 = (unsigned short *)dest; + float value; + int u, v; + + // get source row + + memcpy(src, srcPR + row*width*bytes, width*bytes); + memcpy(dest, destPR + row*width*bytes, width*bytes); + + // get dest row and combine the two + + float t = m_csmooth; + float t2 = m_lsmooth; + + // Values are squared, so that sliders get a nonlinear chracteristic + // for better adjustment accuracy when values are small. + t*=t; + t2*=t2; + + for (u = 0 ; !m_cancel && (u < width) ; u++) + { + float dpix[3], spix[3]; + float lum, red, green, blue; + float lum2, red2, green2, blue2; + + if (m_orgImage.sixteenBit()) // 16 bits image + { + red = (float) s16[2]/(float)m_clampMax; + green = (float) s16[1]/(float)m_clampMax; + blue = (float) s16[0]/(float)m_clampMax; + } + else // 8 bits image + { + red = (float) s[2]/(float)m_clampMax; + green = (float) s[1]/(float)m_clampMax; + blue = (float) s[0]/(float)m_clampMax; + } + + spix[2] = red; + spix[1] = green; + spix[0] = blue; + + lum = (FR*red + FG*green + FB*blue); + + if (m_orgImage.sixteenBit()) // 16 bits image + { + red2 = (float) d16[2]/(float)m_clampMax; + green2 = (float) d16[1]/(float)m_clampMax; + blue2 = (float) d16[0]/(float)m_clampMax; + } + else // 8 bits image + { + red2 = (float) d[2]/(float)m_clampMax; + green2 = (float) d[1]/(float)m_clampMax; + blue2 = (float) d[0]/(float)m_clampMax; + } + + lum2 = (FR*red2 + FG*green2 + FB*blue2); + + // Calculate luminance error (contrast error) for filtered template. + // This error is biggest, where edges are. Edges anyway cannot be filtered. + // Therefore we can correct luminance error in edges without increasing noise. + // Should be adjusted carefully, or not so carefully if you intentionally want to add noise. + // Noise, if not colorized, /can/ look good, so this makes sense. + + float dl = lum - lum2; + + // Multiply dl with first derivative of gamma curve divided by derivative value for midtone 0.5 + // So bright tones will be corrected more (get more luminance noise and -information) than + // darker values because bright parts of image generally are less noisy, this is what we want. + + dl *= pow(lum2/0.5, m_gamma-1.0); + + if (t2 > 0.0) + dl *= (1.0 - exp(-dl*dl/(2.0*t2*t2))); + + // NOTE: commented from original implementation + // if (dl > p) dl = p; + // if (dl < -p) dl = -p; + + dpix[2] = red2 + dl; + dpix[1] = green2 + dl; + dpix[0] = blue2 + dl; + + for (v = 0 ; !m_cancel && (v < 3) ; v++) + { + float value = spix[v]; + float fvalue = dpix[v]; + float mvalue = (value + fvalue)/2.0; + float diff = (value) - (fvalue); + + // Multiply diff with first derivative of gamma curve divided by derivative value for midtone 0.5 + // So midtones will stay unchanged, darker values get more blur and brighter values get less blur + // when we increase gamma. + + diff *= pow(mvalue/0.5, m_gamma-1.0); + + // Calculate noise probability for pixel + // TODO : probably it is not probability but an arbitrary curve. + // Probably we should provide a GUI-interface for this!!! + + if (t > 0.0) + prob = exp(-diff*diff/(2.0*t*t)); + else + prob = 0.0; + + // Allow viewing of raw filter output + + if (t >= 0.99) + prob = 1.0; + + dpix[v] = value = fvalue * prob + value * (1.0 - prob); + } + + if (m_orgImage.sixteenBit()) // 16 bits image + { + value = dpix[0]*(float)m_clampMax+0.5; + d16[0] = (unsigned short)CLAMP(value, 0, m_clampMax); + value = dpix[1]*(float)m_clampMax+0.5; + d16[1] = (unsigned short)CLAMP(value, 0, m_clampMax); + value = dpix[2]*(float)m_clampMax+0.5; + d16[2] = (unsigned short)CLAMP(value, 0, m_clampMax); + + d16 += 4; + s16 += 4; + } + else // 8 bits image + { + value = dpix[0]*(float)m_clampMax+0.5; + d[0] = (uchar)CLAMP(value, 0, m_clampMax); + value = dpix[1]*(float)m_clampMax+0.5; + d[1] = (uchar)CLAMP(value, 0, m_clampMax); + value = dpix[2]*(float)m_clampMax+0.5; + d[2] = (uchar)CLAMP(value, 0, m_clampMax); + + d += 4; + s += 4; + } + } + + memcpy(destPR + row*width*bytes, dest, width*bytes); + + progress = (int)(40.0 + ((double)row * 60.0) / height); + if ( progress%2 == 0 ) + postProgress( progress ); + } + + delete [] data; + delete [] data2; + delete [] buffer; + delete [] rbuf; + delete [] tbuf; + delete [] dest; + delete [] src; +} + +// This function is written as if it is blurring a column at a time, +// even though it can operate on rows, too. There is no difference +// in the processing of the lines, at least to the blur_line function. +// 'len' is the length of src and dest + +void NoiseReduction::blur_line(float* const data, float* const data2, float* const buffer, + float* rbuf, float* tbuf, const uchar *src, uchar *dest, int len) +{ + int b; + int row; + int idx; + + unsigned short *src16 = (unsigned short *)src; + unsigned short *dest16 = (unsigned short *)dest; + + // Calculate radius factors + + for (row = 0, idx = 0 ; !m_cancel && (idx < len) ; row += 4, idx++) + { + // Color weigths are chosen proportional to Bayer Sensor pixel count + + if (m_orgImage.sixteenBit()) // 16 bits image + { + data[idx] = (float) dest16[row+2] / (float)m_clampMax * 0.25; // Red color + data[idx] += (float) dest16[row+1] / (float)m_clampMax * 0.5; // Green color + data[idx] += (float) dest16[row] / (float)m_clampMax * 0.25; // Blue color + data[idx] = mypow(data[idx], m_gamma); + } + else // 8 bits image + { + data[idx] = (float) dest[row+2] / (float)m_clampMax * 0.25; // Red color + data[idx] += (float) dest[row+1] / (float)m_clampMax * 0.5; // Green color + data[idx] += (float) dest[row] / (float)m_clampMax * 0.25; // Blue color + data[idx] = mypow(data[idx], m_gamma); + } + } + + filter(data, data2, buffer, rbuf, tbuf, len, -1); + + // Do actual filtering + + for (b = 0 ; !m_cancel && (b < 3) ; b++) + { + for (row = b, idx = 0 ; !m_cancel && (idx < len) ; row += 4, idx++) + { + if (m_orgImage.sixteenBit()) // 16 bits image + data[idx] = (float)src16[row] / (float)m_clampMax; + else // 8 bits image + data[idx] = (float)src[row] / (float)m_clampMax; + } + + filter(data, data2, buffer, rbuf, tbuf, len, b); + + for (row = b, idx = 0 ; !m_cancel && (idx < len) ; row += 4, idx++) + { + int value = (int)(data[idx] * (float)m_clampMax + 0.5); + + if (m_orgImage.sixteenBit()) // 16 bits image + dest16[row] = (unsigned short)CLAMP( value, 0, m_clampMax); + else // 8 bits image + dest[row] = (uchar)CLAMP( value, 0, m_clampMax); + } + } +} + +void NoiseReduction::iir_init(double r) +{ + if (m_iir.r == r) + return; + + // damping settings; + m_iir.r = r; + + double q; + + if ( r >= 2.5) + q = 0.98711 * r - 0.96330; + else + q = 3.97156 - 4.14554 * sqrt(1.0 - 0.26891 * r); + + m_iir.q = q; + m_iir.b0 = 1.57825 + ((0.422205 * q + 1.4281) * q + 2.44413) * q; + m_iir.b1 = ((1.26661 * q +2.85619) * q + 2.44413) * q / m_iir.b0; + m_iir.b2 = - ((1.26661*q +1.4281) * q * q ) / m_iir.b0; + m_iir.b3 = 0.422205 * q * q * q / m_iir.b0; + m_iir.B = 1.0 - (m_iir.b1 + m_iir.b2 + m_iir.b3); +} + +void NoiseReduction::box_filter(double *src, double *end, double *dest, double radius) +{ + int boxwidth = 1; + float box = (*src); + float fbw = 2.0 * radius; + + if (fbw < 1.0) + fbw = 1.0; + + while(boxwidth+2 <= (int) fbw) boxwidth+=2, box += (src[boxwidth/2]) + (src[-boxwidth/2]); + + double frac = (fbw - (double) boxwidth) / 2.0; + int bh = boxwidth / 2; + int bh1 = boxwidth / 2+1; + + for ( ; src <= end ; src++, dest++) + { + *dest = (box + frac * ((src[bh1])+(src[-bh1]))) / fbw; + box = box - (src[-bh]) + (src[bh1]); + } +} + +// Bidirectional IIR-filter, speed optimized + +void NoiseReduction::iir_filter(float* const start, float* const end, float* dstart, + double radius, const int type) +{ + if (!dstart) + dstart = start; + + int width; + float *src = start; + float *dest = dstart; + float *dend = dstart + (end - start); + + radius = floor((radius + 0.1) / 0.5) * 0.5; + + // NOTE: commented from original implementation + // gfloat boxwidth = radius * 2.0; + // gint bw = (gint) boxwidth; + + int ofs = (int)radius; + if (ofs < 1) ofs = 1; + + double d1, d2, d3; + + width = end - start + 1; + + if (radius < 0.25) + { + if ( start != dest ) + { + memcpy(dest, start, width*sizeof(*dest)); + return; + } + } + + iir_init(radius); + + const double b1 = m_iir.b1; + const double b2 = m_iir.b2 / m_iir.b1; + const double b3 = m_iir.b3 / m_iir.b2; + const double b = m_iir.B / m_iir.b3; + + switch(type) + { + case Gaussian: + + d1 = d2 = d3 = *dest; + dend -= 6; + src--; + dest--; + + while (dest < dend) + { + IIR1(*(++dest), *(++src)); + IIR2(*(++dest), *(++src)); + IIR3(*(++dest), *(++src)); + IIR1(*(++dest), *(++src)); + IIR2(*(++dest), *(++src)); + IIR3(*(++dest), *(++src)); + } + + dend += 6; + + while (1) + { + if (++dest > dend) break; + IIR1(*dest,*(++src)); + if (++dest > dend) break; + IIR2(*dest,*(++src)); + if (++dest > dend) break; + IIR3(*dest,*(++src)); + } + + d1 = d2 = d3 = dest[-1]; + dstart += 6; + + while (dest > dstart) + { + --dest, IIR1(*dest, *dest); + --dest, IIR2(*dest, *dest); + --dest, IIR3(*dest, *dest); + --dest, IIR1(*dest, *dest); + --dest, IIR2(*dest, *dest); + --dest, IIR3(*dest, *dest); + } + + dstart -= 6; + + while (1) + { + if (--dest < dstart) break; + IIR1(*dest, *dest); + if (--dest < dstart) break; + IIR2(*dest, *dest); + if (--dest < dstart) break; + IIR3(*dest, *dest); + } + + break; + + case SecondDerivative: // rectified and filtered second derivative, source and dest may be equal + + d1 = d2 = d3 = 0.0; + dest[0] = dest[ofs] = 0.0; + dend -= 6; + dest--; + src--; + + while (dest < dend) + { + ++src, IIR1(*(++dest), src[ofs]-src[0]); + ++src, IIR2(*(++dest), src[ofs]-src[0]); + ++src, IIR3(*(++dest), src[ofs]-src[0]); + ++src, IIR1(*(++dest), src[ofs]-src[0]); + ++src, IIR2(*(++dest), src[ofs]-src[0]); + ++src, IIR3(*(++dest), src[ofs]-src[0]); + } + + dend += 6; + + while (1) + { + if (++dest > dend) break; + ++src, IIR1(*dest, src[ofs]-src[0]); + if (++dest > dend) break; + ++src, IIR2(*dest, src[ofs]-src[0]); + if (++dest > dend) break; + ++src, IIR3(*dest, src[ofs]-src[0]); + } + + d1 = d2 = d3 = 0.0; + dest[-1] = dest[-ofs-1] = 0.0; + dstart += 6; + + while (dest > dstart) + { + --dest, IIR1A(*dest, dest[0]-dest[-ofs]); + --dest, IIR2A(*dest, dest[0]-dest[-ofs]); + --dest, IIR3A(*dest, dest[0]-dest[-ofs]); + --dest, IIR1A(*dest, dest[0]-dest[-ofs]); + --dest, IIR2A(*dest, dest[0]-dest[-ofs]); + --dest, IIR3A(*dest, dest[0]-dest[-ofs]); + } + + dstart -= 6; + + while (1) + { + if (--dest < dstart) break; + IIR1A(*dest, dest[0]-dest[-ofs]); + if (--dest < dstart) break; + IIR2A(*dest, dest[0]-dest[-ofs]); + if (--dest < dstart) break; + IIR3A(*dest, dest[0]-dest[-ofs]); + } + + break; + } +} + +// A forward-backward box filter is used here and the radius is adapted to luminance jump. +// Radius is calculated fron 1st and 2nd derivative of intensity values. +// (Its not exactly 2nd derivative, but something similar, optimized by experiment) +// The radius variations are filtered. This reduces spatial phase jitter. + +void NoiseReduction::filter(float *buffer, float *data, float *data2, float *rbuf, + float */*tbuf*/, int width, int color) +{ + float *lp = data; + float *rp = data + width-1; + float *lp2 = data2; + float *blp = buffer; + float *brp = buffer + width-1; + float *rbuflp = rbuf; + float *rbufrp = rbuf + width-1; + float fboxwidth = m_radius*2.0; + float fradius = m_radius; + float *p1, *p2; + + if (fboxwidth < 1.0) fboxwidth = 1.0 ; + if (fradius < 0.5) fradius = 0.5; + + int i, pass; + int ofs, ofs2; + float maxrad; + float fbw; + float val; + double rfact = m_effect*m_effect; + double sharp = m_sharp; + + ofs2 = (int)floor(m_damping * 2.0 + 0.1); + ofs = (int)floor(m_lookahead * 2.0 + 0.1); + int w = (int)(fboxwidth + m_damping + m_lookahead + m_phase + 2.0); + + // Mirror image edges + + for (i=1 ; i <= w ; i++) + blp[-i] = blp[i]; + + for (i=1 ; i <= w ; i++) + brp[i] = brp[-i]; + + if (color < 0) // Calc 2nd derivative + { + // boost high frequency in rbuf + + for (p1 = blp, p2 = rbuflp ; p1 <= brp ; p1++, p2++) + { + *p2 = (sharp+1.0) * p1[0] - sharp * 0.5 * (p1[-ofs]+p1[ofs]); + } + + iir_filter(rbuflp-w, rbufrp+w, blp-w, m_lookahead, SecondDerivative); + + // Mirror image edges + + for (i = 1 ; i <= w ; i++) + blp[-i] = blp[i]; + + for (i = 1 ; i <= w ; i++) + brp[i] = brp[-i]; + + // boost high frequency in rbuf + + for (p1 = blp, p2 = rbuflp ; p1 <= brp ; p1++, p2++) + { + *p2 = ((sharp+1.0) * (p1[0]) - sharp * 0.5 * ((p1[-ofs2])+(p1[ofs2]))); + } + + // Mirror rbuf edges + + for (i = 1 ; i <= w ; i++) + rbuflp[-i] = rbuflp[i]; + + for (i = 1 ; i <= w ; i++) + rbufrp[i] = rbufrp[-i]; + + // Lowpass (gauss) filter rbuf, remove phase jitter + + iir_filter(rbuflp-w+5, rbufrp+w-5, rbuflp-w+5, m_damping, Gaussian); + + for (i = -w+5; i < width-1+w-5 ; i++) + { + // NOTE: commented from original implementation + // val = rbuflp[i]; + + val = rbuflp[i]-rfact; + + // Avoid division by zero, clip negative filter overshoot + + if (val < rfact/fradius) val=rfact/fradius; + + val = rfact/val; + + // NOTE: commented from original implementation + // val = pow(val/fradius,m_phase)*fradius; + + if (val < 0.5) val = 0.5; + + rbuflp[i] = val*2.0; + } + + // Mirror rbuf edges + + for (i=1 ; i <= w ; i++) + rbuflp[-i] = rbuflp[i]; + + for (i=1 ; i <= w ; i++) + rbufrp[i] = rbufrp[-i]; + + return; + } + + // Calc lowpass filtered input signal + + iir_filter(blp-w+1, brp+w-1, lp2-w+1, m_radius, Gaussian); + + // Subtract low frequency from input signal (aka original image data) + // and predistort this signal + + val = m_texture + 1.0; + + for (i = -w+1 ; i <= width-1+w-1 ; i++) + { + blp[i] = mypow(blp[i] - lp2[i], val); + } + + float *src, *dest; + val = m_texture + 1.0; + + pass = 2; + + while (pass--) + { + float sum; + int ibw; + src = blp; + dest = lp; + maxrad = 0.0; + + // Mirror left edge + + for (i=1 ; i <= w ; i++) + src[-i] = src[i]; + + sum = (src[-1] += src[-2]); + + // forward pass + + for (rbuf = rbuflp-(int) m_phase ; rbuf <= rbufrp; src++, dest++, rbuf++) + { + // NOTE: commented from original implementation + //fbw = fabs( rbuf[-ofs2]*ll2+rbuf[-ofs2-1]*rl2); + + fbw = *rbuf; + + if (fbw > (maxrad += 1.0)) fbw = maxrad; + else if (fbw < maxrad) maxrad = fbw; + + ibw = (int)fbw; + *src = sum += *src; + *dest = (sum-src[-ibw]+(src[-ibw]-src[-ibw-1])*(fbw-ibw))/fbw; + } + + src = rp; + dest = brp; + maxrad = 0.0; + + // Mirror right edge + + for (i=1 ; i <= w ; i++) + src[i] = src[-i]; + + sum = (src[1] += src[2]); + + // backward pass + + for ( rbuf = rbufrp +(int) m_phase ; rbuf >= rbuflp; src--, dest--, rbuf--) + { + // NOTE: commented from original implementation + //fbw = fabs( rbuf[ofs2]*ll2+rbuf[ofs2+1]*rl2); + + fbw = *rbuf; + + if (fbw > (maxrad +=1.0)) fbw = maxrad; + else if (fbw < maxrad) maxrad = fbw; + + ibw = (int)fbw; + + *src = sum += *src; + *dest = (sum-src[ibw]+(src[ibw]-src[ibw+1])*(fbw-ibw))/fbw; + } + } + + val = 1.0 / (m_texture + 1.0); + + for (i = -w+1 ; i <= width-1+w-1 ; i++) + { + // Undo predistortion + + blp[i]= mypow(blp[i],val); + + // Add in low frequency + + blp[i] += lp2[i]; + + // NOTE: commented from original implementation + // if (blp[i] >= 0.0) blp[i] = pow(blp[i],val); + // else blp[i] = 0.0; + } +} + +} // NameSpace DigikamNoiseReductionImagesPlugin |