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     /*step 1*/
 363     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 364         p[2 * i]     *= F_LFTG_K;
 365     /* step 2*/
 366     for (i = i0 / 2 - 2; i < i1 / 2 + 2; i++)
 367         p[2 * i + 1] *= F_LFTG_X;
 368     /* step 3*/
 369     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 370         p[2 * i]     -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
 371     /* step 4 */
 372     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 373         p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]);
 374     /*step 5*/
 375     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 376         p[2 * i]     += F_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]);
 377     /* step 6 */
 378     for (i = i0 / 2; i < i1 / 2; i++)
 379         p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]);
 380 }
 381
 382 static void dwt_decode97_float(DWTContext *s, float *t)
 383 {
 384     int lev;
 385     int w       = s->linelen[s->ndeclevels - 1][0];
 386     float *line = s->f_linebuf;
 387     float *data = t;
 388     /* position at index O of line range [0-5,w+5] cf. extend function */
 389     line += 5;
 390
 391     for (lev = 0; lev < s->ndeclevels; lev++) {
 392         int lh = s->linelen[lev][0],
 393             lv = s->linelen[lev][1],
 394             mh = s->mod[lev][0],
 395             mv = s->mod[lev][1],
 396             lp;
 397         float *l;
 398         // HOR_SD
 399         l = line + mh;
 400         for (lp = 0; lp < lv; lp++) {
 401             int i, j = 0;
 402             // copy with interleaving
 403             for (i = mh; i < lh; i += 2, j++)
 404                 l[i] = data[w * lp + j];
 405             for (i = 1 - mh; i < lh; i += 2, j++)
 406                 l[i] = data[w * lp + j];
 407
 408             sr_1d97_float(line, mh, mh + lh);
 409
 410             for (i = 0; i < lh; i++)
 411                 data[w * lp + i] = l[i];
 412         }
 413
 414         // VER_SD
 415         l = line + mv;
 416         for (lp = 0; lp < lh; lp++) {
 417             int i, j = 0;
 418             // copy with interleaving
 419             for (i = mv; i < lv; i += 2, j++)
 420                 l[i] = data[w * j + lp];
 421             for (i = 1 - mv; i < lv; i += 2, j++)
 422                 l[i] = data[w * j + lp];
 423
 424             sr_1d97_float(line, mv, mv + lv);
 425
 426             for (i = 0; i < lv; i++)
 427                 data[w * i + lp] = l[i];
 428         }
 429     }
 430 }
 431
 432 static void sr_1d97_int(int32_t *p, int i0, int i1)
 433 {
 434     int i;
 435
 436     if (i1 == i0 + 1)
 437         return;
 438
 439     extend97_int(p, i0, i1);
 440
 441     /*step 1*/
 442     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 443         p[2 * i]      = ((p[2 * i] * I_LFTG_K)     + (1 << 15)) >> 16;
 444     /* step 2*/
 445     for (i = i0 / 2 - 2; i < i1 / 2 + 2; i++)
 446         p[2 * i + 1]  = ((p[2 * i + 1] * I_LFTG_X) + (1 << 15)) >> 16;
 447     /* step 3*/
 448     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 449         p[2 * i]     -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 450     /* step 4 */
 451     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 452         p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 453     /*step 5*/
 454     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 455         p[2 * i]     += (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 456     /* step 6 */
 457     for (i = i0 / 2; i < i1 / 2; i++)
 458         p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 459 }
 460
 461 static void dwt_decode97_int(DWTContext *s, int32_t *t)
 462 {
 463     int lev;
 464     int w       = s->linelen[s->ndeclevels - 1][0];
 465     int32_t *line = s->i_linebuf;
 466     int32_t *data = t;
 467     /* position at index O of line range [0-5,w+5] cf. extend function */
 468     line += 5;
 469
 470     for (lev = 0; lev < s->ndeclevels; lev++) {
 471         int lh = s->linelen[lev][0],
 472             lv = s->linelen[lev][1],
 473             mh = s->mod[lev][0],
 474             mv = s->mod[lev][1],
 475             lp;
 476         int32_t *l;
 477         // HOR_SD
 478         l = line + mh;
 479         for (lp = 0; lp < lv; lp++) {
 480             int i, j = 0;
 481             // copy with interleaving
 482             for (i = mh; i < lh; i += 2, j++)
 483                 l[i] = data[w * lp + j];
 484             for (i = 1 - mh; i < lh; i += 2, j++)
 485                 l[i] = data[w * lp + j];
 486
 487             sr_1d97_int(line, mh, mh + lh);
 488
 489             for (i = 0; i < lh; i++)
 490                 data[w * lp + i] = l[i];
 491         }
 492
 493         // VER_SD
 494         l = line + mv;
 495         for (lp = 0; lp < lh; lp++) {
 496             int i, j = 0;
 497             // copy with interleaving
 498             for (i = mv; i < lv; i += 2, j++)
 499                 l[i] = data[w * j + lp];
 500             for (i = 1 - mv; i < lv; i += 2, j++)
 501                 l[i] = data[w * j + lp];
 502
 503             sr_1d97_int(line, mv, mv + lv);
 504
 505             for (i = 0; i < lv; i++)
 506                 data[w * i + lp] = l[i];
 507         }
 508     }
 509 }
 510
 511 int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
 512                          int decomp_levels, int type)
 513 {
 514     int i, j, lev = decomp_levels, maxlen,
 515         b[2][2];
 516
 517     s->ndeclevels = decomp_levels;
 518     s->type       = type;
 519
 520     for (i = 0; i < 2; i++)
 521         for (j = 0; j < 2; j++)
 522             b[i][j] = border[i][j];
 523
 524     maxlen = FFMAX(b[0][1] - b[0][0],
 525                    b[1][1] - b[1][0]);
 526     while (--lev >= 0)
 527         for (i = 0; i < 2; i++) {
 528             s->linelen[lev][i] = b[i][1] - b[i][0];
 529             s->mod[lev][i]     = b[i][0] & 1;
 530             for (j = 0; j < 2; j++)
 531                 b[i][j] = (b[i][j] + 1) >> 1;
 532         }
 533     switch (type) {
 534     case FF_DWT97:
 535         s->f_linebuf = av_malloc((maxlen + 12) * sizeof(*s->f_linebuf));
 536         if (!s->f_linebuf)
 537             return AVERROR(ENOMEM);
 538         break;
 539      case FF_DWT97_INT:
 540         s->i_linebuf = av_malloc((maxlen + 12) * sizeof(*s->i_linebuf));
 541         if (!s->i_linebuf)
 542             return AVERROR(ENOMEM);
 543         break;
 544     case FF_DWT53:
 545         s->i_linebuf = av_malloc((maxlen +  6) * sizeof(*s->i_linebuf));
 546         if (!s->i_linebuf)
 547             return AVERROR(ENOMEM);
 548         break;
 549     default:
 550         return -1;
 551     }
 552     return 0;
 553 }
 554
 555 int ff_dwt_encode(DWTContext *s, void *t)
 556 {
 557     switch(s->type){
 558         case FF_DWT97:
 559             dwt_encode97_float(s, t); break;
 560         case FF_DWT97_INT:
 561             dwt_encode97_int(s, t); break;
 562         case FF_DWT53:
 563             dwt_encode53(s, t); break;
 564         default:
 565             return -1;
 566     }
 567     return 0;
 568 }
 569
 570 int ff_dwt_decode(DWTContext *s, void *t)
 571 {
 572     switch (s->type) {
 573     case FF_DWT97:
 574         dwt_decode97_float(s, t);
 575         break;
 576     case FF_DWT97_INT:
 577         dwt_decode97_int(s, t);
 578         break;
 579     case FF_DWT53:
 580         dwt_decode53(s, t);
 581         break;
 582     default:
 583         return -1;
 584     }
 585     return 0;
 586 }
 587
 588 void ff_dwt_destroy(DWTContext *s)
 589 {
 590     av_freep(&s->f_linebuf);
 591     av_freep(&s->i_linebuf);
 592 }