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
  55 static inline void extend53(int *p, int i0, int i1)
  56 {
  57     p[i0 - 1] = p[i0 + 1];
  58     p[i1]     = p[i1 - 2];
  59     p[i0 - 2] = p[i0 + 2];
  60     p[i1 + 1] = p[i1 - 3];
  61 }
  62
  63 static inline void extend97_float(float *p, int i0, int i1)
  64 {
  65     int i;
  66
  67     for (i = 1; i <= 4; i++) {
  68         p[i0 - i]     = p[i0 + i];
  69         p[i1 + i - 1] = p[i1 - i - 1];
  70     }
  71 }
  72
  73 static inline void extend97_int(int32_t *p, int i0, int i1)
  74 {
  75     int i;
  76
  77     for (i = 1; i <= 4; i++) {
  78         p[i0 - i]     = p[i0 + i];
  79         p[i1 + i - 1] = p[i1 - i - 1];
  80     }
  81 }
  82
  83 static void sr_1d53(int *p, int i0, int i1)
  84 {
  85     int i;
  86
  87     if (i1 == i0 + 1)
  88         return;
  89
  90     extend53(p, i0, i1);
  91
  92     for (i = i0 / 2; i < i1 / 2 + 1; i++)
  93         p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
  94     for (i = i0 / 2; i < i1 / 2; i++)
  95         p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
  96 }
  97
  98 static void dwt_decode53(DWTContext *s, int *t)
  99 {
 100     int lev;
 101     int w     = s->linelen[s->ndeclevels - 1][0];
 102     int32_t *line = s->i_linebuf;
 103     line += 3;
 104
 105     for (lev = 0; lev < s->ndeclevels; lev++) {
 106         int lh = s->linelen[lev][0],
 107             lv = s->linelen[lev][1],
 108             mh = s->mod[lev][0],
 109             mv = s->mod[lev][1],
 110             lp;
 111         int *l;
 112
 113         // HOR_SD
 114         l = line + mh;
 115         for (lp = 0; lp < lv; lp++) {
 116             int i, j = 0;
 117             // copy with interleaving
 118             for (i = mh; i < lh; i += 2, j++)
 119                 l[i] = t[w * lp + j];
 120             for (i = 1 - mh; i < lh; i += 2, j++)
 121                 l[i] = t[w * lp + j];
 122
 123             sr_1d53(line, mh, mh + lh);
 124
 125             for (i = 0; i < lh; i++)
 126                 t[w * lp + i] = l[i];
 127         }
 128
 129         // VER_SD
 130         l = line + mv;
 131         for (lp = 0; lp < lh; lp++) {
 132             int i, j = 0;
 133             // copy with interleaving
 134             for (i = mv; i < lv; i += 2, j++)
 135                 l[i] = t[w * j + lp];
 136             for (i = 1 - mv; i < lv; i += 2, j++)
 137                 l[i] = t[w * j + lp];
 138
 139             sr_1d53(line, mv, mv + lv);
 140
 141             for (i = 0; i < lv; i++)
 142                 t[w * i + lp] = l[i];
 143         }
 144     }
 145 }
 146
 147 static void sr_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
 156     /*step 1*/
 157     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 158         p[2 * i]     *= F_LFTG_K;
 159     /* step 2*/
 160     for (i = i0 / 2 - 2; i < i1 / 2 + 2; i++)
 161         p[2 * i + 1] *= F_LFTG_X;
 162     /* step 3*/
 163     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 164         p[2 * i]     -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
 165     /* step 4 */
 166     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 167         p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]);
 168     /*step 5*/
 169     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 170         p[2 * i]     += F_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]);
 171     /* step 6 */
 172     for (i = i0 / 2; i < i1 / 2; i++)
 173         p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]);
 174 }
 175
 176 static void dwt_decode97_float(DWTContext *s, float *t)
 177 {
 178     int lev;
 179     int w       = s->linelen[s->ndeclevels - 1][0];
 180     float *line = s->f_linebuf;
 181     float *data = t;
 182     /* position at index O of line range [0-5,w+5] cf. extend function */
 183     line += 5;
 184
 185     for (lev = 0; lev < s->ndeclevels; lev++) {
 186         int lh = s->linelen[lev][0],
 187             lv = s->linelen[lev][1],
 188             mh = s->mod[lev][0],
 189             mv = s->mod[lev][1],
 190             lp;
 191         float *l;
 192         // HOR_SD
 193         l = line + mh;
 194         for (lp = 0; lp < lv; lp++) {
 195             int i, j = 0;
 196             // copy with interleaving
 197             for (i = mh; i < lh; i += 2, j++)
 198                 l[i] = data[w * lp + j];
 199             for (i = 1 - mh; i < lh; i += 2, j++)
 200                 l[i] = data[w * lp + j];
 201
 202             sr_1d97_float(line, mh, mh + lh);
 203
 204             for (i = 0; i < lh; i++)
 205                 data[w * lp + i] = l[i];
 206         }
 207
 208         // VER_SD
 209         l = line + mv;
 210         for (lp = 0; lp < lh; lp++) {
 211             int i, j = 0;
 212             // copy with interleaving
 213             for (i = mv; i < lv; i += 2, j++)
 214                 l[i] = data[w * j + lp];
 215             for (i = 1 - mv; i < lv; i += 2, j++)
 216                 l[i] = data[w * j + lp];
 217
 218             sr_1d97_float(line, mv, mv + lv);
 219
 220             for (i = 0; i < lv; i++)
 221                 data[w * i + lp] = l[i];
 222         }
 223     }
 224 }
 225
 226 static void sr_1d97_int(int32_t *p, int i0, int i1)
 227 {
 228     int i;
 229
 230     if (i1 == i0 + 1)
 231         return;
 232
 233     extend97_int(p, i0, i1);
 234
 235     /*step 1*/
 236     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 237         p[2 * i]      = ((p[2 * i] * I_LFTG_K)     + (1 << 15)) >> 16;
 238     /* step 2*/
 239     for (i = i0 / 2 - 2; i < i1 / 2 + 2; i++)
 240         p[2 * i + 1]  = ((p[2 * i + 1] * I_LFTG_X) + (1 << 15)) >> 16;
 241     /* step 3*/
 242     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 243         p[2 * i]     -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 244     /* step 4 */
 245     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 246         p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 247     /*step 5*/
 248     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 249         p[2 * i]     += (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 250     /* step 6 */
 251     for (i = i0 / 2; i < i1 / 2; i++)
 252         p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 253 }
 254
 255 static void dwt_decode97_int(DWTContext *s, int32_t *t)
 256 {
 257     int lev;
 258     int w       = s->linelen[s->ndeclevels - 1][0];
 259     int32_t *line = s->i_linebuf;
 260     int32_t *data = t;
 261     /* position at index O of line range [0-5,w+5] cf. extend function */
 262     line += 5;
 263
 264     for (lev = 0; lev < s->ndeclevels; lev++) {
 265         int lh = s->linelen[lev][0],
 266             lv = s->linelen[lev][1],
 267             mh = s->mod[lev][0],
 268             mv = s->mod[lev][1],
 269             lp;
 270         int32_t *l;
 271         // HOR_SD
 272         l = line + mh;
 273         for (lp = 0; lp < lv; lp++) {
 274             int i, j = 0;
 275             // copy with interleaving
 276             for (i = mh; i < lh; i += 2, j++)
 277                 l[i] = data[w * lp + j];
 278             for (i = 1 - mh; i < lh; i += 2, j++)
 279                 l[i] = data[w * lp + j];
 280
 281             sr_1d97_int(line, mh, mh + lh);
 282
 283             for (i = 0; i < lh; i++)
 284                 data[w * lp + i] = l[i];
 285         }
 286
 287         // VER_SD
 288         l = line + mv;
 289         for (lp = 0; lp < lh; lp++) {
 290             int i, j = 0;
 291             // copy with interleaving
 292             for (i = mv; i < lv; i += 2, j++)
 293                 l[i] = data[w * j + lp];
 294             for (i = 1 - mv; i < lv; i += 2, j++)
 295                 l[i] = data[w * j + lp];
 296
 297             sr_1d97_int(line, mv, mv + lv);
 298
 299             for (i = 0; i < lv; i++)
 300                 data[w * i + lp] = l[i];
 301         }
 302     }
 303 }
 304
 305 int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
 306                          int decomp_levels, int type)
 307 {
 308     int i, j, lev = decomp_levels, maxlen,
 309         b[2][2];
 310
 311     s->ndeclevels = decomp_levels;
 312     s->type       = type;
 313
 314     for (i = 0; i < 2; i++)
 315         for (j = 0; j < 2; j++)
 316             b[i][j] = border[i][j];
 317
 318     maxlen = FFMAX(b[0][1] - b[0][0],
 319                    b[1][1] - b[1][0]);
 320     while (--lev >= 0)
 321         for (i = 0; i < 2; i++) {
 322             s->linelen[lev][i] = b[i][1] - b[i][0];
 323             s->mod[lev][i]     = b[i][0] & 1;
 324             for (j = 0; j < 2; j++)
 325                 b[i][j] = (b[i][j] + 1) >> 1;
 326         }
 327     switch (type) {
 328     case FF_DWT97:
 329         s->f_linebuf = av_malloc((maxlen + 12) * sizeof(*s->f_linebuf));
 330         if (!s->f_linebuf)
 331             return AVERROR(ENOMEM);
 332         break;
 333      case FF_DWT97_INT:
 334         s->i_linebuf = av_malloc((maxlen + 12) * sizeof(*s->i_linebuf));
 335         if (!s->i_linebuf)
 336             return AVERROR(ENOMEM);
 337         break;
 338     case FF_DWT53:
 339         s->i_linebuf = av_malloc((maxlen +  6) * sizeof(*s->i_linebuf));
 340         if (!s->i_linebuf)
 341             return AVERROR(ENOMEM);
 342         break;
 343     default:
 344         return -1;
 345     }
 346     return 0;
 347 }
 348
 349 int ff_dwt_decode(DWTContext *s, void *t)
 350 {
 351     switch (s->type) {
 352     case FF_DWT97:
 353         dwt_decode97_float(s, t);
 354         break;
 355     case FF_DWT97_INT:
 356         dwt_decode97_int(s, t);
 357         break;
 358     case FF_DWT53:
 359         dwt_decode53(s, t);
 360         break;
 361     default:
 362         return -1;
 363     }
 364     return 0;
 365 }
 366
 367 void ff_dwt_destroy(DWTContext *s)
 368 {
 369     av_freep(&s->f_linebuf);
 370     av_freep(&s->i_linebuf);
 371 }