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