ffmpeg/libavcodec/dwt.c

   1 /*
   2  * Copyright (C) 2004-2010 Michael Niedermayer <michaelni@gmx.at>
   3  * Copyright (C) 2008 David Conrad
   4  *
   5  * This file is part of FFmpeg.
   6  *
   7  * FFmpeg is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU Lesser General Public
   9  * License as published by the Free Software Foundation; either
  10  * version 2.1 of the License, or (at your option) any later version.
  11  *
  12  * FFmpeg is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public
  18  * License along with FFmpeg; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  20  */
  21
  22 #include "libavutil/attributes.h"
  23 #include "dsputil.h"
  24 #include "dwt.h"
  25 #include "libavcodec/x86/dwt.h"
  26
  27 void ff_slice_buffer_init(slice_buffer * buf, int line_count, int max_allocated_lines, int line_width, IDWTELEM * base_buffer)
  28 {
  29     int i;
  30
  31     buf->base_buffer = base_buffer;
  32     buf->line_count  = line_count;
  33     buf->line_width  = line_width;
  34     buf->data_count  = max_allocated_lines;
  35     buf->line        = av_mallocz (sizeof(IDWTELEM *) * line_count);
  36     buf->data_stack  = av_malloc (sizeof(IDWTELEM *) * max_allocated_lines);
  37
  38     for(i = 0; i < max_allocated_lines; i++){
  39         buf->data_stack[i] = av_malloc (sizeof(IDWTELEM) * line_width);
  40     }
  41
  42     buf->data_stack_top = max_allocated_lines - 1;
  43 }
  44
  45 IDWTELEM * ff_slice_buffer_load_line(slice_buffer * buf, int line)
  46 {
  47     IDWTELEM * buffer;
  48
  49     assert(buf->data_stack_top >= 0);
  50     if (buf->line[line])
  51         return buf->line[line];
  52
  53     buffer = buf->data_stack[buf->data_stack_top];
  54     buf->data_stack_top--;
  55     buf->line[line] = buffer;
  56
  57     return buffer;
  58 }
  59
  60 void ff_slice_buffer_release(slice_buffer * buf, int line)
  61 {
  62     IDWTELEM * buffer;
  63
  64     assert(line >= 0 && line < buf->line_count);
  65     assert(buf->line[line]);
  66
  67     buffer = buf->line[line];
  68     buf->data_stack_top++;
  69     buf->data_stack[buf->data_stack_top] = buffer;
  70     buf->line[line] = NULL;
  71 }
  72
  73 void ff_slice_buffer_flush(slice_buffer * buf)
  74 {
  75     int i;
  76     for(i = 0; i < buf->line_count; i++){
  77         if (buf->line[i])
  78             ff_slice_buffer_release(buf, i);
  79     }
  80 }
  81
  82 void ff_slice_buffer_destroy(slice_buffer * buf)
  83 {
  84     int i;
  85     ff_slice_buffer_flush(buf);
  86
  87     for(i = buf->data_count - 1; i >= 0; i--){
  88         av_freep(&buf->data_stack[i]);
  89     }
  90     av_freep(&buf->data_stack);
  91     av_freep(&buf->line);
  92 }
  93
  94 static inline int mirror(int v, int m){
  95     while((unsigned)v > (unsigned)m){
  96         v = -v;
  97         if(v < 0) v+= 2*m;
  98     }
  99     return v;
 100 }
 101
 102 static av_always_inline void
 103 lift(DWTELEM *dst, DWTELEM *src, DWTELEM *ref,
 104      int dst_step, int src_step, int ref_step,
 105      int width, int mul, int add, int shift,
 106      int highpass, int inverse){
 107     const int mirror_left  = !highpass;
 108     const int mirror_right = (width & 1) ^ highpass;
 109     const int w            = (width >> 1) - 1 + (highpass & width);
 110     int i;
 111
 112 #define LIFT(src, ref, inv) ((src) + ((inv) ? - (ref) : + (ref)))
 113     if(mirror_left){
 114         dst[0] = LIFT(src[0], ((mul * 2 * ref[0] + add) >> shift), inverse);
 115         dst   += dst_step;
 116         src   += src_step;
 117     }
 118
 119     for(i=0; i<w; i++){
 120         dst[i * dst_step] =
 121             LIFT(src[i * src_step],
 122                  ((mul * (ref[i * ref_step] + ref[(i + 1) * ref_step]) + add) >> shift),
 123                  inverse);
 124     }
 125
 126     if(mirror_right){
 127         dst[w * dst_step] =
 128             LIFT(src[w * src_step],
 129                  ((mul * 2 * ref[w * ref_step] + add) >> shift),
 130                  inverse);
 131     }
 132 }
 133
 134 static av_always_inline void
 135 inv_lift(IDWTELEM *dst, IDWTELEM *src, IDWTELEM *ref,
 136          int dst_step, int src_step, int ref_step,
 137          int width, int mul, int add, int shift,
 138          int highpass, int inverse){
 139     const int mirror_left  = !highpass;
 140     const int mirror_right = (width&1) ^ highpass;
 141     const int w            = (width >> 1) - 1 + (highpass & width);
 142     int i;
 143
 144 #define LIFT(src, ref, inv) ((src) + ((inv) ? - (ref) : + (ref)))
 145     if(mirror_left){
 146         dst[0] =
 147             LIFT(src[0],
 148                  ((mul * 2 * ref[0] + add) >> shift),
 149                  inverse);
 150         dst   += dst_step;
 151         src   += src_step;
 152     }
 153
 154     for(i = 0; i < w; i++){
 155         dst[i * dst_step] =
 156             LIFT(src[i * src_step],
 157                  ((mul * (ref[i * ref_step] + ref[(i + 1) * ref_step]) + add) >> shift),
 158                  inverse);
 159     }
 160
 161     if(mirror_right){
 162         dst[w * dst_step] =
 163             LIFT(src[w * src_step],
 164                  ((mul * 2 * ref[w * ref_step] + add) >> shift),
 165                  inverse);
 166     }
 167 }
 168
 169 #ifndef liftS
 170 static av_always_inline void
 171 liftS(DWTELEM *dst, DWTELEM *src, DWTELEM *ref,
 172       int dst_step, int src_step, int ref_step,
 173       int width, int mul, int add, int shift,
 174       int highpass, int inverse){
 175     const int mirror_left  = !highpass;
 176     const int mirror_right = (width&1) ^ highpass;
 177     const int w            = (width >> 1) - 1 + (highpass & width);
 178     int i;
 179
 180     assert(shift == 4);
 181 #define LIFTS(src, ref, inv)                                            \
 182     ((inv) ?                                                            \
 183      (src) + (((ref) + 4 * (src)) >> shift):                            \
 184      -((-16 * (src) + (ref) + add / 4 + 1 + (5 << 25)) / (5 * 4) - (1 << 23)))
 185     if(mirror_left){
 186         dst[0] = LIFTS(src[0], mul * 2 * ref[0] + add, inverse);
 187         dst   += dst_step;
 188         src   += src_step;
 189     }
 190
 191     for(i = 0; i < w; i++){
 192         dst[i * dst_step] =
 193             LIFTS(src[i * src_step],
 194                   mul * (ref[i * ref_step] + ref[(i+1) * ref_step]) + add,
 195                   inverse);
 196     }
 197
 198     if(mirror_right){
 199         dst[w * dst_step] =
 200             LIFTS(src[w * src_step], mul * 2 * ref[w * ref_step] + add, inverse);
 201     }
 202 }
 203 static av_always_inline void
 204 inv_liftS(IDWTELEM *dst, IDWTELEM *src, IDWTELEM *ref,
 205           int dst_step, int src_step, int ref_step,
 206           int width, int mul, int add, int shift,
 207           int highpass, int inverse){
 208     const int mirror_left  = !highpass;
 209     const int mirror_right = (width&1) ^ highpass;
 210     const int w            = (width >> 1) - 1 + (highpass & width);
 211     int i;
 212
 213     assert(shift == 4);
 214 #define LIFTS(src, ref, inv)                                            \
 215     ((inv) ?                                                            \
 216      (src) + (((ref) + 4 * (src)) >> shift):                            \
 217      -((-16 * (src) + (ref) + add / 4 + 1 + (5 << 25)) / (5 * 4) - (1 << 23)))
 218     if(mirror_left){
 219         dst[0] = LIFTS(src[0], mul * 2 * ref[0] + add, inverse);
 220         dst += dst_step;
 221         src += src_step;
 222     }
 223
 224     for(i = 0; i < w; i++){
 225         dst[i * dst_step] =
 226             LIFTS(src[i * src_step],
 227                   mul * (ref[i * ref_step] + ref[(i+1) * ref_step]) + add,
 228                   inverse);
 229     }
 230
 231     if(mirror_right){
 232         dst[w * dst_step] =
 233             LIFTS(src[w * src_step], mul * 2 * ref[w * ref_step] + add, inverse);
 234     }
 235 }
 236 #endif /* ! liftS */
 237
 238 static void horizontal_decompose53i(DWTELEM *b, int width){
 239     DWTELEM temp[width];
 240     const int width2 = width>>1;
 241     const int w2     = (width+1)>>1;
 242     int x;
 243
 244     for(x = 0; x < width2; x++){
 245         temp[x   ] = b[2 * x    ];
 246         temp[x+w2] = b[2 * x + 1];
 247     }
 248     if(width & 1)
 249         temp[x   ] = b[2 * x    ];
 250 #if 0
 251     {
 252     int A1,A2,A3,A4;
 253     A2  = temp[1       ];
 254     A4  = temp[0       ];
 255     A1  = temp[0+width2];
 256     A1 -= (A2 + A4) >> 1;
 257     A4 += (A1 + 1) >> 1;
 258     b[0 + width2] = A1;
 259     b[0       ]   = A4;
 260     for(x = 1; x + 1 < width2; x += 2){
 261         A3  = temp[x + width2];
 262         A4  = temp[x + 1     ];
 263         A3 -= (A2 + A4) >> 1;
 264         A2 += (A1 + A3 + 2) >> 2;
 265         b[x + width2] = A3;
 266         b[x         ] = A2;
 267
 268         A1  = temp[x + 1 + width2];
 269         A2  = temp[x + 2         ];
 270         A1 -= (A2 + A4) >> 1;
 271         A4 += (A1 + A3 + 2) >> 2;
 272         b[x + 1 + width2] = A1;
 273         b[x + 1         ] = A4;
 274     }
 275     A3  = temp[width - 1];
 276     A3 -= A2;
 277     A2 += (A1 + A3 + 2) >> 2;
 278     b[width -1] = A3;
 279     b[width2-1] = A2;
 280     }
 281 #else
 282     lift(b + w2, temp + w2, temp, 1, 1, 1, width, -1, 0, 1, 1, 0);
 283     lift(b   , temp   , b + w2, 1, 1, 1, width,  1, 2, 2, 0, 0);
 284 #endif /* 0 */
 285 }
 286
 287 static void vertical_decompose53iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
 288     int i;
 289
 290     for(i = 0; i < width; i++){
 291         b1[i] -= (b0[i] + b2[i]) >> 1;
 292     }
 293 }
 294
 295 static void vertical_decompose53iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
 296     int i;
 297
 298     for(i = 0; i < width; i++){
 299         b1[i] += (b0[i] + b2[i] + 2) >> 2;
 300     }
 301 }
 302
 303 static void spatial_decompose53i(DWTELEM *buffer, int width, int height, int stride){
 304     int y;
 305     DWTELEM *b0 = buffer + mirror(-2 - 1, height-1)*stride;
 306     DWTELEM *b1 = buffer + mirror(-2    , height-1)*stride;
 307
 308     for(y = -2; y < height; y += 2){
 309         DWTELEM *b2 = buffer + mirror(y + 1, height-1) * stride;
 310         DWTELEM *b3 = buffer + mirror(y + 2, height-1) * stride;
 311
 312         if(y + 1 < (unsigned)height) horizontal_decompose53i(b2, width);
 313         if(y + 2 < (unsigned)height) horizontal_decompose53i(b3, width);
 314
 315         if(y + 1 < (unsigned)height) vertical_decompose53iH0(b1, b2, b3, width);
 316         if(y + 0 < (unsigned)height) vertical_decompose53iL0(b0, b1, b2, width);
 317
 318         b0 = b2;
 319         b1 = b3;
 320     }
 321 }
 322
 323 static void horizontal_decompose97i(DWTELEM *b, int width){
 324     DWTELEM temp[width];
 325     const int w2 = (width+1)>>1;
 326
 327     lift (temp + w2, b    + 1 , b        , 1, 2, 2, width,  W_AM, W_AO, W_AS, 1, 1);
 328     liftS(temp     , b        , temp + w2, 1, 2, 1, width,  W_BM, W_BO, W_BS, 0, 0);
 329     lift (b    + w2, temp + w2, temp     , 1, 1, 1, width,  W_CM, W_CO, W_CS, 1, 0);
 330     lift (b        , temp     , b    + w2, 1, 1, 1, width,  W_DM, W_DO, W_DS, 0, 0);
 331 }
 332
 333
 334 static void vertical_decompose97iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
 335     int i;
 336
 337     for(i = 0; i < width; i++){
 338         b1[i] -= (W_AM * (b0[i] + b2[i]) + W_AO) >> W_AS;
 339     }
 340 }
 341
 342 static void vertical_decompose97iH1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
 343     int i;
 344
 345     for(i=0; i < width; i++){
 346         b1[i] += (W_CM * (b0[i] + b2[i]) + W_CO) >> W_CS;
 347     }
 348 }
 349
 350 static void vertical_decompose97iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
 351     int i;
 352
 353     for(i = 0; i < width; i++){
 354 #ifdef liftS
 355         b1[i] -= (W_BM * (b0[i] + b2[i]) + W_BO) >> W_BS;
 356 #else
 357         b1[i] = (16 * 4 * b1[i] - 4 * (b0[i] + b2[i]) + W_BO * 5 + (5 << 27)) / (5 * 16) - (1 << 23);
 358 #endif
 359     }
 360 }
 361
 362 static void vertical_decompose97iL1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2, int width){
 363     int i;
 364
 365     for(i = 0; i < width; i++){
 366         b1[i] += (W_DM * (b0[i] + b2[i]) + W_DO) >> W_DS;
 367     }
 368 }
 369
 370 static void spatial_decompose97i(DWTELEM *buffer, int width, int height, int stride){
 371     int y;
 372     DWTELEM *b0 = buffer + mirror(-4 - 1, height-1) * stride;
 373     DWTELEM *b1 = buffer + mirror(-4    , height-1) * stride;
 374     DWTELEM *b2 = buffer + mirror(-4 + 1, height-1) * stride;
 375     DWTELEM *b3 = buffer + mirror(-4 + 2, height-1) * stride;
 376
 377     for(y = -4; y < height; y += 2){
 378         DWTELEM *b4 = buffer + mirror(y + 3, height-1) * stride;
 379         DWTELEM *b5 = buffer + mirror(y + 4, height-1) * stride;
 380
 381         if(y + 3 < (unsigned)height) horizontal_decompose97i(b4, width);
 382         if(y + 4 < (unsigned)height) horizontal_decompose97i(b5, width);
 383
 384         if(y + 3 < (unsigned)height) vertical_decompose97iH0(b3, b4, b5, width);
 385         if(y + 2 < (unsigned)height) vertical_decompose97iL0(b2, b3, b4, width);
 386         if(y + 1 < (unsigned)height) vertical_decompose97iH1(b1, b2, b3, width);
 387         if(y + 0 < (unsigned)height) vertical_decompose97iL1(b0, b1, b2, width);
 388
 389         b0 = b2;
 390         b1 = b3;
 391         b2 = b4;
 392         b3 = b5;
 393     }
 394 }
 395
 396 void ff_spatial_dwt(DWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){
 397     int level;
 398
 399     for(level = 0; level < decomposition_count; level++){
 400         switch(type){
 401         case DWT_97: spatial_decompose97i(buffer, width >> level, height >> level, stride << level); break;
 402         case DWT_53: spatial_decompose53i(buffer, width >> level, height >> level, stride << level); break;
 403         }
 404     }
 405 }
 406
 407 static void horizontal_compose53i(IDWTELEM *b, int width){
 408     IDWTELEM temp[width];
 409     const int width2 = width       >> 1;
 410     const int w2     = (width + 1) >> 1;
 411     int x;
 412
 413     for(x = 0; x < width2; x++){
 414         temp[2 * x    ] = b[x     ];
 415         temp[2 * x + 1] = b[x + w2];
 416     }
 417     if(width & 1)
 418         temp[2 * x    ] = b[x   ];
 419
 420     b[0] = temp[0] - ((temp[1] + 1) >> 1);
 421     for(x = 2; x < width - 1; x += 2){
 422         b[x    ] = temp[x    ] - ((temp[x - 1] + temp[x + 1] + 2) >> 2);
 423         b[x - 1] = temp[x - 1] + ((b   [x - 2] + b   [x  ] + 1) >> 1);
 424     }
 425     if(width & 1){
 426         b[x    ] = temp[x    ] - ((temp[x - 1] + 1) >> 1);
 427         b[x - 1] = temp[x - 1] + ((b   [x - 2] + b  [x  ] + 1) >> 1);
 428     }else
 429         b[x - 1] = temp[x - 1] + b[x - 2];
 430 }
 431
 432 static void vertical_compose53iH0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
 433     int i;
 434
 435     for(i = 0; i < width; i++){
 436         b1[i] += (b0[i] + b2[i]) >> 1;
 437     }
 438 }
 439
 440 static void vertical_compose53iL0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
 441     int i;
 442
 443     for(i = 0; i < width; i++){
 444         b1[i] -= (b0[i] + b2[i] + 2) >> 2;
 445     }
 446 }
 447
 448 static void spatial_compose53i_buffered_init(DWTCompose *cs, slice_buffer * sb, int height, int stride_line){
 449     cs->b0 = slice_buffer_get_line(sb, mirror(-1-1, height - 1) * stride_line);
 450     cs->b1 = slice_buffer_get_line(sb, mirror(-1  , height - 1) * stride_line);
 451     cs->y  = -1;
 452 }
 453
 454 static void spatial_compose53i_init(DWTCompose *cs, IDWTELEM *buffer, int height, int stride){
 455     cs->b0 = buffer + mirror(-1-1, height - 1) * stride;
 456     cs->b1 = buffer + mirror(-1  , height - 1) * stride;
 457     cs->y  = -1;
 458 }
 459
 460 static void spatial_compose53i_dy_buffered(DWTCompose *cs, slice_buffer * sb, int width, int height, int stride_line){
 461     int y        = cs->y;
 462     IDWTELEM *b0 = cs->b0;
 463     IDWTELEM *b1 = cs->b1;
 464     IDWTELEM *b2 = slice_buffer_get_line(sb, mirror(y + 1, height-1) * stride_line);
 465     IDWTELEM *b3 = slice_buffer_get_line(sb, mirror(y + 2, height-1) * stride_line);
 466
 467     if(y + 1 < (unsigned)height && y < (unsigned)height){
 468         int x;
 469
 470         for(x = 0; x < width; x++){
 471             b2[x] -= (b1[x] + b3[x] + 2) >> 2;
 472             b1[x] += (b0[x] + b2[x]) >> 1;
 473         }
 474     }else{
 475         if(y + 1 < (unsigned)height) vertical_compose53iL0(b1, b2, b3, width);
 476         if(y + 0 < (unsigned)height) vertical_compose53iH0(b0, b1, b2, width);
 477     }
 478
 479         if(y - 1 <(unsigned)height) horizontal_compose53i(b0, width);
 480         if(y + 0 <(unsigned)height) horizontal_compose53i(b1, width);
 481
 482     cs->b0 = b2;
 483     cs->b1 = b3;
 484     cs->y += 2;
 485 }
 486
 487 static void spatial_compose53i_dy(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride){
 488     int y        = cs->y;
 489     IDWTELEM *b0 = cs->b0;
 490     IDWTELEM *b1 = cs->b1;
 491     IDWTELEM *b2 = buffer + mirror(y + 1, height - 1) * stride;
 492     IDWTELEM *b3 = buffer + mirror(y + 2, height - 1) * stride;
 493
 494         if(y + 1 < (unsigned)height) vertical_compose53iL0(b1, b2, b3, width);
 495         if(y + 0 < (unsigned)height) vertical_compose53iH0(b0, b1, b2, width);
 496
 497         if(y - 1 < (unsigned)height) horizontal_compose53i(b0, width);
 498         if(y + 0 < (unsigned)height) horizontal_compose53i(b1, width);
 499
 500     cs->b0 = b2;
 501     cs->b1 = b3;
 502     cs->y += 2;
 503 }
 504
 505 static void av_unused spatial_compose53i(IDWTELEM *buffer, int width, int height, int stride){
 506     DWTCompose cs;
 507     spatial_compose53i_init(&cs, buffer, height, stride);
 508     while(cs.y <= height)
 509         spatial_compose53i_dy(&cs, buffer, width, height, stride);
 510 }
 511
 512
 513 void ff_snow_horizontal_compose97i(IDWTELEM *b, int width){
 514     IDWTELEM temp[width];
 515     const int w2 = (width + 1) >> 1;
 516
 517 #if 0 //maybe more understadable but slower
 518     inv_lift (temp    , b      , b  + w2, 2, 1, 1, width,  W_DM, W_DO, W_DS, 0, 1);
 519     inv_lift (temp + 1, b  + w2, temp   , 2, 1, 2, width,  W_CM, W_CO, W_CS, 1, 1);
 520
 521     inv_liftS(b      ,temp    ,temp + 1, 2, 2, 2, width,  W_BM, W_BO, W_BS, 0, 1);
 522     inv_lift (b + 1  ,temp + 1,b       , 2, 2, 2, width,  W_AM, W_AO, W_AS, 1, 0);
 523 #else
 524     int x;
 525     temp[0] = b[0] - ((3 * b[w2] + 2) >> 2);
 526     for(x = 1; x < (width >> 1); x++){
 527         temp[2 * x    ] = b[x         ] - (( 3 * (b   [ x + w2 - 1] + b[x + w2]) + 4) >> 3);
 528         temp[2 * x - 1] = b[x + w2 - 1] - temp[2 * x - 2] - temp[2 * x];
 529     }
 530     if(width & 1){
 531         temp[2*x  ] = b[x     ] - ((3*b   [x+w2-1]+2)>>2);
 532         temp[2*x-1] = b[x+w2-1] - temp[2*x-2] - temp[2*x];
 533     }else
 534         temp[2*x-1] = b[x+w2-1] - 2*temp[2*x-2];
 535
 536     b[0] = temp[0] + ((2*temp[0] + temp[1]+4)>>3);
 537     for(x=2; x<width-1; x+=2){
 538         b[x  ] = temp[x  ] + ((4*temp[x  ] + temp[x-1] + temp[x+1]+8)>>4);
 539         b[x-1] = temp[x-1] + ((3*(b  [x-2] + b   [x  ] ))>>1);
 540     }
 541     if(width&1){
 542         b[x  ] = temp[x  ] + ((2*temp[x  ] + temp[x-1]+4)>>3);
 543         b[x-1] = temp[x-1] + ((3*(b  [x-2] + b   [x  ] ))>>1);
 544     }else
 545         b[x-1] = temp[x-1] + 3*b [x-2];
 546 #endif
 547 }
 548
 549 static void vertical_compose97iH0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
 550     int i;
 551
 552     for(i=0; i<width; i++){
 553         b1[i] += (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
 554     }
 555 }
 556
 557 static void vertical_compose97iH1(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
 558     int i;
 559
 560     for(i=0; i<width; i++){
 561         b1[i] -= (W_CM*(b0[i] + b2[i])+W_CO)>>W_CS;
 562     }
 563 }
 564
 565 static void vertical_compose97iL0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
 566     int i;
 567
 568     for(i=0; i<width; i++){
 569 #ifdef liftS
 570         b1[i] += (W_BM*(b0[i] + b2[i])+W_BO)>>W_BS;
 571 #else
 572         b1[i] += (W_BM*(b0[i] + b2[i])+4*b1[i]+W_BO)>>W_BS;
 573 #endif
 574     }
 575 }
 576
 577 static void vertical_compose97iL1(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width){
 578     int i;
 579
 580     for(i=0; i<width; i++){
 581         b1[i] -= (W_DM*(b0[i] + b2[i])+W_DO)>>W_DS;
 582     }
 583 }
 584
 585 void ff_snow_vertical_compose97i(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, IDWTELEM *b3, IDWTELEM *b4, IDWTELEM *b5, int width){
 586     int i;
 587
 588     for(i=0; i<width; i++){
 589         b4[i] -= (W_DM*(b3[i] + b5[i])+W_DO)>>W_DS;
 590         b3[i] -= (W_CM*(b2[i] + b4[i])+W_CO)>>W_CS;
 591 #ifdef liftS
 592         b2[i] += (W_BM*(b1[i] + b3[i])+W_BO)>>W_BS;
 593 #else
 594         b2[i] += (W_BM*(b1[i] + b3[i])+4*b2[i]+W_BO)>>W_BS;
 595 #endif
 596         b1[i] += (W_AM*(b0[i] + b2[i])+W_AO)>>W_AS;
 597     }
 598 }
 599
 600 static void spatial_compose97i_buffered_init(DWTCompose *cs, slice_buffer * sb, int height, int stride_line){
 601     cs->b0 = slice_buffer_get_line(sb, mirror(-3-1, height-1) * stride_line);
 602     cs->b1 = slice_buffer_get_line(sb, mirror(-3  , height-1) * stride_line);
 603     cs->b2 = slice_buffer_get_line(sb, mirror(-3+1, height-1) * stride_line);
 604     cs->b3 = slice_buffer_get_line(sb, mirror(-3+2, height-1) * stride_line);
 605     cs->y = -3;
 606 }
 607
 608 static void spatial_compose97i_init(DWTCompose *cs, IDWTELEM *buffer, int height, int stride){
 609     cs->b0 = buffer + mirror(-3-1, height-1)*stride;
 610     cs->b1 = buffer + mirror(-3  , height-1)*stride;
 611     cs->b2 = buffer + mirror(-3+1, height-1)*stride;
 612     cs->b3 = buffer + mirror(-3+2, height-1)*stride;
 613     cs->y = -3;
 614 }
 615
 616 static void spatial_compose97i_dy_buffered(DWTContext *dsp, DWTCompose *cs, slice_buffer * sb, int width, int height, int stride_line){
 617     int y = cs->y;
 618
 619     IDWTELEM *b0= cs->b0;
 620     IDWTELEM *b1= cs->b1;
 621     IDWTELEM *b2= cs->b2;
 622     IDWTELEM *b3= cs->b3;
 623     IDWTELEM *b4= slice_buffer_get_line(sb, mirror(y + 3, height - 1) * stride_line);
 624     IDWTELEM *b5= slice_buffer_get_line(sb, mirror(y + 4, height - 1) * stride_line);
 625
 626     if(y>0 && y+4<height){
 627         dsp->vertical_compose97i(b0, b1, b2, b3, b4, b5, width);
 628     }else{
 629         if(y+3<(unsigned)height) vertical_compose97iL1(b3, b4, b5, width);
 630         if(y+2<(unsigned)height) vertical_compose97iH1(b2, b3, b4, width);
 631         if(y+1<(unsigned)height) vertical_compose97iL0(b1, b2, b3, width);
 632         if(y+0<(unsigned)height) vertical_compose97iH0(b0, b1, b2, width);
 633     }
 634
 635     if(y-1<(unsigned)height) dsp->horizontal_compose97i(b0, width);
 636     if(y+0<(unsigned)height) dsp->horizontal_compose97i(b1, width);
 637
 638     cs->b0=b2;
 639     cs->b1=b3;
 640     cs->b2=b4;
 641     cs->b3=b5;
 642     cs->y += 2;
 643 }
 644
 645 static void spatial_compose97i_dy(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride){
 646     int y = cs->y;
 647     IDWTELEM *b0= cs->b0;
 648     IDWTELEM *b1= cs->b1;
 649     IDWTELEM *b2= cs->b2;
 650     IDWTELEM *b3= cs->b3;
 651     IDWTELEM *b4= buffer + mirror(y+3, height-1)*stride;
 652     IDWTELEM *b5= buffer + mirror(y+4, height-1)*stride;
 653
 654     if(y+3<(unsigned)height) vertical_compose97iL1(b3, b4, b5, width);
 655     if(y+2<(unsigned)height) vertical_compose97iH1(b2, b3, b4, width);
 656     if(y+1<(unsigned)height) vertical_compose97iL0(b1, b2, b3, width);
 657     if(y+0<(unsigned)height) vertical_compose97iH0(b0, b1, b2, width);
 658
 659     if(y-1<(unsigned)height) ff_snow_horizontal_compose97i(b0, width);
 660     if(y+0<(unsigned)height) ff_snow_horizontal_compose97i(b1, width);
 661
 662     cs->b0=b2;
 663     cs->b1=b3;
 664     cs->b2=b4;
 665     cs->b3=b5;
 666     cs->y += 2;
 667 }
 668
 669 static void av_unused spatial_compose97i(IDWTELEM *buffer, int width, int height, int stride){
 670     DWTCompose cs;
 671     spatial_compose97i_init(&cs, buffer, height, stride);
 672     while(cs.y <= height)
 673         spatial_compose97i_dy(&cs, buffer, width, height, stride);
 674 }
 675
 676 void ff_spatial_idwt_buffered_init(DWTCompose *cs, slice_buffer * sb, int width, int height, int stride_line, int type, int decomposition_count){
 677     int level;
 678     for(level=decomposition_count-1; level>=0; level--){
 679         switch(type){
 680         case DWT_97: spatial_compose97i_buffered_init(cs+level, sb, height>>level, stride_line<<level); break;
 681         case DWT_53: spatial_compose53i_buffered_init(cs+level, sb, height>>level, stride_line<<level); break;
 682         }
 683     }
 684 }
 685
 686 void ff_spatial_idwt_buffered_slice(DWTContext *dsp, DWTCompose *cs, slice_buffer * slice_buf, int width, int height, int stride_line, int type, int decomposition_count, int y){
 687     const int support = type==1 ? 3 : 5;
 688     int level;
 689     if(type==2) return;
 690
 691     for(level=decomposition_count-1; level>=0; level--){
 692         while(cs[level].y <= FFMIN((y>>level)+support, height>>level)){
 693             switch(type){
 694             case DWT_97: spatial_compose97i_dy_buffered(dsp, cs+level, slice_buf, width>>level, height>>level, stride_line<<level);
 695                 break;
 696             case DWT_53: spatial_compose53i_dy_buffered(cs+level, slice_buf, width>>level, height>>level, stride_line<<level);
 697                 break;
 698             }
 699         }
 700     }
 701 }
 702
 703 static void ff_spatial_idwt_init(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){
 704     int level;
 705     for(level=decomposition_count-1; level>=0; level--){
 706         switch(type){
 707         case DWT_97: spatial_compose97i_init(cs+level, buffer, height>>level, stride<<level); break;
 708         case DWT_53: spatial_compose53i_init(cs+level, buffer, height>>level, stride<<level); break;
 709         }
 710     }
 711 }
 712
 713 static void ff_spatial_idwt_slice(DWTCompose *cs, IDWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count, int y){
 714     const int support = type==1 ? 3 : 5;
 715     int level;
 716     if(type==2) return;
 717
 718     for(level=decomposition_count-1; level>=0; level--){
 719         while(cs[level].y <= FFMIN((y>>level)+support, height>>level)){
 720             switch(type){
 721             case DWT_97: spatial_compose97i_dy(cs+level, buffer, width>>level, height>>level, stride<<level);
 722                 break;
 723             case DWT_53: spatial_compose53i_dy(cs+level, buffer, width>>level, height>>level, stride<<level);
 724                 break;
 725             }
 726         }
 727     }
 728 }
 729
 730 void ff_spatial_idwt(IDWTELEM *buffer, int width, int height, int stride, int type, int decomposition_count){
 731         DWTCompose cs[MAX_DECOMPOSITIONS];
 732         int y;
 733         ff_spatial_idwt_init(cs, buffer, width, height, stride, type, decomposition_count);
 734         for(y=0; y<height; y+=4)
 735             ff_spatial_idwt_slice(cs, buffer, width, height, stride, type, decomposition_count, y);
 736 }
 737
 738 static inline int w_c(void *v, uint8_t * pix1, uint8_t * pix2, int line_size, int w, int h, int type){
 739     int s, i, j;
 740     const int dec_count= w==8 ? 3 : 4;
 741     int tmp[32*32];
 742     int level, ori;
 743     static const int scale[2][2][4][4]={
 744       {
 745         {
 746             // 9/7 8x8 dec=3
 747             {268, 239, 239, 213},
 748             {  0, 224, 224, 152},
 749             {  0, 135, 135, 110},
 750         },{
 751             // 9/7 16x16 or 32x32 dec=4
 752             {344, 310, 310, 280},
 753             {  0, 320, 320, 228},
 754             {  0, 175, 175, 136},
 755             {  0, 129, 129, 102},
 756         }
 757       },{
 758         {
 759             // 5/3 8x8 dec=3
 760             {275, 245, 245, 218},
 761             {  0, 230, 230, 156},
 762             {  0, 138, 138, 113},
 763         },{
 764             // 5/3 16x16 or 32x32 dec=4
 765             {352, 317, 317, 286},
 766             {  0, 328, 328, 233},
 767             {  0, 180, 180, 140},
 768             {  0, 132, 132, 105},
 769         }
 770       }
 771     };
 772
 773     for (i = 0; i < h; i++) {
 774         for (j = 0; j < w; j+=4) {
 775             tmp[32*i+j+0] = (pix1[j+0] - pix2[j+0])<<4;
 776             tmp[32*i+j+1] = (pix1[j+1] - pix2[j+1])<<4;
 777             tmp[32*i+j+2] = (pix1[j+2] - pix2[j+2])<<4;
 778             tmp[32*i+j+3] = (pix1[j+3] - pix2[j+3])<<4;
 779         }
 780         pix1 += line_size;
 781         pix2 += line_size;
 782     }
 783
 784     ff_spatial_dwt(tmp, w, h, 32, type, dec_count);
 785
 786     s=0;
 787     assert(w==h);
 788     for(level=0; level<dec_count; level++){
 789         for(ori= level ? 1 : 0; ori<4; ori++){
 790             int size= w>>(dec_count-level);
 791             int sx= (ori&1) ? size : 0;
 792             int stride= 32<<(dec_count-level);
 793             int sy= (ori&2) ? stride>>1 : 0;
 794
 795             for(i=0; i<size; i++){
 796                 for(j=0; j<size; j++){
 797                     int v= tmp[sx + sy + i*stride + j] * scale[type][dec_count-3][level][ori];
 798                     s += FFABS(v);
 799                 }
 800             }
 801         }
 802     }
 803     assert(s>=0);
 804     return s>>9;
 805 }
 806
 807 static int w53_8_c(void *v, uint8_t * pix1, uint8_t * pix2, int line_size, int h){
 808     return w_c(v, pix1, pix2, line_size,  8, h, 1);
 809 }
 810
 811 static int w97_8_c(void *v, uint8_t * pix1, uint8_t * pix2, int line_size, int h){
 812     return w_c(v, pix1, pix2, line_size,  8, h, 0);
 813 }
 814
 815 static int w53_16_c(void *v, uint8_t * pix1, uint8_t * pix2, int line_size, int h){
 816     return w_c(v, pix1, pix2, line_size, 16, h, 1);
 817 }
 818
 819 static int w97_16_c(void *v, uint8_t * pix1, uint8_t * pix2, int line_size, int h){
 820     return w_c(v, pix1, pix2, line_size, 16, h, 0);
 821 }
 822
 823 int ff_w53_32_c(void *v, uint8_t * pix1, uint8_t * pix2, int line_size, int h){
 824     return w_c(v, pix1, pix2, line_size, 32, h, 1);
 825 }
 826
 827 int ff_w97_32_c(void *v, uint8_t * pix1, uint8_t * pix2, int line_size, int h){
 828     return w_c(v, pix1, pix2, line_size, 32, h, 0);
 829 }
 830
 831 void ff_dsputil_init_dwt(DSPContext *c)
 832 {
 833     c->w53[0]= w53_16_c;
 834     c->w53[1]= w53_8_c;
 835     c->w97[0]= w97_16_c;
 836     c->w97[1]= w97_8_c;
 837 }
 838
 839 void ff_dwt_init(DWTContext *c)
 840 {
 841     c->vertical_compose97i = ff_snow_vertical_compose97i;
 842     c->horizontal_compose97i = ff_snow_horizontal_compose97i;
 843     c->inner_add_yblock = ff_snow_inner_add_yblock;
 844
 845     if (HAVE_MMX) ff_dwt_init_x86(c);
 846 }
 847
 848
 849 static av_always_inline
 850 void interleave(IDWTELEM *dst, IDWTELEM *src0, IDWTELEM *src1, int w2, int add, int shift)
 851 {
 852     int i;
 853     for (i = 0; i < w2; i++) {
 854         dst[2*i  ] = (src0[i] + add) >> shift;
 855         dst[2*i+1] = (src1[i] + add) >> shift;
 856     }
 857 }
 858
 859 static void horizontal_compose_dirac53i(IDWTELEM *b, IDWTELEM *temp, int w)
 860 {
 861     const int w2 = w >> 1;
 862     int x;
 863
 864     temp[0] = COMPOSE_53iL0(b[w2], b[0], b[w2]);
 865     for (x = 1; x < w2; x++) {
 866         temp[x     ] = COMPOSE_53iL0     (b[x+w2-1], b[x     ], b[x+w2]);
 867         temp[x+w2-1] = COMPOSE_DIRAC53iH0(temp[x-1], b[x+w2-1], temp[x]);
 868     }
 869     temp[w-1] = COMPOSE_DIRAC53iH0(temp[w2-1], b[w-1], temp[w2-1]);
 870
 871     interleave(b, temp, temp+w2, w2, 1, 1);
 872 }
 873
 874 static void horizontal_compose_dd97i(IDWTELEM *b, IDWTELEM *tmp, int w)
 875 {
 876     const int w2 = w >> 1;
 877     int x;
 878
 879     tmp[0] = COMPOSE_53iL0(b[w2], b[0], b[w2]);
 880     for (x = 1; x < w2; x++)
 881         tmp[x] = COMPOSE_53iL0(b[x+w2-1], b[x], b[x+w2]);
 882
 883     // extend the edges
 884     tmp[-1]   = tmp[0];
 885     tmp[w2+1] = tmp[w2] = tmp[w2-1];
 886
 887     for (x = 0; x < w2; x++) {
 888         b[2*x  ] = (tmp[x] + 1)>>1;
 889         b[2*x+1] = (COMPOSE_DD97iH0(tmp[x-1], tmp[x], b[x+w2], tmp[x+1], tmp[x+2]) + 1)>>1;
 890     }
 891 }
 892
 893 static void horizontal_compose_dd137i(IDWTELEM *b, IDWTELEM *tmp, int w)
 894 {
 895     const int w2 = w >> 1;
 896     int x;
 897
 898     tmp[0] = COMPOSE_DD137iL0(b[w2], b[w2], b[0], b[w2  ], b[w2+1]);
 899     tmp[1] = COMPOSE_DD137iL0(b[w2], b[w2], b[1], b[w2+1], b[w2+2]);
 900     for (x = 2; x < w2-1; x++)
 901         tmp[x] = COMPOSE_DD137iL0(b[x+w2-2], b[x+w2-1], b[x], b[x+w2], b[x+w2+1]);
 902     tmp[w2-1] = COMPOSE_DD137iL0(b[w-3], b[w-2], b[w2-1], b[w-1], b[w-1]);
 903
 904     // extend the edges
 905     tmp[-1]   = tmp[0];
 906     tmp[w2+1] = tmp[w2] = tmp[w2-1];
 907
 908     for (x = 0; x < w2; x++) {
 909         b[2*x  ] = (tmp[x] + 1)>>1;
 910         b[2*x+1] = (COMPOSE_DD97iH0(tmp[x-1], tmp[x], b[x+w2], tmp[x+1], tmp[x+2]) + 1)>>1;
 911     }
 912 }
 913
 914 static av_always_inline
 915 void horizontal_compose_haari(IDWTELEM *b, IDWTELEM *temp, int w, int shift)
 916 {
 917     const int w2 = w >> 1;
 918     int x;
 919
 920     for (x = 0; x < w2; x++) {
 921         temp[x   ] = COMPOSE_HAARiL0(b[x   ], b[x+w2]);
 922         temp[x+w2] = COMPOSE_HAARiH0(b[x+w2], temp[x]);
 923     }
 924
 925     interleave(b, temp, temp+w2, w2, shift, shift);
 926 }
 927
 928 static void horizontal_compose_haar0i(IDWTELEM *b, IDWTELEM *temp, int w)
 929 {
 930     horizontal_compose_haari(b, temp, w, 0);
 931 }
 932
 933 static void horizontal_compose_haar1i(IDWTELEM *b, IDWTELEM *temp, int w)
 934 {
 935     horizontal_compose_haari(b, temp, w, 1);
 936 }
 937
 938 static void horizontal_compose_fidelityi(IDWTELEM *b, IDWTELEM *tmp, int w)
 939 {
 940     const int w2 = w >> 1;
 941     int i, x;
 942     IDWTELEM v[8];
 943
 944     for (x = 0; x < w2; x++) {
 945         for (i = 0; i < 8; i++)
 946             v[i] = b[av_clip(x-3+i, 0, w2-1)];
 947         tmp[x] = COMPOSE_FIDELITYiH0(v[0], v[1], v[2], v[3], b[x+w2], v[4], v[5], v[6], v[7]);
 948     }
 949
 950     for (x = 0; x < w2; x++) {
 951         for (i = 0; i < 8; i++)
 952             v[i] = tmp[av_clip(x-4+i, 0, w2-1)];
 953         tmp[x+w2] = COMPOSE_FIDELITYiL0(v[0], v[1], v[2], v[3], b[x], v[4], v[5], v[6], v[7]);
 954     }
 955
 956     interleave(b, tmp+w2, tmp, w2, 0, 0);
 957 }
 958
 959 static void horizontal_compose_daub97i(IDWTELEM *b, IDWTELEM *temp, int w)
 960 {
 961     const int w2 = w >> 1;
 962     int x, b0, b1, b2;
 963
 964     temp[0] = COMPOSE_DAUB97iL1(b[w2], b[0], b[w2]);
 965     for (x = 1; x < w2; x++) {
 966         temp[x     ] = COMPOSE_DAUB97iL1(b[x+w2-1], b[x     ], b[x+w2]);
 967         temp[x+w2-1] = COMPOSE_DAUB97iH1(temp[x-1], b[x+w2-1], temp[x]);
 968     }
 969     temp[w-1] = COMPOSE_DAUB97iH1(temp[w2-1], b[w-1], temp[w2-1]);
 970
 971     // second stage combined with interleave and shift
 972     b0 = b2 = COMPOSE_DAUB97iL0(temp[w2], temp[0], temp[w2]);
 973     b[0] = (b0 + 1) >> 1;
 974     for (x = 1; x < w2; x++) {
 975         b2 = COMPOSE_DAUB97iL0(temp[x+w2-1], temp[x     ], temp[x+w2]);
 976         b1 = COMPOSE_DAUB97iH0(          b0, temp[x+w2-1], b2        );
 977         b[2*x-1] = (b1 + 1) >> 1;
 978         b[2*x  ] = (b2 + 1) >> 1;
 979         b0 = b2;
 980     }
 981     b[w-1] = (COMPOSE_DAUB97iH0(b2, temp[w-1], b2) + 1) >> 1;
 982 }
 983
 984 static void vertical_compose_dirac53iH0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width)
 985 {
 986     int i;
 987
 988     for(i=0; i<width; i++){
 989         b1[i] = COMPOSE_DIRAC53iH0(b0[i], b1[i], b2[i]);
 990     }
 991 }
 992
 993 static void vertical_compose_dd97iH0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2,
 994                                   IDWTELEM *b3, IDWTELEM *b4, int width)
 995 {
 996     int i;
 997
 998     for(i=0; i<width; i++){
 999         b2[i] = COMPOSE_DD97iH0(b0[i], b1[i], b2[i], b3[i], b4[i]);
1000     }
1001 }
1002
1003 static void vertical_compose_dd137iL0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2,
1004                                       IDWTELEM *b3, IDWTELEM *b4, int width)
1005 {
1006     int i;
1007
1008     for(i=0; i<width; i++){
1009         b2[i] = COMPOSE_DD137iL0(b0[i], b1[i], b2[i], b3[i], b4[i]);
1010     }
1011 }
1012
1013 static void vertical_compose_haar(IDWTELEM *b0, IDWTELEM *b1, int width)
1014 {
1015     int i;
1016
1017     for (i = 0; i < width; i++) {
1018         b0[i] = COMPOSE_HAARiL0(b0[i], b1[i]);
1019         b1[i] = COMPOSE_HAARiH0(b1[i], b0[i]);
1020     }
1021 }
1022
1023 static void vertical_compose_fidelityiH0(IDWTELEM *dst, IDWTELEM *b[8], int width)
1024 {
1025     int i;
1026
1027     for(i=0; i<width; i++){
1028         dst[i] = COMPOSE_FIDELITYiH0(b[0][i], b[1][i], b[2][i], b[3][i], dst[i], b[4][i], b[5][i], b[6][i], b[7][i]);
1029     }
1030 }
1031
1032 static void vertical_compose_fidelityiL0(IDWTELEM *dst, IDWTELEM *b[8], int width)
1033 {
1034     int i;
1035
1036     for(i=0; i<width; i++){
1037         dst[i] = COMPOSE_FIDELITYiL0(b[0][i], b[1][i], b[2][i], b[3][i], dst[i], b[4][i], b[5][i], b[6][i], b[7][i]);
1038     }
1039 }
1040
1041 static void vertical_compose_daub97iH0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width)
1042 {
1043     int i;
1044
1045     for(i=0; i<width; i++){
1046         b1[i] = COMPOSE_DAUB97iH0(b0[i], b1[i], b2[i]);
1047     }
1048 }
1049
1050 static void vertical_compose_daub97iH1(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width)
1051 {
1052     int i;
1053
1054     for(i=0; i<width; i++){
1055         b1[i] = COMPOSE_DAUB97iH1(b0[i], b1[i], b2[i]);
1056     }
1057 }
1058
1059 static void vertical_compose_daub97iL0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width)
1060 {
1061     int i;
1062
1063     for(i=0; i<width; i++){
1064         b1[i] = COMPOSE_DAUB97iL0(b0[i], b1[i], b2[i]);
1065     }
1066 }
1067
1068 static void vertical_compose_daub97iL1(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2, int width)
1069 {
1070     int i;
1071
1072     for(i=0; i<width; i++){
1073         b1[i] = COMPOSE_DAUB97iL1(b0[i], b1[i], b2[i]);
1074     }
1075 }
1076
1077
1078 static void spatial_compose_dd97i_dy(DWTContext *d, int level, int width, int height, int stride)
1079 {
1080     vertical_compose_3tap vertical_compose_l0 = d->vertical_compose_l0;
1081     vertical_compose_5tap vertical_compose_h0 = d->vertical_compose_h0;
1082     DWTCompose *cs = d->cs + level;
1083
1084     int i, y = cs->y;
1085     IDWTELEM *b[8];
1086     for (i = 0; i < 6; i++)
1087         b[i] = cs->b[i];
1088     b[6] = d->buffer + av_clip(y+5, 0, height-2)*stride;
1089     b[7] = d->buffer + av_clip(y+6, 1, height-1)*stride;
1090
1091         if(y+5<(unsigned)height) vertical_compose_l0(      b[5], b[6], b[7],       width);
1092         if(y+1<(unsigned)height) vertical_compose_h0(b[0], b[2], b[3], b[4], b[6], width);
1093
1094         if(y-1<(unsigned)height) d->horizontal_compose(b[0], d->temp, width);
1095         if(y+0<(unsigned)height) d->horizontal_compose(b[1], d->temp, width);
1096
1097     for (i = 0; i < 6; i++)
1098         cs->b[i] = b[i+2];
1099     cs->y += 2;
1100 }
1101
1102 static void spatial_compose_dirac53i_dy(DWTContext *d, int level, int width, int height, int stride)
1103 {
1104     vertical_compose_3tap vertical_compose_l0 = d->vertical_compose_l0;
1105     vertical_compose_3tap vertical_compose_h0 = d->vertical_compose_h0;
1106     DWTCompose *cs = d->cs + level;
1107
1108     int y= cs->y;
1109     IDWTELEM *b[4] = { cs->b[0], cs->b[1] };
1110     b[2] = d->buffer + mirror(y+1, height-1)*stride;
1111     b[3] = d->buffer + mirror(y+2, height-1)*stride;
1112
1113         if(y+1<(unsigned)height) vertical_compose_l0(b[1], b[2], b[3], width);
1114         if(y+0<(unsigned)height) vertical_compose_h0(b[0], b[1], b[2], width);
1115
1116         if(y-1<(unsigned)height) d->horizontal_compose(b[0], d->temp, width);
1117         if(y+0<(unsigned)height) d->horizontal_compose(b[1], d->temp, width);
1118
1119     cs->b[0] = b[2];
1120     cs->b[1] = b[3];
1121     cs->y += 2;
1122 }
1123
1124
1125 static void spatial_compose_dd137i_dy(DWTContext *d, int level, int width, int height, int stride)
1126 {
1127     vertical_compose_5tap vertical_compose_l0 = d->vertical_compose_l0;
1128     vertical_compose_5tap vertical_compose_h0 = d->vertical_compose_h0;
1129     DWTCompose *cs = d->cs + level;
1130
1131     int i, y = cs->y;
1132     IDWTELEM *b[10];
1133     for (i = 0; i < 8; i++)
1134         b[i] = cs->b[i];
1135     b[8] = d->buffer + av_clip(y+7, 0, height-2)*stride;
1136     b[9] = d->buffer + av_clip(y+8, 1, height-1)*stride;
1137
1138         if(y+5<(unsigned)height) vertical_compose_l0(b[3], b[5], b[6], b[7], b[9], width);
1139         if(y+1<(unsigned)height) vertical_compose_h0(b[0], b[2], b[3], b[4], b[6], width);
1140
1141         if(y-1<(unsigned)height) d->horizontal_compose(b[0], d->temp, width);
1142         if(y+0<(unsigned)height) d->horizontal_compose(b[1], d->temp, width);
1143
1144     for (i = 0; i < 8; i++)
1145         cs->b[i] = b[i+2];
1146     cs->y += 2;
1147 }
1148
1149 // haar makes the assumption that height is even (always true for dirac)
1150 static void spatial_compose_haari_dy(DWTContext *d, int level, int width, int height, int stride)
1151 {
1152     vertical_compose_2tap vertical_compose = d->vertical_compose;
1153     int y = d->cs[level].y;
1154     IDWTELEM *b0 = d->buffer + (y-1)*stride;
1155     IDWTELEM *b1 = d->buffer + (y  )*stride;
1156
1157     vertical_compose(b0, b1, width);
1158     d->horizontal_compose(b0, d->temp, width);
1159     d->horizontal_compose(b1, d->temp, width);
1160
1161     d->cs[level].y += 2;
1162 }
1163
1164 // Don't do sliced idwt for fidelity; the 9 tap filter makes it a bit annoying
1165 // Fortunately, this filter isn't used in practice.
1166 static void spatial_compose_fidelity(DWTContext *d, int level, int width, int height, int stride)
1167 {
1168     vertical_compose_9tap vertical_compose_l0 = d->vertical_compose_l0;
1169     vertical_compose_9tap vertical_compose_h0 = d->vertical_compose_h0;
1170     int i, y;
1171     IDWTELEM *b[8];
1172
1173     for (y = 1; y < height; y += 2) {
1174         for (i = 0; i < 8; i++)
1175             b[i] = d->buffer + av_clip((y-7 + 2*i), 0, height-2)*stride;
1176         vertical_compose_h0(d->buffer + y*stride, b, width);
1177     }
1178
1179     for (y = 0; y < height; y += 2) {
1180         for (i = 0; i < 8; i++)
1181             b[i] = d->buffer + av_clip((y-7 + 2*i), 1, height-1)*stride;
1182         vertical_compose_l0(d->buffer + y*stride, b, width);
1183     }
1184
1185     for (y = 0; y < height; y++)
1186         d->horizontal_compose(d->buffer + y*stride, d->temp, width);
1187
1188     d->cs[level].y = height+1;
1189 }
1190
1191 static void spatial_compose_daub97i_dy(DWTContext *d, int level, int width, int height, int stride)
1192 {
1193     vertical_compose_3tap vertical_compose_l0 = d->vertical_compose_l0;
1194     vertical_compose_3tap vertical_compose_h0 = d->vertical_compose_h0;
1195     vertical_compose_3tap vertical_compose_l1 = d->vertical_compose_l1;
1196     vertical_compose_3tap vertical_compose_h1 = d->vertical_compose_h1;
1197     DWTCompose *cs = d->cs + level;
1198
1199     int i, y = cs->y;
1200     IDWTELEM *b[6];
1201     for (i = 0; i < 4; i++)
1202         b[i] = cs->b[i];
1203     b[4] = d->buffer + mirror(y+3, height-1)*stride;
1204     b[5] = d->buffer + mirror(y+4, height-1)*stride;
1205
1206         if(y+3<(unsigned)height) vertical_compose_l1(b[3], b[4], b[5], width);
1207         if(y+2<(unsigned)height) vertical_compose_h1(b[2], b[3], b[4], width);
1208         if(y+1<(unsigned)height) vertical_compose_l0(b[1], b[2], b[3], width);
1209         if(y+0<(unsigned)height) vertical_compose_h0(b[0], b[1], b[2], width);
1210
1211         if(y-1<(unsigned)height) d->horizontal_compose(b[0], d->temp, width);
1212         if(y+0<(unsigned)height) d->horizontal_compose(b[1], d->temp, width);
1213
1214     for (i = 0; i < 4; i++)
1215         cs->b[i] = b[i+2];
1216     cs->y += 2;
1217 }
1218
1219
1220 static void spatial_compose97i_init2(DWTCompose *cs, IDWTELEM *buffer, int height, int stride)
1221 {
1222     cs->b[0] = buffer + mirror(-3-1, height-1)*stride;
1223     cs->b[1] = buffer + mirror(-3  , height-1)*stride;
1224     cs->b[2] = buffer + mirror(-3+1, height-1)*stride;
1225     cs->b[3] = buffer + mirror(-3+2, height-1)*stride;
1226     cs->y = -3;
1227 }
1228
1229 static void spatial_compose53i_init2(DWTCompose *cs, IDWTELEM *buffer, int height, int stride)
1230 {
1231     cs->b[0] = buffer + mirror(-1-1, height-1)*stride;
1232     cs->b[1] = buffer + mirror(-1  , height-1)*stride;
1233     cs->y = -1;
1234 }
1235
1236 static void spatial_compose_dd97i_init(DWTCompose *cs, IDWTELEM *buffer, int height, int stride)
1237 {
1238     cs->b[0] = buffer + av_clip(-5-1, 0, height-2)*stride;
1239     cs->b[1] = buffer + av_clip(-5  , 1, height-1)*stride;
1240     cs->b[2] = buffer + av_clip(-5+1, 0, height-2)*stride;
1241     cs->b[3] = buffer + av_clip(-5+2, 1, height-1)*stride;
1242     cs->b[4] = buffer + av_clip(-5+3, 0, height-2)*stride;
1243     cs->b[5] = buffer + av_clip(-5+4, 1, height-1)*stride;
1244     cs->y = -5;
1245 }
1246
1247 static void spatial_compose_dd137i_init(DWTCompose *cs, IDWTELEM *buffer, int height, int stride)
1248 {
1249     cs->b[0] = buffer + av_clip(-5-1, 0, height-2)*stride;
1250     cs->b[1] = buffer + av_clip(-5  , 1, height-1)*stride;
1251     cs->b[2] = buffer + av_clip(-5+1, 0, height-2)*stride;
1252     cs->b[3] = buffer + av_clip(-5+2, 1, height-1)*stride;
1253     cs->b[4] = buffer + av_clip(-5+3, 0, height-2)*stride;
1254     cs->b[5] = buffer + av_clip(-5+4, 1, height-1)*stride;
1255     cs->b[6] = buffer + av_clip(-5+5, 0, height-2)*stride;
1256     cs->b[7] = buffer + av_clip(-5+6, 1, height-1)*stride;
1257     cs->y = -5;
1258 }
1259
1260 int ff_spatial_idwt_init2(DWTContext *d, IDWTELEM *buffer, int width, int height,
1261                           int stride, enum dwt_type type, int decomposition_count,
1262                           IDWTELEM *temp)
1263 {
1264     int level;
1265
1266     d->buffer = buffer;
1267     d->width = width;
1268     d->height = height;
1269     d->stride = stride;
1270     d->decomposition_count = decomposition_count;
1271     d->temp = temp + 8;
1272
1273     for(level=decomposition_count-1; level>=0; level--){
1274         int hl = height >> level;
1275         int stride_l = stride << level;
1276
1277         switch(type){
1278         case DWT_DIRAC_DD9_7:
1279             spatial_compose_dd97i_init(d->cs+level, buffer, hl, stride_l);
1280             break;
1281         case DWT_DIRAC_LEGALL5_3:
1282             spatial_compose53i_init2(d->cs+level, buffer, hl, stride_l);
1283             break;
1284         case DWT_DIRAC_DD13_7:
1285             spatial_compose_dd137i_init(d->cs+level, buffer, hl, stride_l);
1286             break;
1287         case DWT_DIRAC_HAAR0:
1288         case DWT_DIRAC_HAAR1:
1289             d->cs[level].y = 1;
1290             break;
1291         case DWT_DIRAC_DAUB9_7:
1292             spatial_compose97i_init2(d->cs+level, buffer, hl, stride_l);
1293             break;
1294         default:
1295             d->cs[level].y = 0;
1296             break;
1297         }
1298     }
1299
1300     switch (type) {
1301     case DWT_DIRAC_DD9_7:
1302         d->spatial_compose = spatial_compose_dd97i_dy;
1303         d->vertical_compose_l0 = vertical_compose53iL0;
1304         d->vertical_compose_h0 = vertical_compose_dd97iH0;
1305         d->horizontal_compose = horizontal_compose_dd97i;
1306         d->support = 7;
1307         break;
1308     case DWT_DIRAC_LEGALL5_3:
1309         d->spatial_compose = spatial_compose_dirac53i_dy;
1310         d->vertical_compose_l0 = vertical_compose53iL0;
1311         d->vertical_compose_h0 = vertical_compose_dirac53iH0;
1312         d->horizontal_compose = horizontal_compose_dirac53i;
1313         d->support = 3;
1314         break;
1315     case DWT_DIRAC_DD13_7:
1316         d->spatial_compose = spatial_compose_dd137i_dy;
1317         d->vertical_compose_l0 = vertical_compose_dd137iL0;
1318         d->vertical_compose_h0 = vertical_compose_dd97iH0;
1319         d->horizontal_compose = horizontal_compose_dd137i;
1320         d->support = 7;
1321         break;
1322     case DWT_DIRAC_HAAR0:
1323     case DWT_DIRAC_HAAR1:
1324         d->spatial_compose = spatial_compose_haari_dy;
1325         d->vertical_compose = vertical_compose_haar;
1326         if (type == DWT_DIRAC_HAAR0)
1327             d->horizontal_compose = horizontal_compose_haar0i;
1328         else
1329             d->horizontal_compose = horizontal_compose_haar1i;
1330         d->support = 1;
1331         break;
1332     case DWT_DIRAC_FIDELITY:
1333         d->spatial_compose = spatial_compose_fidelity;
1334         d->vertical_compose_l0 = vertical_compose_fidelityiL0;
1335         d->vertical_compose_h0 = vertical_compose_fidelityiH0;
1336         d->horizontal_compose = horizontal_compose_fidelityi;
1337         break;
1338     case DWT_DIRAC_DAUB9_7:
1339         d->spatial_compose = spatial_compose_daub97i_dy;
1340         d->vertical_compose_l0 = vertical_compose_daub97iL0;
1341         d->vertical_compose_h0 = vertical_compose_daub97iH0;
1342         d->vertical_compose_l1 = vertical_compose_daub97iL1;
1343         d->vertical_compose_h1 = vertical_compose_daub97iH1;
1344         d->horizontal_compose = horizontal_compose_daub97i;
1345         d->support = 5;
1346         break;
1347     default:
1348         av_log(NULL, AV_LOG_ERROR, "Unknown wavelet type %d\n", type);
1349         return -1;
1350     }
1351
1352     if (HAVE_MMX) ff_spatial_idwt_init_mmx(d, type);
1353
1354     return 0;
1355 }
1356
1357 void ff_spatial_idwt_slice2(DWTContext *d, int y)
1358 {
1359     int level, support = d->support;
1360
1361     for (level = d->decomposition_count-1; level >= 0; level--) {
1362         int wl = d->width  >> level;
1363         int hl = d->height >> level;
1364         int stride_l = d->stride << level;
1365
1366         while (d->cs[level].y <= FFMIN((y>>level)+support, hl))
1367             d->spatial_compose(d, level, wl, hl, stride_l);
1368     }
1369 }
1370
1371 int ff_spatial_idwt2(IDWTELEM *buffer, int width, int height, int stride,
1372                      enum dwt_type type, int decomposition_count, IDWTELEM *temp)
1373 {
1374     DWTContext d;
1375     int y;
1376
1377     if (ff_spatial_idwt_init2(&d, buffer, width, height, stride, type, decomposition_count, temp))
1378         return -1;
1379
1380     for (y = 0; y < d.height; y += 4)
1381         ff_spatial_idwt_slice2(&d, y);
1382
1383     return 0;
1384 }