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, ¶ms);
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