1/* 2 * The copyright in this software is being made available under the 2-clauses 3 * BSD License, included below. This software may be subject to other third 4 * party and contributor rights, including patent rights, and no such rights 5 * are granted under this license. 6 * 7 * Copyright (c) 2002-2014, Universite catholique de Louvain (UCL), Belgium 8 * Copyright (c) 2002-2014, Professor Benoit Macq 9 * Copyright (c) 2001-2003, David Janssens 10 * Copyright (c) 2002-2003, Yannick Verschueren 11 * Copyright (c) 2003-2007, Francois-Olivier Devaux 12 * Copyright (c) 2003-2014, Antonin Descampe 13 * Copyright (c) 2005, Herve Drolon, FreeImage Team 14 * Copyright (c) 2007, Jonathan Ballard <dzonatas@dzonux.net> 15 * Copyright (c) 2007, Callum Lerwick <seg@haxxed.com> 16 * All rights reserved. 17 * 18 * Redistribution and use in source and binary forms, with or without 19 * modification, are permitted provided that the following conditions 20 * are met: 21 * 1. Redistributions of source code must retain the above copyright 22 * notice, this list of conditions and the following disclaimer. 23 * 2. Redistributions in binary form must reproduce the above copyright 24 * notice, this list of conditions and the following disclaimer in the 25 * documentation and/or other materials provided with the distribution. 26 * 27 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS' 28 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 29 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 30 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 31 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 32 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 33 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 34 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 35 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 36 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 37 * POSSIBILITY OF SUCH DAMAGE. 38 */ 39 40#ifdef __SSE__ 41#include <xmmintrin.h> 42#endif 43 44#include "opj_includes.h" 45 46/** @defgroup DWT DWT - Implementation of a discrete wavelet transform */ 47/*@{*/ 48 49#define OPJ_WS(i) v->mem[(i)*2] 50#define OPJ_WD(i) v->mem[(1+(i)*2)] 51 52/** @name Local data structures */ 53/*@{*/ 54 55typedef struct dwt_local { 56 OPJ_INT32* mem; 57 OPJ_INT32 dn; 58 OPJ_INT32 sn; 59 OPJ_INT32 cas; 60} opj_dwt_t; 61 62typedef union { 63 OPJ_FLOAT32 f[4]; 64} opj_v4_t; 65 66typedef struct v4dwt_local { 67 opj_v4_t* wavelet ; 68 OPJ_INT32 dn ; 69 OPJ_INT32 sn ; 70 OPJ_INT32 cas ; 71} opj_v4dwt_t ; 72 73static const OPJ_FLOAT32 opj_dwt_alpha = 1.586134342f; /* 12994 */ 74static const OPJ_FLOAT32 opj_dwt_beta = 0.052980118f; /* 434 */ 75static const OPJ_FLOAT32 opj_dwt_gamma = -0.882911075f; /* -7233 */ 76static const OPJ_FLOAT32 opj_dwt_delta = -0.443506852f; /* -3633 */ 77 78static const OPJ_FLOAT32 opj_K = 1.230174105f; /* 10078 */ 79static const OPJ_FLOAT32 opj_c13318 = 1.625732422f; 80 81/*@}*/ 82 83/** 84Virtual function type for wavelet transform in 1-D 85*/ 86typedef void (*DWT1DFN)(opj_dwt_t* v); 87 88/** @name Local static functions */ 89/*@{*/ 90 91/** 92Forward lazy transform (horizontal) 93*/ 94static void opj_dwt_deinterleave_h(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn, OPJ_INT32 sn, OPJ_INT32 cas); 95/** 96Forward lazy transform (vertical) 97*/ 98static void opj_dwt_deinterleave_v(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn, OPJ_INT32 sn, OPJ_INT32 x, OPJ_INT32 cas); 99/** 100Inverse lazy transform (horizontal) 101*/ 102static void opj_dwt_interleave_h(opj_dwt_t* h, OPJ_INT32 *a); 103/** 104Inverse lazy transform (vertical) 105*/ 106static void opj_dwt_interleave_v(opj_dwt_t* v, OPJ_INT32 *a, OPJ_INT32 x); 107/** 108Forward 5-3 wavelet transform in 1-D 109*/ 110static void opj_dwt_encode_1(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn, OPJ_INT32 cas); 111/** 112Inverse 5-3 wavelet transform in 1-D 113*/ 114static void opj_dwt_decode_1(opj_dwt_t *v); 115static void opj_dwt_decode_1_(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn, OPJ_INT32 cas); 116/** 117Forward 9-7 wavelet transform in 1-D 118*/ 119static void opj_dwt_encode_1_real(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn, OPJ_INT32 cas); 120/** 121Explicit calculation of the Quantization Stepsizes 122*/ 123static void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps, opj_stepsize_t *bandno_stepsize); 124/** 125Inverse wavelet transform in 2-D. 126*/ 127static OPJ_BOOL opj_dwt_decode_tile(opj_tcd_tilecomp_t* tilec, OPJ_UINT32 i, DWT1DFN fn); 128 129static OPJ_BOOL opj_dwt_encode_procedure( opj_tcd_tilecomp_t * tilec, 130 void (*p_function)(OPJ_INT32 *, OPJ_INT32,OPJ_INT32,OPJ_INT32) ); 131 132static OPJ_UINT32 opj_dwt_max_resolution(opj_tcd_resolution_t* restrict r, OPJ_UINT32 i); 133 134/* <summary> */ 135/* Inverse 9-7 wavelet transform in 1-D. */ 136/* </summary> */ 137static void opj_v4dwt_decode(opj_v4dwt_t* restrict dwt); 138 139static void opj_v4dwt_interleave_h(opj_v4dwt_t* restrict w, OPJ_FLOAT32* restrict a, OPJ_INT32 x, OPJ_INT32 size); 140 141static void opj_v4dwt_interleave_v(opj_v4dwt_t* restrict v , OPJ_FLOAT32* restrict a , OPJ_INT32 x, OPJ_INT32 nb_elts_read); 142 143#ifdef __SSE__ 144static void opj_v4dwt_decode_step1_sse(opj_v4_t* w, OPJ_INT32 count, const __m128 c); 145 146static void opj_v4dwt_decode_step2_sse(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k, OPJ_INT32 m, __m128 c); 147 148#else 149static void opj_v4dwt_decode_step1(opj_v4_t* w, OPJ_INT32 count, const OPJ_FLOAT32 c); 150 151static void opj_v4dwt_decode_step2(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k, OPJ_INT32 m, OPJ_FLOAT32 c); 152 153#endif 154 155/*@}*/ 156 157/*@}*/ 158 159#define OPJ_S(i) a[(i)*2] 160#define OPJ_D(i) a[(1+(i)*2)] 161#define OPJ_S_(i) ((i)<0?OPJ_S(0):((i)>=sn?OPJ_S(sn-1):OPJ_S(i))) 162#define OPJ_D_(i) ((i)<0?OPJ_D(0):((i)>=dn?OPJ_D(dn-1):OPJ_D(i))) 163/* new */ 164#define OPJ_SS_(i) ((i)<0?OPJ_S(0):((i)>=dn?OPJ_S(dn-1):OPJ_S(i))) 165#define OPJ_DD_(i) ((i)<0?OPJ_D(0):((i)>=sn?OPJ_D(sn-1):OPJ_D(i))) 166 167/* <summary> */ 168/* This table contains the norms of the 5-3 wavelets for different bands. */ 169/* </summary> */ 170static const OPJ_FLOAT64 opj_dwt_norms[4][10] = { 171 {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3}, 172 {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9}, 173 {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9}, 174 {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93} 175}; 176 177/* <summary> */ 178/* This table contains the norms of the 9-7 wavelets for different bands. */ 179/* </summary> */ 180static const OPJ_FLOAT64 opj_dwt_norms_real[4][10] = { 181 {1.000, 1.965, 4.177, 8.403, 16.90, 33.84, 67.69, 135.3, 270.6, 540.9}, 182 {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0}, 183 {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0}, 184 {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2} 185}; 186 187/* 188========================================================== 189 local functions 190========================================================== 191*/ 192 193/* <summary> */ 194/* Forward lazy transform (horizontal). */ 195/* </summary> */ 196void opj_dwt_deinterleave_h(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn, OPJ_INT32 sn, OPJ_INT32 cas) { 197 OPJ_INT32 i; 198 OPJ_INT32 * l_dest = b; 199 OPJ_INT32 * l_src = a+cas; 200 201 for (i=0; i<sn; ++i) { 202 *l_dest++ = *l_src; 203 l_src += 2; 204 } 205 206 l_dest = b + sn; 207 l_src = a + 1 - cas; 208 209 for (i=0; i<dn; ++i) { 210 *l_dest++=*l_src; 211 l_src += 2; 212 } 213} 214 215/* <summary> */ 216/* Forward lazy transform (vertical). */ 217/* </summary> */ 218void opj_dwt_deinterleave_v(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn, OPJ_INT32 sn, OPJ_INT32 x, OPJ_INT32 cas) { 219 OPJ_INT32 i = sn; 220 OPJ_INT32 * l_dest = b; 221 OPJ_INT32 * l_src = a+cas; 222 223 while (i--) { 224 *l_dest = *l_src; 225 l_dest += x; 226 l_src += 2; 227 } /* b[i*x]=a[2*i+cas]; */ 228 229 l_dest = b + sn * x; 230 l_src = a + 1 - cas; 231 232 i = dn; 233 while (i--) { 234 *l_dest = *l_src; 235 l_dest += x; 236 l_src += 2; 237 } /*b[(sn+i)*x]=a[(2*i+1-cas)];*/ 238} 239 240/* <summary> */ 241/* Inverse lazy transform (horizontal). */ 242/* </summary> */ 243void opj_dwt_interleave_h(opj_dwt_t* h, OPJ_INT32 *a) { 244 OPJ_INT32 *ai = a; 245 OPJ_INT32 *bi = h->mem + h->cas; 246 OPJ_INT32 i = h->sn; 247 while( i-- ) { 248 *bi = *(ai++); 249 bi += 2; 250 } 251 ai = a + h->sn; 252 bi = h->mem + 1 - h->cas; 253 i = h->dn ; 254 while( i-- ) { 255 *bi = *(ai++); 256 bi += 2; 257 } 258} 259 260/* <summary> */ 261/* Inverse lazy transform (vertical). */ 262/* </summary> */ 263void opj_dwt_interleave_v(opj_dwt_t* v, OPJ_INT32 *a, OPJ_INT32 x) { 264 OPJ_INT32 *ai = a; 265 OPJ_INT32 *bi = v->mem + v->cas; 266 OPJ_INT32 i = v->sn; 267 while( i-- ) { 268 *bi = *ai; 269 bi += 2; 270 ai += x; 271 } 272 ai = a + (v->sn * x); 273 bi = v->mem + 1 - v->cas; 274 i = v->dn ; 275 while( i-- ) { 276 *bi = *ai; 277 bi += 2; 278 ai += x; 279 } 280} 281 282 283/* <summary> */ 284/* Forward 5-3 wavelet transform in 1-D. */ 285/* </summary> */ 286void opj_dwt_encode_1(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn, OPJ_INT32 cas) { 287 OPJ_INT32 i; 288 289 if (!cas) { 290 if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */ 291 for (i = 0; i < dn; i++) OPJ_D(i) -= (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1; 292 for (i = 0; i < sn; i++) OPJ_S(i) += (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2; 293 } 294 } else { 295 if (!sn && dn == 1) /* NEW : CASE ONE ELEMENT */ 296 OPJ_S(0) *= 2; 297 else { 298 for (i = 0; i < dn; i++) OPJ_S(i) -= (OPJ_DD_(i) + OPJ_DD_(i - 1)) >> 1; 299 for (i = 0; i < sn; i++) OPJ_D(i) += (OPJ_SS_(i) + OPJ_SS_(i + 1) + 2) >> 2; 300 } 301 } 302} 303 304/* <summary> */ 305/* Inverse 5-3 wavelet transform in 1-D. */ 306/* </summary> */ 307void opj_dwt_decode_1_(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn, OPJ_INT32 cas) { 308 OPJ_INT32 i; 309 310 if (!cas) { 311 if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */ 312 for (i = 0; i < sn; i++) OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2; 313 for (i = 0; i < dn; i++) OPJ_D(i) += (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1; 314 } 315 } else { 316 if (!sn && dn == 1) /* NEW : CASE ONE ELEMENT */ 317 OPJ_S(0) /= 2; 318 else { 319 for (i = 0; i < sn; i++) OPJ_D(i) -= (OPJ_SS_(i) + OPJ_SS_(i + 1) + 2) >> 2; 320 for (i = 0; i < dn; i++) OPJ_S(i) += (OPJ_DD_(i) + OPJ_DD_(i - 1)) >> 1; 321 } 322 } 323} 324 325/* <summary> */ 326/* Inverse 5-3 wavelet transform in 1-D. */ 327/* </summary> */ 328void opj_dwt_decode_1(opj_dwt_t *v) { 329 opj_dwt_decode_1_(v->mem, v->dn, v->sn, v->cas); 330} 331 332/* <summary> */ 333/* Forward 9-7 wavelet transform in 1-D. */ 334/* </summary> */ 335void opj_dwt_encode_1_real(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn, OPJ_INT32 cas) { 336 OPJ_INT32 i; 337 if (!cas) { 338 if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */ 339 for (i = 0; i < dn; i++) 340 OPJ_D(i) -= opj_int_fix_mul(OPJ_S_(i) + OPJ_S_(i + 1), 12993); 341 for (i = 0; i < sn; i++) 342 OPJ_S(i) -= opj_int_fix_mul(OPJ_D_(i - 1) + OPJ_D_(i), 434); 343 for (i = 0; i < dn; i++) 344 OPJ_D(i) += opj_int_fix_mul(OPJ_S_(i) + OPJ_S_(i + 1), 7233); 345 for (i = 0; i < sn; i++) 346 OPJ_S(i) += opj_int_fix_mul(OPJ_D_(i - 1) + OPJ_D_(i), 3633); 347 for (i = 0; i < dn; i++) 348 OPJ_D(i) = opj_int_fix_mul(OPJ_D(i), 5038); /*5038 */ 349 for (i = 0; i < sn; i++) 350 OPJ_S(i) = opj_int_fix_mul(OPJ_S(i), 6659); /*6660 */ 351 } 352 } else { 353 if ((sn > 0) || (dn > 1)) { /* NEW : CASE ONE ELEMENT */ 354 for (i = 0; i < dn; i++) 355 OPJ_S(i) -= opj_int_fix_mul(OPJ_DD_(i) + OPJ_DD_(i - 1), 12993); 356 for (i = 0; i < sn; i++) 357 OPJ_D(i) -= opj_int_fix_mul(OPJ_SS_(i) + OPJ_SS_(i + 1), 434); 358 for (i = 0; i < dn; i++) 359 OPJ_S(i) += opj_int_fix_mul(OPJ_DD_(i) + OPJ_DD_(i - 1), 7233); 360 for (i = 0; i < sn; i++) 361 OPJ_D(i) += opj_int_fix_mul(OPJ_SS_(i) + OPJ_SS_(i + 1), 3633); 362 for (i = 0; i < dn; i++) 363 OPJ_S(i) = opj_int_fix_mul(OPJ_S(i), 5038); /*5038 */ 364 for (i = 0; i < sn; i++) 365 OPJ_D(i) = opj_int_fix_mul(OPJ_D(i), 6659); /*6660 */ 366 } 367 } 368} 369 370void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps, opj_stepsize_t *bandno_stepsize) { 371 OPJ_INT32 p, n; 372 p = opj_int_floorlog2(stepsize) - 13; 373 n = 11 - opj_int_floorlog2(stepsize); 374 bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff; 375 bandno_stepsize->expn = numbps - p; 376} 377 378/* 379========================================================== 380 DWT interface 381========================================================== 382*/ 383 384 385/* <summary> */ 386/* Forward 5-3 wavelet transform in 2-D. */ 387/* </summary> */ 388INLINE OPJ_BOOL opj_dwt_encode_procedure(opj_tcd_tilecomp_t * tilec,void (*p_function)(OPJ_INT32 *, OPJ_INT32,OPJ_INT32,OPJ_INT32) ) 389{ 390 OPJ_INT32 i, j, k; 391 OPJ_INT32 *a = 00; 392 OPJ_INT32 *aj = 00; 393 OPJ_INT32 *bj = 00; 394 OPJ_INT32 w, l; 395 396 OPJ_INT32 rw; /* width of the resolution level computed */ 397 OPJ_INT32 rh; /* height of the resolution level computed */ 398 OPJ_UINT32 l_data_size; 399 400 opj_tcd_resolution_t * l_cur_res = 0; 401 opj_tcd_resolution_t * l_last_res = 0; 402 403 w = tilec->x1-tilec->x0; 404 l = (OPJ_INT32)tilec->numresolutions-1; 405 a = tilec->data; 406 407 l_cur_res = tilec->resolutions + l; 408 l_last_res = l_cur_res - 1; 409 410 l_data_size = opj_dwt_max_resolution( tilec->resolutions,tilec->numresolutions) * (OPJ_UINT32)sizeof(OPJ_INT32); 411 bj = (OPJ_INT32*)opj_malloc((size_t)l_data_size); 412 if (! bj) { 413 return OPJ_FALSE; 414 } 415 i = l; 416 417 while (i--) { 418 OPJ_INT32 rw1; /* width of the resolution level once lower than computed one */ 419 OPJ_INT32 rh1; /* height of the resolution level once lower than computed one */ 420 OPJ_INT32 cas_col; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */ 421 OPJ_INT32 cas_row; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */ 422 OPJ_INT32 dn, sn; 423 424 rw = l_cur_res->x1 - l_cur_res->x0; 425 rh = l_cur_res->y1 - l_cur_res->y0; 426 rw1 = l_last_res->x1 - l_last_res->x0; 427 rh1 = l_last_res->y1 - l_last_res->y0; 428 429 cas_row = l_cur_res->x0 & 1; 430 cas_col = l_cur_res->y0 & 1; 431 432 sn = rh1; 433 dn = rh - rh1; 434 for (j = 0; j < rw; ++j) { 435 aj = a + j; 436 for (k = 0; k < rh; ++k) { 437 bj[k] = aj[k*w]; 438 } 439 440 (*p_function) (bj, dn, sn, cas_col); 441 442 opj_dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col); 443 } 444 445 sn = rw1; 446 dn = rw - rw1; 447 448 for (j = 0; j < rh; j++) { 449 aj = a + j * w; 450 for (k = 0; k < rw; k++) bj[k] = aj[k]; 451 (*p_function) (bj, dn, sn, cas_row); 452 opj_dwt_deinterleave_h(bj, aj, dn, sn, cas_row); 453 } 454 455 l_cur_res = l_last_res; 456 457 --l_last_res; 458 } 459 460 opj_free(bj); 461 return OPJ_TRUE; 462} 463 464/* Forward 5-3 wavelet transform in 2-D. */ 465/* </summary> */ 466OPJ_BOOL opj_dwt_encode(opj_tcd_tilecomp_t * tilec) 467{ 468 return opj_dwt_encode_procedure(tilec,opj_dwt_encode_1); 469} 470 471/* <summary> */ 472/* Inverse 5-3 wavelet transform in 2-D. */ 473/* </summary> */ 474OPJ_BOOL opj_dwt_decode(opj_tcd_tilecomp_t* tilec, OPJ_UINT32 numres) { 475 return opj_dwt_decode_tile(tilec, numres, &opj_dwt_decode_1); 476} 477 478 479/* <summary> */ 480/* Get gain of 5-3 wavelet transform. */ 481/* </summary> */ 482OPJ_UINT32 opj_dwt_getgain(OPJ_UINT32 orient) { 483 if (orient == 0) 484 return 0; 485 if (orient == 1 || orient == 2) 486 return 1; 487 return 2; 488} 489 490/* <summary> */ 491/* Get norm of 5-3 wavelet. */ 492/* </summary> */ 493OPJ_FLOAT64 opj_dwt_getnorm(OPJ_UINT32 level, OPJ_UINT32 orient) { 494 return opj_dwt_norms[orient][level]; 495} 496 497/* <summary> */ 498/* Forward 9-7 wavelet transform in 2-D. */ 499/* </summary> */ 500OPJ_BOOL opj_dwt_encode_real(opj_tcd_tilecomp_t * tilec) 501{ 502 return opj_dwt_encode_procedure(tilec,opj_dwt_encode_1_real); 503} 504 505/* <summary> */ 506/* Get gain of 9-7 wavelet transform. */ 507/* </summary> */ 508OPJ_UINT32 opj_dwt_getgain_real(OPJ_UINT32 orient) { 509 (void)orient; 510 return 0; 511} 512 513/* <summary> */ 514/* Get norm of 9-7 wavelet. */ 515/* </summary> */ 516OPJ_FLOAT64 opj_dwt_getnorm_real(OPJ_UINT32 level, OPJ_UINT32 orient) { 517 return opj_dwt_norms_real[orient][level]; 518} 519 520void opj_dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, OPJ_UINT32 prec) { 521 OPJ_UINT32 numbands, bandno; 522 numbands = 3 * tccp->numresolutions - 2; 523 for (bandno = 0; bandno < numbands; bandno++) { 524 OPJ_FLOAT64 stepsize; 525 OPJ_UINT32 resno, level, orient, gain; 526 527 resno = (bandno == 0) ? 0 : ((bandno - 1) / 3 + 1); 528 orient = (bandno == 0) ? 0 : ((bandno - 1) % 3 + 1); 529 level = tccp->numresolutions - 1 - resno; 530 gain = (tccp->qmfbid == 0) ? 0 : ((orient == 0) ? 0 : (((orient == 1) || (orient == 2)) ? 1 : 2)); 531 if (tccp->qntsty == J2K_CCP_QNTSTY_NOQNT) { 532 stepsize = 1.0; 533 } else { 534 OPJ_FLOAT64 norm = opj_dwt_norms_real[orient][level]; 535 stepsize = (1 << (gain)) / norm; 536 } 537 opj_dwt_encode_stepsize((OPJ_INT32) floor(stepsize * 8192.0), (OPJ_INT32)(prec + gain), &tccp->stepsizes[bandno]); 538 } 539} 540 541/* <summary> */ 542/* Determine maximum computed resolution level for inverse wavelet transform */ 543/* </summary> */ 544OPJ_UINT32 opj_dwt_max_resolution(opj_tcd_resolution_t* restrict r, OPJ_UINT32 i) { 545 OPJ_UINT32 mr = 0; 546 OPJ_UINT32 w; 547 while( --i ) { 548 ++r; 549 if( mr < ( w = (OPJ_UINT32)(r->x1 - r->x0) ) ) 550 mr = w ; 551 if( mr < ( w = (OPJ_UINT32)(r->y1 - r->y0) ) ) 552 mr = w ; 553 } 554 return mr ; 555} 556 557/* <summary> */ 558/* Inverse wavelet transform in 2-D. */ 559/* </summary> */ 560OPJ_BOOL opj_dwt_decode_tile(opj_tcd_tilecomp_t* tilec, OPJ_UINT32 numres, DWT1DFN dwt_1D) { 561 opj_dwt_t h; 562 opj_dwt_t v; 563 564 opj_tcd_resolution_t* tr = tilec->resolutions; 565 566 OPJ_UINT32 rw = (OPJ_UINT32)(tr->x1 - tr->x0); /* width of the resolution level computed */ 567 OPJ_UINT32 rh = (OPJ_UINT32)(tr->y1 - tr->y0); /* height of the resolution level computed */ 568 569 OPJ_UINT32 w = (OPJ_UINT32)(tilec->x1 - tilec->x0); 570 571 h.mem = (OPJ_INT32*) 572 opj_aligned_malloc(opj_dwt_max_resolution(tr, numres) * sizeof(OPJ_INT32)); 573 if (! h.mem){ 574 return OPJ_FALSE; 575 } 576 577 v.mem = h.mem; 578 579 while( --numres) { 580 OPJ_INT32 * restrict tiledp = tilec->data; 581 OPJ_UINT32 j; 582 583 ++tr; 584 h.sn = (OPJ_INT32)rw; 585 v.sn = (OPJ_INT32)rh; 586 587 rw = (OPJ_UINT32)(tr->x1 - tr->x0); 588 rh = (OPJ_UINT32)(tr->y1 - tr->y0); 589 590 h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn); 591 h.cas = tr->x0 % 2; 592 593 for(j = 0; j < rh; ++j) { 594 opj_dwt_interleave_h(&h, &tiledp[j*w]); 595 (dwt_1D)(&h); 596 memcpy(&tiledp[j*w], h.mem, rw * sizeof(OPJ_INT32)); 597 } 598 599 v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn); 600 v.cas = tr->y0 % 2; 601 602 for(j = 0; j < rw; ++j){ 603 OPJ_UINT32 k; 604 opj_dwt_interleave_v(&v, &tiledp[j], (OPJ_INT32)w); 605 (dwt_1D)(&v); 606 for(k = 0; k < rh; ++k) { 607 tiledp[k * w + j] = v.mem[k]; 608 } 609 } 610 } 611 opj_aligned_free(h.mem); 612 return OPJ_TRUE; 613} 614 615void opj_v4dwt_interleave_h(opj_v4dwt_t* restrict w, OPJ_FLOAT32* restrict a, OPJ_INT32 x, OPJ_INT32 size){ 616 OPJ_FLOAT32* restrict bi = (OPJ_FLOAT32*) (w->wavelet + w->cas); 617 OPJ_INT32 count = w->sn; 618 OPJ_INT32 i, k; 619 620 for(k = 0; k < 2; ++k){ 621 if ( count + 3 * x < size && ((size_t) a & 0x0f) == 0 && ((size_t) bi & 0x0f) == 0 && (x & 0x0f) == 0 ) { 622 /* Fast code path */ 623 for(i = 0; i < count; ++i){ 624 OPJ_INT32 j = i; 625 bi[i*8 ] = a[j]; 626 j += x; 627 bi[i*8 + 1] = a[j]; 628 j += x; 629 bi[i*8 + 2] = a[j]; 630 j += x; 631 bi[i*8 + 3] = a[j]; 632 } 633 } 634 else { 635 /* Slow code path */ 636 for(i = 0; i < count; ++i){ 637 OPJ_INT32 j = i; 638 bi[i*8 ] = a[j]; 639 j += x; 640 if(j >= size) continue; 641 bi[i*8 + 1] = a[j]; 642 j += x; 643 if(j >= size) continue; 644 bi[i*8 + 2] = a[j]; 645 j += x; 646 if(j >= size) continue; 647 bi[i*8 + 3] = a[j]; /* This one*/ 648 } 649 } 650 651 bi = (OPJ_FLOAT32*) (w->wavelet + 1 - w->cas); 652 a += w->sn; 653 size -= w->sn; 654 count = w->dn; 655 } 656} 657 658void opj_v4dwt_interleave_v(opj_v4dwt_t* restrict v , OPJ_FLOAT32* restrict a , OPJ_INT32 x, OPJ_INT32 nb_elts_read){ 659 opj_v4_t* restrict bi = v->wavelet + v->cas; 660 OPJ_INT32 i; 661 662 for(i = 0; i < v->sn; ++i){ 663 memcpy(&bi[i*2], &a[i*x], (size_t)nb_elts_read * sizeof(OPJ_FLOAT32)); 664 } 665 666 a += v->sn * x; 667 bi = v->wavelet + 1 - v->cas; 668 669 for(i = 0; i < v->dn; ++i){ 670 memcpy(&bi[i*2], &a[i*x], (size_t)nb_elts_read * sizeof(OPJ_FLOAT32)); 671 } 672} 673 674#ifdef __SSE__ 675 676void opj_v4dwt_decode_step1_sse(opj_v4_t* w, OPJ_INT32 count, const __m128 c){ 677 __m128* restrict vw = (__m128*) w; 678 OPJ_INT32 i; 679 /* 4x unrolled loop */ 680 for(i = 0; i < count >> 2; ++i){ 681 *vw = _mm_mul_ps(*vw, c); 682 vw += 2; 683 *vw = _mm_mul_ps(*vw, c); 684 vw += 2; 685 *vw = _mm_mul_ps(*vw, c); 686 vw += 2; 687 *vw = _mm_mul_ps(*vw, c); 688 vw += 2; 689 } 690 count &= 3; 691 for(i = 0; i < count; ++i){ 692 *vw = _mm_mul_ps(*vw, c); 693 vw += 2; 694 } 695} 696 697void opj_v4dwt_decode_step2_sse(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k, OPJ_INT32 m, __m128 c){ 698 __m128* restrict vl = (__m128*) l; 699 __m128* restrict vw = (__m128*) w; 700 OPJ_INT32 i; 701 __m128 tmp1, tmp2, tmp3; 702 tmp1 = vl[0]; 703 for(i = 0; i < m; ++i){ 704 tmp2 = vw[-1]; 705 tmp3 = vw[ 0]; 706 vw[-1] = _mm_add_ps(tmp2, _mm_mul_ps(_mm_add_ps(tmp1, tmp3), c)); 707 tmp1 = tmp3; 708 vw += 2; 709 } 710 vl = vw - 2; 711 if(m >= k){ 712 return; 713 } 714 c = _mm_add_ps(c, c); 715 c = _mm_mul_ps(c, vl[0]); 716 for(; m < k; ++m){ 717 __m128 tmp = vw[-1]; 718 vw[-1] = _mm_add_ps(tmp, c); 719 vw += 2; 720 } 721} 722 723#else 724 725void opj_v4dwt_decode_step1(opj_v4_t* w, OPJ_INT32 count, const OPJ_FLOAT32 c) 726{ 727 OPJ_FLOAT32* restrict fw = (OPJ_FLOAT32*) w; 728 OPJ_INT32 i; 729 for(i = 0; i < count; ++i){ 730 OPJ_FLOAT32 tmp1 = fw[i*8 ]; 731 OPJ_FLOAT32 tmp2 = fw[i*8 + 1]; 732 OPJ_FLOAT32 tmp3 = fw[i*8 + 2]; 733 OPJ_FLOAT32 tmp4 = fw[i*8 + 3]; 734 fw[i*8 ] = tmp1 * c; 735 fw[i*8 + 1] = tmp2 * c; 736 fw[i*8 + 2] = tmp3 * c; 737 fw[i*8 + 3] = tmp4 * c; 738 } 739} 740 741void opj_v4dwt_decode_step2(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k, OPJ_INT32 m, OPJ_FLOAT32 c) 742{ 743 OPJ_FLOAT32* restrict fl = (OPJ_FLOAT32*) l; 744 OPJ_FLOAT32* restrict fw = (OPJ_FLOAT32*) w; 745 OPJ_INT32 i; 746 for(i = 0; i < m; ++i){ 747 OPJ_FLOAT32 tmp1_1 = fl[0]; 748 OPJ_FLOAT32 tmp1_2 = fl[1]; 749 OPJ_FLOAT32 tmp1_3 = fl[2]; 750 OPJ_FLOAT32 tmp1_4 = fl[3]; 751 OPJ_FLOAT32 tmp2_1 = fw[-4]; 752 OPJ_FLOAT32 tmp2_2 = fw[-3]; 753 OPJ_FLOAT32 tmp2_3 = fw[-2]; 754 OPJ_FLOAT32 tmp2_4 = fw[-1]; 755 OPJ_FLOAT32 tmp3_1 = fw[0]; 756 OPJ_FLOAT32 tmp3_2 = fw[1]; 757 OPJ_FLOAT32 tmp3_3 = fw[2]; 758 OPJ_FLOAT32 tmp3_4 = fw[3]; 759 fw[-4] = tmp2_1 + ((tmp1_1 + tmp3_1) * c); 760 fw[-3] = tmp2_2 + ((tmp1_2 + tmp3_2) * c); 761 fw[-2] = tmp2_3 + ((tmp1_3 + tmp3_3) * c); 762 fw[-1] = tmp2_4 + ((tmp1_4 + tmp3_4) * c); 763 fl = fw; 764 fw += 8; 765 } 766 if(m < k){ 767 OPJ_FLOAT32 c1; 768 OPJ_FLOAT32 c2; 769 OPJ_FLOAT32 c3; 770 OPJ_FLOAT32 c4; 771 c += c; 772 c1 = fl[0] * c; 773 c2 = fl[1] * c; 774 c3 = fl[2] * c; 775 c4 = fl[3] * c; 776 for(; m < k; ++m){ 777 OPJ_FLOAT32 tmp1 = fw[-4]; 778 OPJ_FLOAT32 tmp2 = fw[-3]; 779 OPJ_FLOAT32 tmp3 = fw[-2]; 780 OPJ_FLOAT32 tmp4 = fw[-1]; 781 fw[-4] = tmp1 + c1; 782 fw[-3] = tmp2 + c2; 783 fw[-2] = tmp3 + c3; 784 fw[-1] = tmp4 + c4; 785 fw += 8; 786 } 787 } 788} 789 790#endif 791 792/* <summary> */ 793/* Inverse 9-7 wavelet transform in 1-D. */ 794/* </summary> */ 795void opj_v4dwt_decode(opj_v4dwt_t* restrict dwt) 796{ 797 OPJ_INT32 a, b; 798 if(dwt->cas == 0) { 799 if(!((dwt->dn > 0) || (dwt->sn > 1))){ 800 return; 801 } 802 a = 0; 803 b = 1; 804 }else{ 805 if(!((dwt->sn > 0) || (dwt->dn > 1))) { 806 return; 807 } 808 a = 1; 809 b = 0; 810 } 811#ifdef __SSE__ 812 opj_v4dwt_decode_step1_sse(dwt->wavelet+a, dwt->sn, _mm_set1_ps(opj_K)); 813 opj_v4dwt_decode_step1_sse(dwt->wavelet+b, dwt->dn, _mm_set1_ps(opj_c13318)); 814 opj_v4dwt_decode_step2_sse(dwt->wavelet+b, dwt->wavelet+a+1, dwt->sn, opj_int_min(dwt->sn, dwt->dn-a), _mm_set1_ps(opj_dwt_delta)); 815 opj_v4dwt_decode_step2_sse(dwt->wavelet+a, dwt->wavelet+b+1, dwt->dn, opj_int_min(dwt->dn, dwt->sn-b), _mm_set1_ps(opj_dwt_gamma)); 816 opj_v4dwt_decode_step2_sse(dwt->wavelet+b, dwt->wavelet+a+1, dwt->sn, opj_int_min(dwt->sn, dwt->dn-a), _mm_set1_ps(opj_dwt_beta)); 817 opj_v4dwt_decode_step2_sse(dwt->wavelet+a, dwt->wavelet+b+1, dwt->dn, opj_int_min(dwt->dn, dwt->sn-b), _mm_set1_ps(opj_dwt_alpha)); 818#else 819 opj_v4dwt_decode_step1(dwt->wavelet+a, dwt->sn, opj_K); 820 opj_v4dwt_decode_step1(dwt->wavelet+b, dwt->dn, opj_c13318); 821 opj_v4dwt_decode_step2(dwt->wavelet+b, dwt->wavelet+a+1, dwt->sn, opj_int_min(dwt->sn, dwt->dn-a), opj_dwt_delta); 822 opj_v4dwt_decode_step2(dwt->wavelet+a, dwt->wavelet+b+1, dwt->dn, opj_int_min(dwt->dn, dwt->sn-b), opj_dwt_gamma); 823 opj_v4dwt_decode_step2(dwt->wavelet+b, dwt->wavelet+a+1, dwt->sn, opj_int_min(dwt->sn, dwt->dn-a), opj_dwt_beta); 824 opj_v4dwt_decode_step2(dwt->wavelet+a, dwt->wavelet+b+1, dwt->dn, opj_int_min(dwt->dn, dwt->sn-b), opj_dwt_alpha); 825#endif 826} 827 828 829/* <summary> */ 830/* Inverse 9-7 wavelet transform in 2-D. */ 831/* </summary> */ 832OPJ_BOOL opj_dwt_decode_real(opj_tcd_tilecomp_t* restrict tilec, OPJ_UINT32 numres) 833{ 834 opj_v4dwt_t h; 835 opj_v4dwt_t v; 836 837 opj_tcd_resolution_t* res = tilec->resolutions; 838 839 OPJ_UINT32 rw = (OPJ_UINT32)(res->x1 - res->x0); /* width of the resolution level computed */ 840 OPJ_UINT32 rh = (OPJ_UINT32)(res->y1 - res->y0); /* height of the resolution level computed */ 841 842 OPJ_UINT32 w = (OPJ_UINT32)(tilec->x1 - tilec->x0); 843 844 h.wavelet = (opj_v4_t*) opj_aligned_malloc((opj_dwt_max_resolution(res, numres)+5) * sizeof(opj_v4_t)); 845 v.wavelet = h.wavelet; 846 847 while( --numres) { 848 OPJ_FLOAT32 * restrict aj = (OPJ_FLOAT32*) tilec->data; 849 OPJ_UINT32 bufsize = (OPJ_UINT32)((tilec->x1 - tilec->x0) * (tilec->y1 - tilec->y0)); 850 OPJ_INT32 j; 851 852 h.sn = (OPJ_INT32)rw; 853 v.sn = (OPJ_INT32)rh; 854 855 ++res; 856 857 rw = (OPJ_UINT32)(res->x1 - res->x0); /* width of the resolution level computed */ 858 rh = (OPJ_UINT32)(res->y1 - res->y0); /* height of the resolution level computed */ 859 860 h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn); 861 h.cas = res->x0 % 2; 862 863 for(j = (OPJ_INT32)rh; j > 3; j -= 4) { 864 OPJ_INT32 k; 865 opj_v4dwt_interleave_h(&h, aj, (OPJ_INT32)w, (OPJ_INT32)bufsize); 866 opj_v4dwt_decode(&h); 867 868 for(k = (OPJ_INT32)rw; --k >= 0;){ 869 aj[k ] = h.wavelet[k].f[0]; 870 aj[k+(OPJ_INT32)w ] = h.wavelet[k].f[1]; 871 aj[k+(OPJ_INT32)w*2] = h.wavelet[k].f[2]; 872 aj[k+(OPJ_INT32)w*3] = h.wavelet[k].f[3]; 873 } 874 875 aj += w*4; 876 bufsize -= w*4; 877 } 878 879 if (rh & 0x03) { 880 OPJ_INT32 k; 881 j = rh & 0x03; 882 opj_v4dwt_interleave_h(&h, aj, (OPJ_INT32)w, (OPJ_INT32)bufsize); 883 opj_v4dwt_decode(&h); 884 for(k = (OPJ_INT32)rw; --k >= 0;){ 885 switch(j) { 886 case 3: aj[k+(OPJ_INT32)w*2] = h.wavelet[k].f[2]; 887 case 2: aj[k+(OPJ_INT32)w ] = h.wavelet[k].f[1]; 888 case 1: aj[k ] = h.wavelet[k].f[0]; 889 } 890 } 891 } 892 893 v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn); 894 v.cas = res->y0 % 2; 895 896 aj = (OPJ_FLOAT32*) tilec->data; 897 for(j = (OPJ_INT32)rw; j > 3; j -= 4){ 898 OPJ_UINT32 k; 899 900 opj_v4dwt_interleave_v(&v, aj, (OPJ_INT32)w, 4); 901 opj_v4dwt_decode(&v); 902 903 for(k = 0; k < rh; ++k){ 904 memcpy(&aj[k*w], &v.wavelet[k], 4 * sizeof(OPJ_FLOAT32)); 905 } 906 aj += 4; 907 } 908 909 if (rw & 0x03){ 910 OPJ_UINT32 k; 911 912 j = rw & 0x03; 913 914 opj_v4dwt_interleave_v(&v, aj, (OPJ_INT32)w, j); 915 opj_v4dwt_decode(&v); 916 917 for(k = 0; k < rh; ++k){ 918 memcpy(&aj[k*w], &v.wavelet[k], (size_t)j * sizeof(OPJ_FLOAT32)); 919 } 920 } 921 } 922 923 opj_aligned_free(h.wavelet); 924 return OPJ_TRUE; 925} 926