1 ////////////////////////////////////////////////////////////////
2 //
3 // CFA line denoise by DCT filtering
4 //
5 // copyright (c) 2008-2010 Emil Martinec <ejmartin@uchicago.edu>
6 // parallelized 2013 by Ingo Weyrich
7 //
8 // code dated: June 7, 2010
9 //
10 // cfa_linedn_RT.cc is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // This program is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with this program. If not, see <https://www.gnu.org/licenses/>.
22 //
23 ////////////////////////////////////////////////////////////////
24
25 #include <cmath>
26
27 #include "rawimagesource.h"
28 #include "rt_math.h"
29
30 #define TS 224 // Tile size of 224 instead of 512 speeds up processing
31
32 #define CLASS
33
34 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35
36 using namespace std;
37 using namespace rtengine;
38
39 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40
cfa_linedn(float noise,bool horizontal,bool vertical,const CFALineDenoiseRowBlender & rowblender)41 void RawImageSource::CLASS cfa_linedn(float noise, bool horizontal, bool vertical, const CFALineDenoiseRowBlender &rowblender)
42 {
43 // local variables
44 int height = H, width = W;
45
46 const float clip_pt = 0.8 * initialGain * 65535.0;
47
48 const float eps = 1e-5; //tolerance to avoid dividing by zero
49
50 const float gauss[5] = {0.20416368871516755, 0.18017382291138087, 0.1238315368057753, 0.0662822452863612, 0.02763055063889883};
51 const float rolloff[8] = {0, 0.135335, 0.249352, 0.411112, 0.606531, 0.800737, 0.945959, 1}; //gaussian with sigma=3
52 const float window[8] = {0, .25, .75, 1, 1, .75, .25, 0}; //sine squared
53
54 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
55
56 double progress = 0.0;
57 if (plistener) {
58 plistener->setProgressStr("PROGRESSBAR_LINEDENOISE");
59 plistener->setProgress(progress);
60 }
61
62 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63 float noisevar = SQR(3 * noise * 65535); // _noise_ (as a fraction of saturation) is input to the algorithm
64 float noisevarm4 = 4.0f * noisevar;
65 float* RawDataTmp = (float*)malloc(static_cast<unsigned long>(width) * height * sizeof(float));
66 #ifdef _OPENMP
67 #pragma omp parallel
68 #endif
69 {
70
71 // allocate memory and assure the arrays don't have same 64 byte boundary to avoid L1 conflict misses
72 float *cfain = (float*)malloc(3 * TS * TS * sizeof(float) + 2 * 16 * sizeof(float));
73 float *cfadiff = (cfain + (1 * TS * TS) + 1 * 16);
74 float *cfadn = (cfain + (2 * TS * TS) + 2 * 16);
75 float cfablur[TS];
76
77 float linehvar[4], linevvar[4], noisefactor[4][8][2], coeffsq;
78 float dctblock[4][8][8];
79
80 #ifdef _OPENMP
81 #pragma omp for
82 #endif
83
84 for(int i = 0; i < height; i++)
85 for(int j = 0; j < width; j++) {
86 RawDataTmp[i * width + j] = rawData[i][j];
87 }
88
89 // Main algorithm: Tile loop
90 #ifdef _OPENMP
91 #pragma omp for schedule(dynamic) collapse(2)
92 #endif
93
94 for (int top = 0; top < height - 16; top += TS - 32)
95 for (int left = 0; left < width - 16; left += TS - 32) {
96
97 int bottom = min(top + TS, height);
98 int right = min(left + TS, width);
99 int numrows = bottom - top;
100 int numcols = right - left;
101 int indx1;
102
103 // load CFA data; data should be in linear gamma space, before white balance multipliers are applied
104 for (int rr = top; rr < top + numrows; rr++)
105 for (int cc = left, indx = (rr - top) * TS; cc < left + numcols; cc++, indx++) {
106 cfain[indx] = rawData[rr][cc];
107 }
108
109 //pad the block to a multiple of 16 on both sides
110
111 if (numcols < TS) {
112 indx1 = numcols % 16;
113
114 for (int i = 0; i < (16 - indx1); i++)
115 for (int rr = 0; rr < numrows; rr++) {
116 cfain[(rr)*TS + numcols + i] = cfain[(rr) * TS + numcols - i - 1];
117 }
118
119 numcols += 16 - indx1;
120 }
121
122 if (numrows < TS) {
123 indx1 = numrows % 16;
124
125 for (int i = 0; i < (16 - indx1); i++)
126 for (int cc = 0; cc < numcols; cc++) {
127 cfain[(numrows + i)*TS + cc] = cfain[(numrows - i - 1) * TS + cc];
128 }
129
130 numrows += 16 - indx1;
131 }
132
133 //The cleaning algorithm starts here
134
135 //gaussian blur of CFA data
136 for (int rr = 8; rr < numrows - 8; rr++) {
137 for (int indx = rr * TS, indxb = 0; indx < rr * TS + numcols; indx++, indxb++) {
138 cfablur[indxb] = gauss[0] * cfain[indx];
139
140 for (int i = 1; i < 5; i++) {
141 cfablur[indxb] += gauss[i] * (cfain[indx - (2 * i) * TS] + cfain[indx + (2 * i) * TS]);
142 }
143 }
144
145 for (int indx = rr * TS + 8, indxb = 8; indx < rr * TS + numcols - 8; indx++, indxb++) {
146 cfadn[indx] = gauss[0] * cfablur[indxb];
147
148 for (int i = 1; i < 5; i++) {
149 cfadn[indx] += gauss[i] * (cfablur[indxb - 2 * i] + cfablur[indxb + 2 * i]);
150 }
151
152 cfadiff[indx] = cfain[indx] - cfadn[indx]; // hipass cfa data
153 }
154 }
155
156 //begin block DCT
157 for (int rr = 8; rr < numrows - 22; rr += 8) // (rr,cc) shift by 8 to overlap blocks
158 for (int cc = 8; cc < numcols - 22; cc += 8) {
159 for (int ey = 0; ey < 2; ey++) // (ex,ey) specify RGGB subarray
160 for (int ex = 0; ex < 2; ex++) {
161 //grab an 8x8 block of a given RGGB channel
162 for (int i = 0; i < 8; i++)
163 for (int j = 0; j < 8; j++) {
164 dctblock[2 * ey + ex][i][j] = cfadiff[(rr + 2 * i + ey) * TS + cc + 2 * j + ex];
165 }
166
167 ddct8x8s(-1, dctblock[2 * ey + ex]); //forward DCT
168 }
169
170 for (int ey = 0; ey < 2; ey++) // (ex,ey) specify RGGB subarray
171 for (int ex = 0; ex < 2; ex++) {
172 linehvar[2 * ey + ex] = linevvar[2 * ey + ex] = 0;
173
174 for (int i = 4; i < 8; i++) {
175 linehvar[2 * ey + ex] += SQR(dctblock[2 * ey + ex][0][i]);
176 linevvar[2 * ey + ex] += SQR(dctblock[2 * ey + ex][i][0]);
177 }
178
179 //Wiener filter for line denoising; roll off low frequencies
180 for (int i = 1; i < 8; i++) {
181 coeffsq = SQR(dctblock[2 * ey + ex][i][0]); //vertical
182 noisefactor[2 * ey + ex][i][0] = coeffsq / (coeffsq + rolloff[i] * noisevar + eps);
183 coeffsq = SQR(dctblock[2 * ey + ex][0][i]); //horizontal
184 noisefactor[2 * ey + ex][i][1] = coeffsq / (coeffsq + rolloff[i] * noisevar + eps);
185 // noisefactor labels are [RGGB subarray][row/col position][0=vert,1=hor]
186 }
187 }
188
189 //horizontal lines
190 if (horizontal && noisevarm4 > (linehvar[0] + linehvar[1])) { //horizontal lines
191 for (int i = 1; i < 8; i++) {
192 dctblock[0][0][i] *= 0.5f * (noisefactor[0][i][1] + noisefactor[1][i][1]); //or should we use MIN???
193 dctblock[1][0][i] *= 0.5f * (noisefactor[0][i][1] + noisefactor[1][i][1]); //or should we use MIN???
194 }
195 }
196
197 if (horizontal && noisevarm4 > (linehvar[2] + linehvar[3])) { //horizontal lines
198 for (int i = 1; i < 8; i++) {
199 dctblock[2][0][i] *= 0.5f * (noisefactor[2][i][1] + noisefactor[3][i][1]); //or should we use MIN???
200 dctblock[3][0][i] *= 0.5f * (noisefactor[2][i][1] + noisefactor[3][i][1]); //or should we use MIN???
201 }
202 }
203
204 //vertical lines
205 if (vertical && noisevarm4 > (linevvar[0] + linevvar[2])) { //vertical lines
206 for (int i = 1; i < 8; i++) {
207 dctblock[0][i][0] *= 0.5f * (noisefactor[0][i][0] + noisefactor[2][i][0]); //or should we use MIN???
208 dctblock[2][i][0] *= 0.5f * (noisefactor[0][i][0] + noisefactor[2][i][0]); //or should we use MIN???
209 }
210 }
211
212 if (vertical && noisevarm4 > (linevvar[1] + linevvar[3])) { //vertical lines
213 for (int i = 1; i < 8; i++) {
214 dctblock[1][i][0] *= 0.5f * (noisefactor[1][i][0] + noisefactor[3][i][0]); //or should we use MIN???
215 dctblock[3][i][0] *= 0.5f * (noisefactor[1][i][0] + noisefactor[3][i][0]); //or should we use MIN???
216 }
217 }
218
219 for (int ey = 0; ey < 2; ey++) // (ex,ey) specify RGGB subarray
220 for (int ex = 0; ex < 2; ex++) {
221 ddct8x8s(1, dctblock[2 * ey + ex]); //inverse DCT
222
223 //multiply by window fn and add to output (cfadn)
224 for (int i = 0; i < 8; i++)
225 for (int j = 0; j < 8; j++) {
226 cfadn[(rr + 2 * i + ey)*TS + cc + 2 * j + ex] += window[i] * window[j] * dctblock[2 * ey + ex][i][j];
227 }
228 }
229 }
230
231 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
232 // copy smoothed results to temporary buffer
233 for (int rr = 16; rr < numrows - 16; rr++) {
234 int row = rr + top;
235
236 for (int col = 16 + left, indx = rr * TS + 16; indx < rr * TS + numcols - 16; indx++, col++) {
237 if (rawData[row][col] < clip_pt && cfadn[indx] < clip_pt) {
238 RawDataTmp[row * width + col] = CLIP(cfadn[indx]);
239 }
240 }
241 }
242
243 if(plistener) {
244 progress += (double)((TS - 32) * (TS - 32)) / (height * width);
245
246 if (progress > 1.0) {
247 progress = 1.0;
248 }
249
250 plistener->setProgress(progress);
251 }
252
253 }
254
255 // clean up
256 free(cfain);
257
258 // copy temporary buffer back to image matrix
259 #ifdef _OPENMP
260 #pragma omp for schedule(dynamic,16)
261 #endif
262
263 for(int i = 0; i < height; i++) {
264 float f = rowblender(i);
265 if (f > 0.f) {
266 float f2 = 1.f - f;
267 for(int j = 0; j < width; j++) {
268 rawData[i][j] = f * RawDataTmp[i * width + j] + f2 * rawData[i][j];
269 }
270 }
271 }
272
273 } // end of parallel processing
274
275 free(RawDataTmp);
276 }
277 #undef TS
278
279
280 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
281 /*
282 Discrete Cosine Transform Code
283
284 Copyright(C) 1997 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
285 You may use, copy, modify this code for any purpose and
286 without fee. You may distribute this ORIGINAL package.
287 */
288
289
290 /*
291 Short Discrete Cosine Transform
292 data length :8x8
293 method :row-column, radix 4 FFT
294 functions
295 ddct8x8s : 8x8 DCT
296 function prototypes
297 void ddct8x8s(int isgn, float **a);
298 */
299
300
301 /*
302 -------- 8x8 DCT (Discrete Cosine Transform) / Inverse of DCT --------
303 [definition]
304 <case1> Normalized 8x8 IDCT
305 C[k1][k2] = (1/4) * sum_j1=0^7 sum_j2=0^7
306 a[j1][j2] * s[j1] * s[j2] *
307 cos(pi*j1*(k1+1/2)/8) *
308 cos(pi*j2*(k2+1/2)/8), 0<=k1<8, 0<=k2<8
309 (s[0] = 1/sqrt(2), s[j] = 1, j > 0)
310 <case2> Normalized 8x8 DCT
311 C[k1][k2] = (1/4) * s[k1] * s[k2] * sum_j1=0^7 sum_j2=0^7
312 a[j1][j2] *
313 cos(pi*(j1+1/2)*k1/8) *
314 cos(pi*(j2+1/2)*k2/8), 0<=k1<8, 0<=k2<8
315 (s[0] = 1/sqrt(2), s[j] = 1, j > 0)
316 [usage]
317 <case1>
318 ddct8x8s(1, a);
319 <case2>
320 ddct8x8s(-1, a);
321 [parameters]
322 a[0...7][0...7] :input/output data (double **)
323 output data
324 a[k1][k2] = C[k1][k2], 0<=k1<8, 0<=k2<8
325 */
326
327
328 /* Cn_kR = sqrt(2.0/n) * cos(pi/2*k/n) */
329 /* Cn_kI = sqrt(2.0/n) * sin(pi/2*k/n) */
330 /* Wn_kR = cos(pi/2*k/n) */
331 /* Wn_kI = sin(pi/2*k/n) */
332 #define C8_1R 0.49039264020161522456f
333 #define C8_1I 0.09754516100806413392f
334 #define C8_2R 0.46193976625564337806f
335 #define C8_2I 0.19134171618254488586f
336 #define C8_3R 0.41573480615127261854f
337 #define C8_3I 0.27778511650980111237f
338 #define C8_4R 0.35355339059327376220f
339 #define W8_4R 0.70710678118654752440f
340
341
ddct8x8s(int isgn,float a[8][8])342 void RawImageSource::ddct8x8s(int isgn, float a[8][8])
343 {
344 int j;
345 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
346 float xr, xi;
347
348 if (isgn < 0) {
349 for (j = 0; j <= 7; j++) {
350 x0r = a[0][j] + a[7][j];
351 x1r = a[0][j] - a[7][j];
352 x0i = a[2][j] + a[5][j];
353 x1i = a[2][j] - a[5][j];
354 x2r = a[4][j] + a[3][j];
355 x3r = a[4][j] - a[3][j];
356 x2i = a[6][j] + a[1][j];
357 x3i = a[6][j] - a[1][j];
358 xr = x0r + x2r;
359 xi = x0i + x2i;
360 a[0][j] = C8_4R * (xr + xi);
361 a[4][j] = C8_4R * (xr - xi);
362 xr = x0r - x2r;
363 xi = x0i - x2i;
364 a[2][j] = C8_2R * xr - C8_2I * xi;
365 a[6][j] = C8_2R * xi + C8_2I * xr;
366 xr = W8_4R * (x1i - x3i);
367 x1i = W8_4R * (x1i + x3i);
368 x3i = x1i - x3r;
369 x1i += x3r;
370 x3r = x1r - xr;
371 x1r += xr;
372 a[1][j] = C8_1R * x1r - C8_1I * x1i;
373 a[7][j] = C8_1R * x1i + C8_1I * x1r;
374 a[3][j] = C8_3R * x3r - C8_3I * x3i;
375 a[5][j] = C8_3R * x3i + C8_3I * x3r;
376 }
377
378 for (j = 0; j <= 7; j++) {
379 x0r = a[j][0] + a[j][7];
380 x1r = a[j][0] - a[j][7];
381 x0i = a[j][2] + a[j][5];
382 x1i = a[j][2] - a[j][5];
383 x2r = a[j][4] + a[j][3];
384 x3r = a[j][4] - a[j][3];
385 x2i = a[j][6] + a[j][1];
386 x3i = a[j][6] - a[j][1];
387 xr = x0r + x2r;
388 xi = x0i + x2i;
389 a[j][0] = C8_4R * (xr + xi);
390 a[j][4] = C8_4R * (xr - xi);
391 xr = x0r - x2r;
392 xi = x0i - x2i;
393 a[j][2] = C8_2R * xr - C8_2I * xi;
394 a[j][6] = C8_2R * xi + C8_2I * xr;
395 xr = W8_4R * (x1i - x3i);
396 x1i = W8_4R * (x1i + x3i);
397 x3i = x1i - x3r;
398 x1i += x3r;
399 x3r = x1r - xr;
400 x1r += xr;
401 a[j][1] = C8_1R * x1r - C8_1I * x1i;
402 a[j][7] = C8_1R * x1i + C8_1I * x1r;
403 a[j][3] = C8_3R * x3r - C8_3I * x3i;
404 a[j][5] = C8_3R * x3i + C8_3I * x3r;
405 }
406 } else {
407 for (j = 0; j <= 7; j++) {
408 x1r = C8_1R * a[1][j] + C8_1I * a[7][j];
409 x1i = C8_1R * a[7][j] - C8_1I * a[1][j];
410 x3r = C8_3R * a[3][j] + C8_3I * a[5][j];
411 x3i = C8_3R * a[5][j] - C8_3I * a[3][j];
412 xr = x1r - x3r;
413 xi = x1i + x3i;
414 x1r += x3r;
415 x3i -= x1i;
416 x1i = W8_4R * (xr + xi);
417 x3r = W8_4R * (xr - xi);
418 xr = C8_2R * a[2][j] + C8_2I * a[6][j];
419 xi = C8_2R * a[6][j] - C8_2I * a[2][j];
420 x0r = C8_4R * (a[0][j] + a[4][j]);
421 x0i = C8_4R * (a[0][j] - a[4][j]);
422 x2r = x0r - xr;
423 x2i = x0i - xi;
424 x0r += xr;
425 x0i += xi;
426 a[0][j] = x0r + x1r;
427 a[7][j] = x0r - x1r;
428 a[2][j] = x0i + x1i;
429 a[5][j] = x0i - x1i;
430 a[4][j] = x2r - x3i;
431 a[3][j] = x2r + x3i;
432 a[6][j] = x2i - x3r;
433 a[1][j] = x2i + x3r;
434 }
435
436 for (j = 0; j <= 7; j++) {
437 x1r = C8_1R * a[j][1] + C8_1I * a[j][7];
438 x1i = C8_1R * a[j][7] - C8_1I * a[j][1];
439 x3r = C8_3R * a[j][3] + C8_3I * a[j][5];
440 x3i = C8_3R * a[j][5] - C8_3I * a[j][3];
441 xr = x1r - x3r;
442 xi = x1i + x3i;
443 x1r += x3r;
444 x3i -= x1i;
445 x1i = W8_4R * (xr + xi);
446 x3r = W8_4R * (xr - xi);
447 xr = C8_2R * a[j][2] + C8_2I * a[j][6];
448 xi = C8_2R * a[j][6] - C8_2I * a[j][2];
449 x0r = C8_4R * (a[j][0] + a[j][4]);
450 x0i = C8_4R * (a[j][0] - a[j][4]);
451 x2r = x0r - xr;
452 x2i = x0i - xi;
453 x0r += xr;
454 x0i += xi;
455 a[j][0] = x0r + x1r;
456 a[j][7] = x0r - x1r;
457 a[j][2] = x0i + x1i;
458 a[j][5] = x0i - x1i;
459 a[j][4] = x2r - x3i;
460 a[j][3] = x2r + x3i;
461 a[j][6] = x2i - x3r;
462 a[j][1] = x2i + x3r;
463 }
464 }
465 }
466
467