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/common.h"
  29 #include "libavutil/mem.h"
  30 #include "jpeg2000dwt.h"
  31 #include "internal.h"
  32
  33 /* Defines for 9/7 DWT lifting parameters.
  34  * Parameters are in float. */
  35 #define F_LFTG_ALPHA  1.586134342059924f
  36 #define F_LFTG_BETA   0.052980118572961f
  37 #define F_LFTG_GAMMA  0.882911075530934f
  38 #define F_LFTG_DELTA  0.443506852043971f
  39 #define F_LFTG_K      1.230174104914001f
  40 #define F_LFTG_X      1.625732422f
  41 /* FIXME: Why use 1.625732422 instead of 1/F_LFTG_K?
  42  * Incorrect value in JPEG2000 norm.
  43  * see (ISO/IEC 15444:1 (version 2002) F.3.8.2 */
  44
  45 /* Lifting parameters in integer format.
  46  * Computed as param = (float param) * (1 << 16) */
  47 #define I_LFTG_ALPHA  103949
  48 #define I_LFTG_BETA     3472
  49 #define I_LFTG_GAMMA   57862
  50 #define I_LFTG_DELTA   29066
  51 #define I_LFTG_K       80621
  52 #define I_LFTG_X      106544
  53
  54 static inline void extend53(int *p, int i0, int i1)
  55 {
  56     p[i0 - 1] = p[i0 + 1];
  57     p[i1]     = p[i1 - 2];
  58     p[i0 - 2] = p[i0 + 2];
  59     p[i1 + 1] = p[i1 - 3];
  60 }
  61
  62 static inline void extend97_float(float *p, int i0, int i1)
  63 {
  64     int i;
  65
  66     for (i = 1; i <= 4; i++) {
  67         p[i0 - i]     = p[i0 + i];
  68         p[i1 + i - 1] = p[i1 - i - 1];
  69     }
  70 }
  71
  72 static inline void extend97_int(int32_t *p, int i0, int i1)
  73 {
  74     int i;
  75
  76     for (i = 1; i <= 4; i++) {
  77         p[i0 - i]     = p[i0 + i];
  78         p[i1 + i - 1] = p[i1 - i - 1];
  79     }
  80 }
  81
  82 static void sd_1d53(int *p, int i0, int i1)
  83 {
  84     int i;
  85
  86     if (i1 == i0 + 1)
  87         return;
  88
  89     extend53(p, i0, i1);
  90
  91     for (i = (i0+1)/2 - 1; i < (i1+1)/2; i++)
  92         p[2*i+1] -= (p[2*i] + p[2*i+2]) >> 1;
  93     for (i = (i0+1)/2; i < (i1+1)/2; 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         // HOR_SD
 113         l = line + mh;
 114         for (lp = 0; lp < lv; lp++){
 115             int i, j = 0;
 116
 117             for (i = 0; i < lh; i++)
 118                 l[i] = t[w*lp + i];
 119
 120             sd_1d53(line, mh, mh + lh);
 121
 122             // copy back and deinterleave
 123             for (i =   mh; i < lh; i+=2, j++)
 124                 t[w*lp + j] = l[i];
 125             for (i = 1-mh; i < lh; i+=2, j++)
 126                 t[w*lp + j] = l[i];
 127         }
 128
 129         // VER_SD
 130         l = line + mv;
 131         for (lp = 0; lp < lh; lp++) {
 132             int i, j = 0;
 133
 134             for (i = 0; i < lv; i++)
 135                 l[i] = t[w*i + lp];
 136
 137             sd_1d53(line, mv, mv + lv);
 138
 139             // copy back and deinterleave
 140             for (i =   mv; i < lv; i+=2, j++)
 141                 t[w*j + lp] = l[i];
 142             for (i = 1-mv; i < lv; i+=2, j++)
 143                 t[w*j + lp] = 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         return;
 153
 154     extend97_float(p, i0, i1);
 155     i0++; i1++;
 156
 157     for (i = i0/2 - 2; i < i1/2 + 1; i++)
 158         p[2*i+1] -= 1.586134 * (p[2*i] + p[2*i+2]);
 159     for (i = i0/2 - 1; i < i1/2 + 1; i++)
 160         p[2*i] -= 0.052980 * (p[2*i-1] + p[2*i+1]);
 161     for (i = i0/2 - 1; i < i1/2; i++)
 162         p[2*i+1] += 0.882911 * (p[2*i] + p[2*i+2]);
 163     for (i = i0/2; i < i1/2; i++)
 164         p[2*i] += 0.443506 * (p[2*i-1] + p[2*i+1]);
 165 }
 166
 167 static void dwt_encode97_float(DWTContext *s, float *t)
 168 {
 169     int lev,
 170         w = s->linelen[s->ndeclevels-1][0];
 171     float *line = s->f_linebuf;
 172     line += 5;
 173
 174     for (lev = s->ndeclevels-1; lev >= 0; lev--){
 175         int lh = s->linelen[lev][0],
 176             lv = s->linelen[lev][1],
 177             mh = s->mod[lev][0],
 178             mv = s->mod[lev][1],
 179             lp;
 180         float *l;
 181
 182         // HOR_SD
 183         l = line + mh;
 184         for (lp = 0; lp < lv; lp++){
 185             int i, j = 0;
 186
 187             for (i = 0; i < lh; i++)
 188                 l[i] = t[w*lp + i];
 189
 190             sd_1d97_float(line, mh, mh + lh);
 191
 192             // copy back and deinterleave
 193             for (i =   mh; i < lh; i+=2, j++)
 194                 t[w*lp + j] = F_LFTG_X * l[i] / 2;
 195             for (i = 1-mh; i < lh; i+=2, j++)
 196                 t[w*lp + j] = F_LFTG_K * l[i] / 2;
 197         }
 198
 199         // VER_SD
 200         l = line + mv;
 201         for (lp = 0; lp < lh; lp++) {
 202             int i, j = 0;
 203
 204             for (i = 0; i < lv; i++)
 205                 l[i] = t[w*i + lp];
 206
 207             sd_1d97_float(line, mv, mv + lv);
 208
 209             // copy back and deinterleave
 210             for (i =   mv; i < lv; i+=2, j++)
 211                 t[w*j + lp] = F_LFTG_X * l[i] / 2;
 212             for (i = 1-mv; i < lv; i+=2, j++)
 213                 t[w*j + lp] = F_LFTG_K * l[i] / 2;
 214         }
 215     }
 216 }
 217
 218 static void sd_1d97_int(int *p, int i0, int i1)
 219 {
 220     int i;
 221
 222     if (i1 == i0 + 1)
 223         return;
 224
 225     extend97_int(p, i0, i1);
 226     i0++; i1++;
 227
 228     for (i = i0/2 - 2; i < i1/2 + 1; i++)
 229         p[2 * i + 1] -= (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 230     for (i = i0/2 - 1; i < i1/2 + 1; i++)
 231         p[2 * i]     -= (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 232     for (i = i0/2 - 1; i < i1/2; i++)
 233         p[2 * i + 1] += (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 234     for (i = i0/2; i < i1/2; i++)
 235         p[2 * i]     += (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 236 }
 237
 238 static void dwt_encode97_int(DWTContext *s, int *t)
 239 {
 240     int lev,
 241         w = s->linelen[s->ndeclevels-1][0];
 242     int *line = s->i_linebuf;
 243     line += 5;
 244
 245     for (lev = s->ndeclevels-1; lev >= 0; lev--){
 246         int lh = s->linelen[lev][0],
 247             lv = s->linelen[lev][1],
 248             mh = s->mod[lev][0],
 249             mv = s->mod[lev][1],
 250             lp;
 251         int *l;
 252
 253         // HOR_SD
 254         l = line + mh;
 255         for (lp = 0; lp < lv; lp++){
 256             int i, j = 0;
 257
 258             for (i = 0; i < lh; i++)
 259                 l[i] = t[w*lp + i];
 260
 261             sd_1d97_int(line, mh, mh + lh);
 262
 263             // copy back and deinterleave
 264             for (i =   mh; i < lh; i+=2, j++)
 265                 t[w*lp + j] = ((l[i] * I_LFTG_X) + (1 << 16)) >> 17;
 266             for (i = 1-mh; i < lh; i+=2, j++)
 267                 t[w*lp + j] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
 268         }
 269
 270         // VER_SD
 271         l = line + mv;
 272         for (lp = 0; lp < lh; lp++) {
 273             int i, j = 0;
 274
 275             for (i = 0; i < lv; i++)
 276                 l[i] = t[w*i + lp];
 277
 278             sd_1d97_int(line, mv, mv + lv);
 279
 280             // copy back and deinterleave
 281             for (i =   mv; i < lv; i+=2, j++)
 282                 t[w*j + lp] = ((l[i] * I_LFTG_X) + (1 << 16)) >> 17;
 283             for (i = 1-mv; i < lv; i+=2, j++)
 284                 t[w*j + lp] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
 285         }
 286     }
 287 }
 288
 289 static void sr_1d53(int *p, int i0, int i1)
 290 {
 291     int i;
 292
 293     if (i1 <= i0 + 1) {
 294         if (i0 == 1)
 295             p[1] >>= 1;
 296         return;
 297     }
 298
 299     extend53(p, i0, i1);
 300
 301     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 302         p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
 303     for (i = i0 / 2; i < i1 / 2; i++)
 304         p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
 305 }
 306
 307 static void dwt_decode53(DWTContext *s, int *t)
 308 {
 309     int lev;
 310     int w     = s->linelen[s->ndeclevels - 1][0];
 311     int32_t *line = s->i_linebuf;
 312     line += 3;
 313
 314     for (lev = 0; lev < s->ndeclevels; lev++) {
 315         int lh = s->linelen[lev][0],
 316             lv = s->linelen[lev][1],
 317             mh = s->mod[lev][0],
 318             mv = s->mod[lev][1],
 319             lp;
 320         int *l;
 321
 322         // HOR_SD
 323         l = line + mh;
 324         for (lp = 0; lp < lv; lp++) {
 325             int i, j = 0;
 326             // copy with interleaving
 327             for (i = mh; i < lh; i += 2, j++)
 328                 l[i] = t[w * lp + j];
 329             for (i = 1 - mh; i < lh; i += 2, j++)
 330                 l[i] = t[w * lp + j];
 331
 332             sr_1d53(line, mh, mh + lh);
 333
 334             for (i = 0; i < lh; i++)
 335                 t[w * lp + i] = l[i];
 336         }
 337
 338         // VER_SD
 339         l = line + mv;
 340         for (lp = 0; lp < lh; lp++) {
 341             int i, j = 0;
 342             // copy with interleaving
 343             for (i = mv; i < lv; i += 2, j++)
 344                 l[i] = t[w * j + lp];
 345             for (i = 1 - mv; i < lv; i += 2, j++)
 346                 l[i] = t[w * j + lp];
 347
 348             sr_1d53(line, mv, mv + lv);
 349
 350             for (i = 0; i < lv; i++)
 351                 t[w * i + lp] = l[i];
 352         }
 353     }
 354 }
 355
 356 static void sr_1d97_float(float *p, int i0, int i1)
 357 {
 358     int i;
 359
 360     if (i1 <= i0 + 1) {
 361         if (i0 == 1)
 362             p[1] *= F_LFTG_K/2;
 363         return;
 364     }
 365
 366     extend97_float(p, i0, i1);
 367
 368     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 369         p[2 * i]     -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
 370     /* step 4 */
 371     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 372         p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]);
 373     /*step 5*/
 374     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 375         p[2 * i]     += F_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]);
 376     /* step 6 */
 377     for (i = i0 / 2; i < i1 / 2; i++)
 378         p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]);
 379 }
 380
 381 static void dwt_decode97_float(DWTContext *s, float *t)
 382 {
 383     int lev;
 384     int w       = s->linelen[s->ndeclevels - 1][0];
 385     float *line = s->f_linebuf;
 386     float *data = t;
 387     /* position at index O of line range [0-5,w+5] cf. extend function */
 388     line += 5;
 389
 390     for (lev = 0; lev < s->ndeclevels; lev++) {
 391         int lh = s->linelen[lev][0],
 392             lv = s->linelen[lev][1],
 393             mh = s->mod[lev][0],
 394             mv = s->mod[lev][1],
 395             lp;
 396         float *l;
 397         // HOR_SD
 398         l = line + mh;
 399         for (lp = 0; lp < lv; lp++) {
 400             int i, j = 0;
 401             // copy with interleaving
 402             for (i = mh; i < lh; i += 2, j++)
 403                 l[i] = data[w * lp + j] * F_LFTG_K;
 404             for (i = 1 - mh; i < lh; i += 2, j++)
 405                 l[i] = data[w * lp + j] * F_LFTG_X;
 406
 407             sr_1d97_float(line, mh, mh + lh);
 408
 409             for (i = 0; i < lh; i++)
 410                 data[w * lp + i] = l[i];
 411         }
 412
 413         // VER_SD
 414         l = line + mv;
 415         for (lp = 0; lp < lh; lp++) {
 416             int i, j = 0;
 417             // copy with interleaving
 418             for (i = mv; i < lv; i += 2, j++)
 419                 l[i] = data[w * j + lp] * F_LFTG_K;
 420             for (i = 1 - mv; i < lv; i += 2, j++)
 421                 l[i] = data[w * j + lp] * F_LFTG_X;
 422
 423             sr_1d97_float(line, mv, mv + lv);
 424
 425             for (i = 0; i < lv; i++)
 426                 data[w * i + lp] = l[i];
 427         }
 428     }
 429 }
 430
 431 static void sr_1d97_int(int32_t *p, int i0, int i1)
 432 {
 433     int i;
 434
 435     if (i1 <= i0 + 1) {
 436         if (i0 == 1)
 437             p[1] = (p[1] * I_LFTG_K + (1<<16)) >> 17;
 438         return;
 439     }
 440
 441     extend97_int(p, i0, i1);
 442
 443     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 444         p[2 * i]     -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 445     /* step 4 */
 446     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 447         p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 448     /*step 5*/
 449     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 450         p[2 * i]     += (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 451     /* step 6 */
 452     for (i = i0 / 2; i < i1 / 2; i++)
 453         p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 454 }
 455
 456 static void dwt_decode97_int(DWTContext *s, int32_t *t)
 457 {
 458     int lev;
 459     int w       = s->linelen[s->ndeclevels - 1][0];
 460     int32_t *line = s->i_linebuf;
 461     int32_t *data = t;
 462     /* position at index O of line range [0-5,w+5] cf. extend function */
 463     line += 5;
 464
 465     for (lev = 0; lev < s->ndeclevels; lev++) {
 466         int lh = s->linelen[lev][0],
 467             lv = s->linelen[lev][1],
 468             mh = s->mod[lev][0],
 469             mv = s->mod[lev][1],
 470             lp;
 471         int32_t *l;
 472         // HOR_SD
 473         l = line + mh;
 474         for (lp = 0; lp < lv; lp++) {
 475             int i, j = 0;
 476             // rescale with interleaving
 477             for (i = mh; i < lh; i += 2, j++)
 478                 l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
 479             for (i = 1 - mh; i < lh; i += 2, j++)
 480                 l[i] = ((data[w * lp + j] * I_LFTG_X) + (1 << 15)) >> 16;
 481
 482             sr_1d97_int(line, mh, mh + lh);
 483
 484             for (i = 0; i < lh; i++)
 485                 data[w * lp + i] = l[i];
 486         }
 487
 488         // VER_SD
 489         l = line + mv;
 490         for (lp = 0; lp < lh; lp++) {
 491             int i, j = 0;
 492             // rescale with interleaving
 493             for (i = mv; i < lv; i += 2, j++)
 494                 l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
 495             for (i = 1 - mv; i < lv; i += 2, j++)
 496                 l[i] = ((data[w * j + lp] * I_LFTG_X) + (1 << 15)) >> 16;
 497
 498             sr_1d97_int(line, mv, mv + lv);
 499
 500             for (i = 0; i < lv; i++)
 501                 data[w * i + lp] = l[i];
 502         }
 503     }
 504 }
 505
 506 int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
 507                          int decomp_levels, int type)
 508 {
 509     int i, j, lev = decomp_levels, maxlen,
 510         b[2][2];
 511
 512     s->ndeclevels = decomp_levels;
 513     s->type       = type;
 514
 515     for (i = 0; i < 2; i++)
 516         for (j = 0; j < 2; j++)
 517             b[i][j] = border[i][j];
 518
 519     maxlen = FFMAX(b[0][1] - b[0][0],
 520                    b[1][1] - b[1][0]);
 521     while (--lev >= 0)
 522         for (i = 0; i < 2; i++) {
 523             s->linelen[lev][i] = b[i][1] - b[i][0];
 524             s->mod[lev][i]     = b[i][0] & 1;
 525             for (j = 0; j < 2; j++)
 526                 b[i][j] = (b[i][j] + 1) >> 1;
 527         }
 528     switch (type) {
 529     case FF_DWT97:
 530         s->f_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->f_linebuf));
 531         if (!s->f_linebuf)
 532             return AVERROR(ENOMEM);
 533         break;
 534      case FF_DWT97_INT:
 535         s->i_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->i_linebuf));
 536         if (!s->i_linebuf)
 537             return AVERROR(ENOMEM);
 538         break;
 539     case FF_DWT53:
 540         s->i_linebuf = av_malloc_array((maxlen +  6), sizeof(*s->i_linebuf));
 541         if (!s->i_linebuf)
 542             return AVERROR(ENOMEM);
 543         break;
 544     default:
 545         return -1;
 546     }
 547     return 0;
 548 }
 549
 550 int ff_dwt_encode(DWTContext *s, void *t)
 551 {
 552     switch(s->type){
 553         case FF_DWT97:
 554             dwt_encode97_float(s, t); break;
 555         case FF_DWT97_INT:
 556             dwt_encode97_int(s, t); break;
 557         case FF_DWT53:
 558             dwt_encode53(s, t); break;
 559         default:
 560             return -1;
 561     }
 562     return 0;
 563 }
 564
 565 int ff_dwt_decode(DWTContext *s, void *t)
 566 {
 567     switch (s->type) {
 568     case FF_DWT97:
 569         dwt_decode97_float(s, t);
 570         break;
 571     case FF_DWT97_INT:
 572         dwt_decode97_int(s, t);
 573         break;
 574     case FF_DWT53:
 575         dwt_decode53(s, t);
 576         break;
 577     default:
 578         return -1;
 579     }
 580     return 0;
 581 }
 582
 583 void ff_dwt_destroy(DWTContext *s)
 584 {
 585     av_freep(&s->f_linebuf);
 586     av_freep(&s->i_linebuf);
 587 }