libavcodec/j2k_dwt.c

   1 /*
   2  * Discrete wavelet transform
   3  * Copyright (c) 2007 Kamil Nowosad
   4  *
   5  * This file is part of FFmpeg.
   6  *
   7  * FFmpeg is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU Lesser General Public
   9  * License as published by the Free Software Foundation; either
  10  * version 2.1 of the License, or (at your option) any later version.
  11  *
  12  * FFmpeg is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public
  18  * License along with FFmpeg; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  20  */
  21
  22 /**
  23  * Discrete wavelet transform
  24  * @file
  25  * @author Kamil Nowosad
  26  */
  27
  28 #include "j2k_dwt.h"
  29
  30 const static float scale97[] = {1.625786, 1.230174};
  31
  32 static inline void extend53(int *p, int i0, int i1)
  33 {
  34     p[i0 - 1] = p[i0 + 1];
  35     p[i1    ] = p[i1 - 2];
  36     p[i0 - 2] = p[i0 + 2];
  37     p[i1 + 1] = p[i1 - 3];
  38 }
  39
  40 static inline void extend97(float *p, int i0, int i1)
  41 {
  42     int i;
  43
  44     for (i = 1; i <= 4; i++){
  45         p[i0 - i] = p[i0 + i];
  46         p[i1 + i - 1] = p[i1 - i - 1];
  47     }
  48 }
  49
  50 static void sd_1d53(int *p, int i0, int i1)
  51 {
  52     int i;
  53
  54     if (i1 == i0 + 1)
  55         return;
  56
  57     extend53(p, i0, i1);
  58
  59     for (i = (i0+1)/2 - 1; i < (i1+1)/2; i++)
  60         p[2*i+1] -= (p[2*i] + p[2*i+2]) >> 1;
  61     for (i = (i0+1)/2; i < (i1+1)/2; i++)
  62         p[2*i] += (p[2*i-1] + p[2*i+1] + 2) >> 2;
  63 }
  64
  65 static void dwt_encode53(DWTContext *s, int *t)
  66 {
  67     int lev,
  68         w = s->linelen[s->ndeclevels-1][0];
  69     int *line = s->linebuf;
  70     line += 3;
  71
  72     for (lev = s->ndeclevels-1; lev >= 0; lev--){
  73         int lh = s->linelen[lev][0],
  74             lv = s->linelen[lev][1],
  75             mh = s->mod[lev][0],
  76             mv = s->mod[lev][1],
  77             lp;
  78         int *l;
  79
  80         // HOR_SD
  81         l = line + mh;
  82         for (lp = 0; lp < lv; lp++){
  83             int i, j = 0;
  84
  85             for (i = 0; i < lh; i++)
  86                 l[i] = t[w*lp + i];
  87
  88             sd_1d53(line, mh, mh + lh);
  89
  90             // copy back and deinterleave
  91             for (i =   mh; i < lh; i+=2, j++)
  92                 t[w*lp + j] = l[i];
  93             for (i = 1-mh; i < lh; i+=2, j++)
  94                 t[w*lp + j] = l[i];
  95         }
  96
  97         // VER_SD
  98         l = line + mv;
  99         for (lp = 0; lp < lh; lp++) {
 100             int i, j = 0;
 101
 102             for (i = 0; i < lv; i++)
 103                 l[i] = t[w*i + lp];
 104
 105             sd_1d53(line, mv, mv + lv);
 106
 107             // copy back and deinterleave
 108             for (i =   mv; i < lv; i+=2, j++)
 109                 t[w*j + lp] = l[i];
 110             for (i = 1-mv; i < lv; i+=2, j++)
 111                 t[w*j + lp] = l[i];
 112         }
 113     }
 114 }
 115
 116 static void sd_1d97(float *p, int i0, int i1)
 117 {
 118     int i;
 119
 120     if (i1 == i0 + 1)
 121         return;
 122
 123     extend97(p, i0, i1);
 124     i0++; i1++;
 125
 126     for (i = i0/2 - 2; i < i1/2 + 1; i++)
 127         p[2*i+1] -= 1.586134 * (p[2*i] + p[2*i+2]);
 128     for (i = i0/2 - 1; i < i1/2 + 1; i++)
 129         p[2*i] -= 0.052980 * (p[2*i-1] + p[2*i+1]);
 130     for (i = i0/2 - 1; i < i1/2; i++)
 131         p[2*i+1] += 0.882911 * (p[2*i] + p[2*i+2]);
 132     for (i = i0/2; i < i1/2; i++)
 133         p[2*i] += 0.443506 * (p[2*i-1] + p[2*i+1]);
 134 }
 135
 136 static void dwt_encode97(DWTContext *s, int *t)
 137 {
 138     int lev,
 139         w = s->linelen[s->ndeclevels-1][0];
 140     float *line = s->linebuf;
 141     line += 5;
 142
 143     for (lev = s->ndeclevels-1; lev >= 0; lev--){
 144         int lh = s->linelen[lev][0],
 145             lv = s->linelen[lev][1],
 146             mh = s->mod[lev][0],
 147             mv = s->mod[lev][1],
 148             lp;
 149         float *l;
 150
 151         // HOR_SD
 152         l = line + mh;
 153         for (lp = 0; lp < lv; lp++){
 154             int i, j = 0;
 155
 156             for (i = 0; i < lh; i++)
 157                 l[i] = t[w*lp + i];
 158
 159             sd_1d97(line, mh, mh + lh);
 160
 161             // copy back and deinterleave
 162             for (i =   mh; i < lh; i+=2, j++)
 163                 t[w*lp + j] = scale97[mh] * l[i] / 2;
 164             for (i = 1-mh; i < lh; i+=2, j++)
 165                 t[w*lp + j] = scale97[mh] * l[i] / 2;
 166         }
 167
 168         // VER_SD
 169         l = line + mv;
 170         for (lp = 0; lp < lh; lp++) {
 171             int i, j = 0;
 172
 173             for (i = 0; i < lv; i++)
 174                 l[i] = t[w*i + lp];
 175
 176             sd_1d97(line, mv, mv + lv);
 177
 178             // copy back and deinterleave
 179             for (i =   mv; i < lv; i+=2, j++)
 180                 t[w*j + lp] = scale97[mv] * l[i] / 2;
 181             for (i = 1-mv; i < lv; i+=2, j++)
 182                 t[w*j + lp] = scale97[mv] * l[i] / 2;
 183         }
 184     }
 185 }
 186
 187 static void sr_1d53(int *p, int i0, int i1)
 188 {
 189     int i;
 190
 191     if (i1 == i0 + 1)
 192         return;
 193
 194     extend53(p, i0, i1);
 195
 196     for (i = i0/2; i < i1/2 + 1; i++)
 197         p[2*i] -= (p[2*i-1] + p[2*i+1] + 2) >> 2;
 198     for (i = i0/2; i < i1/2; i++)
 199         p[2*i+1] += (p[2*i] + p[2*i+2]) >> 1;
 200 }
 201
 202 static void dwt_decode53(DWTContext *s, int *t)
 203 {
 204     int lev,
 205         w = s->linelen[s->ndeclevels-1][0];
 206     int *line = s->linebuf;
 207     line += 3;
 208
 209     for (lev = 0; lev < s->ndeclevels; lev++){
 210         int lh = s->linelen[lev][0],
 211             lv = s->linelen[lev][1],
 212             mh = s->mod[lev][0],
 213             mv = s->mod[lev][1],
 214             lp;
 215         int *l;
 216
 217         // HOR_SD
 218         l = line + mh;
 219         for (lp = 0; lp < lv; lp++){
 220             int i, j = 0;
 221             // copy with interleaving
 222             for (i =   mh; i < lh; i+=2, j++)
 223                 l[i] = t[w*lp + j];
 224             for (i = 1-mh; i < lh; i+=2, j++)
 225                 l[i] = t[w*lp + j];
 226
 227             sr_1d53(line, mh, mh + lh);
 228
 229             for (i = 0; i < lh; i++)
 230                 t[w*lp + i] = l[i];
 231         }
 232
 233         // VER_SD
 234         l = line + mv;
 235         for (lp = 0; lp < lh; lp++){
 236             int i, j = 0;
 237             // copy with interleaving
 238             for (i =   mv; i < lv; i+=2, j++)
 239                 l[i] = t[w*j + lp];
 240             for (i = 1-mv; i < lv; i+=2, j++)
 241                 l[i] = t[w*j + lp];
 242
 243             sr_1d53(line, mv, mv + lv);
 244
 245             for (i = 0; i < lv; i++)
 246                 t[w*i + lp] = l[i];
 247         }
 248     }
 249 }
 250
 251 static void sr_1d97(float *p, int i0, int i1)
 252 {
 253     int i;
 254
 255     if (i1 == i0 + 1)
 256         return;
 257
 258     extend97(p, i0, i1);
 259
 260     for (i = i0/2 - 1; i < i1/2 + 2; i++)
 261         p[2*i] -= 0.443506 * (p[2*i-1] + p[2*i+1]);
 262     for (i = i0/2 - 1; i < i1/2 + 1; i++)
 263         p[2*i+1] -= 0.882911 * (p[2*i] + p[2*i+2]);
 264     for (i = i0/2; i < i1/2 + 1; i++)
 265         p[2*i] += 0.052980 * (p[2*i-1] + p[2*i+1]);
 266     for (i = i0/2; i < i1/2; i++)
 267         p[2*i+1] += 1.586134 * (p[2*i] + p[2*i+2]);
 268 }
 269
 270 static void dwt_decode97(DWTContext *s, int *t)
 271 {
 272     int lev,
 273         w = s->linelen[s->ndeclevels-1][0];
 274     float *line = s->linebuf;
 275     line += 5;
 276
 277     for (lev = 0; lev < s->ndeclevels; lev++){
 278         int lh = s->linelen[lev][0],
 279             lv = s->linelen[lev][1],
 280             mh = s->mod[lev][0],
 281             mv = s->mod[lev][1],
 282             lp;
 283         float *l;
 284
 285         // HOR_SD
 286         l = line + mh;
 287         for (lp = 0; lp < lv; lp++){
 288             int i, j = 0;
 289             // copy with interleaving
 290             for (i =   mh; i < lh; i+=2, j++)
 291                 l[i] = scale97[1-mh] * t[w*lp + j];
 292             for (i = 1-mh; i < lh; i+=2, j++)
 293                 l[i] = scale97[1-mh] * t[w*lp + j];
 294
 295             sr_1d97(line, mh, mh + lh);
 296
 297             for (i = 0; i < lh; i++)
 298                 t[w*lp + i] = l[i];
 299         }
 300
 301         // VER_SD
 302         l = line + mv;
 303         for (lp = 0; lp < lh; lp++){
 304             int i, j = 0;
 305             // copy with interleaving
 306             for (i =   mv; i < lv; i+=2, j++)
 307                 l[i] = scale97[1-mv] * t[w*j + lp];
 308             for (i = 1-mv; i < lv; i+=2, j++)
 309                 l[i] = scale97[1-mv] * t[w*j + lp];
 310
 311             sr_1d97(line, mv, mv + lv);
 312
 313             for (i = 0; i < lv; i++)
 314                 t[w*i + lp] = l[i];
 315         }
 316     }
 317 }
 318
 319 int ff_j2k_dwt_init(DWTContext *s, uint16_t border[2][2], int decomp_levels, int type)
 320 {
 321     int i, j, lev = decomp_levels, maxlen,
 322         b[2][2];
 323
 324     if (decomp_levels >= FF_DWT_MAX_DECLVLS)
 325         return AVERROR_INVALIDDATA;
 326     s->ndeclevels = decomp_levels;
 327     s->type = type;
 328
 329     for (i = 0; i < 2; i++)
 330         for(j = 0; j < 2; j++)
 331             b[i][j] = border[i][j];
 332
 333     maxlen = FFMAX(b[0][1] - b[0][0],
 334                    b[1][1] - b[1][0]);
 335
 336     while(--lev >= 0){
 337         for (i = 0; i < 2; i++){
 338             s->linelen[lev][i] = b[i][1] - b[i][0];
 339             s->mod[lev][i] = b[i][0] & 1;
 340             for (j = 0; j < 2; j++)
 341                 b[i][j] = (b[i][j] + 1) >> 1;
 342         }
 343     }
 344     if (type == FF_DWT97)
 345         s->linebuf = av_malloc((maxlen + 12) * sizeof(float));
 346     else if (type == FF_DWT53)
 347         s->linebuf = av_malloc((maxlen + 6) * sizeof(int));
 348     else
 349         return -1;
 350
 351     if (!s->linebuf)
 352         return AVERROR(ENOMEM);
 353
 354     return 0;
 355 }
 356
 357 int ff_j2k_dwt_encode(DWTContext *s, int *t)
 358 {
 359     switch(s->type){
 360         case FF_DWT97:
 361             dwt_encode97(s, t); break;
 362         case FF_DWT53:
 363             dwt_encode53(s, t); break;
 364         default:
 365             return -1;
 366     }
 367     return 0;
 368 }
 369
 370 int ff_j2k_dwt_decode(DWTContext *s, int *t)
 371 {
 372     switch(s->type){
 373         case FF_DWT97:
 374             dwt_decode97(s, t); break;
 375         case FF_DWT53:
 376             dwt_decode53(s, t); break;
 377         default:
 378             return -1;
 379     }
 380     return 0;
 381 }
 382
 383 void ff_j2k_dwt_destroy(DWTContext *s)
 384 {
 385     av_freep(&s->linebuf);
 386 }