git.sesse.net Git - ffmpeg/blob - 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         av_assert1(!mh && !mv);
 192         // HOR_SD
 193         l = line + mh;
 194         for (lp = 0; lp < lv; lp++){
 195             int i, j = 0;
 196
 197             for (i = 0; i < lh; i++)
 198                 l[i] = t[w*lp + i];
 199
 200             sd_1d97_float(line, mh, mh + lh);
 201
 202             // copy back and deinterleave
 203             for (i =   mh; i < lh; i+=2, j++)
 204                 t[w*lp + j] = F_LFTG_X * l[i] / 2;
 205             for (i = 1-mh; i < lh; i+=2, j++)
 206                 t[w*lp + j] = F_LFTG_K * l[i] / 2;
 207         }
 208
 209         // VER_SD
 210         l = line + mv;
 211         for (lp = 0; lp < lh; lp++) {
 212             int i, j = 0;
 213
 214             for (i = 0; i < lv; i++)
 215                 l[i] = t[w*i + lp];
 216
 217             sd_1d97_float(line, mv, mv + lv);
 218
 219             // copy back and deinterleave
 220             for (i =   mv; i < lv; i+=2, j++)
 221                 t[w*j + lp] = F_LFTG_X * l[i] / 2;
 222             for (i = 1-mv; i < lv; i+=2, j++)
 223                 t[w*j + lp] = F_LFTG_K * l[i] / 2;
 224         }
 225     }
 226 }
 227
 228 static void sd_1d97_int(int *p, int i0, int i1)
 229 {
 230     int i;
 231
 232     if (i1 <= i0 + 1) {
 233         if (i0 == 1)
 234             p[1] = (p[1] * I_LFTG_X + (1<<15)) >> 16;
 235         else
 236             p[0] = (p[0] * I_LFTG_K + (1<<15)) >> 16;
 237         return;
 238     }
 239
 240     extend97_int(p, i0, i1);
 241     i0++; i1++;
 242
 243     for (i = i0/2 - 2; i < i1/2 + 1; i++)
 244         p[2 * i + 1] -= (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 245     for (i = i0/2 - 1; i < i1/2 + 1; i++)
 246         p[2 * i]     -= (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 247     for (i = i0/2 - 1; i < i1/2; i++)
 248         p[2 * i + 1] += (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 249     for (i = i0/2; i < i1/2; i++)
 250         p[2 * i]     += (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 251 }
 252
 253 static void dwt_encode97_int(DWTContext *s, int *t)
 254 {
 255     int lev,
 256         w = s->linelen[s->ndeclevels-1][0];
 257     int *line = s->i_linebuf;
 258     line += 5;
 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 << 16)) >> 17;
 281             for (i = 1-mv; i < lv; i+=2, j++)
 282                 t[w*j + lp] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
 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 << 16)) >> 17;
 298             for (i = 1-mh; i < lh; i+=2, j++)
 299                 t[w*lp + j] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
 300         }
 301
 302     }
 303 }
 304
 305 static void sr_1d53(int *p, int i0, int i1)
 306 {
 307     int i;
 308
 309     if (i1 <= i0 + 1) {
 310         if (i0 == 1)
 311             p[1] >>= 1;
 312         return;
 313     }
 314
 315     extend53(p, i0, i1);
 316
 317     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 318         p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
 319     for (i = i0 / 2; i < i1 / 2; i++)
 320         p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
 321 }
 322
 323 static void dwt_decode53(DWTContext *s, int *t)
 324 {
 325     int lev;
 326     int w     = s->linelen[s->ndeclevels - 1][0];
 327     int32_t *line = s->i_linebuf;
 328     line += 3;
 329
 330     for (lev = 0; lev < s->ndeclevels; lev++) {
 331         int lh = s->linelen[lev][0],
 332             lv = s->linelen[lev][1],
 333             mh = s->mod[lev][0],
 334             mv = s->mod[lev][1],
 335             lp;
 336         int *l;
 337
 338         // HOR_SD
 339         l = line + mh;
 340         for (lp = 0; lp < lv; lp++) {
 341             int i, j = 0;
 342             // copy with interleaving
 343             for (i = mh; i < lh; i += 2, j++)
 344                 l[i] = t[w * lp + j];
 345             for (i = 1 - mh; i < lh; i += 2, j++)
 346                 l[i] = t[w * lp + j];
 347
 348             sr_1d53(line, mh, mh + lh);
 349
 350             for (i = 0; i < lh; i++)
 351                 t[w * lp + i] = l[i];
 352         }
 353
 354         // VER_SD
 355         l = line + mv;
 356         for (lp = 0; lp < lh; lp++) {
 357             int i, j = 0;
 358             // copy with interleaving
 359             for (i = mv; i < lv; i += 2, j++)
 360                 l[i] = t[w * j + lp];
 361             for (i = 1 - mv; i < lv; i += 2, j++)
 362                 l[i] = t[w * j + lp];
 363
 364             sr_1d53(line, mv, mv + lv);
 365
 366             for (i = 0; i < lv; i++)
 367                 t[w * i + lp] = l[i];
 368         }
 369     }
 370 }
 371
 372 static void sr_1d97_float(float *p, int i0, int i1)
 373 {
 374     int i;
 375
 376     if (i1 <= i0 + 1) {
 377         if (i0 == 1)
 378             p[1] *= F_LFTG_K/2;
 379         else
 380             p[0] *= F_LFTG_X/2;
 381         return;
 382     }
 383
 384     extend97_float(p, i0, i1);
 385
 386     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 387         p[2 * i]     -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
 388     /* step 4 */
 389     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 390         p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]);
 391     /*step 5*/
 392     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 393         p[2 * i]     += F_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]);
 394     /* step 6 */
 395     for (i = i0 / 2; i < i1 / 2; i++)
 396         p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]);
 397 }
 398
 399 static void dwt_decode97_float(DWTContext *s, float *t)
 400 {
 401     int lev;
 402     int w       = s->linelen[s->ndeclevels - 1][0];
 403     float *line = s->f_linebuf;
 404     float *data = t;
 405     /* position at index O of line range [0-5,w+5] cf. extend function */
 406     line += 5;
 407
 408     for (lev = 0; lev < s->ndeclevels; lev++) {
 409         int lh = s->linelen[lev][0],
 410             lv = s->linelen[lev][1],
 411             mh = s->mod[lev][0],
 412             mv = s->mod[lev][1],
 413             lp;
 414         float *l;
 415         // HOR_SD
 416         l = line + mh;
 417         for (lp = 0; lp < lv; lp++) {
 418             int i, j = 0;
 419             // copy with interleaving
 420             for (i = mh; i < lh; i += 2, j++)
 421                 l[i] = data[w * lp + j] * F_LFTG_K;
 422             for (i = 1 - mh; i < lh; i += 2, j++)
 423                 l[i] = data[w * lp + j] * F_LFTG_X;
 424
 425             sr_1d97_float(line, mh, mh + lh);
 426
 427             for (i = 0; i < lh; i++)
 428                 data[w * lp + i] = l[i];
 429         }
 430
 431         // VER_SD
 432         l = line + mv;
 433         for (lp = 0; lp < lh; lp++) {
 434             int i, j = 0;
 435             // copy with interleaving
 436             for (i = mv; i < lv; i += 2, j++)
 437                 l[i] = data[w * j + lp] * F_LFTG_K;
 438             for (i = 1 - mv; i < lv; i += 2, j++)
 439                 l[i] = data[w * j + lp] * F_LFTG_X;
 440
 441             sr_1d97_float(line, mv, mv + lv);
 442
 443             for (i = 0; i < lv; i++)
 444                 data[w * i + lp] = l[i];
 445         }
 446     }
 447 }
 448
 449 static void sr_1d97_int(int32_t *p, int i0, int i1)
 450 {
 451     int i;
 452
 453     if (i1 <= i0 + 1) {
 454         if (i0 == 1)
 455             p[1] = (p[1] * I_LFTG_K + (1<<16)) >> 17;
 456         else
 457             p[0] = (p[0] * I_LFTG_X + (1<<16)) >> 17;
 458         return;
 459     }
 460
 461     extend97_int(p, i0, i1);
 462
 463     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 464         p[2 * i]     -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 465     /* step 4 */
 466     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 467         p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 468     /*step 5*/
 469     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 470         p[2 * i]     += (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 471     /* step 6 */
 472     for (i = i0 / 2; i < i1 / 2; i++)
 473         p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 474 }
 475
 476 static void dwt_decode97_int(DWTContext *s, int32_t *t)
 477 {
 478     int lev;
 479     int w       = s->linelen[s->ndeclevels - 1][0];
 480     int32_t *line = s->i_linebuf;
 481     int32_t *data = t;
 482     /* position at index O of line range [0-5,w+5] cf. extend function */
 483     line += 5;
 484
 485     for (lev = 0; lev < s->ndeclevels; lev++) {
 486         int lh = s->linelen[lev][0],
 487             lv = s->linelen[lev][1],
 488             mh = s->mod[lev][0],
 489             mv = s->mod[lev][1],
 490             lp;
 491         int32_t *l;
 492         // HOR_SD
 493         l = line + mh;
 494         for (lp = 0; lp < lv; lp++) {
 495             int i, j = 0;
 496             // rescale with interleaving
 497             for (i = mh; i < lh; i += 2, j++)
 498                 l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
 499             for (i = 1 - mh; i < lh; i += 2, j++)
 500                 l[i] = ((data[w * lp + j] * I_LFTG_X) + (1 << 15)) >> 16;
 501
 502             sr_1d97_int(line, mh, mh + lh);
 503
 504             for (i = 0; i < lh; i++)
 505                 data[w * lp + i] = l[i];
 506         }
 507
 508         // VER_SD
 509         l = line + mv;
 510         for (lp = 0; lp < lh; lp++) {
 511             int i, j = 0;
 512             // rescale with interleaving
 513             for (i = mv; i < lv; i += 2, j++)
 514                 l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
 515             for (i = 1 - mv; i < lv; i += 2, j++)
 516                 l[i] = ((data[w * j + lp] * I_LFTG_X) + (1 << 15)) >> 16;
 517
 518             sr_1d97_int(line, mv, mv + lv);
 519
 520             for (i = 0; i < lv; i++)
 521                 data[w * i + lp] = l[i];
 522         }
 523     }
 524 }
 525
 526 int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
 527                          int decomp_levels, int type)
 528 {
 529     int i, j, lev = decomp_levels, maxlen,
 530         b[2][2];
 531
 532     s->ndeclevels = decomp_levels;
 533     s->type       = type;
 534
 535     for (i = 0; i < 2; i++)
 536         for (j = 0; j < 2; j++)
 537             b[i][j] = border[i][j];
 538
 539     maxlen = FFMAX(b[0][1] - b[0][0],
 540                    b[1][1] - b[1][0]);
 541     while (--lev >= 0)
 542         for (i = 0; i < 2; i++) {
 543             s->linelen[lev][i] = b[i][1] - b[i][0];
 544             s->mod[lev][i]     = b[i][0] & 1;
 545             for (j = 0; j < 2; j++)
 546                 b[i][j] = (b[i][j] + 1) >> 1;
 547         }
 548     switch (type) {
 549     case FF_DWT97:
 550         s->f_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->f_linebuf));
 551         if (!s->f_linebuf)
 552             return AVERROR(ENOMEM);
 553         break;
 554      case FF_DWT97_INT:
 555         s->i_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->i_linebuf));
 556         if (!s->i_linebuf)
 557             return AVERROR(ENOMEM);
 558         break;
 559     case FF_DWT53:
 560         s->i_linebuf = av_malloc_array((maxlen +  6), sizeof(*s->i_linebuf));
 561         if (!s->i_linebuf)
 562             return AVERROR(ENOMEM);
 563         break;
 564     default:
 565         return -1;
 566     }
 567     return 0;
 568 }
 569
 570 int ff_dwt_encode(DWTContext *s, void *t)
 571 {
 572     switch(s->type){
 573         case FF_DWT97:
 574             dwt_encode97_float(s, t); break;
 575         case FF_DWT97_INT:
 576             dwt_encode97_int(s, t); break;
 577         case FF_DWT53:
 578             dwt_encode53(s, t); break;
 579         default:
 580             return -1;
 581     }
 582     return 0;
 583 }
 584
 585 int ff_dwt_decode(DWTContext *s, void *t)
 586 {
 587     switch (s->type) {
 588     case FF_DWT97:
 589         dwt_decode97_float(s, t);
 590         break;
 591     case FF_DWT97_INT:
 592         dwt_decode97_int(s, t);
 593         break;
 594     case FF_DWT53:
 595         dwt_decode53(s, t);
 596         break;
 597     default:
 598         return -1;
 599     }
 600     return 0;
 601 }
 602
 603 void ff_dwt_destroy(DWTContext *s)
 604 {
 605     av_freep(&s->f_linebuf);
 606     av_freep(&s->i_linebuf);
 607 }
 608
 609 #ifdef TEST
 610
 611 #include "libavutil/lfg.h"
 612
 613 #define MAX_W 256
 614
 615 static int test_dwt(int *array, int *ref, uint16_t border[2][2], int decomp_levels, int type, int max_diff) {
 616     int ret, j;
 617     DWTContext s1={{{0}}}, *s= &s1;
 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         array[j] = ref[j];
 641     }
 642     ff_dwt_destroy(s);
 643
 644     return 0;
 645 }
 646
 647 int main(void) {
 648     int array[MAX_W * MAX_W];
 649     int ref  [MAX_W * MAX_W];
 650     AVLFG prng;
 651     int i,j;
 652     uint16_t border[2][2];
 653     int ret, decomp_levels;
 654
 655     av_lfg_init(&prng, 1);
 656
 657     for (i = 0; i<MAX_W * MAX_W; i++)
 658         array[i] = ref[i] =  av_lfg_get(&prng) % 2048;
 659
 660     for (i = 0; i < 40; i++) {
 661         for (j=0; j<4; j++)
 662             border[j>>1][j&1] = av_lfg_get(&prng) % MAX_W;
 663         if (border[0][0] >= border[0][1] || border[1][0] >= border[1][1])
 664             continue;
 665         decomp_levels = av_lfg_get(&prng) % FF_DWT_MAX_DECLVLS;
 666
 667         ret = test_dwt(array, ref, border, decomp_levels, FF_DWT53, 0);
 668         if (ret)
 669             return ret;
 670         ret = test_dwt(array, ref, border, decomp_levels, FF_DWT97_INT, FFMIN(7+5*decomp_levels, 15+3*decomp_levels));
 671         if (ret)
 672             return ret;
 673     }
 674
 675     return 0;
 676 }
 677
 678 #endif