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