1 
2 #ifdef HAVE_CONFIG_H
3 #include "config.h"
4 #endif
5 
6 #include <schroedinger/schrofilter.h>
7 #include <schroedinger/schrodebug.h>
8 #include <schroedinger/schrowavelet.h>
9 #include <schroedinger/schrotables.h>
10 #include <schroedinger/schrobitstream.h>
11 #include <schroedinger/schrohistogram.h>
12 #include <schroedinger/schroparams.h>
13 #include <schroedinger/schrovirtframe.h>
14 
15 #include <string.h>
16 #include <math.h>
17 
18 static void
sort_u8(uint8_t * d,int n)19 sort_u8 (uint8_t * d, int n)
20 {
21   int start = 0;
22   int end = n;
23   int i;
24   int x;
25 
26   /* OMG bubble sort! */
27   while (start < end) {
28     for (i = start; i < end - 1; i++) {
29       if (d[i] > d[i + 1]) {
30         x = d[i];
31         d[i] = d[i + 1];
32         d[i + 1] = x;
33       }
34     }
35     end--;
36     for (i = end - 2; i >= start; i--) {
37       if (d[i] > d[i + 1]) {
38         x = d[i];
39         d[i] = d[i + 1];
40         d[i + 1] = x;
41       }
42     }
43     start++;
44   }
45 }
46 
47 #ifdef unused
48 /* reference */
49 void
schro_filter_cwmN_ref(uint8_t * d,uint8_t * s1,uint8_t * s2,uint8_t * s3,int n,int weight)50 schro_filter_cwmN_ref (uint8_t * d, uint8_t * s1, uint8_t * s2, uint8_t * s3,
51     int n, int weight)
52 {
53   int i;
54   int j;
55   uint8_t list[8 + 12];
56 
57   for (i = 0; i < n; i++) {
58     list[0] = s1[i + 0];
59     list[1] = s1[i + 1];
60     list[2] = s1[i + 2];
61     list[3] = s2[i + 0];
62     list[4] = s2[i + 2];
63     list[5] = s3[i + 0];
64     list[6] = s3[i + 1];
65     list[7] = s3[i + 2];
66     for (j = 0; j < weight; j++) {
67       list[8 + j] = s2[i + 1];
68     }
69 
70     sort_u8 (list, 8 + weight);
71 
72     d[i] = list[(8 + weight) / 2];
73   }
74 }
75 #endif
76 
77 void
schro_filter_cwmN(uint8_t * d,uint8_t * s1,uint8_t * s2,uint8_t * s3,int n,int weight)78 schro_filter_cwmN (uint8_t * d, uint8_t * s1, uint8_t * s2, uint8_t * s3, int n,
79     int weight)
80 {
81   int i;
82   int j;
83   uint8_t list[8 + 12];
84   int low, hi;
85 
86   for (i = 0; i < n; i++) {
87     list[0] = s1[i + 0];
88     list[1] = s1[i + 1];
89     list[2] = s1[i + 2];
90     list[3] = s2[i + 0];
91     list[4] = s2[i + 2];
92     list[5] = s3[i + 0];
93     list[6] = s3[i + 1];
94     list[7] = s3[i + 2];
95 
96     low = 0;
97     hi = 0;
98     for (j = 0; j < 8; j++) {
99       if (list[j] < s2[i + 1])
100         low++;
101       if (list[j] > s2[i + 1])
102         hi++;
103     }
104 
105     if (low < ((9 - weight) / 2) || hi < ((9 - weight) / 2)) {
106       for (j = 0; j < weight; j++) {
107         list[8 + j] = s2[i + 1];
108       }
109 
110       sort_u8 (list, 8 + weight);
111 
112       d[i] = list[(8 + weight) / 2];
113     } else {
114       d[i] = s2[i + 1];
115     }
116   }
117 }
118 
119 static void
schro_frame_component_filter_cwmN(SchroFrameData * comp,int weight)120 schro_frame_component_filter_cwmN (SchroFrameData * comp, int weight)
121 {
122   int i;
123   uint8_t *tmp;
124   uint8_t *tmp1;
125   uint8_t *tmp2;
126 
127   tmp1 = schro_malloc (comp->width);
128   tmp2 = schro_malloc (comp->width);
129 
130   schro_filter_cwmN (tmp1,
131       OFFSET (comp->data, comp->stride * 0),
132       OFFSET (comp->data, comp->stride * 1),
133       OFFSET (comp->data, comp->stride * 2), comp->width - 2, weight);
134   schro_filter_cwmN (tmp2,
135       OFFSET (comp->data, comp->stride * 1),
136       OFFSET (comp->data, comp->stride * 2),
137       OFFSET (comp->data, comp->stride * 3), comp->width - 2, weight);
138 
139   for (i = 3; i < comp->height - 1; i++) {
140     memcpy (OFFSET (comp->data, comp->stride * (i - 2) + 1),
141         tmp1, comp->width - 2);
142     tmp = tmp1;
143     tmp1 = tmp2;
144     tmp2 = tmp;
145 
146     schro_filter_cwmN (tmp2,
147         OFFSET (comp->data, comp->stride * (i - 1)),
148         OFFSET (comp->data, comp->stride * (i + 0)),
149         OFFSET (comp->data, comp->stride * (i + 1)), comp->width - 2, weight);
150   }
151   memcpy (OFFSET (comp->data, comp->stride * (i - 2) + 1),
152       tmp1, comp->width - 2);
153   memcpy (OFFSET (comp->data, comp->stride * (i - 1) + 1),
154       tmp2, comp->width - 2);
155 
156   schro_free (tmp1);
157   schro_free (tmp2);
158 }
159 
160 void
schro_frame_filter_cwmN(SchroFrame * frame,int weight)161 schro_frame_filter_cwmN (SchroFrame * frame, int weight)
162 {
163   schro_frame_component_filter_cwmN (&frame->components[0], weight);
164   schro_frame_component_filter_cwmN (&frame->components[1], weight);
165   schro_frame_component_filter_cwmN (&frame->components[2], weight);
166 }
167 
168 
169 #ifdef unused
170 static void
schro_frame_component_filter_cwmN_ref(SchroFrameData * comp,int weight)171 schro_frame_component_filter_cwmN_ref (SchroFrameData * comp, int weight)
172 {
173   int i;
174   uint8_t *tmp;
175   uint8_t *tmp1;
176   uint8_t *tmp2;
177 
178   tmp1 = schro_malloc (comp->width);
179   tmp2 = schro_malloc (comp->width);
180 
181   schro_filter_cwmN_ref (tmp1,
182       OFFSET (comp->data, comp->stride * 0),
183       OFFSET (comp->data, comp->stride * 1),
184       OFFSET (comp->data, comp->stride * 2), comp->width - 2, weight);
185   schro_filter_cwmN_ref (tmp2,
186       OFFSET (comp->data, comp->stride * 1),
187       OFFSET (comp->data, comp->stride * 2),
188       OFFSET (comp->data, comp->stride * 3), comp->width - 2, weight);
189 
190   for (i = 3; i < comp->height - 1; i++) {
191     memcpy (OFFSET (comp->data, comp->stride * (i - 2) + 1),
192         tmp1, comp->width - 2);
193     tmp = tmp1;
194     tmp1 = tmp2;
195     tmp2 = tmp;
196 
197     schro_filter_cwmN_ref (tmp2,
198         OFFSET (comp->data, comp->stride * (i - 1)),
199         OFFSET (comp->data, comp->stride * (i + 0)),
200         OFFSET (comp->data, comp->stride * (i + 1)), comp->width - 2, weight);
201   }
202   memcpy (OFFSET (comp->data, comp->stride * (i - 2) + 1),
203       tmp1, comp->width - 2);
204   memcpy (OFFSET (comp->data, comp->stride * (i - 1) + 1),
205       tmp2, comp->width - 2);
206 
207   schro_free (tmp1);
208   schro_free (tmp2);
209 }
210 #endif
211 
212 #ifdef unused
213 void
schro_frame_filter_cwmN_ref(SchroFrame * frame,int weight)214 schro_frame_filter_cwmN_ref (SchroFrame * frame, int weight)
215 {
216   schro_frame_component_filter_cwmN_ref (&frame->components[0], weight);
217   schro_frame_component_filter_cwmN_ref (&frame->components[1], weight);
218   schro_frame_component_filter_cwmN_ref (&frame->components[2], weight);
219 }
220 #endif
221 
222 
223 #if 0
224 /* reference */
225 void
226 schro_filter_cwm7 (uint8_t * d, uint8_t * s1, uint8_t * s2, uint8_t * s3, int n)
227 {
228   int i;
229   int min, max;
230 
231   for (i = 0; i < n; i++) {
232     min = MIN (s1[i + 0], s1[i + 1]);
233     max = MAX (s1[i + 0], s1[i + 1]);
234     min = MIN (min, s1[i + 2]);
235     max = MAX (max, s1[i + 2]);
236     min = MIN (min, s2[i + 0]);
237     max = MAX (max, s2[i + 0]);
238     min = MIN (min, s2[i + 2]);
239     max = MAX (max, s2[i + 2]);
240     min = MIN (min, s3[i + 0]);
241     max = MAX (max, s3[i + 0]);
242     min = MIN (min, s3[i + 1]);
243     max = MAX (max, s3[i + 1]);
244     min = MIN (min, s3[i + 2]);
245     max = MAX (max, s3[i + 2]);
246 
247     d[i] = MIN (max, MAX (min, s2[i + 1]));
248   }
249 }
250 #endif
251 
252 #ifdef unused
253 /* FIXME move to schrooil */
254 void
schro_filter_cwm7(uint8_t * d,uint8_t * s1,uint8_t * s2,uint8_t * s3,int n)255 schro_filter_cwm7 (uint8_t * d, uint8_t * s1, uint8_t * s2, uint8_t * s3, int n)
256 {
257   int i;
258   int min, max;
259 
260   for (i = 0; i < n; i++) {
261     if (s1[i + 0] < s2[i + 1]) {
262       max = MAX (s1[i + 0], s1[i + 1]);
263       max = MAX (max, s1[i + 2]);
264       max = MAX (max, s2[i + 0]);
265       max = MAX (max, s2[i + 2]);
266       max = MAX (max, s3[i + 0]);
267       max = MAX (max, s3[i + 1]);
268       max = MAX (max, s3[i + 2]);
269       d[i] = MIN (max, s2[i + 1]);
270     } else if (s1[i + 0] > s2[i + 1]) {
271       min = MIN (s1[i + 0], s1[i + 1]);
272       min = MIN (min, s1[i + 2]);
273       min = MIN (min, s2[i + 0]);
274       min = MIN (min, s2[i + 2]);
275       min = MIN (min, s3[i + 0]);
276       min = MIN (min, s3[i + 1]);
277       min = MIN (min, s3[i + 2]);
278       d[i] = MAX (min, s2[i + 1]);
279     } else {
280       d[i] = s2[i + 1];
281     }
282   }
283 }
284 #endif
285 
286 #ifdef unused
287 static void
schro_frame_component_filter_cwm7(SchroFrameData * comp)288 schro_frame_component_filter_cwm7 (SchroFrameData * comp)
289 {
290   int i;
291   uint8_t *tmp;
292   uint8_t *tmp1;
293   uint8_t *tmp2;
294 
295   tmp1 = schro_malloc (comp->width);
296   tmp2 = schro_malloc (comp->width);
297 
298   schro_filter_cwm7 (tmp1,
299       OFFSET (comp->data, comp->stride * 0),
300       OFFSET (comp->data, comp->stride * 1),
301       OFFSET (comp->data, comp->stride * 2), comp->width - 2);
302   schro_filter_cwm7 (tmp2,
303       OFFSET (comp->data, comp->stride * 1),
304       OFFSET (comp->data, comp->stride * 2),
305       OFFSET (comp->data, comp->stride * 3), comp->width - 2);
306 
307   for (i = 3; i < comp->height - 1; i++) {
308     memcpy (OFFSET (comp->data, comp->stride * (i - 2) + 1),
309         tmp1, comp->width - 2);
310     tmp = tmp1;
311     tmp1 = tmp2;
312     tmp2 = tmp;
313 
314     schro_filter_cwm7 (tmp2,
315         OFFSET (comp->data, comp->stride * (i - 1)),
316         OFFSET (comp->data, comp->stride * (i + 0)),
317         OFFSET (comp->data, comp->stride * (i + 1)), comp->width - 2);
318   }
319   memcpy (OFFSET (comp->data, comp->stride * (i - 2) + 1),
320       tmp1, comp->width - 2);
321   memcpy (OFFSET (comp->data, comp->stride * (i - 1) + 1),
322       tmp2, comp->width - 2);
323 
324   schro_free (tmp1);
325   schro_free (tmp2);
326 }
327 #endif
328 
329 #ifdef unused
330 void
schro_frame_filter_cwm7(SchroFrame * frame)331 schro_frame_filter_cwm7 (SchroFrame * frame)
332 {
333   schro_frame_component_filter_cwm7 (&frame->components[0]);
334   schro_frame_component_filter_cwm7 (&frame->components[1]);
335   schro_frame_component_filter_cwm7 (&frame->components[2]);
336 }
337 #endif
338 
339 
340 static void
lowpass3_h_u8(SchroFrame * frame,void * _dest,int component,int i)341 lowpass3_h_u8 (SchroFrame *frame, void *_dest, int component, int i)
342 {
343   uint8_t *dest = _dest;
344   uint8_t *src;
345   int tap1, tap2;
346 
347   tap1 = *(int *)frame->virt_priv2;
348   tap2 = 256 - 2*tap1;
349 
350   src = schro_virt_frame_get_line (frame->virt_frame1, component, i);
351 
352   if (component > 0) {
353     memcpy (dest, src, frame->components[component].width);
354     return;
355   }
356 
357   i = 0;
358   dest[i] = (src[i] * tap1 + src[i] * tap2 + src[i+1] * tap1 + 128)>>8;
359 
360   for(i=1;i<frame->width-1;i++) {
361     dest[i] = (src[i-1] * tap1 + src[i] * tap2 + src[i+1] * tap1 + 128)>>8;
362   }
363 
364   i = frame->width - 1;
365   dest[i] = (src[i-1] * tap1 + src[i] * tap2 + src[i] * tap1 + 128)>>8;
366 
367 }
368 
369 static void
lowpass3_v_u8(SchroFrame * frame,void * _dest,int component,int i)370 lowpass3_v_u8 (SchroFrame *frame, void *_dest, int component, int i)
371 {
372   uint8_t *dest = _dest;
373   uint8_t *src1, *src2, *src3;
374   int tap1, tap2;
375 
376   if (component > 0) {
377     src2 = schro_virt_frame_get_line (frame->virt_frame1, component, i);
378     memcpy (dest, src2, frame->components[component].width);
379     return;
380   }
381 
382   tap1 = *(int *)frame->virt_priv2;
383   tap2 = 256 - 2*tap1;
384 
385   src1 = schro_virt_frame_get_line (frame->virt_frame1, component,
386       CLAMP(i-1,0,frame->height));
387   src2 = schro_virt_frame_get_line (frame->virt_frame1, component, i);
388   src3 = schro_virt_frame_get_line (frame->virt_frame1, component,
389       CLAMP(i+1,0,frame->height));
390 
391   for(i=0;i<frame->width;i++) {
392     dest[i] = (src1[i] * tap1 + src2[i] * tap2 + src3[i] * tap1 + 128)>>8;
393   }
394 
395 }
396 
397 void
schro_frame_filter_lowpass(SchroFrame * frame,int tap)398 schro_frame_filter_lowpass (SchroFrame * frame, int tap)
399 {
400   SchroFrame *vf;
401   SchroFrame *vf2;
402   SchroFrame *dup;
403 
404   dup = schro_frame_dup (frame);
405 
406   vf = schro_frame_new_virtual (NULL, frame->format, frame->width, frame->height);
407   vf->virt_frame1 = schro_frame_ref(frame);
408   vf->render_line = lowpass3_h_u8;
409   vf->virt_priv2 = (void *) &tap;
410 
411   vf2 = schro_frame_new_virtual (NULL, frame->format, frame->width, frame->height);
412   vf2->virt_frame1 = vf;
413   vf2->render_line = lowpass3_v_u8;
414   vf2->virt_priv2 = (void *) &tap;
415 
416   schro_virt_frame_render (vf2, dup);
417   schro_frame_convert (frame, dup);
418 
419   schro_frame_unref (vf2);
420   schro_frame_unref (dup);
421 }
422 
423 
424 #ifdef unused
425 static void
lowpass_s16(int16_t * d,int16_t * s,int n)426 lowpass_s16 (int16_t * d, int16_t * s, int n)
427 {
428   int i;
429   int j;
430   int x;
431   const int32_t taps[] = { 2, 9, 28, 55, 68, 55, 28, 9, 2, 0 };
432   const int32_t offsetshift[] = { 128, 8 };
433 
434   for (i = 0; i < 4; i++) {
435     x = 0;
436     for (j = 0; j < 9; j++) {
437       x += s[CLAMP (i + j - 4, 0, n - 1)] * taps[j];
438     }
439     d[i] = (x + 128) >> 8;
440   }
441   schro_mas10_s16 (d + 4, s, taps, offsetshift, n - 9);
442   for (i = n - 6; i < n; i++) {
443     x = 0;
444     for (j = 0; j < 9; j++) {
445       x += s[CLAMP (i + j - 4, 0, n - 1)] * taps[j];
446     }
447     d[i] = (x + 128) >> 8;
448   }
449 }
450 #endif
451 
452 #ifdef unused
453 /* FIXME move to schrooil */
454 static void
lowpass_vert_s16(int16_t * d,int16_t * s,int n)455 lowpass_vert_s16 (int16_t * d, int16_t * s, int n)
456 {
457   int i;
458   int j;
459   int x;
460   static int taps[] = { 2, 9, 28, 55, 68, 55, 28, 9, 2, 0 };
461 
462   for (i = 0; i < n; i++) {
463     x = 0;
464     for (j = 0; j < 9; j++) {
465       x += s[j * n + i] * taps[j];
466     }
467     d[i] = (x + 128) >> 8;
468   }
469 }
470 #endif
471 
472 
473 #ifdef unused
474 static void
schro_frame_component_filter_lowpass_16(SchroFrameData * comp)475 schro_frame_component_filter_lowpass_16 (SchroFrameData * comp)
476 {
477   int i;
478   int16_t *tmp;
479 
480   tmp = schro_malloc (comp->width * 9 * sizeof (int16_t));
481 
482   lowpass_s16 (tmp + 0 * comp->width,
483       OFFSET (comp->data, comp->stride * 0), comp->width);
484   memcpy (tmp + 1 * comp->width, tmp + 0 * comp->width, comp->width * 2);
485   memcpy (tmp + 2 * comp->width, tmp + 0 * comp->width, comp->width * 2);
486   memcpy (tmp + 3 * comp->width, tmp + 0 * comp->width, comp->width * 2);
487   memcpy (tmp + 4 * comp->width, tmp + 0 * comp->width, comp->width * 2);
488   lowpass_s16 (tmp + 5 * comp->width,
489       OFFSET (comp->data, comp->stride * 1), comp->width);
490   lowpass_s16 (tmp + 6 * comp->width,
491       OFFSET (comp->data, comp->stride * 2), comp->width);
492   lowpass_s16 (tmp + 7 * comp->width,
493       OFFSET (comp->data, comp->stride * 3), comp->width);
494   for (i = 0; i < comp->height; i++) {
495     lowpass_s16 (tmp + 8 * comp->width,
496         OFFSET (comp->data, comp->stride * CLAMP (i + 4, 0, comp->height - 1)),
497         comp->width);
498     lowpass_vert_s16 (OFFSET (comp->data, comp->stride * i), tmp, comp->width);
499     memmove (tmp, tmp + comp->width * 1, comp->width * 8 * sizeof (int16_t));
500   }
501 
502   schro_free (tmp);
503 }
504 #endif
505 
506 #ifdef unused
507 void
schro_frame_filter_lowpass_16(SchroFrame * frame)508 schro_frame_filter_lowpass_16 (SchroFrame * frame)
509 {
510   schro_frame_component_filter_lowpass_16 (&frame->components[0]);
511   schro_frame_component_filter_lowpass_16 (&frame->components[1]);
512   schro_frame_component_filter_lowpass_16 (&frame->components[2]);
513 }
514 #endif
515 
516 static void
schro_convert_f64_u8(double * dest,uint8_t * src,int n)517 schro_convert_f64_u8 (double *dest, uint8_t * src, int n)
518 {
519   int i;
520   for (i = 0; i < n; i++) {
521     dest[i] = src[i];
522   }
523 }
524 
525 static void
schro_iir3_s16_f64(int16_t * d,int16_t * s,double * i_3,double * s2_4,int n)526 schro_iir3_s16_f64 (int16_t * d, int16_t * s, double *i_3, double *s2_4, int n)
527 {
528   int i;
529 
530   for (i = 0; i < n; i++) {
531     double x;
532 
533     x = s2_4[0] * s[i] + s2_4[1] * i_3[0] + s2_4[2] * i_3[1] + s2_4[3] * i_3[2];
534     i_3[2] = i_3[1];
535     i_3[1] = i_3[0];
536     i_3[0] = x;
537     d[i] = rint (x);
538   }
539 }
540 
541 static void
schro_iir3_rev_s16_f64(int16_t * d,int16_t * s,double * i_3,double * s2_4,int n)542 schro_iir3_rev_s16_f64 (int16_t * d, int16_t * s, double *i_3, double *s2_4,
543     int n)
544 {
545   int i;
546 
547   for (i = n - 1; i >= 0; i--) {
548     double x;
549 
550     x = s2_4[0] * s[i] + s2_4[1] * i_3[0] + s2_4[2] * i_3[1] + s2_4[3] * i_3[2];
551     i_3[2] = i_3[1];
552     i_3[1] = i_3[0];
553     i_3[0] = x;
554     d[i] = rint (x);
555   }
556 }
557 
558 static void
schro_iir3_across_u8_f64(uint8_t * d,uint8_t * s,double * i1,double * i2,double * i3,double * s2_4,int n)559 schro_iir3_across_u8_f64 (uint8_t * d, uint8_t * s, double *i1, double *i2,
560     double *i3, double *s2_4, int n)
561 {
562   int i;
563 
564   for (i = 0; i < n; i++) {
565     double x;
566 
567     x = s2_4[0] * s[i] + s2_4[1] * i1[i] + s2_4[2] * i2[i] + s2_4[3] * i3[i];
568     i3[i] = i2[i];
569     i2[i] = i1[i];
570     i1[i] = x;
571     d[i] = rint (x);
572   }
573 }
574 
575 static void
schro_iir3_across_s16_f64(int16_t * d,int16_t * s,double * i1,double * i2,double * i3,double * s2_4,int n)576 schro_iir3_across_s16_f64 (int16_t * d, int16_t * s, double *i1, double *i2,
577     double *i3, double *s2_4, int n)
578 {
579   int i;
580 
581   for (i = 0; i < n; i++) {
582     double x;
583 
584     x = s2_4[0] * s[i] + s2_4[1] * i1[i] + s2_4[2] * i2[i] + s2_4[3] * i3[i];
585     i3[i] = i2[i];
586     i2[i] = i1[i];
587     i1[i] = x;
588     d[i] = rint (x);
589   }
590 }
591 
592 static void
schro_convert_f64_s16(double * dest,int16_t * src,int n)593 schro_convert_f64_s16 (double *dest, int16_t * src, int n)
594 {
595   int i;
596   for (i = 0; i < n; i++) {
597     dest[i] = src[i];
598   }
599 }
600 
601 static void
schro_iir3_u8_f64(uint8_t * d,uint8_t * s,double * i_3,double * s2_4,int n)602 schro_iir3_u8_f64 (uint8_t * d, uint8_t * s, double *i_3, double *s2_4, int n)
603 {
604   int i;
605 
606   for (i = 0; i < n; i++) {
607     double x;
608 
609     x = s2_4[0] * s[i] + s2_4[1] * i_3[0] + s2_4[2] * i_3[1] + s2_4[3] * i_3[2];
610     i_3[2] = i_3[1];
611     i_3[1] = i_3[0];
612     i_3[0] = x;
613     d[i] = rint (x);
614   }
615 }
616 
617 static void
schro_iir3_rev_u8_f64(uint8_t * d,uint8_t * s,double * i_3,double * s2_4,int n)618 schro_iir3_rev_u8_f64 (uint8_t * d, uint8_t * s, double *i_3, double *s2_4,
619     int n)
620 {
621   int i;
622 
623   for (i = n - 1; i >= 0; i--) {
624     double x;
625 
626     x = s2_4[0] * s[i] + s2_4[1] * i_3[0] + s2_4[2] * i_3[1] + s2_4[3] * i_3[2];
627     i_3[2] = i_3[1];
628     i_3[1] = i_3[0];
629     i_3[0] = x;
630     d[i] = rint (x);
631   }
632 }
633 
634 static void
lowpass2_u8(uint8_t * d,uint8_t * s,double * coeff,int n)635 lowpass2_u8 (uint8_t * d, uint8_t * s, double *coeff, int n)
636 {
637   double state[3];
638 
639   state[0] = s[0];
640   state[1] = s[0];
641   state[2] = s[0];
642   schro_iir3_u8_f64 (d, s, state, coeff, n);
643 
644   state[0] = d[n - 1];
645   state[1] = d[n - 1];
646   state[2] = d[n - 1];
647   schro_iir3_rev_u8_f64 (d, s, state, coeff, n);
648 }
649 
650 static void
lowpass2_s16(int16_t * d,int16_t * s,double * coeff,int n)651 lowpass2_s16 (int16_t * d, int16_t * s, double *coeff, int n)
652 {
653   double state[3];
654 
655   state[0] = s[0];
656   state[1] = s[0];
657   state[2] = s[0];
658   schro_iir3_s16_f64 (d, s, state, coeff, n);
659 
660   state[0] = d[n - 1];
661   state[1] = d[n - 1];
662   state[2] = d[n - 1];
663   schro_iir3_rev_s16_f64 (d, s, state, coeff, n);
664 }
665 
666 static void
generate_coeff(double * coeff,double sigma)667 generate_coeff (double *coeff, double sigma)
668 {
669   double q;
670   double b0, b0inv, b1, b2, b3, B;
671 
672   if (sigma >= 2.5) {
673     q = 0.98711 * sigma - 0.96330;
674   } else {
675     q = 3.97156 - 4.41554 * sqrt (1 - 0.26891 * sigma);
676   }
677 
678   b0 = 1.57825 + 2.44413 * q + 1.4281 * q * q + 0.422205 * q * q * q;
679   b0inv = 1.0 / b0;
680   b1 = 2.44413 * q + 2.85619 * q * q + 1.26661 * q * q * q;
681   b2 = -1.4281 * q * q - 1.26661 * q * q * q;
682   b3 = 0.422205 * q * q * q;
683   B = 1 - (b1 + b2 + b3) / b0;
684 
685   coeff[0] = B;
686   coeff[1] = b1 * b0inv;
687   coeff[2] = b2 * b0inv;
688   coeff[3] = b3 * b0inv;
689 }
690 
691 static void
schro_frame_component_filter_lowpass2_u8(SchroFrameData * comp,double h_sigma,double v_sigma)692 schro_frame_component_filter_lowpass2_u8 (SchroFrameData * comp,
693     double h_sigma, double v_sigma)
694 {
695   int i;
696   double h_coeff[4];
697   double v_coeff[4];
698   double *i1, *i2, *i3;
699 
700   generate_coeff (h_coeff, h_sigma);
701   generate_coeff (v_coeff, v_sigma);
702 
703   i1 = schro_malloc (sizeof (double) * comp->width);
704   i2 = schro_malloc (sizeof (double) * comp->width);
705   i3 = schro_malloc (sizeof (double) * comp->width);
706 
707   for (i = 0; i < comp->height; i++) {
708     lowpass2_u8 (OFFSET (comp->data, comp->stride * i),
709         OFFSET (comp->data, comp->stride * i), h_coeff, comp->width);
710   }
711 
712   schro_convert_f64_u8 (i1, OFFSET (comp->data, comp->stride * 0), comp->width);
713   memcpy (i2, i1, sizeof (double) * comp->width);
714   memcpy (i3, i1, sizeof (double) * comp->width);
715   for (i = 0; i < comp->height; i++) {
716     schro_iir3_across_u8_f64 (OFFSET (comp->data, comp->stride * i),
717         OFFSET (comp->data, comp->stride * i),
718         i1, i2, i3, v_coeff, comp->width);
719   }
720 
721   schro_convert_f64_u8 (i1, OFFSET (comp->data,
722           comp->stride * (comp->height - 1)), comp->width);
723   memcpy (i2, i1, sizeof (double) * comp->width);
724   memcpy (i3, i1, sizeof (double) * comp->width);
725   for (i = comp->height - 1; i >= 0; i--) {
726     schro_iir3_across_u8_f64 (OFFSET (comp->data, comp->stride * i),
727         OFFSET (comp->data, comp->stride * i),
728         i1, i2, i3, v_coeff, comp->width);
729   }
730 
731   schro_free (i1);
732   schro_free (i2);
733   schro_free (i3);
734 }
735 
736 static void
schro_frame_component_filter_lowpass2_s16(SchroFrameData * comp,double h_sigma,double v_sigma)737 schro_frame_component_filter_lowpass2_s16 (SchroFrameData * comp,
738     double h_sigma, double v_sigma)
739 {
740   int i;
741   double h_coeff[4];
742   double v_coeff[4];
743   double *i1, *i2, *i3;
744 
745   generate_coeff (h_coeff, h_sigma);
746   generate_coeff (v_coeff, v_sigma);
747 
748   i1 = schro_malloc (sizeof (double) * comp->width);
749   i2 = schro_malloc (sizeof (double) * comp->width);
750   i3 = schro_malloc (sizeof (double) * comp->width);
751 
752   for (i = 0; i < comp->height; i++) {
753     lowpass2_s16 (OFFSET (comp->data, comp->stride * i),
754         OFFSET (comp->data, comp->stride * i), h_coeff, comp->width);
755   }
756 
757   schro_convert_f64_s16 (i1, OFFSET (comp->data, comp->stride * 0),
758       comp->width);
759   memcpy (i2, i1, sizeof (double) * comp->width);
760   memcpy (i3, i1, sizeof (double) * comp->width);
761   for (i = 0; i < comp->height; i++) {
762     schro_iir3_across_s16_f64 (OFFSET (comp->data, comp->stride * i),
763         OFFSET (comp->data, comp->stride * i),
764         i1, i2, i3, v_coeff, comp->width);
765   }
766 
767   schro_convert_f64_s16 (i1, OFFSET (comp->data,
768           comp->stride * (comp->height - 1)), comp->width);
769   memcpy (i2, i1, sizeof (double) * comp->width);
770   memcpy (i3, i1, sizeof (double) * comp->width);
771   for (i = comp->height - 1; i >= 0; i--) {
772     schro_iir3_across_s16_f64 (OFFSET (comp->data, comp->stride * i),
773         OFFSET (comp->data, comp->stride * i),
774         i1, i2, i3, v_coeff, comp->width);
775   }
776 
777 
778 
779   schro_free (i1);
780   schro_free (i2);
781   schro_free (i3);
782 }
783 
784 void
schro_frame_filter_lowpass2(SchroFrame * frame,double sigma)785 schro_frame_filter_lowpass2 (SchroFrame * frame, double sigma)
786 {
787   double chroma_sigma_h;
788   double chroma_sigma_v;
789 
790   chroma_sigma_h = sigma / (1 << SCHRO_FRAME_FORMAT_H_SHIFT (frame->format));
791   chroma_sigma_v = sigma / (1 << SCHRO_FRAME_FORMAT_V_SHIFT (frame->format));
792 
793   switch (SCHRO_FRAME_FORMAT_DEPTH (frame->format)) {
794     case SCHRO_FRAME_FORMAT_DEPTH_U8:
795       schro_frame_component_filter_lowpass2_u8 (&frame->components[0], sigma,
796           sigma);
797       schro_frame_component_filter_lowpass2_u8 (&frame->components[1],
798           chroma_sigma_h, chroma_sigma_v);
799       schro_frame_component_filter_lowpass2_u8 (&frame->components[2],
800           chroma_sigma_h, chroma_sigma_v);
801       break;
802     case SCHRO_FRAME_FORMAT_DEPTH_S16:
803       schro_frame_component_filter_lowpass2_s16 (&frame->components[0], sigma,
804           sigma);
805       schro_frame_component_filter_lowpass2_s16 (&frame->components[1],
806           chroma_sigma_h, chroma_sigma_v);
807       schro_frame_component_filter_lowpass2_s16 (&frame->components[2],
808           chroma_sigma_h, chroma_sigma_v);
809       break;
810     default:
811       SCHRO_ASSERT (0);
812       break;
813   }
814 }
815 
816 
817 #ifdef unused
818 void
schro_frame_filter_wavelet(SchroFrame * frame)819 schro_frame_filter_wavelet (SchroFrame * frame)
820 {
821   SchroFrame *tmpframe;
822   SchroFrameData *comp;
823   SchroHistogram hist;
824   int component;
825   int16_t *tmp;
826   SchroParams params;
827   int i;
828 
829   tmp = schro_malloc (2 * frame->width * sizeof (int16_t));
830 
831   tmpframe = schro_frame_new_and_alloc (NULL,
832       SCHRO_FRAME_FORMAT_S16_444 | frame->format,
833       ROUND_UP_POW2 (frame->width, 5), ROUND_UP_POW2 (frame->height, 5));
834   schro_frame_convert (tmpframe, frame);
835 
836   params.transform_depth = 1;
837   params.iwt_luma_width = frame->width;
838   params.iwt_luma_height = frame->height;
839   params.iwt_chroma_width = frame->components[1].width;
840   params.iwt_chroma_height = frame->components[1].height;
841 
842   for (component = 0; component < 3; component++) {
843     comp = &tmpframe->components[component];
844 
845     schro_wavelet_transform_2d (comp, SCHRO_WAVELET_LE_GALL_5_3, tmp);
846 
847     for (i = 1; i < 4; i++) {
848       SchroFrameData fd;
849       int y;
850       int cutoff;
851 
852       schro_subband_get_frame_data (&fd, tmpframe, component, i, &params);
853       schro_histogram_init (&hist);
854 
855       for (y = 0; y < fd.height; y++) {
856         schro_histogram_add_array_s16 (&hist, OFFSET (fd.data, y * fd.stride),
857             fd.width);
858       }
859 
860       cutoff = 100;
861       for (y = 0; y < fd.height; y++) {
862         int16_t *line = OFFSET (fd.data, fd.stride * y);
863         int x;
864         for (x = 0; x < fd.width; x++) {
865           if (line[x] > -cutoff && line[x] < cutoff)
866             line[x] = 0;
867         }
868       }
869     }
870 
871     schro_wavelet_inverse_transform_2d (comp, SCHRO_WAVELET_LE_GALL_5_3, tmp);
872   }
873 
874   schro_frame_convert (frame, tmpframe);
875   schro_frame_unref (tmpframe);
876 }
877 #endif
878 
879 
880 static double
random_std(void)881 random_std (void)
882 {
883   double x;
884   double y;
885 
886   while (1) {
887     x = -5.0 + rand () * (1.0 / RAND_MAX) * 10;
888     y = rand () * (1.0 / RAND_MAX);
889 
890     if (y < exp (-x * x * 0.5))
891       return x;
892   }
893 }
894 
895 static void
addnoise_u8(uint8_t * dest,int n,double sigma)896 addnoise_u8 (uint8_t * dest, int n, double sigma)
897 {
898   int i;
899   int x;
900 
901   for (i = 0; i < n; i++) {
902     x = rint (random_std () * sigma) + dest[i];
903     dest[i] = CLAMP (x, 0, 255);
904   }
905 }
906 
907 static void
schro_frame_component_filter_addnoise(SchroFrameData * comp,double sigma)908 schro_frame_component_filter_addnoise (SchroFrameData * comp, double sigma)
909 {
910   int i;
911 
912   for (i = 0; i < comp->height; i++) {
913     addnoise_u8 (OFFSET (comp->data, comp->stride * i), comp->width, sigma);
914   }
915 }
916 
917 void
schro_frame_filter_addnoise(SchroFrame * frame,double sigma)918 schro_frame_filter_addnoise (SchroFrame * frame, double sigma)
919 {
920   schro_frame_component_filter_addnoise (&frame->components[0], sigma);
921   schro_frame_component_filter_addnoise (&frame->components[1], sigma);
922   schro_frame_component_filter_addnoise (&frame->components[2], sigma);
923 }
924 
925 
926 static int
ilogx_size(int i)927 ilogx_size (int i)
928 {
929   if (i < (1 << SCHRO_HISTOGRAM_SHIFT))
930     return 1;
931   return 1 << ((i >> SCHRO_HISTOGRAM_SHIFT) - 1);
932 }
933 
934 static int
iexpx(int x)935 iexpx (int x)
936 {
937   if (x < (1 << SCHRO_HISTOGRAM_SHIFT))
938     return x;
939 
940   return ((1 << SCHRO_HISTOGRAM_SHIFT) | (x & ((1 << SCHRO_HISTOGRAM_SHIFT) -
941               1))) << ((x >> SCHRO_HISTOGRAM_SHIFT) - 1);
942 }
943 
944 
945 void
schro_frame_filter_adaptive_lowpass(SchroFrame * frame)946 schro_frame_filter_adaptive_lowpass (SchroFrame * frame)
947 {
948   SchroHistogram hist;
949   double slope;
950   SchroFrame *tmp;
951   int16_t tmpdata[2048];
952   double sigma;
953   int j;
954   int i;
955 
956   tmp = schro_frame_new_and_alloc (NULL,
957       SCHRO_FRAME_FORMAT_S16_444 | frame->format, frame->width, frame->height);
958   schro_frame_convert (tmp, frame);
959 
960   schro_wavelet_transform_2d (&tmp->components[0], SCHRO_WAVELET_LE_GALL_5_3,
961       tmpdata);
962 
963   schro_histogram_init (&hist);
964   for (j = 0; j < tmp->height / 2; j++) {
965     schro_histogram_add_array_s16 (&hist,
966         OFFSET (tmp->components[0].data,
967             tmp->components[0].stride * (2 * j + 1)), tmp->width / 2);
968   }
969 
970   schro_frame_unref (tmp);
971   tmp = NULL;
972 
973   slope = schro_histogram_estimate_slope (&hist);
974 
975   for (i = 0; i < SCHRO_HISTOGRAM_SIZE; i++) {
976     schro_dump (SCHRO_DUMP_HIST_TEST, "%d %g\n",
977         iexpx (i), hist.bins[i] / ilogx_size (i));
978   }
979 
980   /* good for 2 Mb DVD intra-only rip */
981   sigma = -1.0 / slope;
982   if (sigma > 1.0) {
983     SCHRO_DEBUG ("enabling filtering (slope %g)", slope);
984 
985     schro_frame_filter_lowpass2 (frame, sigma);
986   }
987 }
988