jpeg2000dwt.c 16.4 KB
Newer Older
1 2 3 4 5
/*
 * Discrete wavelet transform
 * Copyright (c) 2007 Kamil Nowosad
 * Copyright (c) 2013 Nicolas Bertrand <nicoinattendu@gmail.com>
 *
6
 * This file is part of FFmpeg.
7
 *
8
 * FFmpeg is free software; you can redistribute it and/or
9 10 11 12
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
13
 * FFmpeg is distributed in the hope that it will be useful,
14 15 16 17 18
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * 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 21 22 23 24 25 26 27
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */

/**
 * @file
 * Discrete wavelet transform
 */

28
#include "libavutil/avassert.h"
29 30 31 32 33 34 35 36 37 38 39 40 41 42
#include "libavutil/common.h"
#include "libavutil/mem.h"
#include "jpeg2000dwt.h"
#include "internal.h"

/* Defines for 9/7 DWT lifting parameters.
 * Parameters are in float. */
#define F_LFTG_ALPHA  1.586134342059924f
#define F_LFTG_BETA   0.052980118572961f
#define F_LFTG_GAMMA  0.882911075530934f
#define F_LFTG_DELTA  0.443506852043971f

/* Lifting parameters in integer format.
 * Computed as param = (float param) * (1 << 16) */
43 44 45 46 47 48 49
#define I_LFTG_ALPHA  103949ll
#define I_LFTG_BETA     3472ll
#define I_LFTG_GAMMA   57862ll
#define I_LFTG_DELTA   29066ll
#define I_LFTG_K       80621ll
#define I_LFTG_X       53274ll
#define I_PRESHIFT 8
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78

static inline void extend53(int *p, int i0, int i1)
{
    p[i0 - 1] = p[i0 + 1];
    p[i1]     = p[i1 - 2];
    p[i0 - 2] = p[i0 + 2];
    p[i1 + 1] = p[i1 - 3];
}

static inline void extend97_float(float *p, int i0, int i1)
{
    int i;

    for (i = 1; i <= 4; i++) {
        p[i0 - i]     = p[i0 + i];
        p[i1 + i - 1] = p[i1 - i - 1];
    }
}

static inline void extend97_int(int32_t *p, int i0, int i1)
{
    int i;

    for (i = 1; i <= 4; i++) {
        p[i0 - i]     = p[i0 + i];
        p[i1 + i - 1] = p[i1 - i - 1];
    }
}

79 80 81 82
static void sd_1d53(int *p, int i0, int i1)
{
    int i;

83 84 85
    if (i1 <= i0 + 1) {
        if (i0 == 1)
            p[1] <<= 1;
86
        return;
87
    }
88 89 90

    extend53(p, i0, i1);

91
    for (i = ((i0+1)>>1) - 1; i < (i1+1)>>1; i++)
92
        p[2*i+1] -= (p[2*i] + p[2*i+2]) >> 1;
93
    for (i = ((i0+1)>>1); i < (i1+1)>>1; i++)
94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
        p[2*i] += (p[2*i-1] + p[2*i+1] + 2) >> 2;
}

