r972: Fix aspect ratio of YUV4MPEG streams.
[cinelerra_cv/ct.git] / plugins / denoise / denoise.C
blob970f6f34ced7f917d3749da5a3e075f1aa024f7f
1 #include "bcdisplayinfo.h"
2 #include "clip.h"
3 #include "bchash.h"
4 #include "filexml.h"
5 #include "denoise.h"
6 #include "picon_png.h"
7 #include "units.h"
8 #include "vframe.h"
10 #include <math.h>
11 #include <string.h>
13 #include <libintl.h>
14 #define _(String) gettext(String)
15 #define gettext_noop(String) String
16 #define N_(String) gettext_noop (String)
21 #define WINDOW_BORDER (window_size / 2)
22 #define SGN(x) (x<0 ? -1: 1)
25 REGISTER_PLUGIN(DenoiseEffect)
31 DenoiseEffect::DenoiseEffect(PluginServer *server)
32  : PluginAClient(server)
34         reset();
35         PLUGIN_CONSTRUCTOR_MACRO
38 DenoiseEffect::~DenoiseEffect()
40         PLUGIN_DESTRUCTOR_MACRO
41         delete_dsp();
44 NEW_PICON_MACRO(DenoiseEffect)
46 LOAD_CONFIGURATION_MACRO(DenoiseEffect, DenoiseConfig)
48 SHOW_GUI_MACRO(DenoiseEffect, DenoiseThread)
50 RAISE_WINDOW_MACRO(DenoiseEffect)
52 SET_STRING_MACRO(DenoiseEffect)
54 void DenoiseEffect::delete_dsp()
56         if(ex_coeff_d) delete ex_coeff_d;
57         if(ex_coeff_r) delete ex_coeff_r;
58         if(ex_coeff_rn) delete ex_coeff_rn;
59         if(wave_coeff_d) delete wave_coeff_d;
60         if(wave_coeff_r) delete wave_coeff_r;
61         if(decomp_filter) delete decomp_filter;
62         if(recon_filter) delete recon_filter;
63         if(input_buffer) delete [] input_buffer;
64         if(output_buffer) delete [] output_buffer;
65         if(dsp_in) delete [] dsp_in;
66         if(dsp_out) delete [] dsp_out;
67         if(dsp_iteration) delete [] dsp_iteration;
69         ex_coeff_d = 0;
70         ex_coeff_r = 0;
71         ex_coeff_rn = 0;
72         wave_coeff_d = 0;
73         wave_coeff_r = 0;
74         decomp_filter = 0;
75         recon_filter = 0;
76         input_buffer = 0;
77         output_buffer = 0;
78         dsp_in = 0;
79         dsp_out = 0;
80         dsp_iteration = 0;
84 void DenoiseEffect::reset()
86         first_window = 1;
87         thread = 0;
88         ex_coeff_d = 0;
89         ex_coeff_r = 0;
90         ex_coeff_rn = 0;
91         wave_coeff_d = 0;
92         wave_coeff_r = 0;
93         decomp_filter = 0;
94         recon_filter = 0;
95         input_buffer = 0;
96         output_buffer = 0;
97         input_size = 0;
98         output_size = 0;
99         input_allocation = 0;
100         output_allocation = 0;
101         dsp_iteration = 0;
102         in_scale = 0;
103         out_scale = 0;
104         dsp_in = 0;
105         dsp_out = 0;
106         initialized = 0;
109         alpha = 1.359803732;
110         beta = -0.782106385;
111         window_size = 4096;
112         output_level = 1.0;
113         levels = 1;
114         iterations = 1;
117 char* DenoiseEffect::plugin_title() { return N_("Denoise"); }
118 int DenoiseEffect::is_realtime() { return 1; }
122 void DenoiseEffect::read_data(KeyFrame *keyframe)
124         FileXML input;
125         input.set_shared_string(keyframe->data, strlen(keyframe->data));
127         int result = 0;
128         while(!result)
129         {
130                 result = input.read_tag();
132                 if(!result)
133                 {
134                         if(input.tag.title_is("DENOISE"))
135                         {
136                                 config.level = input.tag.get_property("LEVEL", config.level);
137                         }
138                 }
139         }
142 void DenoiseEffect::save_data(KeyFrame *keyframe)
144         FileXML output;
145         output.set_shared_string(keyframe->data, MESSAGESIZE);
147         output.tag.set_title("DENOISE");
148         output.tag.set_property("LEVEL", config.level);
149         output.append_tag();
150         output.append_newline();
152         output.terminate_string();
155 int DenoiseEffect::load_defaults()
157         char directory[BCTEXTLEN], string[BCTEXTLEN];
158         sprintf(directory, "%sdenoise.rc", BCASTDIR);
159         defaults = new BC_Hash(directory);
160         defaults->load();
161         
162         config.level = defaults->get("LEVEL", config.level);
163         return 0;
166 int DenoiseEffect::save_defaults()
168         char string[BCTEXTLEN];
170         defaults->update("LEVEL", config.level);
171         defaults->save();
173         return 0;
176 void DenoiseEffect::update_gui()
178         if(thread)
179         {
180                 thread->window->lock_window();
181                 thread->window->update();
182                 thread->window->unlock_window();
183         }
188 double DenoiseEffect::dot_product(double *data, double *filter, char filtlen)
190         static int i;
191         static double sum;
193         sum = 0.0;
194         for(i = 0; i < filtlen; i++) sum += *data-- * *filter++;
195         return sum;
198 int DenoiseEffect::convolve_dec_2(double *input_sequence, 
199         int64_t length,
200         double *filter, 
201         int filtlen, 
202         double *output_sequence)
204 // convolve the input sequence with the filter and decimate by two
205         int i, shortlen, offset;
206         int64_t lengthp4 = length + 4;
207         int64_t lengthm4 = length - 4;
208         int64_t lengthp5 = length + 5;
209         int64_t lengthp8 = length + 8;
211         for(i = 0; (i <= lengthp8) && ((i - filtlen) <= lengthp8); i += 2)
212         {
213                 if(i < filtlen)
214                         *output_sequence++ = dot_product(input_sequence + i, filter, i + 1);
215                 else 
216                 if(i > lengthp5)
217                 {
218                         offset = i - lengthm4;
219                         shortlen = filtlen - offset;
220                         *output_sequence++ = dot_product(input_sequence + lengthp4,
221                                                                 filter + offset, shortlen);
222                 }
223                 else
224                         *output_sequence++ = dot_product(input_sequence + i, filter, filtlen);
225         }
226         return 0;
229 int64_t DenoiseEffect::decompose_branches(double *in_data, 
230         int64_t length, 
231         WaveletFilters *decomp_filter, 
232         double *out_low, 
233         double *out_high)
235 // Take input data and filters and form two branches of half the
236 // original length. Length of branches is returned.
237         convolve_dec_2(in_data, length, decomp_filter->h, decomp_filter->length, out_low);
238         convolve_dec_2(in_data, length, decomp_filter->g, decomp_filter->length, out_high);
239         return (length / 2);
242 int DenoiseEffect::wavelet_decomposition(double *in_data, 
243         int64_t in_length, 
244         double **out_data)
246         for(int i = 0; i < levels; i++)
247         {
248                 in_length = decompose_branches(in_data, 
249                         in_length, 
250                         decomp_filter, 
251                         out_data[2 * i], 
252                         out_data[(2 * i) + 1]);
254                 in_data = out_data[2 * i];
255         }
256         return 0;
259 int DenoiseEffect::tree_copy(double **output, 
260         double **input, 
261         int length, 
262         int levels)
264         register int i, j, k, l, m;
266         for(i = 0, k = 1; k < levels; i++, k++)
267         {
268                 length /= 2;
269                 l = 2 * i;
270                 m = l + 1;
272                 for(j = 0; j < length + 5; j++)
273                 {
274                         output[l][j] = 0.0;
275                         output[m][j] = input[m][j];
276                 }
277         }
279         length /= 2;
280         l = 2 * i;
281         m = l + 1;
283         for(j = 0; j < length + 5; j++)
284         {
285                 output[l][j] = input[l][j];
286                 output[m][j] = input[m][j];
287         }
288         return 0;
291 int DenoiseEffect::threshold(int window_size, double gammas, int levels)
293         int i, j;
294         double threshold, cv, cvb, abs_coeff_r;
295         double *coeff_r, *coeff_l;
296         int length;
298         for(i = 0; i < levels; i++) 
299         {
300                 length = (window_size >> (i + 1)) + 5;
301                 threshold = sqrt(2 * log(length) / log(2)) * gammas / sqrt(length);
303                 for(j = 0; j < length; j++) 
304                 {
305                         coeff_r = &(ex_coeff_r->values[(2 * i) + 1][j]);
306                         coeff_l = &(ex_coeff_rn->values[(2 * i) + 1][j]);
308                         cv = SGN(*coeff_r);
309                         abs_coeff_r = fabs(*coeff_r);
310                         cvb = abs_coeff_r - threshold;
311                         cv *= cvb;
313                         if(abs_coeff_r > threshold) 
314                         {
315                                 *coeff_r = cv;
316                                 *coeff_l = 0.0;
317                         }
318                         else 
319                         {
320                                 *coeff_l = *coeff_r;
321                                 *coeff_r = 0.0;
322                         }
323                 }
324         }
328 double DenoiseEffect::dot_product_even(double *data, double *filter, int filtlen)
330         static int i;
331         static double sum;
333         sum = 0.0;
334         for(i = 0; i < filtlen; i += 2) sum += *data-- * filter[i];
335         return sum;
339 double DenoiseEffect::dot_product_odd(double *data, double *filter, int filtlen)
341         static int i;
342         static double sum;
344         sum = 0.0;
345         for(i = 1; i < filtlen; i += 2) sum += *data-- * filter[i];
346         return sum;
349 int DenoiseEffect::convolve_int_2(double *input_sequence, 
350         int64_t length, 
351         double *filter, 
352         int filtlen, 
353         int sum_output, 
354         double *output_sequence)
355 // insert zeros between each element of the input sequence and
356 // convolve with the filter to interpolate the data
358         register int i, j;
359         int endpoint = length + filtlen - 2;
361         if (sum_output)
362         {
363 // summation with previous convolution
364 // every other dot product interpolates the data
365                 for(i = (filtlen / 2) - 1, j = (filtlen / 2); i < endpoint; i++, j++)
366                 {
367                         *output_sequence++ += dot_product_odd(input_sequence + i, filter, filtlen);
368                         *output_sequence++ += dot_product_even(input_sequence + j, filter, filtlen);
369                 }
371                 *output_sequence++ += dot_product_odd(input_sequence + i, filter, filtlen);
372         }
373         else
374         {
375 // first convolution of pair
376 // every other dot product interpolates the data
377                 for(i = (filtlen / 2) - 1, j = (filtlen / 2); i < endpoint; i++, j++)
378                 {
379                         *output_sequence++ = dot_product_odd(input_sequence + i, filter, filtlen);
380                         *output_sequence++ = dot_product_even(input_sequence + j, filter, filtlen);
381                 }
383                 *output_sequence++ = dot_product_odd(input_sequence + i, filter, filtlen);
384         }
385         return 0;
389 int64_t DenoiseEffect::reconstruct_branches(double *in_low, 
390         double *in_high, 
391         int64_t in_length,
392         WaveletFilters *recon_filter, 
393         double *output)
395 // take input data and filters and form two branches of half the
396 // original length. length of branches is returned
397         convolve_int_2(in_low, in_length, recon_filter->h, 
398                                         recon_filter->length, 0, output);
399         convolve_int_2(in_high, in_length, recon_filter->g, 
400                                         recon_filter->length, 1, output);
401         return in_length * 2;
404 int DenoiseEffect::wavelet_reconstruction(double **in_data, 
405         int64_t in_length, 
406         double *out_data)
408         double *output;
409         int i;
411         in_length = in_length >> levels;
412 // destination of all but last branch reconstruction is the next
413 // higher intermediate approximation
414         for(i = levels - 1; i > 0; i--)
415         {
416                 output = in_data[2 * (i - 1)];
417                 in_length = reconstruct_branches(in_data[2 * i], 
418                         in_data[(2 * i) + 1],
419                         in_length, 
420                         recon_filter, 
421                         output);
422         }
424 // destination of the last branch reconstruction is the output data
425         reconstruct_branches(in_data[0], 
426                 in_data[1], 
427                 in_length, 
428                 recon_filter, 
429                 out_data);
431         return 0;
434 void DenoiseEffect::process_window()
436         int i, j;
437         for(j = 0; j < iterations; j++)
438         {
439                 wavelet_decomposition(dsp_in, window_size, ex_coeff_d->values);
441                 tree_copy(ex_coeff_r->values, ex_coeff_d->values, window_size, levels);
442                 tree_copy(ex_coeff_rn->values, ex_coeff_d->values, window_size, levels);
444 // qualify coeffs
445 //printf("DenoiseEffect::process_window %f\n", config.level);
446                 threshold(window_size, config.level * 10.0, levels);
448                 wavelet_reconstruction(ex_coeff_r->values, window_size, dsp_iteration);
449                 wavelet_reconstruction(ex_coeff_rn->values, window_size, dsp_in);
451                 for(i = 0; i < window_size; i++)
452                         dsp_out[i] += dsp_iteration[i];
453         }
459 int DenoiseEffect::process_realtime(int64_t size, double *input_ptr, double *output_ptr)
461         load_configuration();
463         if(!initialized)
464         {
465                 int64_t size_factor = (int)(pow(2, levels));
466                 dsp_in = new double[window_size * size_factor];
467                 dsp_out = new double[window_size * 2];
468                 dsp_iteration = new double[window_size * 2];
471                 ex_coeff_d = new Tree(window_size, levels);
472                 ex_coeff_r = new Tree(window_size, levels);
473                 ex_coeff_rn = new Tree(window_size, levels);
474                 wave_coeff_d = new WaveletCoeffs(alpha, beta);
475                 wave_coeff_r = new WaveletCoeffs(alpha, beta);
476                 decomp_filter = new WaveletFilters(wave_coeff_d, DECOMP);
477                 recon_filter = new WaveletFilters(wave_coeff_r, RECON);
478                 in_scale = 65535 / sqrt(window_size) / iterations;
479                 out_scale = output_level / 65535 * sqrt(window_size);
480                 initialized = 1;
481         }
482         
483 // Append input buffer
484         if(input_size + size > input_allocation)
485         {
486                 double *new_input = new double[input_size + size];
487                 if(input_buffer)
488                 {
489                         memcpy(new_input, input_buffer, sizeof(double) * input_size);
490                         delete [] input_buffer;
491                 }
492                 input_buffer = new_input;
493                 input_allocation = input_size + size;
494         }
495         memcpy(input_buffer + input_size, 
496                 input_ptr, 
497                 size * sizeof(double));
498         input_size += size;
501 // Have enough to do some windows
502         while(input_size >= window_size)
503         {
504 // Load dsp_in
505                 for(int i = 0; i < window_size; i++)
506                 {
507                         dsp_in[i] = input_buffer[i] * in_scale;
508                 }
509                 bzero(dsp_out, sizeof(double) * window_size);
516 // First window produces garbage
517                 if(!first_window)
518                         process_window();
519                 first_window = 0;
526 // Crossfade into the output buffer
527                 int64_t new_allocation = output_size + window_size;
528                 if(new_allocation > output_allocation)
529                 {
530                         double *new_output = new double[new_allocation];
532                         if(output_buffer)
533                         {
534                                 memcpy(new_output, output_buffer, sizeof(double) * output_size);
535 //printf("CrossfadeFFT::process_fifo 1 %p\n", output_buffer);
536                                 delete [] output_buffer;
537 //printf("CrossfadeFFT::process_fifo 2\n");
538                         }
539                         output_buffer = new_output;
540                         output_allocation = new_allocation;
541                 }
543                 if(output_size >= WINDOW_BORDER)
544                 {
545                         for(int i = 0, j = output_size - WINDOW_BORDER; 
546                                 i < WINDOW_BORDER; 
547                                 i++, j++)
548                         {
549                                 double src_level = (double)i / WINDOW_BORDER;
550                                 double dst_level = (double)(WINDOW_BORDER - i) / WINDOW_BORDER;
551                                 output_buffer[j] = output_buffer[j] * dst_level + out_scale * dsp_out[i] * src_level;
552                         }
554                         for(int i = 0; i < window_size - WINDOW_BORDER; i++)
555                                 output_buffer[output_size + i] = dsp_out[WINDOW_BORDER + i] * out_scale;
556                         output_size += window_size - WINDOW_BORDER;
557                 }
558                 else
559                 {
560 // First buffer has no crossfade
561                         memcpy(output_buffer + output_size, 
562                                 dsp_out, 
563                                 sizeof(double) * window_size);
564                         output_size += window_size;
565                 }
568 // Shift input buffer forward
569                 for(int i = window_size - WINDOW_BORDER, j = 0; 
570                         i < input_size; 
571                         i++, j++)
572                         input_buffer[j] = input_buffer[i];
573                 input_size -= window_size - WINDOW_BORDER;
574         }
577 // Have enough to send to output
578         if(output_size - WINDOW_BORDER >= size)
579         {
580                 memcpy(output_ptr, output_buffer, sizeof(double) * size);
581                 for(int i = size, j = 0; i < output_size; i++, j++)
582                         output_buffer[j] = output_buffer[i];
583                 output_size -= size;
584         }
585         else
586         {
587 //printf("DenoiseEffect::process_realtime 1\n");
588                 bzero(output_ptr, sizeof(double) * size);
589         }
591         return 0;
600 Tree::Tree(int input_length, int levels)
602         this->input_length = input_length;
603         this->levels = levels;
604         int i, j;
606 // create decomposition tree
607         values = new double*[2 * levels];
608         j = input_length;
609         for (i = 0; i < levels; i++)
610         {
611                 j /= 2;
612                 if (j == 0)
613                 {
614                         levels = i;
615                         continue;
616                 }
617                 values[2 * i] = new double[j + 5];
618                 values[2 * i + 1] = new double[j + 5];
619         }
622 Tree::~Tree()
624         int i;
626         for (i = 2 * levels - 1; i >= 0; i--)
627                 delete values[i];
629         delete values;
632 WaveletCoeffs::WaveletCoeffs(double alpha, double beta)
634         int i;
635         double tcosa = cos(alpha);
636         double tcosb = cos(beta);
637         double tsina = sin(alpha);
638         double tsinb = sin(beta);
640 // calculate first two wavelet coefficients  a = a(-2) and b = a(-1)
641         values[0] = ((1.0 + tcosa + tsina) * (1.0 - tcosb - tsinb)
642                                         + 2.0 * tsinb * tcosa) / 4.0;
643         values[1] = ((1.0 - tcosa + tsina) * (1.0 + tcosb - tsinb)
644                                         - 2.0 * tsinb * tcosa) / 4.0;
646         tcosa = cos(alpha - beta);
647         tsina = sin(alpha - beta);
649 // calculate last four wavelet coefficients  c = a(0), d = a(1), 
650 // e = a(2), and f = a(3)
651         values[2]  = (1.0 + tcosa + tsina) / 2.0;
652         values[3]  = (1.0 + tcosa - tsina) / 2.0;
653         values[4]  = 1 - values[0] - values[2];
654         values[5]  = 1 - values[1] - values[3];
656 // zero out very small coefficient values caused by truncation error
657         for (i = 0; i < 6; i++)
658         {
659                 if (fabs(values[i]) < 1.0e-15) values[i] = 0.0;
660         }
663 WaveletCoeffs::~WaveletCoeffs()
668 WaveletFilters::WaveletFilters(WaveletCoeffs *wave_coeffs, wavetype transform)
670         int i, j, k;
672 // find the first non-zero wavelet coefficient
673         i = 0;
674         while(wave_coeffs->values[i] == 0.0) i++;
676 // find the last non-zero wavelet coefficient
677         j = 5;
678         while(wave_coeffs->values[j] == 0.0) j--;
680 // Form the decomposition filters h~ and g~ or the reconstruction
681 // filters h and g.  The division by 2 in the construction
682 // of the decomposition filters is for normalization.
683         length = j - i + 1;
684         for(k = 0; k < length; k++)
685         {
686                 if (transform == DECOMP)
687                 {
688                         h[k] = wave_coeffs->values[j--] / 2.0;
689                         g[k] = (double) (((i++ & 0x01) * 2) - 1) * wave_coeffs->values[i] / 2.0;
690                 }
691                 else
692                 {
693                         h[k] = wave_coeffs->values[i++];
694                         g[k] = (double) (((j-- & 0x01) * 2) - 1) * wave_coeffs->values[j];
695                 }
696         }
698 // clear out the additional array locations, if any
699         while (k < 6)
700         {
701                 h[k] = 0.0;
702                 g[k++] = 0.0;
703         }
706 WaveletFilters::~WaveletFilters()
718 DenoiseConfig::DenoiseConfig()
720         level = 1.0;
723 void DenoiseConfig::copy_from(DenoiseConfig &that)
725         level = that.level;
728 int DenoiseConfig::equivalent(DenoiseConfig &that)
730         return EQUIV(level, that.level);
733 void DenoiseConfig::interpolate(DenoiseConfig &prev, 
734         DenoiseConfig &next, 
735         int64_t prev_frame, 
736         int64_t next_frame, 
737         int64_t current_frame)
739         double next_scale = (double)(current_frame - prev_frame) / (next_frame - prev_frame);
740         double prev_scale = (double)(next_frame - current_frame) / (next_frame - prev_frame);
741         this->level = prev.level * prev_scale + next.level * next_scale;
753 PLUGIN_THREAD_OBJECT(DenoiseEffect, DenoiseThread, DenoiseWindow)
764 DenoiseWindow::DenoiseWindow(DenoiseEffect *plugin, int x, int y)
765  : BC_Window(plugin->gui_string, 
766         x, 
767         y, 
768         150, 
769         50, 
770         150, 
771         50,
772         0, 
773         0,
774         1)
776         this->plugin = plugin;
779 void DenoiseWindow::create_objects()
781         int x = 10, y = 10;
782         
783         add_subwindow(new BC_Title(x, y, _("Level:")));
784         x += 70;
785         add_subwindow(scale = new DenoiseLevel(plugin, x, y));
786         show_window();
787         flush();
790 int DenoiseWindow::close_event()
792 // Set result to 1 to indicate a client side close
793         set_done(1);
794         return 1;
797 void DenoiseWindow::update()
799         scale->update(plugin->config.level);
813 DenoiseLevel::DenoiseLevel(DenoiseEffect *plugin, int x, int y)
814  : BC_FPot(x, y, (float)plugin->config.level, 0, 1.0)
816         this->plugin = plugin;
817         set_precision(0.01);
820 int DenoiseLevel::handle_event()
822         plugin->config.level = get_value();
823         plugin->send_configure_change();
824         return 1;