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  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] *= 1 << 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, int 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     if (s->ndeclevels == 0)
 584         return 0;
 585
 586     switch(s->type){
 587         case FF_DWT97:
 588             dwt_encode97_float(s, t); break;
 589         case FF_DWT97_INT:
 590             dwt_encode97_int(s, t); break;
 591         case FF_DWT53:
 592             dwt_encode53(s, t); break;
 593         default:
 594             return -1;
 595     }
 596     return 0;
 597 }
 598
 599 int ff_dwt_decode(DWTContext *s, void *t)
 600 {
 601     if (s->ndeclevels == 0)
 602         return 0;
 603
 604     switch (s->type) {
 605     case FF_DWT97:
 606         dwt_decode97_float(s, t);
 607         break;
 608     case FF_DWT97_INT:
 609         dwt_decode97_int(s, t);
 610         break;
 611     case FF_DWT53:
 612         dwt_decode53(s, t);
 613         break;
 614     default:
 615         return -1;
 616     }
 617     return 0;
 618 }
 619
 620 void ff_dwt_destroy(DWTContext *s)
 621 {
 622     av_freep(&s->f_linebuf);
 623     av_freep(&s->i_linebuf);
 624 }