static void dwt_encode53(DWTContext *s, int *t)
{
    int lev,
        w = s->linelen[s->ndeclevels-1][0];
    int *line = s->i_linebuf;
    line += 3;

    for (lev = s->ndeclevels-1; lev >= 0; lev--){
        int lh = s->linelen[lev][0],
            lv = s->linelen[lev][1],
            mh = s->mod[lev][0],
            mv = s->mod[lev][1],
            lp;
        int *l;

        // VER_SD
        l = line + mv;
        for (lp = 0; lp < lh; lp++) {
            int i, j = 0;

            for (i = 0; i < lv; i++)
                l[i] = t[w*i + lp];

            sd_1d53(line, mv, mv + lv);

            // copy back and deinterleave
            for (i =   mv; i < lv; i+=2, j++)
                t[w*j + lp] = l[i];
            for (i = 1-mv; i < lv; i+=2, j++)
                t[w*j + lp] = l[i];
        }
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144

        // HOR_SD
        l = line + mh;
        for (lp = 0; lp < lv; lp++){
            int i, j = 0;

            for (i = 0; i < lh; i++)
                l[i] = t[w*lp + i];

            sd_1d53(line, mh, mh + lh);

            // copy back and deinterleave
            for (i =   mh; i < lh; i+=2, j++)
                t[w*lp + j] = l[i];
            for (i = 1-mh; i < lh; i+=2, j++)
                t[w*lp + j] = l[i];
        }
145 146
    }
}
147 148 149 150
static void sd_1d97_float(float *p, int i0, int i1)
{
    int i;

151 152
    if (i1 <= i0 + 1) {
        if (i0 == 1)
153
            p[1] *= F_LFTG_X * 2;
154 155
        else
            p[0] *= F_LFTG_K;
156
        return;
157
    }
158 159 160 161

    extend97_float(p, i0, i1);
    i0++; i1++;

162
    for (i = (i0>>1) - 2; i < (i1>>1) + 1; i++)
163
        p[2*i+1] -= 1.586134 * (p[2*i] + p[2*i+2]);
164
    for (i = (i0>>1) - 1; i < (i1>>1) + 1; i++)
165
        p[2*i] -= 0.052980 * (p[2*i-1] + p[2*i+1]);
166
    for (i = (i0>>1) - 1; i < (i1>>1); i++)
167
        p[2*i+1] += 0.882911 * (p[2*i] + p[2*i+2]);
168
    for (i = (i0>>1); i < (i1>>1); i++)
169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
        p[2*i] += 0.443506 * (p[2*i-1] + p[2*i+1]);
}

static void dwt_encode97_float(DWTContext *s, float *t)
{
    int lev,
        w = s->linelen[s->ndeclevels-1][0];
    float *line = s->f_linebuf;
    line += 5;

    for (lev = s->ndeclevels-1; lev >= 0; lev--){
        int lh = s->linelen[lev][0],
            lv = s->linelen[lev][1],
            mh = s->mod[lev][0],
            mv = s->mod[lev][1],
            lp;
        float *l;

        // HOR_SD
        l = line + mh;
        for (lp = 0; lp < lv; lp++){
            int i, j = 0;

            for (i = 0; i < lh; i++)
                l[i] = t[w*lp + i];

            sd_1d97_float(line, mh, mh + lh);

            // copy back and deinterleave
            for (i =   mh; i < lh; i+=2, j++)
199
                t[w*lp + j] = l[i];
200
            for (i = 1-mh; i < lh; i+=2, j++)
201
                t[w*lp + j] = l[i];
202 203 204 205 206 207 208 209 210 211 212 213 214 215
        }

        // VER_SD
        l = line + mv;
        for (lp = 0; lp < lh; lp++) {
            int i, j = 0;

            for (i = 0; i < lv; i++)
                l[i] = t[w*i + lp];

            sd_1d97_float(line, mv, mv + lv);

            // copy back and deinterleave
            for (i =   mv; i < lv; i+=2, j++)
216
                t[w*j + lp] = l[i];
217
            for (i = 1-mv; i < lv; i+=2, j++)
218
                t[w*j + lp] = l[i];
219 220 221
        }
    }
}
222

