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 Libav.
   7  *
   8  * Libav 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  * Libav 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 Libav; 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     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 157         p[2 * i]     -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
 158     /* step 4 */
 159     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 160         p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]);
 161     /*step 5*/
 162     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 163         p[2 * i]     += F_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]);
 164     /* step 6 */
 165     for (i = i0 / 2; i < i1 / 2; i++)
 166         p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]);
 167 }
 168
 169 static void dwt_decode97_float(DWTContext *s, float *t)
 170 {
 171     int lev;
 172     int w       = s->linelen[s->ndeclevels - 1][0];
 173     float *line = s->f_linebuf;
 174     float *data = t;
 175     /* position at index O of line range [0-5,w+5] cf. extend function */
 176     line += 5;
 177
 178     for (lev = 0; lev < s->ndeclevels; lev++) {
 179         int lh = s->linelen[lev][0],
 180             lv = s->linelen[lev][1],
 181             mh = s->mod[lev][0],
 182             mv = s->mod[lev][1],
 183             lp;
 184         float *l;
 185         // HOR_SD
 186         l = line + mh;
 187         for (lp = 0; lp < lv; lp++) {
 188             int i, j = 0;
 189             // copy with interleaving
 190             for (i = mh; i < lh; i += 2, j++)
 191                 l[i] = data[w * lp + j] * F_LFTG_K;
 192             for (i = 1 - mh; i < lh; i += 2, j++)
 193                 l[i] = data[w * lp + j] * F_LFTG_X;
 194
 195             sr_1d97_float(line, mh, mh + lh);
 196
 197             for (i = 0; i < lh; i++)
 198                 data[w * lp + i] = l[i];
 199         }
 200
 201         // VER_SD
 202         l = line + mv;
 203         for (lp = 0; lp < lh; lp++) {
 204             int i, j = 0;
 205             // copy with interleaving
 206             for (i = mv; i < lv; i += 2, j++)
 207                 l[i] = data[w * j + lp] * F_LFTG_K;
 208             for (i = 1 - mv; i < lv; i += 2, j++)
 209                 l[i] = data[w * j + lp] * F_LFTG_X;
 210
 211             sr_1d97_float(line, mv, mv + lv);
 212
 213             for (i = 0; i < lv; i++)
 214                 data[w * i + lp] = l[i];
 215         }
 216     }
 217 }
 218
 219 static void sr_1d97_int(int32_t *p, int i0, int i1)
 220 {
 221     int i;
 222
 223     if (i1 == i0 + 1)
 224         return;
 225
 226     extend97_int(p, i0, i1);
 227
 228     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
 229         p[2 * i]     -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 230     /* step 4 */
 231     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
 232         p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 233     /*step 5*/
 234     for (i = i0 / 2; i < i1 / 2 + 1; i++)
 235         p[2 * i]     += (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
 236     /* step 6 */
 237     for (i = i0 / 2; i < i1 / 2; i++)
 238         p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
 239 }
 240
 241 static void dwt_decode97_int(DWTContext *s, int32_t *t)
 242 {
 243     int lev;
 244     int w       = s->linelen[s->ndeclevels - 1][0];
 245     int32_t *line = s->i_linebuf;
 246     int32_t *data = t;
 247     /* position at index O of line range [0-5,w+5] cf. extend function */
 248     line += 5;
 249
 250     for (lev = 0; lev < s->ndeclevels; lev++) {
 251         int lh = s->linelen[lev][0],
 252             lv = s->linelen[lev][1],
 253             mh = s->mod[lev][0],
 254             mv = s->mod[lev][1],
 255             lp;
 256         int32_t *l;
 257         // HOR_SD
 258         l = line + mh;
 259         for (lp = 0; lp < lv; lp++) {
 260             int i, j = 0;
 261             // rescale with interleaving
 262             for (i = mh; i < lh; i += 2, j++)
 263                 l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
 264             for (i = 1 - mh; i < lh; i += 2, j++)
 265                 l[i] = ((data[w * lp + j] * I_LFTG_X) + (1 << 15)) >> 16;
 266
 267             sr_1d97_int(line, mh, mh + lh);
 268
 269             for (i = 0; i < lh; i++)
 270                 data[w * lp + i] = l[i];
 271         }
 272
 273         // VER_SD
 274         l = line + mv;
 275         for (lp = 0; lp < lh; lp++) {
 276             int i, j = 0;
 277             // rescale with interleaving
 278             for (i = mv; i < lv; i += 2, j++)
 279                 l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
 280             for (i = 1 - mv; i < lv; i += 2, j++)
 281                 l[i] = ((data[w * j + lp] * I_LFTG_X) + (1 << 15)) >> 16;
 282
 283             sr_1d97_int(line, mv, mv + lv);
 284
 285             for (i = 0; i < lv; i++)
 286                 data[w * i + lp] = l[i];
 287         }
 288     }
 289 }
 290
 291 int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
 292                          int decomp_levels, int type)
 293 {
 294     int i, j, lev = decomp_levels, maxlen,
 295         b[2][2];
 296
 297     s->ndeclevels = decomp_levels;
 298     s->type       = type;
 299
 300     for (i = 0; i < 2; i++)
 301         for (j = 0; j < 2; j++)
 302             b[i][j] = border[i][j];
 303
 304     maxlen = FFMAX(b[0][1] - b[0][0],
 305                    b[1][1] - b[1][0]);
 306     while (--lev >= 0)
 307         for (i = 0; i < 2; i++) {
 308             s->linelen[lev][i] = b[i][1] - b[i][0];
 309             s->mod[lev][i]     = b[i][0] & 1;
 310             for (j = 0; j < 2; j++)
 311                 b[i][j] = (b[i][j] + 1) >> 1;
 312         }
 313     switch (type) {
 314     case FF_DWT97:
 315         s->f_linebuf = av_malloc((maxlen + 12) * sizeof(*s->f_linebuf));
 316         if (!s->f_linebuf)
 317             return AVERROR(ENOMEM);
 318         break;
 319      case FF_DWT97_INT:
 320         s->i_linebuf = av_malloc((maxlen + 12) * sizeof(*s->i_linebuf));
 321         if (!s->i_linebuf)
 322             return AVERROR(ENOMEM);
 323         break;
 324     case FF_DWT53:
 325         s->i_linebuf = av_malloc((maxlen +  6) * sizeof(*s->i_linebuf));
 326         if (!s->i_linebuf)
 327             return AVERROR(ENOMEM);
 328         break;
 329     default:
 330         return -1;
 331     }
 332     return 0;
 333 }
 334
 335 int ff_dwt_decode(DWTContext *s, void *t)
 336 {
 337     switch (s->type) {
 338     case FF_DWT97:
 339         dwt_decode97_float(s, t);
 340         break;
 341     case FF_DWT97_INT:
 342         dwt_decode97_int(s, t);
 343         break;
 344     case FF_DWT53:
 345         dwt_decode53(s, t);
 346         break;
 347     default:
 348         return -1;
 349     }
 350     return 0;
 351 }
 352
 353 void ff_dwt_destroy(DWTContext *s)
 354 {
 355     av_freep(&s->f_linebuf);
 356     av_freep(&s->i_linebuf);
 357 }