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