Merge branch 'ct' of git.pipapo.org:cinelerra-ct into ct
[cinelerra_cv/ct.git] / plugins / denoise / denoise.C
blob6bdc8c57e2e9b013435df602db04af2603d65a6e
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.tag.set_title("/DENOISE");
151         output.append_tag();
152         output.append_newline();
154         output.terminate_string();
157 int DenoiseEffect::load_defaults()
159         char directory[BCTEXTLEN], string[BCTEXTLEN];
160         sprintf(directory, "%sdenoise.rc", BCASTDIR);
161         defaults = new BC_Hash(directory);
162         defaults->load();
163         
164         config.level = defaults->get("LEVEL", config.level);
165         return 0;
168 int DenoiseEffect::save_defaults()
170         char string[BCTEXTLEN];
172         defaults->update("LEVEL", config.level);
173         defaults->save();
175         return 0;
178 void DenoiseEffect::update_gui()
180         if(thread)
181         {
182                 thread->window->lock_window();
183                 thread->window->update();
184                 thread->window->unlock_window();
185         }
190 double DenoiseEffect::dot_product(double *data, double *filter, char filtlen)
192         static int i;
193         static double sum;
195         sum = 0.0;
196         for(i = 0; i < filtlen; i++) sum += *data-- * *filter++;
197         return sum;
200 int DenoiseEffect::convolve_dec_2(double *input_sequence, 
201         int64_t length,
202         double *filter, 
203         int filtlen, 
204         double *output_sequence)
206 // convolve the input sequence with the filter and decimate by two
207         int i, shortlen, offset;
208         int64_t lengthp4 = length + 4;
209         int64_t lengthm4 = length - 4;
210         int64_t lengthp5 = length + 5;
211         int64_t lengthp8 = length + 8;
213         for(i = 0; (i <= lengthp8) && ((i - filtlen) <= lengthp8); i += 2)
214         {
215                 if(i < filtlen)
216                         *output_sequence++ = dot_product(input_sequence + i, filter, i + 1);
217                 else 
218                 if(i > lengthp5)
219                 {
220                         offset = i - lengthm4;
221                         shortlen = filtlen - offset;
222                         *output_sequence++ = dot_product(input_sequence + lengthp4,
223                                                                 filter + offset, shortlen);
224                 }
225                 else
226                         *output_sequence++ = dot_product(input_sequence + i, filter, filtlen);
227         }
228         return 0;
231 int64_t DenoiseEffect::decompose_branches(double *in_data, 
232         int64_t length, 
233         WaveletFilters *decomp_filter, 
234         double *out_low, 
235         double *out_high)
237 // Take input data and filters and form two branches of half the
238 // original length. Length of branches is returned.
239         convolve_dec_2(in_data, length, decomp_filter->h, decomp_filter->length, out_low);
240         convolve_dec_2(in_data, length, decomp_filter->g, decomp_filter->length, out_high);
241         return (length / 2);
244 int DenoiseEffect::wavelet_decomposition(double *in_data, 
245         int64_t in_length, 
246         double **out_data)
248         for(int i = 0; i < levels; i++)
249         {
250                 in_length = decompose_branches(in_data, 
251                         in_length, 
252                         decomp_filter, 
253                         out_data[2 * i], 
254                         out_data[(2 * i) + 1]);
256                 in_data = out_data[2 * i];
257         }
258         return 0;
261 int DenoiseEffect::tree_copy(double **output, 
262         double **input, 
263         int length, 
264         int levels)
266         register int i, j, k, l, m;
268         for(i = 0, k = 1; k < levels; i++, k++)
269         {
270                 length /= 2;
271                 l = 2 * i;
272                 m = l + 1;
274                 for(j = 0; j < length + 5; j++)
275                 {
276                         output[l][j] = 0.0;
277                         output[m][j] = input[m][j];
278                 }
279         }
281         length /= 2;
282         l = 2 * i;
283         m = l + 1;
285         for(j = 0; j < length + 5; j++)
286         {
287                 output[l][j] = input[l][j];
288                 output[m][j] = input[m][j];
289         }
290         return 0;
293 int DenoiseEffect::threshold(int window_size, double gammas, int levels)
295         int i, j;
296         double threshold, cv, cvb, abs_coeff_r;
297         double *coeff_r, *coeff_l;
298         int length;
300         for(i = 0; i < levels; i++) 
301         {
302                 length = (window_size >> (i + 1)) + 5;
303                 threshold = sqrt(2 * log(length) / log(2)) * gammas / sqrt(length);
305                 for(j = 0; j < length; j++) 
306                 {
307                         coeff_r = &(ex_coeff_r->values[(2 * i) + 1][j]);
308                         coeff_l = &(ex_coeff_rn->values[(2 * i) + 1][j]);
310                         cv = SGN(*coeff_r);
311                         abs_coeff_r = fabs(*coeff_r);
312                         cvb = abs_coeff_r - threshold;
313                         cv *= cvb;
315                         if(abs_coeff_r > threshold) 
316                         {
317                                 *coeff_r = cv;
318                                 *coeff_l = 0.0;
319                         }
320                         else 
321                         {
322                                 *coeff_l = *coeff_r;
323                                 *coeff_r = 0.0;
324                         }
325                 }
326         }
330 double DenoiseEffect::dot_product_even(double *data, double *filter, int filtlen)
332         static int i;
333         static double sum;
335         sum = 0.0;
336         for(i = 0; i < filtlen; i += 2) sum += *data-- * filter[i];
337         return sum;
341 double DenoiseEffect::dot_product_odd(double *data, double *filter, int filtlen)
343         static int i;
344         static double sum;
346         sum = 0.0;
347         for(i = 1; i < filtlen; i += 2) sum += *data-- * filter[i];
348         return sum;
351 int DenoiseEffect::convolve_int_2(double *input_sequence, 
352         int64_t length, 
353         double *filter, 
354         int filtlen, 
355         int sum_output, 
356         double *output_sequence)
357 // insert zeros between each element of the input sequence and
358 // convolve with the filter to interpolate the data
360         register int i, j;
361         int endpoint = length + filtlen - 2;
363         if (sum_output)
364         {
365 // summation with previous convolution
366 // every other dot product interpolates the data
367                 for(i = (filtlen / 2) - 1, j = (filtlen / 2); i < endpoint; i++, j++)
368                 {
369                         *output_sequence++ += dot_product_odd(input_sequence + i, filter, filtlen);
370                         *output_sequence++ += dot_product_even(input_sequence + j, filter, filtlen);
371                 }
373                 *output_sequence++ += dot_product_odd(input_sequence + i, filter, filtlen);
374         }
375         else
376         {
377 // first convolution of pair
378 // every other dot product interpolates the data
379                 for(i = (filtlen / 2) - 1, j = (filtlen / 2); i < endpoint; i++, j++)
380                 {
381                         *output_sequence++ = dot_product_odd(input_sequence + i, filter, filtlen);
382                         *output_sequence++ = dot_product_even(input_sequence + j, filter, filtlen);
383                 }
385                 *output_sequence++ = dot_product_odd(input_sequence + i, filter, filtlen);
386         }
387         return 0;
391 int64_t DenoiseEffect::reconstruct_branches(double *in_low, 
392         double *in_high, 
393         int64_t in_length,
394         WaveletFilters *recon_filter, 
395         double *output)
397 // take input data and filters and form two branches of half the
398 // original length. length of branches is returned
399         convolve_int_2(in_low, in_length, recon_filter->h, 
400                                         recon_filter->length, 0, output);
401         convolve_int_2(in_high, in_length, recon_filter->g, 
402                                         recon_filter->length, 1, output);
403         return in_length * 2;
406 int DenoiseEffect::wavelet_reconstruction(double **in_data, 
407         int64_t in_length, 
408         double *out_data)
410         double *output;
411         int i;
413         in_length = in_length >> levels;
414 // destination of all but last branch reconstruction is the next
415 // higher intermediate approximation
416         for(i = levels - 1; i > 0; i--)
417         {
418                 output = in_data[2 * (i - 1)];
419                 in_length = reconstruct_branches(in_data[2 * i], 
420                         in_data[(2 * i) + 1],
421                         in_length, 
422                         recon_filter, 
423                         output);
424         }
426 // destination of the last branch reconstruction is the output data
427         reconstruct_branches(in_data[0], 
428                 in_data[1], 
429                 in_length, 
430                 recon_filter, 
431                 out_data);
433         return 0;
436 void DenoiseEffect::process_window()
438         int i, j;
439         for(j = 0; j < iterations; j++)
440         {
441                 wavelet_decomposition(dsp_in, window_size, ex_coeff_d->values);
443                 tree_copy(ex_coeff_r->values, ex_coeff_d->values, window_size, levels);
444                 tree_copy(ex_coeff_rn->values, ex_coeff_d->values, window_size, levels);
446 // qualify coeffs
447 //printf("DenoiseEffect::process_window %f\n", config.level);
448                 threshold(window_size, config.level * 10.0, levels);
450                 wavelet_reconstruction(ex_coeff_r->values, window_size, dsp_iteration);
451                 wavelet_reconstruction(ex_coeff_rn->values, window_size, dsp_in);
453                 for(i = 0; i < window_size; i++)
454                         dsp_out[i] += dsp_iteration[i];
455         }
461 int DenoiseEffect::process_realtime(int64_t size, double *input_ptr, double *output_ptr)
463         load_configuration();
465         if(!initialized)
466         {
467                 int64_t size_factor = (int)(pow(2, levels));
468                 dsp_in = new double[window_size * size_factor];
469                 dsp_out = new double[window_size * 2];
470                 dsp_iteration = new double[window_size * 2];
473                 ex_coeff_d = new Tree(window_size, levels);
474                 ex_coeff_r = new Tree(window_size, levels);
475                 ex_coeff_rn = new Tree(window_size, levels);
476                 wave_coeff_d = new WaveletCoeffs(alpha, beta);
477                 wave_coeff_r = new WaveletCoeffs(alpha, beta);
478                 decomp_filter = new WaveletFilters(wave_coeff_d, DECOMP);
479                 recon_filter = new WaveletFilters(wave_coeff_r, RECON);
480                 in_scale = 65535 / sqrt(window_size) / iterations;
481                 out_scale = output_level / 65535 * sqrt(window_size);
482                 initialized = 1;
483         }
484         
485 // Append input buffer
486         if(input_size + size > input_allocation)
487         {
488                 double *new_input = new double[input_size + size];
489                 if(input_buffer)
490                 {
491                         memcpy(new_input, input_buffer, sizeof(double) * input_size);
492                         delete [] input_buffer;
493                 }
494                 input_buffer = new_input;
495                 input_allocation = input_size + size;
496         }
497         memcpy(input_buffer + input_size, 
498                 input_ptr, 
499                 size * sizeof(double));
500         input_size += size;
503 // Have enough to do some windows
504         while(input_size >= window_size)
505         {
506 // Load dsp_in
507                 for(int i = 0; i < window_size; i++)
508                 {
509                         dsp_in[i] = input_buffer[i] * in_scale;
510                 }
511                 bzero(dsp_out, sizeof(double) * window_size);
518 // First window produces garbage
519                 if(!first_window)
520                         process_window();
521                 first_window = 0;
528 // Crossfade into the output buffer
529                 int64_t new_allocation = output_size + window_size;
530                 if(new_allocation > output_allocation)
531                 {
532                         double *new_output = new double[new_allocation];
534                         if(output_buffer)
535                         {
536                                 memcpy(new_output, output_buffer, sizeof(double) * output_size);
537 //printf("CrossfadeFFT::process_fifo 1 %p\n", output_buffer);
538                                 delete [] output_buffer;
539 //printf("CrossfadeFFT::process_fifo 2\n");
540                         }
541                         output_buffer = new_output;
542                         output_allocation = new_allocation;
543                 }
545                 if(output_size >= WINDOW_BORDER)
546                 {
547                         for(int i = 0, j = output_size - WINDOW_BORDER; 
548                                 i < WINDOW_BORDER; 
549                                 i++, j++)
550                         {
551                                 double src_level = (double)i / WINDOW_BORDER;
552                                 double dst_level = (double)(WINDOW_BORDER - i) / WINDOW_BORDER;
553                                 output_buffer[j] = output_buffer[j] * dst_level + out_scale * dsp_out[i] * src_level;
554                         }
556                         for(int i = 0; i < window_size - WINDOW_BORDER; i++)
557                                 output_buffer[output_size + i] = dsp_out[WINDOW_BORDER + i] * out_scale;
558                         output_size += window_size - WINDOW_BORDER;
559                 }
560                 else
561                 {
562 // First buffer has no crossfade
563                         memcpy(output_buffer + output_size, 
564                                 dsp_out, 
565                                 sizeof(double) * window_size);
566                         output_size += window_size;
567                 }
570 // Shift input buffer forward
571                 for(int i = window_size - WINDOW_BORDER, j = 0; 
572                         i < input_size; 
573                         i++, j++)
574                         input_buffer[j] = input_buffer[i];
575                 input_size -= window_size - WINDOW_BORDER;
576         }
579 // Have enough to send to output
580         if(output_size - WINDOW_BORDER >= size)
581         {
582                 memcpy(output_ptr, output_buffer, sizeof(double) * size);
583                 for(int i = size, j = 0; i < output_size; i++, j++)
584                         output_buffer[j] = output_buffer[i];
585                 output_size -= size;
586         }
587         else
588         {
589 //printf("DenoiseEffect::process_realtime 1\n");
590                 bzero(output_ptr, sizeof(double) * size);
591         }
593         return 0;
602 Tree::Tree(int input_length, int levels)
604         this->input_length = input_length;
605         this->levels = levels;
606         int i, j;
608 // create decomposition tree
609         values = new double*[2 * levels];
610         j = input_length;
611         for (i = 0; i < levels; i++)
612         {
613                 j /= 2;
614                 if (j == 0)
615                 {
616                         levels = i;
617                         continue;
618                 }
619                 values[2 * i] = new double[j + 5];
620                 values[2 * i + 1] = new double[j + 5];
621         }
624 Tree::~Tree()
626         int i;
628         for (i = 2 * levels - 1; i >= 0; i--)
629                 delete values[i];
631         delete values;
634 WaveletCoeffs::WaveletCoeffs(double alpha, double beta)
636         int i;
637         double tcosa = cos(alpha);
638         double tcosb = cos(beta);
639         double tsina = sin(alpha);
640         double tsinb = sin(beta);
642 // calculate first two wavelet coefficients  a = a(-2) and b = a(-1)
643         values[0] = ((1.0 + tcosa + tsina) * (1.0 - tcosb - tsinb)
644                                         + 2.0 * tsinb * tcosa) / 4.0;
645         values[1] = ((1.0 - tcosa + tsina) * (1.0 + tcosb - tsinb)
646                                         - 2.0 * tsinb * tcosa) / 4.0;
648         tcosa = cos(alpha - beta);
649         tsina = sin(alpha - beta);
651 // calculate last four wavelet coefficients  c = a(0), d = a(1), 
652 // e = a(2), and f = a(3)
653         values[2]  = (1.0 + tcosa + tsina) / 2.0;
654         values[3]  = (1.0 + tcosa - tsina) / 2.0;
655         values[4]  = 1 - values[0] - values[2];
656         values[5]  = 1 - values[1] - values[3];
658 // zero out very small coefficient values caused by truncation error
659         for (i = 0; i < 6; i++)
660         {
661                 if (fabs(values[i]) < 1.0e-15) values[i] = 0.0;
662         }
665 WaveletCoeffs::~WaveletCoeffs()
670 WaveletFilters::WaveletFilters(WaveletCoeffs *wave_coeffs, wavetype transform)
672         int i, j, k;
674 // find the first non-zero wavelet coefficient
675         i = 0;
676         while(wave_coeffs->values[i] == 0.0) i++;
678 // find the last non-zero wavelet coefficient
679         j = 5;
680         while(wave_coeffs->values[j] == 0.0) j--;
682 // Form the decomposition filters h~ and g~ or the reconstruction
683 // filters h and g.  The division by 2 in the construction
684 // of the decomposition filters is for normalization.
685         length = j - i + 1;
686         for(k = 0; k < length; k++)
687         {
688                 if (transform == DECOMP)
689                 {
690                         h[k] = wave_coeffs->values[j--] / 2.0;
691                         g[k] = (double) (((i++ & 0x01) * 2) - 1) * wave_coeffs->values[i] / 2.0;
692                 }
693                 else
694                 {
695                         h[k] = wave_coeffs->values[i++];
696                         g[k] = (double) (((j-- & 0x01) * 2) - 1) * wave_coeffs->values[j];
697                 }
698         }
700 // clear out the additional array locations, if any
701         while (k < 6)
702         {
703                 h[k] = 0.0;
704                 g[k++] = 0.0;
705         }
708 WaveletFilters::~WaveletFilters()
720 DenoiseConfig::DenoiseConfig()
722         level = 1.0;
725 void DenoiseConfig::copy_from(DenoiseConfig &that)
727         level = that.level;
730 int DenoiseConfig::equivalent(DenoiseConfig &that)
732         return EQUIV(level, that.level);
735 void DenoiseConfig::interpolate(DenoiseConfig &prev, 
736         DenoiseConfig &next, 
737         int64_t prev_frame, 
738         int64_t next_frame, 
739         int64_t current_frame)
741         double next_scale = (double)(current_frame - prev_frame) / (next_frame - prev_frame);
742         double prev_scale = (double)(next_frame - current_frame) / (next_frame - prev_frame);
743         this->level = prev.level * prev_scale + next.level * next_scale;
755 PLUGIN_THREAD_OBJECT(DenoiseEffect, DenoiseThread, DenoiseWindow)
766 DenoiseWindow::DenoiseWindow(DenoiseEffect *plugin, int x, int y)
767  : BC_Window(plugin->gui_string, 
768         x, 
769         y, 
770         150, 
771         50, 
772         150, 
773         50,
774         0, 
775         0,
776         1)
778         this->plugin = plugin;
781 void DenoiseWindow::create_objects()
783         int x = 10, y = 10;
784         
785         add_subwindow(new BC_Title(x, y, _("Level:")));
786         x += 70;
787         add_subwindow(scale = new DenoiseLevel(plugin, x, y));
788         show_window();
789         flush();
792 int DenoiseWindow::close_event()
794 // Set result to 1 to indicate a client side close
795         set_done(1);
796         return 1;
799 void DenoiseWindow::update()
801         scale->update(plugin->config.level);
815 DenoiseLevel::DenoiseLevel(DenoiseEffect *plugin, int x, int y)
816  : BC_FPot(x, y, (float)plugin->config.level, 0, 1.0)
818         this->plugin = plugin;
819         set_precision(0.01);
822 int DenoiseLevel::handle_event()
824         plugin->config.level = get_value();
825         plugin->send_configure_change();
826         return 1;