223
static void sd_1d97_int(int *p, int i0, int i1)
224 225 226
{
    int i;

227 228
    if (i1 <= i0 + 1) {
        if (i0 == 1)
229
            p[1] = (p[1] * I_LFTG_X + (1<<14)) >> 15;
230 231
        else
            p[0] = (p[0] * I_LFTG_K + (1<<15)) >> 16;
232
        return;
233
    }
234

235
    extend97_int(p, i0, i1);
236 237
    i0++; i1++;

238
    for (i = (i0>>1) - 2; i < (i1>>1) + 1; i++)
239
        p[2 * i + 1] -= (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
240
    for (i = (i0>>1) - 1; i < (i1>>1) + 1; i++)
241
        p[2 * i]     -= (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
242
    for (i = (i0>>1) - 1; i < (i1>>1); i++)
243
        p[2 * i + 1] += (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
244
    for (i = (i0>>1); i < (i1>>1); i++)
245
        p[2 * i]     += (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
246 247 248 249
}

static void dwt_encode97_int(DWTContext *s, int *t)
{
250 251 252 253
    int lev;
    int w = s->linelen[s->ndeclevels-1][0];
    int h = s->linelen[s->ndeclevels-1][1];
    int i;
254
    int *line = s->i_linebuf;
255 256
    line += 5;

257
    for (i = 0; i < w * h; i++)
258
        t[i] *= 1 << I_PRESHIFT;
259

260 261 262 263 264 265
    for (lev = s->ndeclevels-1; lev >= 0; lev--){
        int lh = s->linelen[lev][0],
            lv = s->linelen[lev][1],
            mh = s->mod[lev][0],
            mv = s->mod[lev][1],
            lp;
266
        int *l;
267

268 269 270 271 272 273 274 275 276 277 278 279
        // VER_SD
        l = line + mv;
        for (lp = 0; lp < lh; lp++) {
            int i, j = 0;

            for (i = 0; i < lv; i++)
                l[i] = t[w*i + lp];

            sd_1d97_int(line, mv, mv + lv);

            // copy back and deinterleave
            for (i =   mv; i < lv; i+=2, j++)
280
                t[w*j + lp] = ((l[i] * I_LFTG_X) + (1 << 15)) >> 16;
281
            for (i = 1-mv; i < lv; i+=2, j++)
282
                t[w*j + lp] = l[i];
283 284
        }

285 286 287 288 289 290 291 292
        // HOR_SD
        l = line + mh;
        for (lp = 0; lp < lv; lp++){
            int i, j = 0;

            for (i = 0; i < lh; i++)
                l[i] = t[w*lp + i];

293
            sd_1d97_int(line, mh, mh + lh);
294 295 296

            // copy back and deinterleave
            for (i =   mh; i < lh; i+=2, j++)
297
                t[w*lp + j] = ((l[i] * I_LFTG_X) + (1 << 15)) >> 16;
298
            for (i = 1-mh; i < lh; i+=2, j++)
299
                t[w*lp + j] = l[i];
300 301 302
        }

    }
303 304 305

    for (i = 0; i < w * h; i++)
        t[i] = (t[i] + ((1<<I_PRESHIFT)>>1)) >> I_PRESHIFT;
306 307
}

308
static void sr_1d53(unsigned *p, int i0, int i1)
309 310 311
{
    int i;

312 313
    if (i1 <= i0 + 1) {
        if (i0 == 1)
314
            p[1] = (int)p[1] >> 1;
315
        return;
316
    }
317 318 319

    extend53(p, i0, i1);

320
    for (i = (i0 >> 1); i < (i1 >> 1) + 1; i++)
321
        p[2 * i] -= (int)(p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
322
    for (i = (i0 >> 1); i < (i1 >> 1); i++)
323
        p[2 * i + 1] += (int)(p[2 * i] + p[2 * i + 2]) >> 1;
324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378
}

static void dwt_decode53(DWTContext *s, int *t)
{
    int lev;
    int w     = s->linelen[s->ndeclevels - 1][0];
    int32_t *line = s->i_linebuf;
    line += 3;

    for (lev = 0; lev < s->ndeclevels; lev++) {
        int lh = s->linelen[lev][0],
            lv = s->linelen[lev][1],
            mh = s->mod[lev][0],
            mv = s->mod[lev][1],
            lp;
        int *l;

        // HOR_SD
        l = line + mh;
        for (lp = 0; lp < lv; lp++) {
            int i, j = 0;
            // copy with interleaving
            for (i = mh; i < lh; i += 2, j++)
                l[i] = t[w * lp + j];
            for (i = 1 - mh; i < lh; i += 2, j++)
                l[i] = t[w * lp + j];

            sr_1d53(line, mh, mh + lh);

            for (i = 0; i < lh; i++)
                t[w * lp + i] = l[i];
        }

        // VER_SD
        l = line + mv;
        for (lp = 0; lp < lh; lp++) {
            int i, j = 0;
            // copy with interleaving
            for (i = mv; i < lv; i += 2, j++)
                l[i] = t[w * j + lp];
            for (i = 1 - mv; i < lv; i += 2, j++)
                l[i] = t[w * j + lp];

            sr_1d53(line, mv, mv + lv);

            for (i = 0; i < lv; i++)
                t[w * i + lp] = l[i];
        }
    }
}

static void sr_1d97_float(float *p, int i0, int i1)
{
    int i;

379 380 381
    if (i1 <= i0 + 1) {
        if (i0 == 1)
            p[1] *= F_LFTG_K/2;
382
        else
383
            p[0] *= F_LFTG_X;
384
        return;
385
    }
386 387 388

    extend97_float(p, i0, i1);

389
    for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 2; i++)
390 391
        p[2 * i]     -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
    /* step 4 */
392
    for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 1; i++)
393 394
        p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]);
    /*step 5*/
395
    for (i = (i0 >> 1); i < (i1 >> 1) + 1; i++)
396 397
        p[2 * i]     += F_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]);
    /* step 6 */
