1 /*
2 * * Copyright (C) 2006-2011 Anders Brander <anders@brander.dk>,
3 * * Anders Kvist <akv@lnxbx.dk> and Klaus Post <klauspost@gmail.com>
4 *
5 * This program is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU General Public License
7 * as published by the Free Software Foundation; either version 2
8 * of the License, or (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 *
19 * Based on FFT3DFilter plugin for Avisynth 2.5 - 3D Frequency Domain filter
20 * Copyright(C)2004-2005 A.G.Balakhnin aka Fizick, email: bag@hotmail.ru, web: http://bag.hotmail.ru
21 */
22
23 #include "complexfilter.h"
24 #include <math.h>
25 #include "fftwindow.h"
26
27 /*
28 * These classes define the processing that must be done
29 * to each block.
30 * Note that the process() function must be re-entrant,
31 * as all threads will use the same instance for each plane.
32 *
33 */
34 #ifndef MAX
35 #define MAX(x,y) ((x) > (y) ? (x) : (y))
36 #endif
37
38 namespace RawStudio {
39 namespace FFTFilter {
40
41
42 /**** BASE CLASS *****/
43
ComplexFilter(int block_width,int block_height)44 ComplexFilter::ComplexFilter( int block_width, int block_height ) :
45 bw(block_width), bh(block_height), norm(1.0f/(block_width*block_height)),
46 sharpen(0), sigmaSquaredSharpenMin(0), sigmaSquaredSharpenMax(0), sharpenWindow(0)
47 {}
48
~ComplexFilter(void)49 ComplexFilter::~ComplexFilter(void)
50 {
51 if (sharpenWindow)
52 delete sharpenWindow;
53 sharpenWindow = 0;
54 }
55
setSharpen(float _sharpen,float sigmaSharpenMin,float sigmaSharpenMax,float scutoff)56 void ComplexFilter::setSharpen( float _sharpen, float sigmaSharpenMin, float sigmaSharpenMax, float scutoff )
57 {
58 if (ABS(_sharpen) <0.001f)
59 return;
60 sharpen = _sharpen;
61 sigmaSquaredSharpenMin = sigmaSharpenMin*sigmaSharpenMin/norm;
62 sigmaSquaredSharpenMax = sigmaSharpenMax*sigmaSharpenMax/norm;
63 // window for sharpen
64 float svr = 1.0f; // Horizontal to vertical ratio
65 if (!sharpenWindow) {
66 sharpenWindow = new FloatImagePlane(bw, bh);
67 sharpenWindow->allocateImage();
68 }
69
70 for (int j=0; j<bh; j++)
71 {
72 int dj = j;
73 if (j>=bh/2)
74 dj = bh-j;
75 float d2v = float(dj*dj)*(svr*svr)/((bh/2)*(bh/2));
76 float *wsharpen = sharpenWindow->getLine(j);
77 for (int i=0; i<bw; i++)
78 {
79 float d2 = d2v + float(i*i)/((bw/2)*(bw/2)); // distance_2 - v1.7
80 wsharpen[i] = sharpen * (1 - exp(-d2/(2*scutoff*scutoff)));
81 }
82 }
83 /* In sharpen function, remember: Sharpen factor is already applied to wsharpen*/
84 }
85
process(ComplexBlock * block)86 void ComplexFilter::process( ComplexBlock* block )
87 {
88 if (ABS(sharpen) >0.001f)
89 processSharpen(block);
90 else
91 processNoSharpen(block);
92 }
93
skipBlock()94 gboolean ComplexFilter::skipBlock() {
95 if (ABS(sharpen) >0.001f)
96 return false;
97 return true;
98 }
99
100 /** DeGridComplexFilter **/
DeGridComplexFilter(int block_width,int block_height,float _degrid,FFTWindow * _window,fftwf_plan plan_forward)101 DeGridComplexFilter::DeGridComplexFilter(int block_width, int block_height, float _degrid, FFTWindow *_window, fftwf_plan plan_forward) :
102 ComplexFilter(block_width, block_height),
103 degrid(_degrid),
104 window(_window)
105 {
106 grid = new ComplexBlock(bw, bh);
107 FloatImagePlane realGrid(bw, bh);
108 realGrid.allocateImage();
109 float* f = realGrid.data;
110 int count = bh*realGrid.pitch;
111
112 for (int i = 0 ; i < count; i++){
113 f[i] = 65535.0f;
114 }
115 window->applyAnalysisWindow(&realGrid,&realGrid);
116 fftwf_execute_dft_r2c(plan_forward, f, grid->complex);
117 }
118
~DeGridComplexFilter(void)119 DeGridComplexFilter::~DeGridComplexFilter( void )
120 {
121 delete grid;
122 }
123
processSharpenOnly(ComplexBlock * block)124 void DeGridComplexFilter::processSharpenOnly(ComplexBlock* block) {
125
126 #if defined (__i386__) || defined (__x86_64__)
127 guint cpu = rs_detect_cpu_features();
128 if (cpu & RS_CPU_FLAG_SSE3)
129 return processSharpenOnlySSE3(block);
130 else if (cpu & RS_CPU_FLAG_SSE)
131 return processSharpenOnlySSE(block);
132 #endif
133
134 int x,y;
135 fftwf_complex* outcur = block->complex;
136 fftwf_complex* gridsample = grid->complex;
137
138 float gridfraction = degrid*outcur[0][0]/gridsample[0][0];
139 for (y=0; y<bh; y++) {
140 float *wsharpen = sharpenWindow->getLine(y);
141 for (x=0; x<bw; x++) {
142 float gridcorrection0 = gridfraction*gridsample[x][0];
143 float re = outcur[x][0] - gridcorrection0;
144 float gridcorrection1 = gridfraction*gridsample[x][1];
145 float im = outcur[x][1] - gridcorrection1;
146 float psd = (re*re + im*im) + 1e-15f;// power spectrum density
147 //improved sharpen mode to prevent grid artifactes and to limit sharpening both fo low and high amplitudes
148 float sfact = (1 + wsharpen[x]*sqrt( psd*sigmaSquaredSharpenMax/((psd + sigmaSquaredSharpenMin)*(psd + sigmaSquaredSharpenMax)) )) ;
149 re *= sfact; // apply filter on real part
150 im *= sfact; // apply filter on imaginary part
151 outcur[x][0] = re + gridcorrection0;
152 outcur[x][1] = im + gridcorrection1;
153 }
154 gridsample += bw;
155 outcur += bw;
156 wsharpen += bw;
157 }
158 }
159
160
161 /**** Basic Wiener Filter *****/
162
163
ComplexWienerFilter(int block_width,int block_height,float _beta,float _sigma)164 ComplexWienerFilter::ComplexWienerFilter( int block_width, int block_height,float _beta, float _sigma ) :
165 ComplexFilter(block_width, block_height)
166 {
167 lowlimit = (_beta-1)/_beta;
168 sigmaSquaredNoiseNormed = _sigma*_sigma/norm;
169 }
170
171
~ComplexWienerFilter(void)172 ComplexWienerFilter::~ComplexWienerFilter( void ){}
173
skipBlock()174 gboolean ComplexWienerFilter::skipBlock() {
175 if (ABS(sharpen) >0.001f)
176 return false;
177 if (sigmaSquaredNoiseNormed > 1e-15f)
178 return false;
179 return true;
180 }
181
processNoSharpen(ComplexBlock * block)182 void ComplexWienerFilter::processNoSharpen( ComplexBlock* block )
183 {
184 int x,y;
185 float psd;
186 float WienerFactor;
187 fftwf_complex* outcur = block->complex;
188 g_assert(bw == block->w);
189 g_assert(bh == block->h);
190
191 for (y=0; y<bh; y++) {
192 for (x=0; x<bw; x++) {
193 psd = (outcur[x][0]*outcur[x][0] + outcur[x][1]*outcur[x][1]) + 1e-15f;// power spectrum density
194 WienerFactor = MAX((psd - sigmaSquaredNoiseNormed)/psd, lowlimit); // limited Wiener filter
195 outcur[x][0] *= WienerFactor; // apply filter on real part
196 outcur[x][1] *= WienerFactor; // apply filter on imaginary part
197 }
198 outcur += bw;
199 }
200
201 }
202
processSharpen(ComplexBlock * block)203 void ComplexWienerFilter::processSharpen( ComplexBlock* block )
204 {
205 int x,y;
206 float psd;
207 float WienerFactor;
208 fftwf_complex* outcur = block->complex;
209 g_assert(bw == block->w);
210 g_assert(bh == block->h);
211 for (y=0; y<bh; y++) {
212 float *wsharpen = sharpenWindow->getLine(y);
213 for (x=0; x<bw; x++) {
214 psd = (outcur[x][0]*outcur[x][0] + outcur[x][1]*outcur[x][1]) + 1e-15f;// power spectrum density
215 WienerFactor = MAX((psd - sigmaSquaredNoiseNormed)/psd, lowlimit); // limited Wiener filter
216 WienerFactor *= 1 + wsharpen[x]*sqrt( psd*sigmaSquaredSharpenMax/((psd + sigmaSquaredSharpenMin)*(psd + sigmaSquaredSharpenMax)) );
217 outcur[x][0] *= WienerFactor; // apply filter on real part
218 outcur[x][1] *= WienerFactor; // apply filter on imaginary part
219 }
220 outcur += bw;
221 wsharpen += bw;
222 }
223
224 }
225
226
227 /**** Apply Pattern Filter *****/
228
ComplexPatternFilter(int block_width,int block_height,float _beta,FloatImagePlane * _pattern,float pattern_strength)229 ComplexPatternFilter::ComplexPatternFilter( int block_width, int block_height, float _beta,
230 FloatImagePlane* _pattern, float pattern_strength ) :
231 ComplexFilter(block_width, block_height),
232 pfactor(pattern_strength)
233 {
234 lowlimit = (_beta-1)/_beta;
235 pattern = _pattern;
236 }
237
~ComplexPatternFilter(void)238 ComplexPatternFilter::~ComplexPatternFilter( void )
239 {
240 if (pattern)
241 delete pattern;
242 }
243
244
processNoSharpen(ComplexBlock * block)245 void ComplexPatternFilter::processNoSharpen( ComplexBlock* block )
246 {
247 g_assert(bw == block->w);
248 g_assert(bh == block->h);
249 int x,y;
250 float psd;
251 fftwf_complex* outcur = block->complex;
252 float* pattern2d = pattern->data;
253 float patternfactor;
254
255 for (y=0; y<bh; y++) {
256 for (x=0; x<bw; x++) {
257 psd = (outcur[x][0]*outcur[x][0] + outcur[x][1]*outcur[x][1]) + 1e-15f;
258 patternfactor = MAX((psd - pfactor*pattern2d[x])/psd, lowlimit);
259 outcur[x][0] *= patternfactor;
260 outcur[x][1] *= patternfactor;
261 }
262 outcur += bw;
263 pattern2d += pattern->pitch;
264 }
265 }
266
processSharpen(ComplexBlock * block)267 void ComplexPatternFilter::processSharpen( ComplexBlock* block )
268 {
269 g_assert(!"Not implemented");
270 }
271
skipBlock()272 gboolean ComplexPatternFilter::skipBlock() {
273 if (ABS(sharpen) >0.001f)
274 return false;
275 if (pfactor > 1e-15f)
276 return false;
277 return true;
278 }
279
ComplexWienerFilterDeGrid(int block_width,int block_height,float _beta,float _sigma,float _degrid,fftwf_plan plan_forward,FFTWindow * _window)280 ComplexWienerFilterDeGrid::ComplexWienerFilterDeGrid( int block_width, int block_height,
281 float _beta, float _sigma, float _degrid,
282 fftwf_plan plan_forward, FFTWindow *_window)
283 : DeGridComplexFilter(block_width, block_height, _degrid, _window, plan_forward)
284 {
285 lowlimit = (_beta-1)/_beta;
286 sigmaSquaredNoiseNormed = _sigma*_sigma/norm;
287 }
288
~ComplexWienerFilterDeGrid(void)289 ComplexWienerFilterDeGrid::~ComplexWienerFilterDeGrid( void )
290 {
291 }
292
skipBlock()293 gboolean ComplexWienerFilterDeGrid::skipBlock() {
294 if (ABS(sharpen) >0.001f)
295 return false;
296 if (sigmaSquaredNoiseNormed > 1e-15f)
297 return false;
298 return true;
299 }
300
processNoSharpen(ComplexBlock * block)301 void ComplexWienerFilterDeGrid::processNoSharpen( ComplexBlock* block )
302 {
303 if (sigmaSquaredNoiseNormed <= 1e-15f)
304 return;
305
306 #if defined (__i386__) || defined (__x86_64__)
307 guint cpu = rs_detect_cpu_features();
308 if (cpu & RS_CPU_FLAG_SSE3)
309 return processNoSharpen_SSE3(block);
310 else if (cpu & RS_CPU_FLAG_SSE)
311 return processNoSharpen_SSE(block);
312 #endif
313
314 int x,y;
315 float psd;
316 float WienerFactor;
317 fftwf_complex* outcur = block->complex;
318 fftwf_complex* gridsample = grid->complex;
319
320 float gridfraction = degrid*outcur[0][0]/gridsample[0][0];
321 for (y=0; y<bh; y++) {
322 for (x=0; x<bw; x++) {
323 float gridcorrection0 = gridfraction*gridsample[x][0];
324 float corrected0 = outcur[x][0] - gridcorrection0;
325 float gridcorrection1 = gridfraction*gridsample[x][1];
326 float corrected1 = outcur[x][1] - gridcorrection1;
327 psd = (corrected0*corrected0 + corrected1*corrected1 ) + 1e-15f;// power spectrum density
328 WienerFactor = MAX((psd - sigmaSquaredNoiseNormed)/psd, lowlimit); // limited Wiener filter
329 corrected0 *= WienerFactor; // apply filter on real part
330 corrected1 *= WienerFactor; // apply filter on imaginary part
331 outcur[x][0] = corrected0 + gridcorrection0;
332 outcur[x][1] = corrected1 + gridcorrection1;
333 }
334 outcur += bw;
335 gridsample += bw;
336 }
337
338 }
339
340
processSharpen(ComplexBlock * block)341 void ComplexWienerFilterDeGrid::processSharpen( ComplexBlock* block )
342 {
343 if (sigmaSquaredNoiseNormed <= 1e-15f)
344 return processSharpenOnly(block);
345
346 #if defined (__i386__) || defined (__x86_64__)
347 guint cpu = rs_detect_cpu_features();
348 if (cpu & RS_CPU_FLAG_SSE3)
349 return processSharpen_SSE3(block);
350 else if (cpu & RS_CPU_FLAG_SSE)
351 return processSharpen_SSE(block);
352 #endif
353
354 int x,y;
355 float psd;
356 float WienerFactor;
357 fftwf_complex* outcur = block->complex;
358 fftwf_complex* gridsample = grid->complex;
359
360 float gridfraction = degrid*outcur[0][0]/gridsample[0][0];
361 for (y=0; y<bh; y++) {
362 float *wsharpen = sharpenWindow->getLine(y);
363 for (x=0; x<bw; x++) {
364 float gridcorrection0 = gridfraction*gridsample[x][0];
365 float corrected0 = outcur[x][0] - gridcorrection0;
366 float gridcorrection1 = gridfraction*gridsample[x][1];
367 float corrected1 = outcur[x][1] - gridcorrection1;
368 psd = (corrected0*corrected0 + corrected1*corrected1 ) + 1e-15f;// power spectrum density
369 WienerFactor = MAX((psd - sigmaSquaredNoiseNormed)/psd, lowlimit); // limited Wiener filter
370 WienerFactor *= 1 + wsharpen[x]*sqrt( psd*sigmaSquaredSharpenMax/((psd + sigmaSquaredSharpenMin)*(psd + sigmaSquaredSharpenMax)) );
371 corrected0 *= WienerFactor; // apply filter on real part
372 corrected1 *= WienerFactor; // apply filter on imaginary part
373 outcur[x][0] = corrected0 + gridcorrection0;
374 outcur[x][1] = corrected1 + gridcorrection1;
375
376 }
377 outcur += bw;
378 gridsample += bw;
379 }
380 }
381
ComplexFilterPatternDeGrid(int block_width,int block_height,float _beta,float _sigma,float _degrid,fftwf_plan plan_forward,FFTWindow * _window,FloatImagePlane * _pattern)382 ComplexFilterPatternDeGrid::ComplexFilterPatternDeGrid( int block_width, int block_height,
383 float _beta, float _sigma, float _degrid,
384 fftwf_plan plan_forward, FFTWindow *_window,
385 FloatImagePlane *_pattern)
386 :
387 DeGridComplexFilter(block_width, block_height, _degrid, _window, plan_forward),
388 pattern(_pattern)
389 {
390 lowlimit = (_beta-1)/_beta;
391 sigmaSquaredNoiseNormed = _sigma*_sigma/norm;
392 }
393
~ComplexFilterPatternDeGrid(void)394 ComplexFilterPatternDeGrid::~ComplexFilterPatternDeGrid( void )
395 {
396 }
397
skipBlock()398 gboolean ComplexFilterPatternDeGrid::skipBlock() {
399 return false;
400 }
401
402
processNoSharpen(ComplexBlock * block)403 void ComplexFilterPatternDeGrid::processNoSharpen( ComplexBlock* block )
404 {
405 int x,y;
406 float psd;
407 float WienerFactor;
408 fftwf_complex* outcur = block->complex;
409 fftwf_complex* gridsample = grid->complex;
410
411 float gridfraction = degrid*outcur[0][0]/gridsample[0][0];
412 for (y=0; y<bh; y++) {
413 float *pattern2d = pattern->getLine(y);
414 for (x=0; x<bw; x++) {
415 float gridcorrection0 = gridfraction*gridsample[x][0];
416 float corrected0 = outcur[x][0] - gridcorrection0;
417 float gridcorrection1 = gridfraction*gridsample[x][1];
418 float corrected1 = outcur[x][1] - gridcorrection1;
419 psd = (corrected0*corrected0 + corrected1*corrected1 ) + 1e-15f;// power spectrum density
420 WienerFactor = MAX((psd - pattern2d[x])/psd, lowlimit); // limited Wiener filter
421 corrected0 *= WienerFactor; // apply filter on real part
422 corrected1 *= WienerFactor; // apply filter on imaginary part
423 outcur[x][0] = corrected0 + gridcorrection0;
424 outcur[x][1] = corrected1 + gridcorrection1;
425 }
426 outcur += bw;
427 gridsample += bw;
428 }
429 }
430
431 /* Naiive cut together from ApplyPattern2D_degrid_C and Sharpen_degrid_C,
432 it may be possible to factor out some grid correction */
433
processSharpen(ComplexBlock * block)434 void ComplexFilterPatternDeGrid::processSharpen( ComplexBlock* block )
435 {
436 if (sigmaSquaredNoiseNormed <= 1e-15f)
437 return processSharpenOnly(block);
438
439 int x,y;
440 float psd;
441 float WienerFactor;
442 fftwf_complex* outcur = block->complex;
443 fftwf_complex* gridsample = grid->complex;
444
445 float gridfraction = degrid*outcur[0][0]/gridsample[0][0];
446 for (y=0; y<bh; y++) {
447 float *pattern2d = pattern->getLine(y);
448 float *wsharpen = sharpenWindow->getLine(y);
449 for (x=0; x<bw; x++) {
450 float gridcorrection0 = gridfraction*gridsample[x][0];
451 float corrected0 = outcur[x][0] - gridcorrection0;
452 float gridcorrection1 = gridfraction*gridsample[x][1];
453 float corrected1 = outcur[x][1] - gridcorrection1;
454 psd = (corrected0*corrected0 + corrected1*corrected1 ) + 1e-15f;// power spectrum density
455 WienerFactor = MAX((psd - pattern2d[x])/psd, lowlimit); // limited Wiener filter
456
457 corrected0 *= WienerFactor; // apply filter on real part
458 corrected1 *= WienerFactor; // apply filter on imaginary part
459 corrected0 += gridcorrection0; // apply filter on real part
460 corrected1 += gridcorrection1; // apply filter on imaginary part
461
462 gridcorrection0 = gridfraction*corrected0;
463 float re = corrected0 - gridcorrection0;
464 gridcorrection1 = gridfraction*corrected0;
465 float im = corrected1 - gridcorrection1;
466 psd = (re*re + im*im) + 1e-15f;// power spectrum density
467
468 float sfact = (1 + wsharpen[x]*sqrt( psd*sigmaSquaredSharpenMax/((psd + sigmaSquaredSharpenMin)*(psd + sigmaSquaredSharpenMax)) )) ;
469
470 corrected0 *= sfact; // apply filter on real part
471 corrected1 *= sfact; // apply filter on imaginary part
472
473 outcur[x][0] = corrected0 + gridcorrection0;
474 outcur[x][1] = corrected1 + gridcorrection1;
475 }
476 outcur += bw;
477 gridsample += bw;
478 }
479 }
480
481 }}// namespace RawStudio::FFTFilter
482