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