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