dwt: K&R prettyprinting cosmetics
[ffmpeg.git] / libavcodec / dwt.c
1 /*
2  * Copyright (C) 2004-2010 Michael Niedermayer <michaelni@gmx.at>
3  *
4  * This file is part of Libav.
5  *
6  * Libav is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU Lesser General Public
8  * License as published by the Free Software Foundation; either
9  * version 2.1 of the License, or (at your option) any later version.
10  *
11  * Libav is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with Libav; if not, write to the Free Software
18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19  */
20
21 #include "libavutil/attributes.h"
22 #include "dsputil.h"
23 #include "dwt.h"
24
25 void ff_slice_buffer_init(slice_buffer *buf, int line_count,
26                           int max_allocated_lines, int line_width,
27                           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     buf->data_stack_top = max_allocated_lines - 1;
42 }
43
44 IDWTELEM *ff_slice_buffer_load_line(slice_buffer *buf, int line)
45 {
46     IDWTELEM *buffer;
47
48     assert(buf->data_stack_top >= 0);
49 //  assert(!buf->line[line]);
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 void ff_slice_buffer_destroy(slice_buffer *buf)
82 {
83     int i;
84     ff_slice_buffer_flush(buf);
85
86     for (i = buf->data_count - 1; i >= 0; i--)
87         av_freep(&buf->data_stack[i]);
88     av_freep(&buf->data_stack);
89     av_freep(&buf->line);
90 }
91
92 static inline int mirror(int v, int m)
93 {
94     while ((unsigned)v > (unsigned)m) {
95         v = -v;
96         if (v < 0)
97             v += 2 * m;
98     }
99     return v;
100 }
101
102 static av_always_inline void lift(DWTELEM *dst, DWTELEM *src, DWTELEM *ref,
103                                   int dst_step, int src_step, int ref_step,
104                                   int width, int mul, int add, int shift,
105                                   int highpass, int inverse)
106 {
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] = LIFT(src[i * src_step],
121                                  ((mul * (ref[i * ref_step] +
122                                           ref[(i + 1) * ref_step]) +
123                                    add) >> shift),
124                                  inverse);
125
126     if (mirror_right)
127         dst[w * dst_step] = LIFT(src[w * src_step],
128                                  ((mul * 2 * ref[w * ref_step] + add) >> shift),
129                                  inverse);
130 }
131
132 static av_always_inline void inv_lift(IDWTELEM *dst, IDWTELEM *src, IDWTELEM *ref,
133                                       int dst_step, int src_step, int ref_step,
134                                       int width, int mul, int add, int shift,
135                                       int highpass, int inverse)
136 {
137     const int mirror_left  = !highpass;
138     const int mirror_right = (width & 1) ^ highpass;
139     const int w            = (width >> 1) - 1 + (highpass & width);
140     int i;
141
142 #define LIFT(src, ref, inv) ((src) + ((inv) ? -(ref) : +(ref)))
143     if (mirror_left) {
144         dst[0] = LIFT(src[0], ((mul * 2 * ref[0] + add) >> shift), inverse);
145         dst   += dst_step;
146         src   += src_step;
147     }
148
149     for (i = 0; i < w; i++)
150         dst[i * dst_step] = LIFT(src[i * src_step],
151                                  ((mul * (ref[i * ref_step] +
152                                           ref[(i + 1) * ref_step]) +
153                                    add) >> shift),
154                                  inverse);
155
156     if (mirror_right) {
157         dst[w * dst_step] = LIFT(src[w * src_step],
158                                  ((mul * 2 * ref[w * ref_step] + add) >> shift),
159                                  inverse);
160     }
161 }
162
163 #ifndef liftS
164 static av_always_inline void liftS(DWTELEM *dst, DWTELEM *src, DWTELEM *ref,
165                                    int dst_step, int src_step, int ref_step,
166                                    int width, int mul, int add, int shift,
167                                    int highpass, int inverse)
168 {
169     const int mirror_left  = !highpass;
170     const int mirror_right = (width & 1) ^ highpass;
171     const int w            = (width >> 1) - 1 + (highpass & width);
172     int i;
173
174     assert(shift == 4);
175 #define LIFTS(src, ref, inv)                                            \
176     ((inv) ? (src) + (((ref) + 4 * (src)) >> shift)                     \
177            : -((-16 * (src) + (ref) + add /                             \
178                 4 + 1 + (5 << 25)) / (5 * 4) - (1 << 23)))
179     if (mirror_left) {
180         dst[0] = LIFTS(src[0], mul * 2 * ref[0] + add, inverse);
181         dst   += dst_step;
182         src   += src_step;
183     }
184
185     for (i = 0; i < w; i++)
186         dst[i * dst_step] = LIFTS(src[i * src_step],
187                                   mul * (ref[i * ref_step] +
188                                          ref[(i + 1) * ref_step]) + add,
189                                   inverse);
190
191     if (mirror_right)
192         dst[w * dst_step] = LIFTS(src[w * src_step],
193                                   mul * 2 * ref[w * ref_step] + add,
194                                   inverse);
195 }
196
197 static av_always_inline void inv_liftS(IDWTELEM *dst, IDWTELEM *src,
198                                        IDWTELEM *ref, int dst_step,
199                                        int src_step, int ref_step,
200                                        int width, int mul, int add, int shift,
201                                        int highpass, int inverse)
202 {
203     const int mirror_left  = !highpass;
204     const int mirror_right = (width & 1) ^ highpass;
205     const int w            = (width >> 1) - 1 + (highpass & width);
206     int i;
207
208     assert(shift == 4);
209 #define LIFTS(src, ref, inv)                                            \
210     ((inv) ? (src) + (((ref) + 4 * (src)) >> shift)                     \
211            : -((-16 * (src) + (ref) + add /                             \
212                 4 + 1 + (5 << 25)) / (5 * 4) - (1 << 23)))
213     if (mirror_left) {
214         dst[0] = LIFTS(src[0], mul * 2 * ref[0] + add, inverse);
215         dst   += dst_step;
216         src   += src_step;
217     }
218
219     for (i = 0; i < w; i++)
220         dst[i * dst_step] = LIFTS(src[i * src_step],
221                                   mul * (ref[i * ref_step] +
222                                          ref[(i + 1) * ref_step]) + add,
223                                   inverse);
224
225     if (mirror_right)
226         dst[w * dst_step] = LIFTS(src[w * src_step],
227                                   mul * 2 * ref[w * ref_step] + add, inverse);
228 }
229 #endif /* ! liftS */
230
231 static void horizontal_decompose53i(DWTELEM *b, int width)
232 {
233     DWTELEM temp[width];
234     const int width2 = width >> 1;
235     int x;
236     const int w2 = (width + 1) >> 1;
237
238     for (x = 0; x < width2; x++) {
239         temp[x]      = b[2 * x];
240         temp[x + w2] = b[2 * x + 1];
241     }
242     if (width & 1)
243         temp[x] = b[2 * x];
244 #if 0
245     {
246         int A1, A2, A3, A4;
247         A2             = temp[1];
248         A4             = temp[0];
249         A1             = temp[0 + width2];
250         A1            -= (A2 + A4) >> 1;
251         A4            += (A1 +  1) >> 1;
252         b[0 + width2]  = A1;
253         b[0]           = A4;
254         for (x = 1; x + 1 < width2; x += 2) {
255             A3                 = temp[x + width2];
256             A4                 = temp[x + 1];
257             A3                -= (A2 + A4)     >> 1;
258             A2                += (A1 + A3 + 2) >> 2;
259             b[x + width2]      = A3;
260             b[x]               = A2;
261
262             A1                 = temp[x + 1 + width2];
263             A2                 = temp[x + 2];
264             A1                -= (A2 + A4)     >> 1;
265             A4                += (A1 + A3 + 2) >> 2;
266             b[x + 1 + width2]  = A1;
267             b[x + 1]           = A4;
268         }
269         A3            = temp[width - 1];
270         A3           -= A2;
271         A2           += (A1 + A3 + 2) >> 2;
272         b[width  - 1] = A3;
273         b[width2 - 1] = A2;
274     }
275 #else
276     lift(b + w2, temp + w2, temp,   1, 1, 1, width, -1, 0, 1, 1, 0);
277     lift(b,      temp,      b + w2, 1, 1, 1, width,  1, 2, 2, 0, 0);
278 #endif /* 0 */
279 }
280
281 static void vertical_decompose53iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2,
282                                     int width)
283 {
284     int i;
285
286     for (i = 0; i < width; i++)
287         b1[i] -= (b0[i] + b2[i]) >> 1;
288 }
289
290 static void vertical_decompose53iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2,
291                                     int width)
292 {
293     int i;
294
295     for (i = 0; i < width; i++)
296         b1[i] += (b0[i] + b2[i] + 2) >> 2;
297 }
298
299 static void spatial_decompose53i(DWTELEM *buffer, int width, int height,
300                                  int stride)
301 {
302     int y;
303     DWTELEM *b0 = buffer + mirror(-2 - 1, height - 1) * stride;
304     DWTELEM *b1 = buffer + mirror(-2,     height - 1) * stride;
305
306     for (y = -2; y < height; y += 2) {
307         DWTELEM *b2 = buffer + mirror(y + 1, height - 1) * stride;
308         DWTELEM *b3 = buffer + mirror(y + 2, height - 1) * stride;
309
310         if (y + 1 < (unsigned)height)
311             horizontal_decompose53i(b2, width);
312         if (y + 2 < (unsigned)height)
313             horizontal_decompose53i(b3, width);
314
315         if (y + 1 < (unsigned)height)
316             vertical_decompose53iH0(b1, b2, b3, width);
317         if (y + 0 < (unsigned)height)
318             vertical_decompose53iL0(b0, b1, b2, width);
319
320         b0 = b2;
321         b1 = b3;
322     }
323 }
324
325 static void horizontal_decompose97i(DWTELEM *b, int width)
326 {
327     DWTELEM temp[width];
328     const int w2 = (width + 1) >> 1;
329
330     lift(temp + w2, b + 1, b,         1, 2, 2, width, W_AM, W_AO, W_AS, 1, 1);
331     liftS(temp,     b,     temp + w2, 1, 2, 1, width, W_BM, W_BO, W_BS, 0, 0);
332     lift(b + w2, temp + w2, temp,     1, 1, 1, width, W_CM, W_CO, W_CS, 1, 0);
333     lift(b,      temp,      b + w2,   1, 1, 1, width, W_DM, W_DO, W_DS, 0, 0);
334 }
335
336 static void vertical_decompose97iH0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2,
337                                     int width)
338 {
339     int i;
340
341     for (i = 0; i < width; i++)
342         b1[i] -= (W_AM * (b0[i] + b2[i]) + W_AO) >> W_AS;
343 }
344
345 static void vertical_decompose97iH1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2,
346                                     int width)
347 {
348     int i;
349
350     for (i = 0; i < width; i++)
351         b1[i] += (W_CM * (b0[i] + b2[i]) + W_CO) >> W_CS;
352 }
353
354 static void vertical_decompose97iL0(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2,
355                                     int width)
356 {
357     int i;
358
359     for (i = 0; i < width; i++)
360 #ifdef liftS
361         b1[i] -= (W_BM * (b0[i] + b2[i]) + W_BO) >> W_BS;
362 #else
363         b1[i] = (16 * 4 * b1[i] - 4 * (b0[i] + b2[i]) + W_BO * 5 + (5 << 27)) /
364                 (5 * 16) - (1 << 23);
365 #endif
366 }
367
368 static void vertical_decompose97iL1(DWTELEM *b0, DWTELEM *b1, DWTELEM *b2,
369                                     int width)
370 {
371     int i;
372
373     for (i = 0; i < width; i++)
374         b1[i] += (W_DM * (b0[i] + b2[i]) + W_DO) >> W_DS;
375 }
376
377 static void spatial_decompose97i(DWTELEM *buffer, int width, int height,
378                                  int stride)
379 {
380     int y;
381     DWTELEM *b0 = buffer + mirror(-4 - 1, height - 1) * stride;
382     DWTELEM *b1 = buffer + mirror(-4,     height - 1) * stride;
383     DWTELEM *b2 = buffer + mirror(-4 + 1, height - 1) * stride;
384     DWTELEM *b3 = buffer + mirror(-4 + 2, height - 1) * stride;
385
386     for (y = -4; y < height; y += 2) {
387         DWTELEM *b4 = buffer + mirror(y + 3, height - 1) * stride;
388         DWTELEM *b5 = buffer + mirror(y + 4, height - 1) * stride;
389
390         if (y + 3 < (unsigned)height)
391             horizontal_decompose97i(b4, width);
392         if (y + 4 < (unsigned)height)
393             horizontal_decompose97i(b5, width);
394
395         if (y + 3 < (unsigned)height)
396             vertical_decompose97iH0(b3, b4, b5, width);
397         if (y + 2 < (unsigned)height)
398             vertical_decompose97iL0(b2, b3, b4, width);
399         if (y + 1 < (unsigned)height)
400             vertical_decompose97iH1(b1, b2, b3, width);
401         if (y + 0 < (unsigned)height)
402             vertical_decompose97iL1(b0, b1, b2, width);
403
404         b0 = b2;
405         b1 = b3;
406         b2 = b4;
407         b3 = b5;
408     }
409 }
410
411 void ff_spatial_dwt(DWTELEM *buffer, int width, int height, int stride,
412                     int type, int decomposition_count)
413 {
414     int level;
415
416     for (level = 0; level < decomposition_count; level++) {
417         switch (type) {
418         case DWT_97:
419             spatial_decompose97i(buffer,
420                                  width >> level, height >> level,
421                                  stride << level);
422             break;
423         case DWT_53:
424             spatial_decompose53i(buffer,
425                                  width >> level, height >> level,
426                                  stride << level);
427             break;
428         }
429     }
430 }
431
432 static void horizontal_compose53i(IDWTELEM *b, int width)
433 {
434     IDWTELEM temp[width];
435     const int width2 = width >> 1;
436     const int w2     = (width + 1) >> 1;
437     int x;
438
439     for (x = 0; x < width2; x++) {
440         temp[2 * x]     = b[x];
441         temp[2 * x + 1] = b[x + w2];
442     }
443     if (width & 1)
444         temp[2 * x] = b[x];
445
446     b[0] = temp[0] - ((temp[1] + 1) >> 1);
447     for (x = 2; x < width - 1; x += 2) {
448         b[x]     = temp[x]     - ((temp[x - 1] + temp[x + 1] + 2) >> 2);
449         b[x - 1] = temp[x - 1] + ((b[x - 2]    + b[x]        + 1) >> 1);
450     }
451     if (width & 1) {
452         b[x]     = temp[x]     - ((temp[x - 1]     + 1) >> 1);
453         b[x - 1] = temp[x - 1] + ((b[x - 2] + b[x] + 1) >> 1);
454     } else
455         b[x - 1] = temp[x - 1] + b[x - 2];
456 }
457
458 static void vertical_compose53iH0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2,
459                                   int width)
460 {
461     int i;
462
463     for (i = 0; i < width; i++)
464         b1[i] += (b0[i] + b2[i]) >> 1;
465 }
466
467 static void vertical_compose53iL0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2,
468                                   int width)
469 {
470     int i;
471
472     for (i = 0; i < width; i++)
473         b1[i] -= (b0[i] + b2[i] + 2) >> 2;
474 }
475
476 static void spatial_compose53i_buffered_init(DWTCompose *cs, slice_buffer *sb,
477                                              int height, int stride_line)
478 {
479     cs->b0 = slice_buffer_get_line(sb,
480                                    mirror(-1 - 1, height - 1) * stride_line);
481     cs->b1 = slice_buffer_get_line(sb, mirror(-1, height - 1) * stride_line);
482     cs->y  = -1;
483 }
484
485 static void spatial_compose53i_init(DWTCompose *cs, IDWTELEM *buffer,
486                                     int height, int stride)
487 {
488     cs->b0 = buffer + mirror(-1 - 1, height - 1) * stride;
489     cs->b1 = buffer + mirror(-1,     height - 1) * stride;
490     cs->y  = -1;
491 }
492
493 static void spatial_compose53i_dy_buffered(DWTCompose *cs, slice_buffer *sb,
494                                            int width, int height,
495                                            int stride_line)
496 {
497     int y = cs->y;
498
499     IDWTELEM *b0 = cs->b0;
500     IDWTELEM *b1 = cs->b1;
501     IDWTELEM *b2 = slice_buffer_get_line(sb,
502                                          mirror(y + 1, height - 1) *
503                                          stride_line);
504     IDWTELEM *b3 = slice_buffer_get_line(sb,
505                                          mirror(y + 2, height - 1) *
506                                          stride_line);
507
508     if (y + 1 < (unsigned)height && y < (unsigned)height) {
509         int x;
510
511         for (x = 0; x < width; x++) {
512             b2[x] -= (b1[x] + b3[x] + 2) >> 2;
513             b1[x] += (b0[x] + b2[x])     >> 1;
514         }
515     } else {
516         if (y + 1 < (unsigned)height)
517             vertical_compose53iL0(b1, b2, b3, width);
518         if (y + 0 < (unsigned)height)
519             vertical_compose53iH0(b0, b1, b2, width);
520     }
521
522     if (y - 1 < (unsigned)height)
523         horizontal_compose53i(b0, width);
524     if (y + 0 < (unsigned)height)
525         horizontal_compose53i(b1, width);
526
527     cs->b0  = b2;
528     cs->b1  = b3;
529     cs->y  += 2;
530 }
531
532 static void spatial_compose53i_dy(DWTCompose *cs, IDWTELEM *buffer, int width,
533                                   int height, int stride)
534 {
535     int y        = cs->y;
536     IDWTELEM *b0 = cs->b0;
537     IDWTELEM *b1 = cs->b1;
538     IDWTELEM *b2 = buffer + mirror(y + 1, height - 1) * stride;
539     IDWTELEM *b3 = buffer + mirror(y + 2, height - 1) * stride;
540
541     if (y + 1 < (unsigned)height)
542         vertical_compose53iL0(b1, b2, b3, width);
543     if (y + 0 < (unsigned)height)
544         vertical_compose53iH0(b0, b1, b2, width);
545
546     if (y - 1 < (unsigned)height)
547         horizontal_compose53i(b0, width);
548     if (y + 0 < (unsigned)height)
549         horizontal_compose53i(b1, width);
550
551     cs->b0  = b2;
552     cs->b1  = b3;
553     cs->y  += 2;
554 }
555
556 static void av_unused spatial_compose53i(IDWTELEM *buffer, int width,
557                                          int height, int stride)
558 {
559     DWTCompose cs;
560     spatial_compose53i_init(&cs, buffer, height, stride);
561     while (cs.y <= height)
562         spatial_compose53i_dy(&cs, buffer, width, height, stride);
563 }
564
565 void ff_snow_horizontal_compose97i(IDWTELEM *b, int width)
566 {
567     IDWTELEM temp[width];
568     const int w2 = (width + 1) >> 1;
569
570 #if 0 //maybe more understadable but slower
571     inv_lift(temp,     b,      b + w2, 2, 1, 1, width, W_DM, W_DO, W_DS, 0, 1);
572     inv_lift(temp + 1, b + w2, temp,   2, 1, 2, width, W_CM, W_CO, W_CS, 1, 1);
573
574     inv_liftS(b,    temp,     temp + 1, 2, 2, 2, width, W_BM, W_BO, W_BS, 0, 1);
575     inv_lift(b + 1, temp + 1, b,        2, 2, 2, width, W_AM, W_AO, W_AS, 1, 0);
576 #else
577     int x;
578     temp[0] = b[0] - ((3 * b[w2] + 2) >> 2);
579     for (x = 1; x < (width >> 1); x++) {
580         temp[2 * x]     = b[x] - ((3 * (b[x + w2 - 1] + b[x + w2]) + 4) >> 3);
581         temp[2 * x - 1] = b[x + w2 - 1] - temp[2 * x - 2] - temp[2 * x];
582     }
583     if (width & 1) {
584         temp[2 * x]     = b[x] - ((3 * b[x + w2 - 1] + 2) >> 2);
585         temp[2 * x - 1] = b[x + w2 - 1] - temp[2 * x - 2] - temp[2 * x];
586     } else
587         temp[2 * x - 1] = b[x + w2 - 1] - 2 * temp[2 * x - 2];
588
589     b[0] = temp[0] + ((2 * temp[0] + temp[1] + 4) >> 3);
590     for (x = 2; x < width - 1; x += 2) {
591         b[x]     = temp[x] + ((4 * temp[x] + temp[x - 1] + temp[x + 1] + 8) >> 4);
592         b[x - 1] = temp[x - 1] + ((3 * (b[x - 2] + b[x])) >> 1);
593     }
594     if (width & 1) {
595         b[x]     = temp[x] + ((2 * temp[x] + temp[x - 1] + 4) >> 3);
596         b[x - 1] = temp[x - 1] + ((3 * (b[x - 2] + b[x])) >> 1);
597     } else
598         b[x - 1] = temp[x - 1] + 3 * b[x - 2];
599 #endif
600 }
601
602 static void vertical_compose97iH0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2,
603                                   int width)
604 {
605     int i;
606
607     for (i = 0; i < width; i++)
608         b1[i] += (W_AM * (b0[i] + b2[i]) + W_AO) >> W_AS;
609 }
610
611 static void vertical_compose97iH1(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2,
612                                   int width)
613 {
614     int i;
615
616     for (i = 0; i < width; i++)
617         b1[i] -= (W_CM * (b0[i] + b2[i]) + W_CO) >> W_CS;
618 }
619
620 static void vertical_compose97iL0(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2,
621                                   int width)
622 {
623     int i;
624
625     for (i = 0; i < width; i++)
626 #ifdef liftS
627         b1[i] += (W_BM * (b0[i] + b2[i]) + W_BO) >> W_BS;
628 #else
629         b1[i] += (W_BM * (b0[i] + b2[i]) + 4 * b1[i] + W_BO) >> W_BS;
630 #endif
631 }
632
633 static void vertical_compose97iL1(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2,
634                                   int width)
635 {
636     int i;
637
638     for (i = 0; i < width; i++)
639         b1[i] -= (W_DM * (b0[i] + b2[i]) + W_DO) >> W_DS;
640 }
641
642 void ff_snow_vertical_compose97i(IDWTELEM *b0, IDWTELEM *b1, IDWTELEM *b2,
643                                  IDWTELEM *b3, IDWTELEM *b4, IDWTELEM *b5,
644                                  int width)
645 {
646     int i;
647
648     for (i = 0; i < width; i++) {
649         b4[i] -= (W_DM * (b3[i] + b5[i]) + W_DO) >> W_DS;
650         b3[i] -= (W_CM * (b2[i] + b4[i]) + W_CO) >> W_CS;
651 #ifdef liftS
652         b2[i] += (W_BM * (b1[i] + b3[i]) + W_BO) >> W_BS;
653 #else
654         b2[i] += (W_BM * (b1[i] + b3[i]) + 4 * b2[i] + W_BO) >> W_BS;
655 #endif
656         b1[i] += (W_AM * (b0[i] + b2[i]) + W_AO) >> W_AS;
657     }
658 }
659
660 static void spatial_compose97i_buffered_init(DWTCompose *cs, slice_buffer *sb,
661                                              int height, int stride_line)
662 {
663     cs->b0 = slice_buffer_get_line(sb, mirror(-3 - 1, height - 1) * stride_line);
664     cs->b1 = slice_buffer_get_line(sb, mirror(-3,     height - 1) * stride_line);
665     cs->b2 = slice_buffer_get_line(sb, mirror(-3 + 1, height - 1) * stride_line);
666     cs->b3 = slice_buffer_get_line(sb, mirror(-3 + 2, height - 1) * stride_line);
667     cs->y  = -3;
668 }
669
670 static void spatial_compose97i_init(DWTCompose *cs, IDWTELEM *buffer, int height,
671                                     int stride)
672 {
673     cs->b0 = buffer + mirror(-3 - 1, height - 1) * stride;
674     cs->b1 = buffer + mirror(-3,     height - 1) * stride;
675     cs->b2 = buffer + mirror(-3 + 1, height - 1) * stride;
676     cs->b3 = buffer + mirror(-3 + 2, height - 1) * stride;
677     cs->y  = -3;
678 }
679
680 static void spatial_compose97i_dy_buffered(DWTContext *dsp, DWTCompose *cs,
681                                            slice_buffer *sb, int width,
682                                            int height, int stride_line)
683 {
684     int y = cs->y;
685
686     IDWTELEM *b0 = cs->b0;
687     IDWTELEM *b1 = cs->b1;
688     IDWTELEM *b2 = cs->b2;
689     IDWTELEM *b3 = cs->b3;
690     IDWTELEM *b4 = slice_buffer_get_line(sb,
691                                          mirror(y + 3, height - 1) *
692                                          stride_line);
693     IDWTELEM *b5 = slice_buffer_get_line(sb,
694                                          mirror(y + 4, height - 1) *
695                                          stride_line);
696
697     if (y > 0 && y + 4 < height) {
698         dsp->vertical_compose97i(b0, b1, b2, b3, b4, b5, width);
699     } else {
700         if (y + 3 < (unsigned)height)
701             vertical_compose97iL1(b3, b4, b5, width);
702         if (y + 2 < (unsigned)height)
703             vertical_compose97iH1(b2, b3, b4, width);
704         if (y + 1 < (unsigned)height)
705             vertical_compose97iL0(b1, b2, b3, width);
706         if (y + 0 < (unsigned)height)
707             vertical_compose97iH0(b0, b1, b2, width);
708     }
709
710     if (y - 1 < (unsigned)height)
711         dsp->horizontal_compose97i(b0, width);
712     if (y + 0 < (unsigned)height)
713         dsp->horizontal_compose97i(b1, width);
714
715     cs->b0  = b2;
716     cs->b1  = b3;
717     cs->b2  = b4;
718     cs->b3  = b5;
719     cs->y  += 2;
720 }
721
722 static void spatial_compose97i_dy(DWTCompose *cs, IDWTELEM *buffer, int width,
723                                   int height, int stride)
724 {
725     int y        = cs->y;
726     IDWTELEM *b0 = cs->b0;
727     IDWTELEM *b1 = cs->b1;
728     IDWTELEM *b2 = cs->b2;
729     IDWTELEM *b3 = cs->b3;
730     IDWTELEM *b4 = buffer + mirror(y + 3, height - 1) * stride;
731     IDWTELEM *b5 = buffer + mirror(y + 4, height - 1) * stride;
732
733     if (y + 3 < (unsigned)height)
734         vertical_compose97iL1(b3, b4, b5, width);
735     if (y + 2 < (unsigned)height)
736         vertical_compose97iH1(b2, b3, b4, width);
737     if (y + 1 < (unsigned)height)
738         vertical_compose97iL0(b1, b2, b3, width);
739     if (y + 0 < (unsigned)height)
740         vertical_compose97iH0(b0, b1, b2, width);
741
742     if (y - 1 < (unsigned)height)
743         ff_snow_horizontal_compose97i(b0, width);
744     if (y + 0 < (unsigned)height)
745         ff_snow_horizontal_compose97i(b1, width);
746
747     cs->b0  = b2;
748     cs->b1  = b3;
749     cs->b2  = b4;
750     cs->b3  = b5;
751     cs->y  += 2;
752 }
753
754 static void av_unused spatial_compose97i(IDWTELEM *buffer, int width,
755                                          int height, int stride)
756 {
757     DWTCompose cs;
758     spatial_compose97i_init(&cs, buffer, height, stride);
759     while (cs.y <= height)
760         spatial_compose97i_dy(&cs, buffer, width, height, stride);
761 }
762
763 void ff_spatial_idwt_buffered_init(DWTCompose *cs, slice_buffer *sb, int width,
764                                    int height, int stride_line, int type,
765                                    int decomposition_count)
766 {
767     int level;
768     for (level = decomposition_count - 1; level >= 0; level--) {
769         switch (type) {
770         case DWT_97:
771             spatial_compose97i_buffered_init(cs + level, sb, height >> level,
772                                              stride_line << level);
773             break;
774         case DWT_53:
775             spatial_compose53i_buffered_init(cs + level, sb, height >> level,
776                                              stride_line << level);
777             break;
778         }
779     }
780 }
781
782 void ff_spatial_idwt_buffered_slice(DWTContext *dsp, DWTCompose *cs,
783                                     slice_buffer *slice_buf, int width,
784                                     int height, int stride_line, int type,
785                                     int decomposition_count, int y)
786 {
787     const int support = type == 1 ? 3 : 5;
788     int level;
789     if (type == 2)
790         return;
791
792     for (level = decomposition_count - 1; level >= 0; level--)
793         while (cs[level].y <= FFMIN((y >> level) + support, height >> level)) {
794             switch (type) {
795             case DWT_97:
796                 spatial_compose97i_dy_buffered(dsp, cs + level, slice_buf,
797                                                width >> level,
798                                                height >> level,
799                                                stride_line << level);
800                 break;
801             case DWT_53:
802                 spatial_compose53i_dy_buffered(cs + level, slice_buf,
803                                                width >> level,
804                                                height >> level,
805                                                stride_line << level);
806                 break;
807             }
808         }
809 }
810
811 static void ff_spatial_idwt_init(DWTCompose *cs, IDWTELEM *buffer, int width,
812                                  int height, int stride, int type,
813                                  int decomposition_count)
814 {
815     int level;
816     for (level = decomposition_count - 1; level >= 0; level--) {
817         switch (type) {
818         case DWT_97:
819             spatial_compose97i_init(cs + level, buffer, height >> level,
820                                     stride << level);
821             break;
822         case DWT_53:
823             spatial_compose53i_init(cs + level, buffer, height >> level,
824                                     stride << level);
825             break;
826         }
827     }
828 }
829
830 static void ff_spatial_idwt_slice(DWTCompose *cs, IDWTELEM *buffer, int width,
831                                   int height, int stride, int type,
832                                   int decomposition_count, int y)
833 {
834     const int support = type == 1 ? 3 : 5;
835     int level;
836     if (type == 2)
837         return;
838
839     for (level = decomposition_count - 1; level >= 0; level--)
840         while (cs[level].y <= FFMIN((y >> level) + support, height >> level)) {
841             switch (type) {
842             case DWT_97:
843                 spatial_compose97i_dy(cs + level, buffer, width >> level,
844                                       height >> level, stride << level);
845                 break;
846             case DWT_53:
847                 spatial_compose53i_dy(cs + level, buffer, width >> level,
848                                       height >> level, stride << level);
849                 break;
850             }
851         }
852 }
853
854 void ff_spatial_idwt(IDWTELEM *buffer, int width, int height, int stride,
855                      int type, int decomposition_count)
856 {
857     DWTCompose cs[MAX_DECOMPOSITIONS];
858     int y;
859     ff_spatial_idwt_init(cs, buffer, width, height, stride, type,
860                          decomposition_count);
861     for (y = 0; y < height; y += 4)
862         ff_spatial_idwt_slice(cs, buffer, width, height, stride, type,
863                               decomposition_count, y);
864 }
865
866 static inline int w_c(void *v, uint8_t *pix1, uint8_t *pix2, int line_size,
867                       int w, int h, int type)
868 {
869     int s, i, j;
870     const int dec_count = w == 8 ? 3 : 4;
871     int tmp[32 * 32];
872     int level, ori;
873     static const int scale[2][2][4][4] = {
874         {
875             { // 9/7 8x8 dec=3
876                 { 268, 239, 239, 213 },
877                 { 0,   224, 224, 152 },
878                 { 0,   135, 135, 110 },
879             },
880             { // 9/7 16x16 or 32x32 dec=4
881                 { 344, 310, 310, 280 },
882                 { 0,   320, 320, 228 },
883                 { 0,   175, 175, 136 },
884                 { 0,   129, 129, 102 },
885             }
886         },
887         {
888             { // 5/3 8x8 dec=3
889                 { 275, 245, 245, 218 },
890                 { 0,   230, 230, 156 },
891                 { 0,   138, 138, 113 },
892             },
893             { // 5/3 16x16 or 32x32 dec=4
894                 { 352, 317, 317, 286 },
895                 { 0,   328, 328, 233 },
896                 { 0,   180, 180, 140 },
897                 { 0,   132, 132, 105 },
898             }
899         }
900     };
901
902     for (i = 0; i < h; i++) {
903         for (j = 0; j < w; j += 4) {
904             tmp[32 * i + j + 0] = (pix1[j + 0] - pix2[j + 0]) << 4;
905             tmp[32 * i + j + 1] = (pix1[j + 1] - pix2[j + 1]) << 4;
906             tmp[32 * i + j + 2] = (pix1[j + 2] - pix2[j + 2]) << 4;
907             tmp[32 * i + j + 3] = (pix1[j + 3] - pix2[j + 3]) << 4;
908         }
909         pix1 += line_size;
910         pix2 += line_size;
911     }
912
913     ff_spatial_dwt(tmp, w, h, 32, type, dec_count);
914
915     s = 0;
916     assert(w == h);
917     for (level = 0; level < dec_count; level++)
918         for (ori = level ? 1 : 0; ori < 4; ori++) {
919             int size   = w >> (dec_count - level);
920             int sx     = (ori & 1) ? size : 0;
921             int stride = 32 << (dec_count - level);
922             int sy     = (ori & 2) ? stride >> 1 : 0;
923
924             for (i = 0; i < size; i++)
925                 for (j = 0; j < size; j++) {
926                     int v = tmp[sx + sy + i * stride + j] *
927                             scale[type][dec_count - 3][level][ori];
928                     s += FFABS(v);
929                 }
930         }
931     assert(s >= 0);
932     return s >> 9;
933 }
934
935 static int w53_8_c(void *v, uint8_t *pix1, uint8_t *pix2, int line_size, int h)
936 {
937     return w_c(v, pix1, pix2, line_size, 8, h, 1);
938 }
939
940 static int w97_8_c(void *v, uint8_t *pix1, uint8_t *pix2, int line_size, int h)
941 {
942     return w_c(v, pix1, pix2, line_size, 8, h, 0);
943 }
944
945 static int w53_16_c(void *v, uint8_t *pix1, uint8_t *pix2, int line_size, int h)
946 {
947     return w_c(v, pix1, pix2, line_size, 16, h, 1);
948 }
949
950 static int w97_16_c(void *v, uint8_t *pix1, uint8_t *pix2, int line_size, int h)
951 {
952     return w_c(v, pix1, pix2, line_size, 16, h, 0);
953 }
954
955 int ff_w53_32_c(void *v, uint8_t *pix1, uint8_t *pix2, int line_size, int h)
956 {
957     return w_c(v, pix1, pix2, line_size, 32, h, 1);
958 }
959
960 int ff_w97_32_c(void *v, uint8_t *pix1, uint8_t *pix2, int line_size, int h)
961 {
962     return w_c(v, pix1, pix2, line_size, 32, h, 0);
963 }
964
965 void ff_dsputil_init_dwt(DSPContext *c)
966 {
967     c->w53[0] = w53_16_c;
968     c->w53[1] = w53_8_c;
969     c->w97[0] = w97_16_c;
970     c->w97[1] = w97_8_c;
971 }
972
973 void ff_dwt_init(DWTContext *c)
974 {
975     c->vertical_compose97i   = ff_snow_vertical_compose97i;
976     c->horizontal_compose97i = ff_snow_horizontal_compose97i;
977     c->inner_add_yblock      = ff_snow_inner_add_yblock;
978
979     if (HAVE_MMX)
980         ff_dwt_init_x86(c);
981 }