398
    for (i = (i0 >> 1); i < (i1 >> 1); i++)
399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423
        p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]);
}

static void dwt_decode97_float(DWTContext *s, float *t)
{
    int lev;
    int w       = s->linelen[s->ndeclevels - 1][0];
    float *line = s->f_linebuf;
    float *data = t;
    /* position at index O of line range [0-5,w+5] cf. extend function */
    line += 5;

    for (lev = 0; lev < s->ndeclevels; lev++) {
        int lh = s->linelen[lev][0],
            lv = s->linelen[lev][1],
            mh = s->mod[lev][0],
            mv = s->mod[lev][1],
            lp;
        float *l;
        // HOR_SD
        l = line + mh;
        for (lp = 0; lp < lv; lp++) {
            int i, j = 0;
            // copy with interleaving
            for (i = mh; i < lh; i += 2, j++)
424
                l[i] = data[w * lp + j];
425
            for (i = 1 - mh; i < lh; i += 2, j++)
426
                l[i] = data[w * lp + j];
427 428 429 430 431 432 433 434 435 436 437 438 439

            sr_1d97_float(line, mh, mh + lh);

            for (i = 0; i < lh; i++)
                data[w * lp + i] = l[i];
        }

        // VER_SD
        l = line + mv;
        for (lp = 0; lp < lh; lp++) {
            int i, j = 0;
            // copy with interleaving
            for (i = mv; i < lv; i += 2, j++)
440
                l[i] = data[w * j + lp];
441
            for (i = 1 - mv; i < lv; i += 2, j++)
442
                l[i] = data[w * j + lp];
443 444 445 446 447 448 449 450 451 452 453 454 455

            sr_1d97_float(line, mv, mv + lv);

            for (i = 0; i < lv; i++)
                data[w * i + lp] = l[i];
        }
    }
}

static void sr_1d97_int(int32_t *p, int i0, int i1)
{
    int i;

456 457 458
    if (i1 <= i0 + 1) {
        if (i0 == 1)
            p[1] = (p[1] * I_LFTG_K + (1<<16)) >> 17;
459
        else
460
            p[0] = (p[0] * I_LFTG_X + (1<<15)) >> 16;
461
        return;
462
    }
463 464 465

    extend97_int(p, i0, i1);

466
    for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 2; i++)
467
        p[2 * i]     -= (I_LFTG_DELTA * (p[2 * i - 1] + (int64_t)p[2 * i + 1]) + (1 << 15)) >> 16;
468
    /* step 4 */
469
    for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 1; i++)
470
        p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i]     + (int64_t)p[2 * i + 2]) + (1 << 15)) >> 16;
471
    /*step 5*/
472
    for (i = (i0 >> 1); i < (i1 >> 1) + 1; i++)
473
        p[2 * i]     += (I_LFTG_BETA  * (p[2 * i - 1] + (int64_t)p[2 * i + 1]) + (1 << 15)) >> 16;
474
    /* step 6 */
475
    for (i = (i0 >> 1); i < (i1 >> 1); i++)
476
        p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i]     + (int64_t)p[2 * i + 2]) + (1 << 15)) >> 16;
477 478 479 480 481 482
}

