1 /*
2 * This file is part of RawTherapee.
3 *
4 * RawTherapee is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * RawTherapee is distributed in the hope that it will be useful,
10 * but widthITheightOUT ANY widthARRANTY; without even the implied warranty of
11 * MERCheightANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with RawTherapee. If not, see <http://www.gnu.org/licenses/>.
16 *
17 * 2010 Emil Martinec <ejmartin@uchicago.edu>
18 *
19 */
20 #include <cstddef>
21 #include "rt_math.h"
22 #include "labimage.h"
23 #include "improcfun.h"
24 #include "sleef.h"
25 #include "opthelper.h"
26 #include "gauss.h"
27 #include "rt_algo.h"
28
29 using namespace std;
30
31 namespace rtengine {
32
33 namespace {
34
impulse_nr(Imagefloat * lab,double thresh,bool multithread)35 void impulse_nr(Imagefloat *lab, double thresh, bool multithread)
36 {
37 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
38 // impulse noise removal
39 // local variables
40
41 const int width = lab->getWidth();
42 const int height = lab->getHeight();
43 float **lab_L = lab->g.ptrs;
44 float **lab_a = lab->r.ptrs;
45 float **lab_b = lab->b.ptrs;
46
47 // buffer for the lowpass image
48 // float * lpf[height] ALIGNED16;
49 // lpf[0] = new float [width * height];
50 // buffer for the highpass image
51 char * impish[height] ALIGNED16;
52 impish[0] = new char [width * height];
53
54 for (int i = 1; i < height; i++) {
55 // lpf[i] = lpf[i - 1] + width;
56 impish[i] = impish[i - 1] + width;
57 }
58
59 markImpulse(width, height, lab_L, impish, thresh);
60
61 // //The cleaning algorithm starts here
62
63 // //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64 // // modified bilateral filter for lowpass image, omitting input pixel; or Gaussian blur
65
66 const float eps = 1.0;
67
68 //now impulsive values have been identified
69
70 // Issue 1671:
71 // often, noise isn't evenly distributed, e.g. only a few noisy pixels in the bright sky, but many in the dark foreground,
72 // so it's better to schedule dynamic and let every thread only process 16 rows, to avoid running big threads out of work
73 // Measured it and in fact gives better performance than without schedule(dynamic,16). Of course, there could be a better
74 // choice for the chunk_size than 16
75 // race conditions are avoided by the array impish
76 #ifdef _OPENMP
77 # pragma omp parallel if (multithread)
78 #endif
79 {
80 int i1, j1, j;
81 float wtdsum[3], dirwt, norm;
82 #ifdef _OPENMP
83 #pragma omp for schedule(dynamic,16)
84 #endif
85
86 for (int i = 0; i < height; i++) {
87 for (j = 0; j < 2; j++) {
88 if (!impish[i][j]) {
89 continue;
90 }
91
92 norm = 0.0;
93 wtdsum[0] = wtdsum[1] = wtdsum[2] = 0.0;
94
95 for (i1 = max(0, i - 2); i1 <= min(i + 2, height - 1); i1++ )
96 for (j1 = 0; j1 <= j + 2; j1++ ) {
97 if (impish[i1][j1]) {
98 continue;
99 }
100
101 dirwt = 1 / (SQR(lab_L[i1][j1] - lab_L[i][j]) + eps); //use more sophisticated rangefn???
102 wtdsum[0] += dirwt * lab_L[i1][j1];
103 wtdsum[1] += dirwt * lab_a[i1][j1];
104 wtdsum[2] += dirwt * lab_b[i1][j1];
105 norm += dirwt;
106 }
107
108 if (norm) {
109 lab_L[i][j] = wtdsum[0] / norm; //low pass filter
110 lab_a[i][j] = wtdsum[1] / norm; //low pass filter
111 lab_b[i][j] = wtdsum[2] / norm; //low pass filter
112 }
113 }
114
115 for (; j < width - 2; j++) {
116 if (!impish[i][j]) {
117 continue;
118 }
119
120 norm = 0.0;
121 wtdsum[0] = wtdsum[1] = wtdsum[2] = 0.0;
122
123 for (i1 = max(0, i - 2); i1 <= min(i + 2, height - 1); i1++ )
124 for (j1 = j - 2; j1 <= j + 2; j1++ ) {
125 if (impish[i1][j1]) {
126 continue;
127 }
128
129 dirwt = 1 / (SQR(lab_L[i1][j1] - lab_L[i][j]) + eps); //use more sophisticated rangefn???
130 wtdsum[0] += dirwt * lab_L[i1][j1];
131 wtdsum[1] += dirwt * lab_a[i1][j1];
132 wtdsum[2] += dirwt * lab_b[i1][j1];
133 norm += dirwt;
134 }
135
136 if (norm) {
137 lab_L[i][j] = wtdsum[0] / norm; //low pass filter
138 lab_a[i][j] = wtdsum[1] / norm; //low pass filter
139 lab_b[i][j] = wtdsum[2] / norm; //low pass filter
140 }
141 }
142
143 for (; j < width; j++) {
144 if (!impish[i][j]) {
145 continue;
146 }
147
148 norm = 0.0;
149 wtdsum[0] = wtdsum[1] = wtdsum[2] = 0.0;
150
151 for (i1 = max(0, i - 2); i1 <= min(i + 2, height - 1); i1++ )
152 for (j1 = j - 2; j1 < width; j1++ ) {
153 if (impish[i1][j1]) {
154 continue;
155 }
156
157 dirwt = 1 / (SQR(lab_L[i1][j1] - lab_L[i][j]) + eps); //use more sophisticated rangefn???
158 wtdsum[0] += dirwt * lab_L[i1][j1];
159 wtdsum[1] += dirwt * lab_a[i1][j1];
160 wtdsum[2] += dirwt * lab_b[i1][j1];
161 norm += dirwt;
162 }
163
164 if (norm) {
165 lab_L[i][j] = wtdsum[0] / norm; //low pass filter
166 lab_a[i][j] = wtdsum[1] / norm; //low pass filter
167 lab_b[i][j] = wtdsum[2] / norm; //low pass filter
168 }
169 }
170 }
171 }
172 //now impulsive values have been corrected
173
174 // delete [] lpf[0];
175 delete [] impish[0];
176 }
177
178 } // namespace
179
180
impulsedenoise(Imagefloat * rgb)181 void ImProcFunctions::impulsedenoise(Imagefloat *rgb)
182 {
183
184 if (params->impulseDenoise.enabled && rgb->getWidth() >= 8 && rgb->getHeight() >= 8 && scale >= 0.5)
185
186 {
187 float t = params->impulseDenoise.thresh / scale;
188 rgb->setMode(Imagefloat::Mode::LAB, multiThread);
189 impulse_nr(rgb, t / 20.0, multiThread);
190 }
191 }
192
193 } // namespace rtengine
194