r125: This commit was manufactured by cvs2svn to create tag 'r1_1_7-last'.
[cinelerra_cv/mob.git] / hvirtual / plugins / denoise / denoise.C
blob316ff2bb24d92814d58f57d0b632ce2423eb13c3
1 #include "bcdisplayinfo.h"
2 #include "clip.h"
3 #include "defaults.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()
119         return _("Denoise");
123 int DenoiseEffect::is_realtime()
125         return 1;
130 void DenoiseEffect::read_data(KeyFrame *keyframe)
132         FileXML input;
133         input.set_shared_string(keyframe->data, strlen(keyframe->data));
135         int result = 0;
136         while(!result)
137         {
138                 result = input.read_tag();
140                 if(!result)
141                 {
142                         if(input.tag.title_is("DENOISE"))
143                         {
144                                 config.level = input.tag.get_property("LEVEL", config.level);
145                         }
146                 }
147         }
150 void DenoiseEffect::save_data(KeyFrame *keyframe)
152         FileXML output;
153         output.set_shared_string(keyframe->data, MESSAGESIZE);
155         output.tag.set_title("DENOISE");
156         output.tag.set_property("LEVEL", config.level);
157         output.append_tag();
158         output.append_newline();
160         output.terminate_string();
163 int DenoiseEffect::load_defaults()
165         char directory[BCTEXTLEN], string[BCTEXTLEN];
166         sprintf(directory, "%sdenoise.rc", BCASTDIR);
167         defaults = new Defaults(directory);
168         defaults->load();
169         
170         config.level = defaults->get("LEVEL", config.level);
171         return 0;
174 int DenoiseEffect::save_defaults()
176         char string[BCTEXTLEN];
178         defaults->update("LEVEL", config.level);
179         defaults->save();
181         return 0;
184 void DenoiseEffect::update_gui()
186         if(thread)
187         {
188                 thread->window->lock_window();
189                 thread->window->update();
190                 thread->window->unlock_window();
191         }
196 double DenoiseEffect::dot_product(double *data, double *filter, char filtlen)
198         static int i;
199         static double sum;
201         sum = 0.0;
202         for(i = 0; i < filtlen; i++) sum += *data-- * *filter++;
203         return sum;
206 int DenoiseEffect::convolve_dec_2(double *input_sequence, 
207         int64_t length,
208         double *filter, 
209         int filtlen, 
210         double *output_sequence)
212 // convolve the input sequence with the filter and decimate by two
213         int i, shortlen, offset;
214         int64_t lengthp4 = length + 4;
215         int64_t lengthm4 = length - 4;
216         int64_t lengthp5 = length + 5;
217         int64_t lengthp8 = length + 8;
219         for(i = 0; (i <= lengthp8) && ((i - filtlen) <= lengthp8); i += 2)
220         {
221                 if(i < filtlen)
222                         *output_sequence++ = dot_product(input_sequence + i, filter, i + 1);
223                 else 
224                 if(i > lengthp5)
225                 {
226                         offset = i - lengthm4;
227                         shortlen = filtlen - offset;
228                         *output_sequence++ = dot_product(input_sequence + lengthp4,
229                                                                 filter + offset, shortlen);
230                 }
231                 else
232                         *output_sequence++ = dot_product(input_sequence + i, filter, filtlen);
233         }
234         return 0;
237 int64_t DenoiseEffect::decompose_branches(double *in_data, 
238         int64_t length, 
239         WaveletFilters *decomp_filter, 
240         double *out_low, 
241         double *out_high)
243 // Take input data and filters and form two branches of half the
244 // original length. Length of branches is returned.
245         convolve_dec_2(in_data, length, decomp_filter->h, decomp_filter->length, out_low);
246         convolve_dec_2(in_data, length, decomp_filter->g, decomp_filter->length, out_high);
247         return (length / 2);
250 int DenoiseEffect::wavelet_decomposition(double *in_data, 
251         int64_t in_length, 
252         double **out_data)
254         for(int i = 0; i < levels; i++)
255         {
256                 in_length = decompose_branches(in_data, 
257                         in_length, 
258                         decomp_filter, 
259                         out_data[2 * i], 
260                         out_data[(2 * i) + 1]);
262                 in_data = out_data[2 * i];
263         }
264         return 0;
267 int DenoiseEffect::tree_copy(double **output, 
268         double **input, 
269         int length, 
270         int levels)
272         register int i, j, k, l, m;
274         for(i = 0, k = 1; k < levels; i++, k++)
275         {
276                 length /= 2;
277                 l = 2 * i;
278                 m = l + 1;
280                 for(j = 0; j < length + 5; j++)
281                 {
282                         output[l][j] = 0.0;
283                         output[m][j] = input[m][j];
284                 }
285         }
287         length /= 2;
288         l = 2 * i;
289         m = l + 1;
291         for(j = 0; j < length + 5; j++)
292         {
293                 output[l][j] = input[l][j];
294                 output[m][j] = input[m][j];
295         }
296         return 0;
299 int DenoiseEffect::threshold(int window_size, double gammas, int levels)
301         int i, j;
302         double threshold, cv, cvb, abs_coeff_r;
303         double *coeff_r, *coeff_l;
304         int length;
306         for(i = 0; i < levels; i++) 
307         {
308                 length = (window_size >> (i + 1)) + 5;
309                 threshold = sqrt(2 * log(length) / log(2)) * gammas / sqrt(length);
311                 for(j = 0; j < length; j++) 
312                 {
313                         coeff_r = &(ex_coeff_r->values[(2 * i) + 1][j]);
314                         coeff_l = &(ex_coeff_rn->values[(2 * i) + 1][j]);
316                         cv = SGN(*coeff_r);
317                         abs_coeff_r = fabs(*coeff_r);
318                         cvb = abs_coeff_r - threshold;
319                         cv *= cvb;
321                         if(abs_coeff_r > threshold) 
322                         {
323                                 *coeff_r = cv;
324                                 *coeff_l = 0.0;
325                         }
326                         else 
327                         {
328                                 *coeff_l = *coeff_r;
329                                 *coeff_r = 0.0;
330                         }
331                 }
332         }
336 double DenoiseEffect::dot_product_even(double *data, double *filter, int filtlen)
338         static int i;
339         static double sum;
341         sum = 0.0;
342         for(i = 0; i < filtlen; i += 2) sum += *data-- * filter[i];
343         return sum;
347 double DenoiseEffect::dot_product_odd(double *data, double *filter, int filtlen)
349         static int i;
350         static double sum;
352         sum = 0.0;
353         for(i = 1; i < filtlen; i += 2) sum += *data-- * filter[i];
354         return sum;
357 int DenoiseEffect::convolve_int_2(double *input_sequence, 
358         int64_t length, 
359         double *filter, 
360         int filtlen, 
361         int sum_output, 
362         double *output_sequence)
363 // insert zeros between each element of the input sequence and
364 // convolve with the filter to interpolate the data
366         register int i, j;
367         int endpoint = length + filtlen - 2;
369         if (sum_output)
370         {
371 // summation with previous convolution
372 // every other dot product interpolates the data
373                 for(i = (filtlen / 2) - 1, j = (filtlen / 2); i < endpoint; i++, j++)
374                 {
375                         *output_sequence++ += dot_product_odd(input_sequence + i, filter, filtlen);
376                         *output_sequence++ += dot_product_even(input_sequence + j, filter, filtlen);
377                 }
379                 *output_sequence++ += dot_product_odd(input_sequence + i, filter, filtlen);
380         }
381         else
382         {
383 // first convolution of pair
384 // every other dot product interpolates the data
385                 for(i = (filtlen / 2) - 1, j = (filtlen / 2); i < endpoint; i++, j++)
386                 {
387                         *output_sequence++ = dot_product_odd(input_sequence + i, filter, filtlen);
388                         *output_sequence++ = dot_product_even(input_sequence + j, filter, filtlen);
389                 }
391                 *output_sequence++ = dot_product_odd(input_sequence + i, filter, filtlen);
392         }
393         return 0;
397 int64_t DenoiseEffect::reconstruct_branches(double *in_low, 
398         double *in_high, 
399         int64_t in_length,
400         WaveletFilters *recon_filter, 
401         double *output)
403 // take input data and filters and form two branches of half the
404 // original length. length of branches is returned
405         convolve_int_2(in_low, in_length, recon_filter->h, 
406                                         recon_filter->length, 0, output);
407         convolve_int_2(in_high, in_length, recon_filter->g, 
408                                         recon_filter->length, 1, output);
409         return in_length * 2;
412 int DenoiseEffect::wavelet_reconstruction(double **in_data, 
413         int64_t in_length, 
414         double *out_data)
416         double *output;
417         int i;
419         in_length = in_length >> levels;
420 // destination of all but last branch reconstruction is the next
421 // higher intermediate approximation
422         for(i = levels - 1; i > 0; i--)
423         {
424                 output = in_data[2 * (i - 1)];
425                 in_length = reconstruct_branches(in_data[2 * i], 
426                         in_data[(2 * i) + 1],
427                         in_length, 
428                         recon_filter, 
429                         output);
430         }
432 // destination of the last branch reconstruction is the output data
433         reconstruct_branches(in_data[0], 
434                 in_data[1], 
435                 in_length, 
436                 recon_filter, 
437                 out_data);
439         return 0;
442 void DenoiseEffect::process_window()
444         int i, j;
445         for(j = 0; j < iterations; j++)
446         {
447                 wavelet_decomposition(dsp_in, window_size, ex_coeff_d->values);
449                 tree_copy(ex_coeff_r->values, ex_coeff_d->values, window_size, levels);
450                 tree_copy(ex_coeff_rn->values, ex_coeff_d->values, window_size, levels);
452 // qualify coeffs
453 //printf("DenoiseEffect::process_window %f\n", config.level);
454                 threshold(window_size, config.level * 10.0, levels);
456                 wavelet_reconstruction(ex_coeff_r->values, window_size, dsp_iteration);
457                 wavelet_reconstruction(ex_coeff_rn->values, window_size, dsp_in);
459                 for(i = 0; i < window_size; i++)
460                         dsp_out[i] += dsp_iteration[i];
461         }
467 int DenoiseEffect::process_realtime(int64_t size, double *input_ptr, double *output_ptr)
469         load_configuration();
471         if(!initialized)
472         {
473                 int64_t size_factor = (int)(pow(2, levels));
474                 dsp_in = new double[window_size * size_factor];
475                 dsp_out = new double[window_size * 2];
476                 dsp_iteration = new double[window_size * 2];
479                 ex_coeff_d = new Tree(window_size, levels);
480                 ex_coeff_r = new Tree(window_size, levels);
481                 ex_coeff_rn = new Tree(window_size, levels);
482                 wave_coeff_d = new WaveletCoeffs(alpha, beta);
483                 wave_coeff_r = new WaveletCoeffs(alpha, beta);
484                 decomp_filter = new WaveletFilters(wave_coeff_d, DECOMP);
485                 recon_filter = new WaveletFilters(wave_coeff_r, RECON);
486                 in_scale = 65535 / sqrt(window_size) / iterations;
487                 out_scale = output_level / 65535 * sqrt(window_size);
488                 initialized = 1;
489         }
490         
491 // Append input buffer
492         if(input_size + size > input_allocation)
493         {
494                 double *new_input = new double[input_size + size];
495                 if(input_buffer)
496                 {
497                         memcpy(new_input, input_buffer, sizeof(double) * input_size);
498                         delete [] input_buffer;
499                 }
500                 input_buffer = new_input;
501                 input_allocation = input_size + size;
502         }
503         memcpy(input_buffer + input_size, 
504                 input_ptr, 
505                 size * sizeof(double));
506         input_size += size;
509 // Have enough to do some windows
510         while(input_size >= window_size)
511         {
512 // Load dsp_in
513                 for(int i = 0; i < window_size; i++)
514                 {
515                         dsp_in[i] = input_buffer[i] * in_scale;
516                 }
517                 bzero(dsp_out, sizeof(double) * window_size);
524 // First window produces garbage
525                 if(!first_window)
526                         process_window();
527                 first_window = 0;
534 // Crossfade into the output buffer
535                 int64_t new_allocation = output_size + window_size;
536                 if(new_allocation > output_allocation)
537                 {
538                         double *new_output = new double[new_allocation];
540                         if(output_buffer)
541                         {
542                                 memcpy(new_output, output_buffer, sizeof(double) * output_size);
543 //printf("CrossfadeFFT::process_fifo 1 %p\n", output_buffer);
544                                 delete [] output_buffer;
545 //printf("CrossfadeFFT::process_fifo 2\n");
546                         }
547                         output_buffer = new_output;
548                         output_allocation = new_allocation;
549                 }
551                 if(output_size >= WINDOW_BORDER)
552                 {
553                         for(int i = 0, j = output_size - WINDOW_BORDER; 
554                                 i < WINDOW_BORDER; 
555                                 i++, j++)
556                         {
557                                 double src_level = (double)i / WINDOW_BORDER;
558                                 double dst_level = (double)(WINDOW_BORDER - i) / WINDOW_BORDER;
559                                 output_buffer[j] = output_buffer[j] * dst_level + out_scale * dsp_out[i] * src_level;
560                         }
562                         for(int i = 0; i < window_size - WINDOW_BORDER; i++)
563                                 output_buffer[output_size + i] = dsp_out[WINDOW_BORDER + i] * out_scale;
564                         output_size += window_size - WINDOW_BORDER;
565                 }
566                 else
567                 {
568 // First buffer has no crossfade
569                         memcpy(output_buffer + output_size, 
570                                 dsp_out, 
571                                 sizeof(double) * window_size);
572                         output_size += window_size;
573                 }
576 // Shift input buffer forward
577                 for(int i = window_size - WINDOW_BORDER, j = 0; 
578                         i < input_size; 
579                         i++, j++)
580                         input_buffer[j] = input_buffer[i];
581                 input_size -= window_size - WINDOW_BORDER;
582         }
585 // Have enough to send to output
586         if(output_size - WINDOW_BORDER >= size)
587         {
588                 memcpy(output_ptr, output_buffer, sizeof(double) * size);
589                 for(int i = size, j = 0; i < output_size; i++, j++)
590                         output_buffer[j] = output_buffer[i];
591                 output_size -= size;
592         }
593         else
594         {
595 //printf("DenoiseEffect::process_realtime 1\n");
596                 bzero(output_ptr, sizeof(double) * size);
597         }
599         return 0;
608 Tree::Tree(int input_length, int levels)
610         this->input_length = input_length;
611         this->levels = levels;
612         int i, j;
614 // create decomposition tree
615         values = new double*[2 * levels];
616         j = input_length;
617         for (i = 0; i < levels; i++)
618         {
619                 j /= 2;
620                 if (j == 0)
621                 {
622                         levels = i;
623                         continue;
624                 }
625                 values[2 * i] = new double[j + 5];
626                 values[2 * i + 1] = new double[j + 5];
627         }
630 Tree::~Tree()
632         int i;
634         for (i = 2 * levels - 1; i >= 0; i--)
635                 delete values[i];
637         delete values;
640 WaveletCoeffs::WaveletCoeffs(double alpha, double beta)
642         int i;
643         double tcosa = cos(alpha);
644         double tcosb = cos(beta);
645         double tsina = sin(alpha);
646         double tsinb = sin(beta);
648 // calculate first two wavelet coefficients  a = a(-2) and b = a(-1)
649         values[0] = ((1.0 + tcosa + tsina) * (1.0 - tcosb - tsinb)
650                                         + 2.0 * tsinb * tcosa) / 4.0;
651         values[1] = ((1.0 - tcosa + tsina) * (1.0 + tcosb - tsinb)
652                                         - 2.0 * tsinb * tcosa) / 4.0;
654         tcosa = cos(alpha - beta);
655         tsina = sin(alpha - beta);
657 // calculate last four wavelet coefficients  c = a(0), d = a(1), 
658 // e = a(2), and f = a(3)
659         values[2]  = (1.0 + tcosa + tsina) / 2.0;
660         values[3]  = (1.0 + tcosa - tsina) / 2.0;
661         values[4]  = 1 - values[0] - values[2];
662         values[5]  = 1 - values[1] - values[3];
664 // zero out very small coefficient values caused by truncation error
665         for (i = 0; i < 6; i++)
666         {
667                 if (fabs(values[i]) < 1.0e-15) values[i] = 0.0;
668         }
671 WaveletCoeffs::~WaveletCoeffs()
676 WaveletFilters::WaveletFilters(WaveletCoeffs *wave_coeffs, wavetype transform)
678         int i, j, k;
680 // find the first non-zero wavelet coefficient
681         i = 0;
682         while(wave_coeffs->values[i] == 0.0) i++;
684 // find the last non-zero wavelet coefficient
685         j = 5;
686         while(wave_coeffs->values[j] == 0.0) j--;
688 // Form the decomposition filters h~ and g~ or the reconstruction
689 // filters h and g.  The division by 2 in the construction
690 // of the decomposition filters is for normalization.
691         length = j - i + 1;
692         for(k = 0; k < length; k++)
693         {
694                 if (transform == DECOMP)
695                 {
696                         h[k] = wave_coeffs->values[j--] / 2.0;
697                         g[k] = (double) (((i++ & 0x01) * 2) - 1) * wave_coeffs->values[i] / 2.0;
698                 }
699                 else
700                 {
701                         h[k] = wave_coeffs->values[i++];
702                         g[k] = (double) (((j-- & 0x01) * 2) - 1) * wave_coeffs->values[j];
703                 }
704         }
706 // clear out the additional array locations, if any
707         while (k < 6)
708         {
709                 h[k] = 0.0;
710                 g[k++] = 0.0;
711         }
714 WaveletFilters::~WaveletFilters()
726 DenoiseConfig::DenoiseConfig()
728         level = 1.0;
731 void DenoiseConfig::copy_from(DenoiseConfig &that)
733         level = that.level;
736 int DenoiseConfig::equivalent(DenoiseConfig &that)
738         return EQUIV(level, that.level);
741 void DenoiseConfig::interpolate(DenoiseConfig &prev, 
742         DenoiseConfig &next, 
743         int64_t prev_frame, 
744         int64_t next_frame, 
745         int64_t current_frame)
747         double next_scale = (double)(current_frame - prev_frame) / (next_frame - prev_frame);
748         double prev_scale = (double)(next_frame - current_frame) / (next_frame - prev_frame);
749         this->level = prev.level * prev_scale + next.level * next_scale;
761 PLUGIN_THREAD_OBJECT(DenoiseEffect, DenoiseThread, DenoiseWindow)
772 DenoiseWindow::DenoiseWindow(DenoiseEffect *plugin, int x, int y)
773  : BC_Window(plugin->gui_string, 
774         x, 
775         y, 
776         150, 
777         50, 
778         150, 
779         50,
780         0, 
781         0,
782         1)
784         this->plugin = plugin;
787 void DenoiseWindow::create_objects()
789         int x = 10, y = 10;
790         
791         add_subwindow(new BC_Title(x, y, _("Level:")));
792         x += 70;
793         add_subwindow(scale = new DenoiseLevel(plugin, x, y));
794         show_window();
795         flush();
798 int DenoiseWindow::close_event()
800 // Set result to 1 to indicate a client side close
801         set_done(1);
802         return 1;
805 void DenoiseWindow::update()
807         scale->update(plugin->config.level);
821 DenoiseLevel::DenoiseLevel(DenoiseEffect *plugin, int x, int y)
822  : BC_FPot(x, y, (float)plugin->config.level, 0, 1.0)
824         this->plugin = plugin;
825         set_precision(0.01);
828 int DenoiseLevel::handle_event()
830         plugin->config.level = get_value();
831         plugin->send_configure_change();
832         return 1;