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