← Back to branch summary

~ubuntu-branches/ubuntu/karmic/digikam/karmic-backports

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Ubuntu Archive Auto-Backport
  • Date: 2009-12-07 19:03:53 UTC
  • mfrom: (54.1.4 lucid)
  • Revision ID: james.westby@ubuntu.com-20091207190353-oara3lenjxymto3i
Tags: 2:1.0.0~rc-1ubuntu1~karmic1
Automated backport upload; no source changes.

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