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         return;
 295
 296     extend53(p, i0, i1);
 297
 298     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 299         p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
 300     for (i = i0 / 2; i < i1 / 2; i++)
 301         p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
 302 }
 303
 304 static void dwt_decode53(DWTContext *s, int *t)
 305 {
 306     int lev;
 307     int w     = s->linelen[s->ndeclevels - 1][0];
 308     int32_t *line = s->i_linebuf;
 309     line += 3;
 310
 311     for (lev = 0; lev < s->ndeclevels; lev++) {
 312         int lh = s->linelen[lev][0],
 313             lv = s->linelen[lev][1],
 314             mh = s->mod[lev][0],
 315             mv = s->mod[lev][1],
 316             lp;
 317         int *l;
 318
 319         // HOR_SD
 320         l = line + mh;
 321         for (lp = 0; lp < lv; lp++) {
 322             int i, j = 0;
 323             // copy with interleaving
 324             for (i = mh; i < lh; i += 2, j++)
 325                 l[i] = t[w * lp + j];
 326             for (i = 1 - mh; i < lh; i += 2, j++)
 327                 l[i] = t[w * lp + j];
 328
 329             sr_1d53(line, mh, mh + lh);
 330
 331             for (i = 0; i < lh; i++)
 332                 t[w * lp + i] = l[i];
 333         }
 334
 335         // VER_SD
 336         l = line + mv;
 337         for (lp = 0; lp < lh; lp++) {
 338             int i, j = 0;
 339             // copy with interleaving
 340             for (i = mv; i < lv; i += 2, j++)
 341                 l[i] = t[w * j + lp];
 342             for (i = 1 - mv; i < lv; i += 2, j++)
 343                 l[i] = t[w * j + lp];
 344
 345             sr_1d53(line, mv, mv + lv);
 346
 347             for (i = 0; i < lv; i++)
 348                 t[w * i + lp] = l[i];
 349         }
 350     }
 351 }
 352
 353 static void sr_1d97_float(float *p, int i0, int i1)
 354 {
 355     int i;
 356
 357     if (i1 == i0 + 1)
 358         return;
 359
 360     extend97_float(p, i0, i1);
 361
 362     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 363         p[2 * i]     -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
 364     /* step 4 */
 365     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 366         p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]);
 367     /*step 5*/
 368     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 369         p[2 * i]     += F_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]);
 370     /* step 6 */
 371     for (i = i0 / 2; i < i1 / 2; i++)
 372         p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]);
 373 }
 374
 375 static void dwt_decode97_float(DWTContext *s, float *t)
 376 {
 377     int lev;
 378     int w       = s->linelen[s->ndeclevels - 1][0];
 379     float *line = s->f_linebuf;
 380     float *data = t;
 381     /* position at index O of line range [0-5,w+5] cf. extend function */
 382     line += 5;
 383
 384     for (lev = 0; lev < s->ndeclevels; lev++) {
 385         int lh = s->linelen[lev][0],
 386             lv = s->linelen[lev][1],
 387             mh = s->mod[lev][0],
 388             mv = s->mod[lev][1],
 389             lp;
 390         float *l;
 391         // HOR_SD
 392         l = line + mh;
 393         for (lp = 0; lp < lv; lp++) {
 394             int i, j = 0;
 395             // copy with interleaving
 396             for (i = mh; i < lh; i += 2, j++)
 397                 l[i] = data[w * lp + j] * F_LFTG_K;
 398             for (i = 1 - mh; i < lh; i += 2, j++)
 399                 l[i] = data[w * lp + j] * F_LFTG_X;
 400
 401             sr_1d97_float(line, mh, mh + lh);
 402
 403             for (i = 0; i < lh; i++)
 404                 data[w * lp + i] = l[i];
 405         }
 406
 407         // VER_SD
 408         l = line + mv;
 409         for (lp = 0; lp < lh; lp++) {
 410             int i, j = 0;
 411             // copy with interleaving
 412             for (i = mv; i < lv; i += 2, j++)
 413                 l[i] = data[w * j + lp] * F_LFTG_K;
 414             for (i = 1 - mv; i < lv; i += 2, j++)
 415                 l[i] = data[w * j + lp] * F_LFTG_X;
 416
 417             sr_1d97_float(line, mv, mv + lv);
 418
 419             for (i = 0; i < lv; i++)
 420                 data[w * i + lp] = l[i];
 421         }
 422     }
 423 }
 424
 425 static void sr_1d97_int(int32_t *p, int i0, int i1)
 426 {
 427     int i;
 428
 429     if (i1 == i0 + 1)
 430         return;
 431
 432     extend97_int(p, i0, i1);
 433
 434     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 435         p[2 * i]     -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 436     /* step 4 */
 437     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 438         p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 439     /*step 5*/
 440     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 441         p[2 * i]     += (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 442     /* step 6 */
 443     for (i = i0 / 2; i < i1 / 2; i++)
 444         p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 445 }
 446
 447 static void dwt_decode97_int(DWTContext *s, int32_t *t)
 448 {
 449     int lev;
 450     int w       = s->linelen[s->ndeclevels - 1][0];
 451     int32_t *line = s->i_linebuf;
 452     int32_t *data = t;
 453     /* position at index O of line range [0-5,w+5] cf. extend function */
 454     line += 5;
 455
 456     for (lev = 0; lev < s->ndeclevels; lev++) {
 457         int lh = s->linelen[lev][0],
 458             lv = s->linelen[lev][1],
 459             mh = s->mod[lev][0],
 460             mv = s->mod[lev][1],
 461             lp;
 462         int32_t *l;
 463         // HOR_SD
 464         l = line + mh;
 465         for (lp = 0; lp < lv; lp++) {
 466             int i, j = 0;
 467             // rescale with interleaving
 468             for (i = mh; i < lh; i += 2, j++)
 469                 l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
 470             for (i = 1 - mh; i < lh; i += 2, j++)
 471                 l[i] = ((data[w * lp + j] * I_LFTG_X) + (1 << 15)) >> 16;
 472
 473             sr_1d97_int(line, mh, mh + lh);
 474
 475             for (i = 0; i < lh; i++)
 476                 data[w * lp + i] = l[i];
 477         }
 478
 479         // VER_SD
 480         l = line + mv;
 481         for (lp = 0; lp < lh; lp++) {
 482             int i, j = 0;
 483             // rescale with interleaving
 484             for (i = mv; i < lv; i += 2, j++)
 485                 l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
 486             for (i = 1 - mv; i < lv; i += 2, j++)
 487                 l[i] = ((data[w * j + lp] * I_LFTG_X) + (1 << 15)) >> 16;
 488
 489             sr_1d97_int(line, mv, mv + lv);
 490
 491             for (i = 0; i < lv; i++)
 492                 data[w * i + lp] = l[i];
 493         }
 494     }
 495 }
 496
 497 int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
 498                          int decomp_levels, int type)
 499 {
 500     int i, j, lev = decomp_levels, maxlen,
 501         b[2][2];
 502
 503     s->ndeclevels = decomp_levels;
 504     s->type       = type;
 505
 506     for (i = 0; i < 2; i++)
 507         for (j = 0; j < 2; j++)
 508             b[i][j] = border[i][j];
 509
 510     maxlen = FFMAX(b[0][1] - b[0][0],
 511                    b[1][1] - b[1][0]);
 512     while (--lev >= 0)
 513         for (i = 0; i < 2; i++) {
 514             s->linelen[lev][i] = b[i][1] - b[i][0];
 515             s->mod[lev][i]     = b[i][0] & 1;
 516             for (j = 0; j < 2; j++)
 517                 b[i][j] = (b[i][j] + 1) >> 1;
 518         }
 519     switch (type) {
 520     case FF_DWT97:
 521         s->f_linebuf = av_malloc((maxlen + 12) * sizeof(*s->f_linebuf));
 522         if (!s->f_linebuf)
 523             return AVERROR(ENOMEM);
 524         break;
 525      case FF_DWT97_INT:
 526         s->i_linebuf = av_malloc((maxlen + 12) * sizeof(*s->i_linebuf));
 527         if (!s->i_linebuf)
 528             return AVERROR(ENOMEM);
 529         break;
 530     case FF_DWT53:
 531         s->i_linebuf = av_malloc((maxlen +  6) * sizeof(*s->i_linebuf));
 532         if (!s->i_linebuf)
 533             return AVERROR(ENOMEM);
 534         break;
 535     default:
 536         return -1;
 537     }
 538     return 0;
 539 }
 540
 541 int ff_dwt_encode(DWTContext *s, void *t)
 542 {
 543     switch(s->type){
 544         case FF_DWT97:
 545             dwt_encode97_float(s, t); break;
 546         case FF_DWT97_INT:
 547             dwt_encode97_int(s, t); break;
 548         case FF_DWT53:
 549             dwt_encode53(s, t); break;
 550         default:
 551             return -1;
 552     }
 553     return 0;
 554 }
 555
 556 int ff_dwt_decode(DWTContext *s, void *t)
 557 {
 558     switch (s->type) {
 559     case FF_DWT97:
 560         dwt_decode97_float(s, t);
 561         break;
 562     case FF_DWT97_INT:
 563         dwt_decode97_int(s, t);
 564         break;
 565     case FF_DWT53:
 566         dwt_decode53(s, t);
 567         break;
 568     default:
 569         return -1;
 570     }
 571     return 0;
 572 }
 573
 574 void ff_dwt_destroy(DWTContext *s)
 575 {
 576     av_freep(&s->f_linebuf);
 577     av_freep(&s->i_linebuf);
 578 }