← Back to branch summary

~ubuntu-branches/ubuntu/saucy/digikam/saucy

« back to all changes in this revision

Viewing changes to libs/dimg/filters/waveletsnr.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christian Mangold
  • Date: 2010-04-09 21:30:01 UTC
  • mfrom: (1.2.28 upstream)
  • Revision ID: james.westby@ubuntu.com-20100409213001-4bfyibrd359rn7o3
Tags: 2:1.2.0-0ubuntu1
https://launchpad.net/bugs/560576
* New upstream release (LP: #560576)
* Remove all patches, fixed upstream
  - Remove quilt build-depend

Show diffs side-by-side

added added

removed removed

Lines of Context:
libs/dimg/filters/waveletsnr.cpp
1
 
/* ============================================================
2
 
 *
3
 
 * This file is a part of digiKam project
4
 
 * http://www.digikam.org
5
 
 *
6
 
 * Date        : 2005-05-25
7
 
 * Description : Wavelets Noise Reduction threaded image filter.
8
 
 *               This filter work in YCrCb color space.
9
 
 *
10
 
 * Copyright (C) 2005-2009 by Gilles Caulier <caulier dot gilles at gmail dot com>
11
 
 * Copyright (C) 2008 by Marco Rossini <marco dot rossini at gmx dot net>
12
 
 *
13
 
 * This program is free software; you can redistribute it
14
 
 * and/or modify it under the terms of the GNU General
15
 
 * Public License as published by the Free Software Foundation;
16
 
 * either version 2, or (at your option)
17
 
 * any later version.
18
 
 *
19
 
 * This program is distributed in the hope that it will be useful,
20
 
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21
 
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
22
 
 * GNU General Public License for more details.
23
 
 *
24
 
 * ============================================================ */
25
 
 
26
 
#include "waveletsnr.h"
27
 
 
28
 
// C++ includes
29
 
 
30
 
#include <cmath>
31
 
 
32
 
// Local includes
33
 
 
34
 
#include "dimg.h"
35
 
#include "dcolor.h"
36
 
#include "dimgimagefilters.h"
37
 
 
38
 
namespace Digikam
39
 
{
40
 
 
41
 
class WaveletsNRPriv
42
 
{
43
 
public:
44
 
 
45
 
    WaveletsNRPriv()
46
 
    {
47
 
        for (int c = 0 ; c < 3; c++)
48
 
        {
49
 
            fimg[c]   = 0;
50
 
            buffer[c] = 0;
51
 
        }
52
 
    }
53
 
 
54
 
    float*              fimg[3];
55
 
    float*              buffer[3];
56
 
 
57
 
    WaveletsNRContainer settings;  
58
 
};
59
 
  
60
 
WaveletsNR::WaveletsNR(DImg *orgImage, QObject *parent, const WaveletsNRContainer& settings)
61
 
          : DImgThreadedFilter(orgImage, parent, "WaveletsNR"),
62
 
            d(new WaveletsNRPriv)
63
 
{
64
 
    d->settings = settings;
65
 
 
66
 
    initFilter();
67
 
}
68
 
 
69
 
WaveletsNR::~WaveletsNR()
70
 
{
71
 
    delete d;
72
 
}
73
 
 
74
 
void WaveletsNR::filterImage()
75
 
{
76
 
    DColor col;
77
 
    int    progress;
78
 
 
79
 
    int width  = m_orgImage.width();
80
 
    int height = m_orgImage.height();
81
 
    float clip = m_orgImage.sixteenBit() ? 65535.0 : 255.0;
82
 
 
83
 
    // Allocate buffers.
84
 
 
85
 
    for (int c = 0; c < 3; c++)
86
 
        d->fimg[c] = new float[width * height];
87
 
 
88
 
    d->buffer[1] = new float[width * height];
89
 
    d->buffer[2] = new float[width * height];
90
 
 
91
 
    // Read the full image and convert pixel values to float [0,1].
92
 
 
93
 
    int j = 0;
94
 
 
95
 
    for (int y = 0; !m_cancel && (y < height); y++)
96
 
    {
97
 
        for (int x = 0; !m_cancel && (x < width); x++)
98
 
        {
99
 
            col           = m_orgImage.getPixelColor(x, y);
100
 
            d->fimg[0][j] = col.red()   / clip;
101
 
            d->fimg[1][j] = col.green() / clip;
102
 
            d->fimg[2][j] = col.blue()  / clip;
103
 
            j++;
104
 
        }
105
 
    }
106
 
 
107
 
    postProgress( 10 );
108
 
 
109
 
    // do colour model conversion sRGB[0,1] -> YCrCb.
110
 
 
111
 
    if (!m_cancel)
112
 
    {
113
 
        srgb2ycbcr(d->fimg, width * height);
114
 
    }
115
 
 
116
 
    postProgress( 20 );
117
 
 
118
 
    // denoise the channels individually
119
 
 
120
 
    for (int c = 0; !m_cancel && (c < 3); c++)
121
 
    {
122
 
        d->buffer[0] = d->fimg[c];
123
 
 
124
 
        if (d->settings.thresholds[c] > 0.0)
125
 
        {
126
 
            waveletDenoise(d->buffer, width, height, d->settings.thresholds[c], d->settings.softness[c]);
127
 
 
128
 
            progress = (int) (30.0 + ((double)c * 60.0) / 4);
129
 
            if ( progress%5 == 0 )
130
 
                postProgress( progress );
131
 
        }
132
 
    }
133
 
 
134
 
    // Retransform the image data to sRGB[0,1].
135
 
 
136
 
    if (!m_cancel)
137
 
    {
138
 
        ycbcr2srgb(d->fimg, width * height);
139
 
    }
140
 
 
141
 
    postProgress( 80 );
142
 
 
143
 
    // clip the values
144
 
 
145
 
    for (int c = 0; !m_cancel && (c < 3); c++)
146
 
    {
147
 
        for (int i = 0; i < width * height; i++)
148
 
        {
149
 
            d->fimg[c][i] = qBound(0.0F, d->fimg[c][i] * clip, clip);
150
 
        }
151
 
    }
152
 
 
153
 
    postProgress( 90 );
154
 
 
155
 
    // Write back the full image and convert pixel values from float [0,1].
156
 
 
157
 
    j = 0;
158
 
 
159
 
    for (int y = 0; !m_cancel && (y < height); y++)
160
 
    {
161
 
        for (int x = 0; x < width; x++)
162
 
        {
163
 
            col.setRed(   (int)(d->fimg[0][j] + 0.5) );
164
 
            col.setGreen( (int)(d->fimg[1][j] + 0.5) );
165
 
            col.setBlue(  (int)(d->fimg[2][j] + 0.5) );
166
 
            col.setAlpha( m_orgImage.getPixelColor(x, y).alpha() );
167
 
            j++;
168
 
 
169
 
            m_destImage.setPixelColor(x, y, col);
170
 
        }
171
 
    }
172
 
 
173
 
    postProgress( 100 );
174
 
 
175
 
    // Free buffers.
176
 
 
177
 
    for (int c = 0; c < 3; c++)
178
 
        delete [] d->fimg[c];
179
 
 
180
 
    delete [] d->buffer[1];
181
 
    delete [] d->buffer[2];
182
 
}
183
 
 
184
 
// -- Wavelets denoise methods -----------------------------------------------------------
185
 
 
186
 
void WaveletsNR::waveletDenoise(float* fimg[3], unsigned int width, unsigned int height,
187
 
                                float threshold, double softness)
188
 
{
189
 
    float        *temp=0, thold;
190
 
    unsigned int i, lev, lpass, hpass, size, col, row;
191
 
    double       stdev[5];
192
 
    unsigned int samples[5];
193
 
 
194
 
    size  = width * height;
195
 
    temp  = new float[qMax(width, height)];
196
 
    hpass = 0;
197
 
 
198
 
    for (lev = 0; !m_cancel && (lev < 5); lev++)
199
 
    {
200
 
        lpass = ((lev & 1) + 1);
201
 
        for (row = 0; !m_cancel && (row < height); row++)
202
 
        {
203
 
            hatTransform(temp, fimg[hpass] + row * width, 1, width, 1 << lev);
204
 
            for (col = 0; col < width; col++)
205
 
            {
206
 
                fimg[lpass][row * width + col] = temp[col] * 0.25;
207
 
            }
208
 
        }
209
 
 
210
 
        for (col = 0; !m_cancel && (col < width); col++)
211
 
        {
212
 
            hatTransform(temp, fimg[lpass] + col, width, height, 1 << lev);
213
 
            for (row = 0; row < height; row++)
214
 
            {
215
 
                fimg[lpass][row * width + col] = temp[row] * 0.25;
216
 
            }
217
 
        }
218
 
 
219
 
        thold = 5.0 / (1 << 6) * exp (-2.6 * sqrt (lev + 1.0)) * 0.8002 / exp (-2.6);
220
 
 
221
 
        // initialize stdev values for all intensities
222
 
 
223
 
        stdev[0]   = stdev[1]   = stdev[2]   = stdev[3]   = stdev[4]   = 0.0;
224
 
        samples[0] = samples[1] = samples[2] = samples[3] = samples[4] = 0;
225
 
 
226
 
        // calculate stdevs for all intensities
227
 
 
228
 
        for (i = 0; !m_cancel && (i < size); i++)
229
 
        {
230
 
            fimg[hpass][i] -= fimg[lpass][i];
231
 
 
232
 
            if (fimg[hpass][i] < thold && fimg[hpass][i] > -thold)
233
 
            {
234
 
                if (fimg[lpass][i] > 0.8)
235
 
                {
236
 
                    stdev[4] += fimg[hpass][i] * fimg[hpass][i];
237
 
                    samples[4]++;
238
 
                }
239
 
                else if (fimg[lpass][i] > 0.6)
240
 
                {
241
 
                    stdev[3] += fimg[hpass][i] * fimg[hpass][i];
242
 
                    samples[3]++;
243
 
                }
244
 
                else if (fimg[lpass][i] > 0.4)
245
 
                {
246
 
                    stdev[2] += fimg[hpass][i] * fimg[hpass][i];
247
 
                    samples[2]++;
248
 
                }
249
 
                else if (fimg[lpass][i] > 0.2)
250
 
                {
251
 
                    stdev[1] += fimg[hpass][i] * fimg[hpass][i];
252
 
                    samples[1]++;
253
 
                }
254
 
                else
255
 
                {
256
 
                    stdev[0] += fimg[hpass][i] * fimg[hpass][i];
257
 
                    samples[0]++;
258
 
                }
259
 
            }
260
 
        }
261
 
 
262
 
        stdev[0] = sqrt(stdev[0] / (samples[0] + 1));
263
 
        stdev[1] = sqrt(stdev[1] / (samples[1] + 1));
264
 
        stdev[2] = sqrt(stdev[2] / (samples[2] + 1));
265
 
        stdev[3] = sqrt(stdev[3] / (samples[3] + 1));
266
 
        stdev[4] = sqrt(stdev[4] / (samples[4] + 1));
267
 
 
268
 
        // do thresholding
269
 
 
270
 
        for (i = 0; !m_cancel && (i < size); i++)
271
 
        {
272
 
            if (fimg[lpass][i] > 0.8)
273
 
            {
274
 
                thold = threshold * stdev[4];
275
 
            }
276
 
            else if (fimg[lpass][i] > 0.6)
277
 
            {
278
 
                thold = threshold * stdev[3];
279
 
            }
280
 
            else if (fimg[lpass][i] > 0.4)
281
 
            {
282
 
                thold = threshold * stdev[2];
283
 
            }
284
 
            else if (fimg[lpass][i] > 0.2)
285
 
            {
286
 
                thold = threshold * stdev[1];
287
 
            }
288
 
            else
289
 
            {
290
 
                thold = threshold * stdev[0];
291
 
            }
292
 
 
293
 
            if (fimg[hpass][i] < -thold)
294
 
                fimg[hpass][i] += thold - thold * softness;
295
 
            else if (fimg[hpass][i] > thold)
296
 
                fimg[hpass][i] -= thold - thold * softness;
297
 
            else
298
 
                fimg[hpass][i] *= softness;
299
 
 
300
 
            if (hpass)
301
 
                fimg[0][i] += fimg[hpass][i];
302
 
        }
303
 
 
304
 
        hpass = lpass;
305
 
    }
306
 
 
307
 
    for (i = 0; !m_cancel && (i < size); i++)
308
 
        fimg[0][i] = fimg[0][i] + fimg[lpass][i];
309
 
 
310
 
    delete [] temp;
311
 
}
312
 
 
313
 
void WaveletsNR::hatTransform (float* temp, float* base, int st, int size, int sc)
314
 
{
315
 
    int i;
316
 
 
317
 
    for (i = 0; i < sc; i++)
318
 
      temp[i] = 2 * base[st * i] + base[st * (sc - i)] + base[st * (i + sc)];
319
 
 
320
 
    for (; i + sc < size; i++)
321
 
      temp[i] = 2 * base[st * i] + base[st * (i - sc)] + base[st * (i + sc)];
322
 
 
323
 
    for (; i < size; i++)
324
 
      temp[i] = 2 * base[st * i] + base[st * (i - sc)] + base[st * (2 * size - 2 - (i + sc))];
325
 
}
326
 
 
327
 
// -- Color Space conversion methods --------------------------------------------------
328
 
 
329
 
void WaveletsNR::srgb2ycbcr(float** fimg, int size)
330
 
{
331
 
    /* using JPEG conversion here - expecting all channels to be in [0:255] range */
332
 
    float y, cb, cr;
333
 
 
334
 
    for (int i = 0; i < size; i++)
335
 
    {
336
 
        y          =  0.2990 * fimg[0][i] + 0.5870 * fimg[1][i] + 0.1140 * fimg[2][i];
337
 
        cb         = -0.1687 * fimg[0][i] - 0.3313 * fimg[1][i] + 0.5000 * fimg[2][i] + 0.5;
338
 
        cr         =  0.5000 * fimg[0][i] - 0.4187 * fimg[1][i] - 0.0813 * fimg[2][i] + 0.5;
339
 
        fimg[0][i] = y;
340
 
        fimg[1][i] = cb;
341
 
        fimg[2][i] = cr;
342
 
    }
343
 
}
344
 
 
345
 
void WaveletsNR::ycbcr2srgb(float** fimg, int size)
346
 
{
347
 
    /* using JPEG conversion here - expecting all channels to be in [0:255] range */
348
 
    float r, g, b;
349
 
 
350
 
    for (int i = 0; i < size; i++)
351
 
    {
352
 
        r          = fimg[0][i] + 1.40200 * (fimg[2][i] - 0.5);
353
 
        g          = fimg[0][i] - 0.34414 * (fimg[1][i] - 0.5) - 0.71414 * (fimg[2][i] - 0.5);
354
 
        b          = fimg[0][i] + 1.77200 * (fimg[1][i] - 0.5);
355
 
        fimg[0][i] = r;
356
 
        fimg[1][i] = g;
357
 
        fimg[2][i] = b;
358
 
    }
359
 
}
360
 
 
361
 
void WaveletsNR::srgb2xyz(float** fimg, int size)
362
 
{
363
 
    /* fimg in [0:1], sRGB */
364
 
    float x, y, z;
365
 
 
366
 
    for (int i = 0; i < size; i++)
367
 
    {
368
 
        /* scaling and gamma correction (approximate) */
369
 
        fimg[0][i] = pow(fimg[0][i], (float)2.2);
370
 
        fimg[1][i] = pow(fimg[1][i], (float)2.2);
371
 
        fimg[2][i] = pow(fimg[2][i], (float)2.2);
372
 
 
373
 
        /* matrix RGB -> XYZ, with D65 reference white (www.brucelindbloom.com) */
374
 
        x = 0.412424  * fimg[0][i] + 0.357579 * fimg[1][i] + 0.180464  * fimg[2][i];
375
 
        y = 0.212656  * fimg[0][i] + 0.715158 * fimg[1][i] + 0.0721856 * fimg[2][i];
376
 
        z = 0.0193324 * fimg[0][i] + 0.119193 * fimg[1][i] + 0.950444  * fimg[2][i];
377
 
 
378
 
        /*
379
 
        x = 0.412424 * fimg[0][i] + 0.212656  * fimg[1][i] + 0.0193324 * fimg[2][i];
380
 
        y = 0.357579 * fimg[0][i] + 0.715158  * fimg[1][i] + 0.119193  * fimg[2][i];
381
 
        z = 0.180464 * fimg[0][i] + 0.0721856 * fimg[1][i] + 0.950444  * fimg[2][i];
382
 
        */
383
 
 
384
 
        fimg[0][i] = x;
385
 
        fimg[1][i] = y;
386
 
        fimg[2][i] = z;
387
 
    }
388
 
}
389
 
 
390
 
void WaveletsNR::xyz2srgb(float** fimg, int size)
391
 
{
392
 
    float r, g, b;
393
 
 
394
 
    for (int i = 0; i < size; i++)
395
 
    {
396
 
        /* matrix RGB -> XYZ, with D65 reference white (www.brucelindbloom.com) */
397
 
        r = 3.24071   * fimg[0][i] - 1.53726  * fimg[1][i] - 0.498571  * fimg[2][i];
398
 
        g = -0.969258 * fimg[0][i] + 1.87599  * fimg[1][i] + 0.0415557 * fimg[2][i];
399
 
        b = 0.0556352 * fimg[0][i] - 0.203996 * fimg[1][i] + 1.05707   * fimg[2][i];
400
 
 
401
 
        /*
402
 
        r =  3.24071  * fimg[0][i] - 0.969258  * fimg[1][i]
403
 
          + 0.0556352 * fimg[2][i];
404
 
        g = -1.53726  * fimg[0][i] + 1.87599   * fimg[1][i]
405
 
          - 0.203996  * fimg[2][i];
406
 
        b = -0.498571 * fimg[0][i] + 0.0415557 * fimg[1][i]
407
 
          + 1.05707   * fimg[2][i];
408
 
        */
409
 
 
410
 
        /* scaling and gamma correction (approximate) */
411
 
        r = r < 0 ? 0 : pow(r, (float)(1.0 / 2.2));
412
 
        g = g < 0 ? 0 : pow(g, (float)(1.0 / 2.2));
413
 
        b = b < 0 ? 0 : pow(b, (float)(1.0 / 2.2));
414
 
 
415
 
        fimg[0][i] = r;
416
 
        fimg[1][i] = g;
417
 
        fimg[2][i] = b;
418
 
    }
419
 
}
420
 
 
421
 
void WaveletsNR::lab2srgb(float** fimg, int size)
422
 
{
423
 
    float x, y, z;
424
 
 
425
 
    for (int i = 0; i < size; i++)
426
 
    {
427
 
        /* convert back to normal LAB */
428
 
        fimg[0][i] = (fimg[0][i] - 0 * 16 * 27 / 24389.0) * 116;
429
 
        fimg[1][i] = (fimg[1][i] - 0.5) * 500 * 2;
430
 
        fimg[2][i] = (fimg[2][i] - 0.5) * 200 * 2.2;
431
 
 
432
 
        /* matrix */
433
 
        y = (fimg[0][i] + 16) / 116;
434
 
        z = y - fimg[2][i] / 200.0;
435
 
        x = fimg[1][i] / 500.0 + y;
436
 
 
437
 
        /* scale */
438
 
        if (x * x * x > 216 / 24389.0)
439
 
            x = x * x * x;
440
 
        else
441
 
            x = (116 * x - 16) * 27 / 24389.0;
442
 
        if (fimg[0][i] > 216 / 27.0)
443
 
            y = y * y * y;
444
 
        else
445
 
            /*y = fimg[0][i] * 27 / 24389.0;*/
446
 
            y = (116 * y - 16) * 27 / 24389.0;
447
 
        if (z * z * z > 216 / 24389.0)
448
 
            z = z * z * z;
449
 
        else
450
 
            z = (116 * z - 16) * 27 / 24389.0;
451
 
 
452
 
        /* white reference */
453
 
        fimg[0][i] = x * 0.95047;
454
 
        fimg[1][i] = y;
455
 
        fimg[2][i] = z * 1.08883;
456
 
    }
457
 
 
458
 
    xyz2srgb(fimg, size);
459
 
}
460
 
 
461
 
void WaveletsNR::srgb2lab(float** fimg, int size)
462
 
{
463
 
    float l, a, b;
464
 
 
465
 
    srgb2xyz(fimg, size);
466
 
 
467
 
    for (int i = 0; i < size; i++)
468
 
    {
469
 
        /* reference white */
470
 
        fimg[0][i] /= 0.95047;
471
 
        /* (just for completeness)
472
 
        fimg[1][i] /= 1.00000; */
473
 
        fimg[2][i] /= 1.08883;
474
 
 
475
 
        /* scale */
476
 
        if (fimg[0][i] > 216.0 / 24389.0)
477
 
        {
478
 
            fimg[0][i] = pow(fimg[0][i], (float)(1.0 / 3.0));
479
 
        }
480
 
        else
481
 
        {
482
 
            fimg[0][i] = (24389.0 * fimg[0][i] / 27.0 + 16.0) / 116.0;
483
 
        }
484
 
 
485
 
        if (fimg[1][i] > 216.0 / 24389.0)
486
 
        {
487
 
            fimg[1][i] = pow(fimg[1][i], (float)(1.0 / 3.0));
488
 
        }
489
 
        else
490
 
        {
491
 
            fimg[1][i] = (24389 * fimg[1][i] / 27.0 + 16.0) / 116.0;
492
 
        }
493
 
 
494
 
        if (fimg[2][i] > 216.0 / 24389.0)
495
 
        {
496
 
            fimg[2][i] = (float)pow(fimg[2][i], (float)(1.0 / 3.0));
497
 
        }
498
 
        else
499
 
        {
500
 
            fimg[2][i] = (24389.0 * fimg[2][i] / 27.0 + 16.0) / 116.0;
501
 
        }
502
 
 
503
 
        l          = 116 * fimg[1][i]  - 16;
504
 
        a          = 500 * (fimg[0][i] - fimg[1][i]);
505
 
        b          = 200 * (fimg[1][i] - fimg[2][i]);
506
 
        fimg[0][i] = l / 116.0; /* + 16 * 27 / 24389.0; */
507
 
        fimg[1][i] = a / 500.0 / 2.0 + 0.5;
508
 
        fimg[2][i] = b / 200.0 / 2.2 + 0.5;
509
 
 
510
 
        if (fimg[0][i] < 0)
511
 
            fimg[0][i] = 0;
512
 
    }
513
 
}
514
 
 
515
 
}  // namespace Digikam