schroedinger/schrofilter.c

   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
  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
  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
  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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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 }