b7363a03580832079d8824c293f2e7ccd0bb654b
[ffmpeg.git] / libavcodec / jpeg2000dwt.c
1 /*
2  * Discrete wavelet transform
3  * Copyright (c) 2007 Kamil Nowosad
4  * Copyright (c) 2013 Nicolas Bertrand <nicoinattendu@gmail.com>
5  *
6  * This file is part of FFmpeg.
7  *
8  * FFmpeg is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU Lesser General Public
10  * License as published by the Free Software Foundation; either
11  * version 2.1 of the License, or (at your option) any later version.
12  *
13  * FFmpeg is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16  * Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with FFmpeg; if not, write to the Free Software
20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21  */
22
23 /**
24  * @file
25  * Discrete wavelet transform
26  */
27
28 #include "libavutil/avassert.h"
29 #include "libavutil/common.h"
30 #include "libavutil/mem.h"
31 #include "jpeg2000dwt.h"
32 #include "internal.h"
33
34 /* Defines for 9/7 DWT lifting parameters.
35  * Parameters are in float. */
36 #define F_LFTG_ALPHA  1.586134342059924f
37 #define F_LFTG_BETA   0.052980118572961f
38 #define F_LFTG_GAMMA  0.882911075530934f
39 #define F_LFTG_DELTA  0.443506852043971f
40 #define F_LFTG_K      1.230174104914001f
41 #define F_LFTG_X      1.625732422f
42 /* FIXME: Why use 1.625732422 instead of 1/F_LFTG_K?
43  * Incorrect value in JPEG2000 norm.
44  * see (ISO/IEC 15444:1 (version 2002) F.3.8.2 */
45
46 /* Lifting parameters in integer format.
47  * Computed as param = (float param) * (1 << 16) */
48 #define I_LFTG_ALPHA  103949
49 #define I_LFTG_BETA     3472
50 #define I_LFTG_GAMMA   57862
51 #define I_LFTG_DELTA   29066
52 #define I_LFTG_K       80621
53 #define I_LFTG_X      106544
54
55 static inline void extend53(int *p, int i0, int i1)
56 {
57     p[i0 - 1] = p[i0 + 1];
58     p[i1]     = p[i1 - 2];
59     p[i0 - 2] = p[i0 + 2];
60     p[i1 + 1] = p[i1 - 3];
61 }
62
63 static inline void extend97_float(float *p, int i0, int i1)
64 {
65     int i;
66
67     for (i = 1; i <= 4; i++) {
68         p[i0 - i]     = p[i0 + i];
69         p[i1 + i - 1] = p[i1 - i - 1];
70     }
71 }
72
73 static inline void extend97_int(int32_t *p, int i0, int i1)
74 {
75     int i;
76
77     for (i = 1; i <= 4; i++) {
78         p[i0 - i]     = p[i0 + i];
79         p[i1 + i - 1] = p[i1 - i - 1];
80     }
81 }
82
83 static void sd_1d53(int *p, int i0, int i1)
84 {
85     int i;
86
87     if (i1 <= i0 + 1) {
88         if (i0 == 1)
89             p[1] <<= 1;
90         return;
91     }
92
93     extend53(p, i0, i1);
94
95     for (i = (i0+1)/2 - 1; i < (i1+1)/2; i++)
96         p[2*i+1] -= (p[2*i] + p[2*i+2]) >> 1;
97     for (i = (i0+1)/2; i < (i1+1)/2; i++)
98         p[2*i] += (p[2*i-1] + p[2*i+1] + 2) >> 2;
99 }
100
101 static void dwt_encode53(DWTContext *s, int *t)
102 {
103     int lev,
104         w = s->linelen[s->ndeclevels-1][0];
105     int *line = s->i_linebuf;
106     line += 3;
107
108     for (lev = s->ndeclevels-1; lev >= 0; lev--){
109         int lh = s->linelen[lev][0],
110             lv = s->linelen[lev][1],
111             mh = s->mod[lev][0],
112             mv = s->mod[lev][1],
113             lp;
114         int *l;
115
116         // VER_SD
117         l = line + mv;
118         for (lp = 0; lp < lh; lp++) {
119             int i, j = 0;
120
121             for (i = 0; i < lv; i++)
122                 l[i] = t[w*i + lp];
123
124             sd_1d53(line, mv, mv + lv);
125
126             // copy back and deinterleave
127             for (i =   mv; i < lv; i+=2, j++)
128                 t[w*j + lp] = l[i];
129             for (i = 1-mv; i < lv; i+=2, j++)
130                 t[w*j + lp] = l[i];
131         }
132
133         // HOR_SD
134         l = line + mh;
135         for (lp = 0; lp < lv; lp++){
136             int i, j = 0;
137
138             for (i = 0; i < lh; i++)
139                 l[i] = t[w*lp + i];
140
141             sd_1d53(line, mh, mh + lh);
142
143             // copy back and deinterleave
144             for (i =   mh; i < lh; i+=2, j++)
145                 t[w*lp + j] = l[i];
146             for (i = 1-mh; i < lh; i+=2, j++)
147                 t[w*lp + j] = l[i];
148         }
149     }
150 }
151 static void sd_1d97_float(float *p, int i0, int i1)
152 {
153     int i;
154
155     if (i1 <= i0 + 1) {
156         if (i0 == 1)
157             p[1] *= F_LFTG_X;
158         else
159             p[0] *= F_LFTG_K;
160         return;
161     }
162
163     extend97_float(p, i0, i1);
164     i0++; i1++;
165
166     for (i = i0/2 - 2; i < i1/2 + 1; i++)
167         p[2*i+1] -= 1.586134 * (p[2*i] + p[2*i+2]);
168     for (i = i0/2 - 1; i < i1/2 + 1; i++)
169         p[2*i] -= 0.052980 * (p[2*i-1] + p[2*i+1]);
170     for (i = i0/2 - 1; i < i1/2; i++)
171         p[2*i+1] += 0.882911 * (p[2*i] + p[2*i+2]);
172     for (i = i0/2; i < i1/2; i++)
173         p[2*i] += 0.443506 * (p[2*i-1] + p[2*i+1]);
174 }
175
176 static void dwt_encode97_float(DWTContext *s, float *t)
177 {
178     int lev,
179         w = s->linelen[s->ndeclevels-1][0];
180     float *line = s->f_linebuf;
181     line += 5;
182
183     for (lev = s->ndeclevels-1; lev >= 0; lev--){
184         int lh = s->linelen[lev][0],
185             lv = s->linelen[lev][1],
186             mh = s->mod[lev][0],
187             mv = s->mod[lev][1],
188             lp;
189         float *l;
190
191         // HOR_SD
192         l = line + mh;
193         for (lp = 0; lp < lv; lp++){
194             int i, j = 0;
195
196             for (i = 0; i < lh; i++)
197                 l[i] = t[w*lp + i];
198
199             sd_1d97_float(line, mh, mh + lh);
200
201             // copy back and deinterleave
202             for (i =   mh; i < lh; i+=2, j++)
203                 t[w*lp + j] = F_LFTG_X * l[i] / 2;
204             for (i = 1-mh; i < lh; i+=2, j++)
205                 t[w*lp + j] = F_LFTG_K * l[i] / 2;
206         }
207
208         // VER_SD
209         l = line + mv;
210         for (lp = 0; lp < lh; lp++) {
211             int i, j = 0;
212
213             for (i = 0; i < lv; i++)
214                 l[i] = t[w*i + lp];
215
216             sd_1d97_float(line, mv, mv + lv);
217
218             // copy back and deinterleave
219             for (i =   mv; i < lv; i+=2, j++)
220                 t[w*j + lp] = F_LFTG_X * l[i] / 2;
221             for (i = 1-mv; i < lv; i+=2, j++)
222                 t[w*j + lp] = F_LFTG_K * l[i] / 2;
223         }
224     }
225 }
226
227 static void sd_1d97_int(int *p, int i0, int i1)
228 {
229     int i;
230
231     if (i1 <= i0 + 1) {
232         if (i0 == 1)
233             p[1] = (p[1] * I_LFTG_X + (1<<15)) >> 16;
234         else
235             p[0] = (p[0] * I_LFTG_K + (1<<15)) >> 16;
236         return;
237     }
238
239     extend97_int(p, i0, i1);
240     i0++; i1++;
241
242     for (i = i0/2 - 2; i < i1/2 + 1; i++)
243         p[2 * i + 1] -= (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
244     for (i = i0/2 - 1; i < i1/2 + 1; i++)
245         p[2 * i]     -= (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
246     for (i = i0/2 - 1; i < i1/2; i++)
247         p[2 * i + 1] += (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
248     for (i = i0/2; i < i1/2; i++)
249         p[2 * i]     += (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
250 }
251
252 static void dwt_encode97_int(DWTContext *s, int *t)
253 {
254     int lev,
255         w = s->linelen[s->ndeclevels-1][0];
256     int *line = s->i_linebuf;
257     line += 5;
258
259     for (lev = s->ndeclevels-1; lev >= 0; lev--){
260         int lh = s->linelen[lev][0],
261             lv = s->linelen[lev][1],
262             mh = s->mod[lev][0],
263             mv = s->mod[lev][1],
264             lp;
265         int *l;
266
267         // VER_SD
268         l = line + mv;
269         for (lp = 0; lp < lh; lp++) {
270             int i, j = 0;
271
272             for (i = 0; i < lv; i++)
273                 l[i] = t[w*i + lp];
274
275             sd_1d97_int(line, mv, mv + lv);
276
277             // copy back and deinterleave
278             for (i =   mv; i < lv; i+=2, j++)
279                 t[w*j + lp] = ((l[i] * I_LFTG_X) + (1 << 16)) >> 17;
280             for (i = 1-mv; i < lv; i+=2, j++)
281                 t[w*j + lp] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
282         }
283
284         // HOR_SD
285         l = line + mh;
286         for (lp = 0; lp < lv; lp++){
287             int i, j = 0;
288
289             for (i = 0; i < lh; i++)
290                 l[i] = t[w*lp + i];
291
292             sd_1d97_int(line, mh, mh + lh);
293
294             // copy back and deinterleave
295             for (i =   mh; i < lh; i+=2, j++)
296                 t[w*lp + j] = ((l[i] * I_LFTG_X) + (1 << 16)) >> 17;
297             for (i = 1-mh; i < lh; i+=2, j++)
298                 t[w*lp + j] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
299         }
300
301     }
302 }
303
304 static void sr_1d53(int *p, int i0, int i1)
305 {
306     int i;
307
308     if (i1 <= i0 + 1) {
309         if (i0 == 1)
310             p[1] >>= 1;
311         return;
312     }
313
314     extend53(p, i0, i1);
315
316     for (i = i0 / 2; i < i1 / 2 + 1; i++)
317         p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
318     for (i = i0 / 2; i < i1 / 2; i++)
319         p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
320 }
321
322 static void dwt_decode53(DWTContext *s, int *t)
323 {
324     int lev;
325     int w     = s->linelen[s->ndeclevels - 1][0];
326     int32_t *line = s->i_linebuf;
327     line += 3;
328
329     for (lev = 0; lev < s->ndeclevels; lev++) {
330         int lh = s->linelen[lev][0],
331             lv = s->linelen[lev][1],
332             mh = s->mod[lev][0],
333             mv = s->mod[lev][1],
334             lp;
335         int *l;
336
337         // HOR_SD
338         l = line + mh;
339         for (lp = 0; lp < lv; lp++) {
340             int i, j = 0;
341             // copy with interleaving
342             for (i = mh; i < lh; i += 2, j++)
343                 l[i] = t[w * lp + j];
344             for (i = 1 - mh; i < lh; i += 2, j++)
345                 l[i] = t[w * lp + j];
346
347             sr_1d53(line, mh, mh + lh);
348
349             for (i = 0; i < lh; i++)
350                 t[w * lp + i] = l[i];
351         }
352
353         // VER_SD
354         l = line + mv;
355         for (lp = 0; lp < lh; lp++) {
356             int i, j = 0;
357             // copy with interleaving
358             for (i = mv; i < lv; i += 2, j++)
359                 l[i] = t[w * j + lp];
360             for (i = 1 - mv; i < lv; i += 2, j++)
361                 l[i] = t[w * j + lp];
362
363             sr_1d53(line, mv, mv + lv);
364
365             for (i = 0; i < lv; i++)
366                 t[w * i + lp] = l[i];
367         }
368     }
369 }
370
371 static void sr_1d97_float(float *p, int i0, int i1)
372 {
373     int i;
374
375     if (i1 <= i0 + 1) {
376         if (i0 == 1)
377             p[1] *= F_LFTG_K/2;
378         else
379             p[0] *= F_LFTG_X/2;
380         return;
381     }
382
383     extend97_float(p, i0, i1);
384
385     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
386         p[2 * i]     -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
387     /* step 4 */
388     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
389         p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]);
390     /*step 5*/
391     for (i = i0 / 2; i < i1 / 2 + 1; i++)
392         p[2 * i]     += F_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]);
393     /* step 6 */
394     for (i = i0 / 2; i < i1 / 2; i++)
395         p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]);
396 }
397
398 static void dwt_decode97_float(DWTContext *s, float *t)
399 {
400     int lev;
401     int w       = s->linelen[s->ndeclevels - 1][0];
402     float *line = s->f_linebuf;
403     float *data = t;
404     /* position at index O of line range [0-5,w+5] cf. extend function */
405     line += 5;
406
407     for (lev = 0; lev < s->ndeclevels; lev++) {
408         int lh = s->linelen[lev][0],
409             lv = s->linelen[lev][1],
410             mh = s->mod[lev][0],
411             mv = s->mod[lev][1],
412             lp;
413         float *l;
414         // HOR_SD
415         l = line + mh;
416         for (lp = 0; lp < lv; lp++) {
417             int i, j = 0;
418             // copy with interleaving
419             for (i = mh; i < lh; i += 2, j++)
420                 l[i] = data[w * lp + j] * F_LFTG_K;
421             for (i = 1 - mh; i < lh; i += 2, j++)
422                 l[i] = data[w * lp + j] * F_LFTG_X;
423
424             sr_1d97_float(line, mh, mh + lh);
425
426             for (i = 0; i < lh; i++)
427                 data[w * lp + i] = l[i];
428         }
429
430         // VER_SD
431         l = line + mv;
432         for (lp = 0; lp < lh; lp++) {
433             int i, j = 0;
434             // copy with interleaving
435             for (i = mv; i < lv; i += 2, j++)
436                 l[i] = data[w * j + lp] * F_LFTG_K;
437             for (i = 1 - mv; i < lv; i += 2, j++)
438                 l[i] = data[w * j + lp] * F_LFTG_X;
439
440             sr_1d97_float(line, mv, mv + lv);
441
442             for (i = 0; i < lv; i++)
443                 data[w * i + lp] = l[i];
444         }
445     }
446 }
447
448 static void sr_1d97_int(int32_t *p, int i0, int i1)
449 {
450     int i;
451
452     if (i1 <= i0 + 1) {
453         if (i0 == 1)
454             p[1] = (p[1] * I_LFTG_K + (1<<16)) >> 17;
455         else
456             p[0] = (p[0] * I_LFTG_X + (1<<16)) >> 17;
457         return;
458     }
459
460     extend97_int(p, i0, i1);
461
462     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
463         p[2 * i]     -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
464     /* step 4 */
465     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
466         p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
467     /*step 5*/
468     for (i = i0 / 2; i < i1 / 2 + 1; i++)
469         p[2 * i]     += (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
470     /* step 6 */
471     for (i = i0 / 2; i < i1 / 2; i++)
472         p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
473 }
474
475 static void dwt_decode97_int(DWTContext *s, int32_t *t)
476 {
477     int lev;
478     int w       = s->linelen[s->ndeclevels - 1][0];
479     int32_t *line = s->i_linebuf;
480     int32_t *data = t;
481     /* position at index O of line range [0-5,w+5] cf. extend function */
482     line += 5;
483
484     for (lev = 0; lev < s->ndeclevels; lev++) {
485         int lh = s->linelen[lev][0],
486             lv = s->linelen[lev][1],
487             mh = s->mod[lev][0],
488             mv = s->mod[lev][1],
489             lp;
490         int32_t *l;
491         // HOR_SD
492         l = line + mh;
493         for (lp = 0; lp < lv; lp++) {
494             int i, j = 0;
495             // rescale with interleaving
496             for (i = mh; i < lh; i += 2, j++)
497                 l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
498             for (i = 1 - mh; i < lh; i += 2, j++)
499                 l[i] = ((data[w * lp + j] * I_LFTG_X) + (1 << 15)) >> 16;
500
501             sr_1d97_int(line, mh, mh + lh);
502
503             for (i = 0; i < lh; i++)
504                 data[w * lp + i] = l[i];
505         }
506
507         // VER_SD
508         l = line + mv;
509         for (lp = 0; lp < lh; lp++) {
510             int i, j = 0;
511             // rescale with interleaving
512             for (i = mv; i < lv; i += 2, j++)
513                 l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
514             for (i = 1 - mv; i < lv; i += 2, j++)
515                 l[i] = ((data[w * j + lp] * I_LFTG_X) + (1 << 15)) >> 16;
516
517             sr_1d97_int(line, mv, mv + lv);
518
519             for (i = 0; i < lv; i++)
520                 data[w * i + lp] = l[i];
521         }
522     }
523 }
524
525 int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
526                          int decomp_levels, int type)
527 {
528     int i, j, lev = decomp_levels, maxlen,
529         b[2][2];
530
531     s->ndeclevels = decomp_levels;
532     s->type       = type;
533
534     for (i = 0; i < 2; i++)
535         for (j = 0; j < 2; j++)
536             b[i][j] = border[i][j];
537
538     maxlen = FFMAX(b[0][1] - b[0][0],
539                    b[1][1] - b[1][0]);
540     while (--lev >= 0)
541         for (i = 0; i < 2; i++) {
542             s->linelen[lev][i] = b[i][1] - b[i][0];
543             s->mod[lev][i]     = b[i][0] & 1;
544             for (j = 0; j < 2; j++)
545                 b[i][j] = (b[i][j] + 1) >> 1;
546         }
547     switch (type) {
548     case FF_DWT97:
549         s->f_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->f_linebuf));
550         if (!s->f_linebuf)
551             return AVERROR(ENOMEM);
552         break;
553      case FF_DWT97_INT:
554         s->i_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->i_linebuf));
555         if (!s->i_linebuf)
556             return AVERROR(ENOMEM);
557         break;
558     case FF_DWT53:
559         s->i_linebuf = av_malloc_array((maxlen +  6), sizeof(*s->i_linebuf));
560         if (!s->i_linebuf)
561             return AVERROR(ENOMEM);
562         break;
563     default:
564         return -1;
565     }
566     return 0;
567 }
568
569 int ff_dwt_encode(DWTContext *s, void *t)
570 {
571     switch(s->type){
572         case FF_DWT97:
573             dwt_encode97_float(s, t); break;
574         case FF_DWT97_INT:
575             dwt_encode97_int(s, t); break;
576         case FF_DWT53:
577             dwt_encode53(s, t); break;
578         default:
579             return -1;
580     }
581     return 0;
582 }
583
584 int ff_dwt_decode(DWTContext *s, void *t)
585 {
586     switch (s->type) {
587     case FF_DWT97:
588         dwt_decode97_float(s, t);
589         break;
590     case FF_DWT97_INT:
591         dwt_decode97_int(s, t);
592         break;
593     case FF_DWT53:
594         dwt_decode53(s, t);
595         break;
596     default:
597         return -1;
598     }
599     return 0;
600 }
601
602 void ff_dwt_destroy(DWTContext *s)
603 {
604     av_freep(&s->f_linebuf);
605     av_freep(&s->i_linebuf);
606 }
607
608 #ifdef TEST
609
610 #include "libavutil/lfg.h"
611
612 #define MAX_W 256
613
614 static int test_dwt(int *array, int *ref, uint16_t border[2][2], int decomp_levels, int type, int max_diff) {
615     int ret, j;
616     DWTContext s1={{{0}}}, *s= &s1;
617     int64_t err2 = 0;
618
619     ret = ff_jpeg2000_dwt_init(s,  border, decomp_levels, type);
620     if (ret < 0) {
621         fprintf(stderr, "ff_jpeg2000_dwt_init failed\n");
622         return 1;
623     }
624     ret = ff_dwt_encode(s, array);
625     if (ret < 0) {
626         fprintf(stderr, "ff_dwt_encode failed\n");
627         return 1;
628     }
629     ret = ff_dwt_decode(s, array);
630     if (ret < 0) {
631         fprintf(stderr, "ff_dwt_encode failed\n");
632         return 1;
633     }
634     for (j = 0; j<MAX_W * MAX_W; j++) {
635         if (FFABS(array[j] - ref[j]) > max_diff) {
636             fprintf(stderr, "missmatch at %d (%d != %d) decomp:%d border %d %d %d %d\n",
637                     j, array[j], ref[j],decomp_levels, border[0][0], border[0][1], border[1][0], border[1][1]);
638             return 2;
639         }
640         err2 += (array[j] - ref[j]) * (array[j] - ref[j]);
641         array[j] = ref[j];
642     }
643     ff_dwt_destroy(s);
644
645     printf("%s, decomp:%2d border %3d %3d %3d %3d milli-err2:%9"PRId64"\n",
646            type == FF_DWT53 ? "5/3i" : "9/7i",
647            decomp_levels, border[0][0], border[0][1], border[1][0], border[1][1],
648            1000*err2 / ((border[0][1] - border[0][0])*(border[1][1] - border[1][0])));
649
650     return 0;
651 }
652
653 static int test_dwtf(float *array, float *ref, uint16_t border[2][2], int decomp_levels, float max_diff) {
654     int ret, j;
655     DWTContext s1={{{0}}}, *s= &s1;
656     double err2 = 0;
657
658     ret = ff_jpeg2000_dwt_init(s,  border, decomp_levels, FF_DWT97);
659     if (ret < 0) {
660         fprintf(stderr, "ff_jpeg2000_dwt_init failed\n");
661         return 1;
662     }
663     ret = ff_dwt_encode(s, array);
664     if (ret < 0) {
665         fprintf(stderr, "ff_dwt_encode failed\n");
666         return 1;
667     }
668     ret = ff_dwt_decode(s, array);
669     if (ret < 0) {
670         fprintf(stderr, "ff_dwt_encode failed\n");
671         return 1;
672     }
673     for (j = 0; j<MAX_W * MAX_W; j++) {
674         if (FFABS(array[j] - ref[j]) > max_diff) {
675             fprintf(stderr, "missmatch at %d (%f != %f) decomp:%d border %d %d %d %d\n",
676                     j, array[j], ref[j],decomp_levels, border[0][0], border[0][1], border[1][0], border[1][1]);
677             return 2;
678         }
679         err2 += (array[j] - ref[j]) * (array[j] - ref[j]);
680         array[j] = ref[j];
681     }
682     ff_dwt_destroy(s);
683
684     printf("9/7f, decomp:%2d border %3d %3d %3d %3d err2:%20.4f\n",
685            decomp_levels, border[0][0], border[0][1], border[1][0], border[1][1],
686            err2 / ((border[0][1] - border[0][0])*(border[1][1] - border[1][0])));
687
688     return 0;
689 }
690
691 int main(void) {
692     int array[MAX_W * MAX_W];
693     int ref  [MAX_W * MAX_W];
694     float arrayf[MAX_W * MAX_W];
695     float reff  [MAX_W * MAX_W];
696     AVLFG prng;
697     int i,j;
698     uint16_t border[2][2];
699     int ret, decomp_levels;
700
701     av_lfg_init(&prng, 1);
702
703     for (i = 0; i<MAX_W * MAX_W; i++)
704         arrayf[i] = reff[i] = array[i] = ref[i] =  av_lfg_get(&prng) % 2048;
705
706     for (i = 0; i < 100; i++) {
707         for (j=0; j<4; j++)
708             border[j>>1][j&1] = av_lfg_get(&prng) % MAX_W;
709         if (border[0][0] >= border[0][1] || border[1][0] >= border[1][1])
710             continue;
711         decomp_levels = av_lfg_get(&prng) % FF_DWT_MAX_DECLVLS;
712
713         ret = test_dwt(array, ref, border, decomp_levels, FF_DWT53, 0);
714         if (ret)
715             return ret;
716         ret = test_dwt(array, ref, border, decomp_levels, FF_DWT97_INT, FFMIN(7+5*decomp_levels, 15+3*decomp_levels));
717         if (ret)
718             return ret;
719         ret = test_dwtf(arrayf, reff, border, decomp_levels, 1.0);
720         if (ret)
721             return ret;
722     }
723
724     return 0;
725 }
726
727 #endif