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
  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] = 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] = 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];
 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];
 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];
 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];
 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];
 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];
 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