static void dwt_decode97_int(DWTContext *s, int32_t *t)
{
    int lev;
    int w       = s->linelen[s->ndeclevels - 1][0];
483 484
    int h       = s->linelen[s->ndeclevels - 1][1];
    int i;
485 486 487 488 489
    int32_t *line = s->i_linebuf;
    int32_t *data = t;
    /* position at index O of line range [0-5,w+5] cf. extend function */
    line += 5;

490
    for (i = 0; i < w * h; i++)
491
        data[i] *= 1LL << I_PRESHIFT;
492

493 494 495 496 497 498 499 500 501 502 503
    for (lev = 0; lev < s->ndeclevels; lev++) {
        int lh = s->linelen[lev][0],
            lv = s->linelen[lev][1],
            mh = s->mod[lev][0],
            mv = s->mod[lev][1],
            lp;
        int32_t *l;
        // HOR_SD
        l = line + mh;
        for (lp = 0; lp < lv; lp++) {
            int i, j = 0;
504
            // rescale with interleaving
505
            for (i = mh; i < lh; i += 2, j++)
506
                l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
507
            for (i = 1 - mh; i < lh; i += 2, j++)
508
                l[i] = data[w * lp + j];
509 510 511 512 513 514 515 516 517 518 519

            sr_1d97_int(line, mh, mh + lh);

            for (i = 0; i < lh; i++)
                data[w * lp + i] = l[i];
        }

        // VER_SD
        l = line + mv;
        for (lp = 0; lp < lh; lp++) {
            int i, j = 0;
520
            // rescale with interleaving
521
            for (i = mv; i < lv; i += 2, j++)
522
                l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
523
            for (i = 1 - mv; i < lv; i += 2, j++)
524
                l[i] = data[w * j + lp];
525 526 527 528 529 530 531

            sr_1d97_int(line, mv, mv + lv);

            for (i = 0; i < lv; i++)
                data[w * i + lp] = l[i];
        }
    }
532 533

    for (i = 0; i < w * h; i++)
534
        data[i] = (data[i] + ((1LL<<I_PRESHIFT)>>1)) >> I_PRESHIFT;
535 536
}

537
int ff_jpeg2000_dwt_init(DWTContext *s, int border[2][2],
538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560
                         int decomp_levels, int type)
{
    int i, j, lev = decomp_levels, maxlen,
        b[2][2];

    s->ndeclevels = decomp_levels;
    s->type       = type;

    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
            b[i][j] = border[i][j];

    maxlen = FFMAX(b[0][1] - b[0][0],
                   b[1][1] - b[1][0]);
    while (--lev >= 0)
        for (i = 0; i < 2; i++) {
            s->linelen[lev][i] = b[i][1] - b[i][0];
            s->mod[lev][i]     = b[i][0] & 1;
            for (j = 0; j < 2; j++)
                b[i][j] = (b[i][j] + 1) >> 1;
        }
    switch (type) {
    case FF_DWT97:
561
        s->f_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->f_linebuf));
562 563 564 565
        if (!s->f_linebuf)
            return AVERROR(ENOMEM);
        break;
     case FF_DWT97_INT:
566
        s->i_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->i_linebuf));
567 568 569 570
        if (!s->i_linebuf)
            return AVERROR(ENOMEM);
        break;
    case FF_DWT53:
571
        s->i_linebuf = av_malloc_array((maxlen +  6), sizeof(*s->i_linebuf));
572 573 574 575 576 577 578 579 580
        if (!s->i_linebuf)
            return AVERROR(ENOMEM);
        break;
    default:
        return -1;
    }
    return 0;
}

581
int ff_dwt_encode(DWTContext *s, void *t)
582
{
583 584 585
    if (s->ndeclevels == 0)
        return 0;

586
    switch(s->type){
587 588
        case FF_DWT97:
            dwt_encode97_float(s, t); break;
589 590 591 592 593 594 595 596 597 598
        case FF_DWT97_INT:
            dwt_encode97_int(s, t); break;
        case FF_DWT53:
            dwt_encode53(s, t); break;
        default:
            return -1;
    }
    return 0;
}

599 600
int ff_dwt_decode(DWTContext *s, void *t)
{
601 602 603
    if (s->ndeclevels == 0)
        return 0;

604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624
    switch (s->type) {
    case FF_DWT97:
        dwt_decode97_float(s, t);
        break;
    case FF_DWT97_INT:
        dwt_decode97_int(s, t);
        break;
    case FF_DWT53:
        dwt_decode53(s, t);
        break;
    default:
        return -1;
    }
    return 0;
}

void ff_dwt_destroy(DWTContext *s)
{
    av_freep(&s->f_linebuf);
    av_freep(&s->i_linebuf);
}