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) 2008, 2011-2012, Centre National d'Etudes Spatiales (CNES), FR 15 * Copyright (c) 2012, CS Systemes d'Information, France 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 #ifdef __SSE2__ 44 #include <emmintrin.h> 45 #endif 46 #ifdef __SSE4_1__ 47 #include <smmintrin.h> 48 #endif 49 50 #include "opj_includes.h" 51 52 /* <summary> */ 53 /* This table contains the norms of the basis function of the reversible MCT. */ 54 /* </summary> */ 55 static const OPJ_FLOAT64 opj_mct_norms[3] = { 1.732, .8292, .8292 }; 56 57 /* <summary> */ 58 /* This table contains the norms of the basis function of the irreversible MCT. */ 59 /* </summary> */ 60 static const OPJ_FLOAT64 opj_mct_norms_real[3] = { 1.732, 1.805, 1.573 }; 61 62 const OPJ_FLOAT64 * opj_mct_get_mct_norms () 63 { 64 return opj_mct_norms; 65 } 66 67 const OPJ_FLOAT64 * opj_mct_get_mct_norms_real () 68 { 69 return opj_mct_norms_real; 70 } 71 72 /* <summary> */ 73 /* Forward reversible MCT. */ 74 /* </summary> */ 75 #ifdef __SSE2__ 76 void opj_mct_encode( 77 OPJ_INT32* restrict c0, 78 OPJ_INT32* restrict c1, 79 OPJ_INT32* restrict c2, 80 OPJ_UINT32 n) 81 { 82 OPJ_SIZE_T i; 83 const OPJ_SIZE_T len = n; 84 85 for(i = 0; i < (len & ~3U); i += 4) { 86 __m128i y, u, v; 87 __m128i r = _mm_load_si128((const __m128i *)&(c0[i])); 88 __m128i g = _mm_load_si128((const __m128i *)&(c1[i])); 89 __m128i b = _mm_load_si128((const __m128i *)&(c2[i])); 90 y = _mm_add_epi32(g, g); 91 y = _mm_add_epi32(y, b); 92 y = _mm_add_epi32(y, r); 93 y = _mm_srai_epi32(y, 2); 94 u = _mm_sub_epi32(b, g); 95 v = _mm_sub_epi32(r, g); 96 _mm_store_si128((__m128i *)&(c0[i]), y); 97 _mm_store_si128((__m128i *)&(c1[i]), u); 98 _mm_store_si128((__m128i *)&(c2[i]), v); 99 } 100 101 for(; i < len; ++i) { 102 OPJ_INT32 r = c0[i]; 103 OPJ_INT32 g = c1[i]; 104 OPJ_INT32 b = c2[i]; 105 OPJ_INT32 y = (r + (g * 2) + b) >> 2; 106 OPJ_INT32 u = b - g; 107 OPJ_INT32 v = r - g; 108 c0[i] = y; 109 c1[i] = u; 110 c2[i] = v; 111 } 112 } 113 #else 114 void opj_mct_encode( 115 OPJ_INT32* restrict c0, 116 OPJ_INT32* restrict c1, 117 OPJ_INT32* restrict c2, 118 OPJ_UINT32 n) 119 { 120 OPJ_SIZE_T i; 121 const OPJ_SIZE_T len = n; 122 123 for(i = 0; i < len; ++i) { 124 OPJ_INT32 r = c0[i]; 125 OPJ_INT32 g = c1[i]; 126 OPJ_INT32 b = c2[i]; 127 OPJ_INT32 y = (r + (g * 2) + b) >> 2; 128 OPJ_INT32 u = b - g; 129 OPJ_INT32 v = r - g; 130 c0[i] = y; 131 c1[i] = u; 132 c2[i] = v; 133 } 134 } 135 #endif 136 137 /* <summary> */ 138 /* Inverse reversible MCT. */ 139 /* </summary> */ 140 #ifdef __SSE2__ 141 void opj_mct_decode( 142 OPJ_INT32* restrict c0, 143 OPJ_INT32* restrict c1, 144 OPJ_INT32* restrict c2, 145 OPJ_UINT32 n) 146 { 147 OPJ_SIZE_T i; 148 const OPJ_SIZE_T len = n; 149 150 for(i = 0; i < (len & ~3U); i += 4) { 151 __m128i r, g, b; 152 __m128i y = _mm_load_si128((const __m128i *)&(c0[i])); 153 __m128i u = _mm_load_si128((const __m128i *)&(c1[i])); 154 __m128i v = _mm_load_si128((const __m128i *)&(c2[i])); 155 g = y; 156 g = _mm_sub_epi32(g, _mm_srai_epi32(_mm_add_epi32(u, v), 2)); 157 r = _mm_add_epi32(v, g); 158 b = _mm_add_epi32(u, g); 159 _mm_store_si128((__m128i *)&(c0[i]), r); 160 _mm_store_si128((__m128i *)&(c1[i]), g); 161 _mm_store_si128((__m128i *)&(c2[i]), b); 162 } 163 for (; i < len; ++i) { 164 OPJ_INT32 y = c0[i]; 165 OPJ_INT32 u = c1[i]; 166 OPJ_INT32 v = c2[i]; 167 OPJ_INT32 g = y - ((u + v) >> 2); 168 OPJ_INT32 r = v + g; 169 OPJ_INT32 b = u + g; 170 c0[i] = r; 171 c1[i] = g; 172 c2[i] = b; 173 } 174 } 175 #else 176 void opj_mct_decode( 177 OPJ_INT32* restrict c0, 178 OPJ_INT32* restrict c1, 179 OPJ_INT32* restrict c2, 180 OPJ_UINT32 n) 181 { 182 OPJ_UINT32 i; 183 for (i = 0; i < n; ++i) { 184 OPJ_INT32 y = c0[i]; 185 OPJ_INT32 u = c1[i]; 186 OPJ_INT32 v = c2[i]; 187 OPJ_INT32 g = y - ((u + v) >> 2); 188 OPJ_INT32 r = v + g; 189 OPJ_INT32 b = u + g; 190 c0[i] = r; 191 c1[i] = g; 192 c2[i] = b; 193 } 194 } 195 #endif 196 197 /* <summary> */ 198 /* Get norm of basis function of reversible MCT. */ 199 /* </summary> */ 200 OPJ_FLOAT64 opj_mct_getnorm(OPJ_UINT32 compno) { 201 return opj_mct_norms[compno]; 202 } 203 204 /* <summary> */ 205 /* Forward irreversible MCT. */ 206 /* </summary> */ 207 #ifdef __SSE4_1__ 208 void opj_mct_encode_real( 209 OPJ_INT32* restrict c0, 210 OPJ_INT32* restrict c1, 211 OPJ_INT32* restrict c2, 212 OPJ_UINT32 n) 213 { 214 OPJ_SIZE_T i; 215 const OPJ_SIZE_T len = n; 216 217 const __m128i ry = _mm_set1_epi32(2449); 218 const __m128i gy = _mm_set1_epi32(4809); 219 const __m128i by = _mm_set1_epi32(934); 220 const __m128i ru = _mm_set1_epi32(1382); 221 const __m128i gu = _mm_set1_epi32(2714); 222 /* const __m128i bu = _mm_set1_epi32(4096); */ 223 /* const __m128i rv = _mm_set1_epi32(4096); */ 224 const __m128i gv = _mm_set1_epi32(3430); 225 const __m128i bv = _mm_set1_epi32(666); 226 const __m128i mulround = _mm_shuffle_epi32(_mm_cvtsi32_si128(4096), _MM_SHUFFLE(1, 0, 1, 0)); 227 228 for(i = 0; i < (len & ~3U); i += 4) { 229 __m128i lo, hi; 230 __m128i y, u, v; 231 __m128i r = _mm_load_si128((const __m128i *)&(c0[i])); 232 __m128i g = _mm_load_si128((const __m128i *)&(c1[i])); 233 __m128i b = _mm_load_si128((const __m128i *)&(c2[i])); 234 235 lo = r; 236 hi = _mm_shuffle_epi32(r, _MM_SHUFFLE(3, 3, 1, 1)); 237 lo = _mm_mul_epi32(lo, ry); 238 hi = _mm_mul_epi32(hi, ry); 239 lo = _mm_add_epi64(lo, mulround); 240 hi = _mm_add_epi64(hi, mulround); 241 lo = _mm_srli_epi64(lo, 13); 242 hi = _mm_slli_epi64(hi, 32-13); 243 y = _mm_blend_epi16(lo, hi, 0xCC); 244 245 lo = g; 246 hi = _mm_shuffle_epi32(g, _MM_SHUFFLE(3, 3, 1, 1)); 247 lo = _mm_mul_epi32(lo, gy); 248 hi = _mm_mul_epi32(hi, gy); 249 lo = _mm_add_epi64(lo, mulround); 250 hi = _mm_add_epi64(hi, mulround); 251 lo = _mm_srli_epi64(lo, 13); 252 hi = _mm_slli_epi64(hi, 32-13); 253 y = _mm_add_epi32(y, _mm_blend_epi16(lo, hi, 0xCC)); 254 255 lo = b; 256 hi = _mm_shuffle_epi32(b, _MM_SHUFFLE(3, 3, 1, 1)); 257 lo = _mm_mul_epi32(lo, by); 258 hi = _mm_mul_epi32(hi, by); 259 lo = _mm_add_epi64(lo, mulround); 260 hi = _mm_add_epi64(hi, mulround); 261 lo = _mm_srli_epi64(lo, 13); 262 hi = _mm_slli_epi64(hi, 32-13); 263 y = _mm_add_epi32(y, _mm_blend_epi16(lo, hi, 0xCC)); 264 _mm_store_si128((__m128i *)&(c0[i]), y); 265 266 /*lo = b; 267 hi = _mm_shuffle_epi32(b, _MM_SHUFFLE(3, 3, 1, 1)); 268 lo = _mm_mul_epi32(lo, mulround); 269 hi = _mm_mul_epi32(hi, mulround);*/ 270 lo = _mm_cvtepi32_epi64(_mm_shuffle_epi32(b, _MM_SHUFFLE(3, 2, 2, 0))); 271 hi = _mm_cvtepi32_epi64(_mm_shuffle_epi32(b, _MM_SHUFFLE(3, 2, 3, 1))); 272 lo = _mm_slli_epi64(lo, 12); 273 hi = _mm_slli_epi64(hi, 12); 274 lo = _mm_add_epi64(lo, mulround); 275 hi = _mm_add_epi64(hi, mulround); 276 lo = _mm_srli_epi64(lo, 13); 277 hi = _mm_slli_epi64(hi, 32-13); 278 u = _mm_blend_epi16(lo, hi, 0xCC); 279 280 lo = r; 281 hi = _mm_shuffle_epi32(r, _MM_SHUFFLE(3, 3, 1, 1)); 282 lo = _mm_mul_epi32(lo, ru); 283 hi = _mm_mul_epi32(hi, ru); 284 lo = _mm_add_epi64(lo, mulround); 285 hi = _mm_add_epi64(hi, mulround); 286 lo = _mm_srli_epi64(lo, 13); 287 hi = _mm_slli_epi64(hi, 32-13); 288 u = _mm_sub_epi32(u, _mm_blend_epi16(lo, hi, 0xCC)); 289 290 lo = g; 291 hi = _mm_shuffle_epi32(g, _MM_SHUFFLE(3, 3, 1, 1)); 292 lo = _mm_mul_epi32(lo, gu); 293 hi = _mm_mul_epi32(hi, gu); 294 lo = _mm_add_epi64(lo, mulround); 295 hi = _mm_add_epi64(hi, mulround); 296 lo = _mm_srli_epi64(lo, 13); 297 hi = _mm_slli_epi64(hi, 32-13); 298 u = _mm_sub_epi32(u, _mm_blend_epi16(lo, hi, 0xCC)); 299 _mm_store_si128((__m128i *)&(c1[i]), u); 300 301 /*lo = r; 302 hi = _mm_shuffle_epi32(r, _MM_SHUFFLE(3, 3, 1, 1)); 303 lo = _mm_mul_epi32(lo, mulround); 304 hi = _mm_mul_epi32(hi, mulround);*/ 305 lo = _mm_cvtepi32_epi64(_mm_shuffle_epi32(r, _MM_SHUFFLE(3, 2, 2, 0))); 306 hi = _mm_cvtepi32_epi64(_mm_shuffle_epi32(r, _MM_SHUFFLE(3, 2, 3, 1))); 307 lo = _mm_slli_epi64(lo, 12); 308 hi = _mm_slli_epi64(hi, 12); 309 lo = _mm_add_epi64(lo, mulround); 310 hi = _mm_add_epi64(hi, mulround); 311 lo = _mm_srli_epi64(lo, 13); 312 hi = _mm_slli_epi64(hi, 32-13); 313 v = _mm_blend_epi16(lo, hi, 0xCC); 314 315 lo = g; 316 hi = _mm_shuffle_epi32(g, _MM_SHUFFLE(3, 3, 1, 1)); 317 lo = _mm_mul_epi32(lo, gv); 318 hi = _mm_mul_epi32(hi, gv); 319 lo = _mm_add_epi64(lo, mulround); 320 hi = _mm_add_epi64(hi, mulround); 321 lo = _mm_srli_epi64(lo, 13); 322 hi = _mm_slli_epi64(hi, 32-13); 323 v = _mm_sub_epi32(v, _mm_blend_epi16(lo, hi, 0xCC)); 324 325 lo = b; 326 hi = _mm_shuffle_epi32(b, _MM_SHUFFLE(3, 3, 1, 1)); 327 lo = _mm_mul_epi32(lo, bv); 328 hi = _mm_mul_epi32(hi, bv); 329 lo = _mm_add_epi64(lo, mulround); 330 hi = _mm_add_epi64(hi, mulround); 331 lo = _mm_srli_epi64(lo, 13); 332 hi = _mm_slli_epi64(hi, 32-13); 333 v = _mm_sub_epi32(v, _mm_blend_epi16(lo, hi, 0xCC)); 334 _mm_store_si128((__m128i *)&(c2[i]), v); 335 } 336 for(; i < len; ++i) { 337 OPJ_INT32 r = c0[i]; 338 OPJ_INT32 g = c1[i]; 339 OPJ_INT32 b = c2[i]; 340 OPJ_INT32 y = opj_int_fix_mul(r, 2449) + opj_int_fix_mul(g, 4809) + opj_int_fix_mul(b, 934); 341 OPJ_INT32 u = -opj_int_fix_mul(r, 1382) - opj_int_fix_mul(g, 2714) + opj_int_fix_mul(b, 4096); 342 OPJ_INT32 v = opj_int_fix_mul(r, 4096) - opj_int_fix_mul(g, 3430) - opj_int_fix_mul(b, 666); 343 c0[i] = y; 344 c1[i] = u; 345 c2[i] = v; 346 } 347 } 348 #else 349 void opj_mct_encode_real( 350 OPJ_INT32* restrict c0, 351 OPJ_INT32* restrict c1, 352 OPJ_INT32* restrict c2, 353 OPJ_UINT32 n) 354 { 355 OPJ_UINT32 i; 356 for(i = 0; i < n; ++i) { 357 OPJ_INT32 r = c0[i]; 358 OPJ_INT32 g = c1[i]; 359 OPJ_INT32 b = c2[i]; 360 OPJ_INT32 y = opj_int_fix_mul(r, 2449) + opj_int_fix_mul(g, 4809) + opj_int_fix_mul(b, 934); 361 OPJ_INT32 u = -opj_int_fix_mul(r, 1382) - opj_int_fix_mul(g, 2714) + opj_int_fix_mul(b, 4096); 362 OPJ_INT32 v = opj_int_fix_mul(r, 4096) - opj_int_fix_mul(g, 3430) - opj_int_fix_mul(b, 666); 363 c0[i] = y; 364 c1[i] = u; 365 c2[i] = v; 366 } 367 } 368 #endif 369 370 /* <summary> */ 371 /* Inverse irreversible MCT. */ 372 /* </summary> */ 373 void opj_mct_decode_real( 374 OPJ_FLOAT32* restrict c0, 375 OPJ_FLOAT32* restrict c1, 376 OPJ_FLOAT32* restrict c2, 377 OPJ_UINT32 n) 378 { 379 OPJ_UINT32 i; 380 #ifdef __SSE__ 381 __m128 vrv, vgu, vgv, vbu; 382 vrv = _mm_set1_ps(1.402f); 383 vgu = _mm_set1_ps(0.34413f); 384 vgv = _mm_set1_ps(0.71414f); 385 vbu = _mm_set1_ps(1.772f); 386 for (i = 0; i < (n >> 3); ++i) { 387 __m128 vy, vu, vv; 388 __m128 vr, vg, vb; 389 390 vy = _mm_load_ps(c0); 391 vu = _mm_load_ps(c1); 392 vv = _mm_load_ps(c2); 393 vr = _mm_add_ps(vy, _mm_mul_ps(vv, vrv)); 394 vg = _mm_sub_ps(_mm_sub_ps(vy, _mm_mul_ps(vu, vgu)), _mm_mul_ps(vv, vgv)); 395 vb = _mm_add_ps(vy, _mm_mul_ps(vu, vbu)); 396 _mm_store_ps(c0, vr); 397 _mm_store_ps(c1, vg); 398 _mm_store_ps(c2, vb); 399 c0 += 4; 400 c1 += 4; 401 c2 += 4; 402 403 vy = _mm_load_ps(c0); 404 vu = _mm_load_ps(c1); 405 vv = _mm_load_ps(c2); 406 vr = _mm_add_ps(vy, _mm_mul_ps(vv, vrv)); 407 vg = _mm_sub_ps(_mm_sub_ps(vy, _mm_mul_ps(vu, vgu)), _mm_mul_ps(vv, vgv)); 408 vb = _mm_add_ps(vy, _mm_mul_ps(vu, vbu)); 409 _mm_store_ps(c0, vr); 410 _mm_store_ps(c1, vg); 411 _mm_store_ps(c2, vb); 412 c0 += 4; 413 c1 += 4; 414 c2 += 4; 415 } 416 n &= 7; 417 #endif 418 for(i = 0; i < n; ++i) { 419 OPJ_FLOAT32 y = c0[i]; 420 OPJ_FLOAT32 u = c1[i]; 421 OPJ_FLOAT32 v = c2[i]; 422 OPJ_FLOAT32 r = y + (v * 1.402f); 423 OPJ_FLOAT32 g = y - (u * 0.34413f) - (v * (0.71414f)); 424 OPJ_FLOAT32 b = y + (u * 1.772f); 425 c0[i] = r; 426 c1[i] = g; 427 c2[i] = b; 428 } 429 } 430 431 /* <summary> */ 432 /* Get norm of basis function of irreversible MCT. */ 433 /* </summary> */ 434 OPJ_FLOAT64 opj_mct_getnorm_real(OPJ_UINT32 compno) { 435 return opj_mct_norms_real[compno]; 436 } 437 438 439 OPJ_BOOL opj_mct_encode_custom( 440 OPJ_BYTE * pCodingdata, 441 OPJ_UINT32 n, 442 OPJ_BYTE ** pData, 443 OPJ_UINT32 pNbComp, 444 OPJ_UINT32 isSigned) 445 { 446 OPJ_FLOAT32 * lMct = (OPJ_FLOAT32 *) pCodingdata; 447 OPJ_UINT32 i; 448 OPJ_UINT32 j; 449 OPJ_UINT32 k; 450 OPJ_UINT32 lNbMatCoeff = pNbComp * pNbComp; 451 OPJ_INT32 * lCurrentData = 00; 452 OPJ_INT32 * lCurrentMatrix = 00; 453 OPJ_INT32 ** lData = (OPJ_INT32 **) pData; 454 OPJ_UINT32 lMultiplicator = 1 << 13; 455 OPJ_INT32 * lMctPtr; 456 457 OPJ_ARG_NOT_USED(isSigned); 458 459 lCurrentData = (OPJ_INT32 *) opj_malloc((pNbComp + lNbMatCoeff) * sizeof(OPJ_INT32)); 460 if (! lCurrentData) { 461 return OPJ_FALSE; 462 } 463 464 lCurrentMatrix = lCurrentData + pNbComp; 465 466 for (i =0;i<lNbMatCoeff;++i) { 467 lCurrentMatrix[i] = (OPJ_INT32) (*(lMct++) * (OPJ_FLOAT32)lMultiplicator); 468 } 469 470 for (i = 0; i < n; ++i) { 471 lMctPtr = lCurrentMatrix; 472 for (j=0;j<pNbComp;++j) { 473 lCurrentData[j] = (*(lData[j])); 474 } 475 476 for (j=0;j<pNbComp;++j) { 477 *(lData[j]) = 0; 478 for (k=0;k<pNbComp;++k) { 479 *(lData[j]) += opj_int_fix_mul(*lMctPtr, lCurrentData[k]); 480 ++lMctPtr; 481 } 482 483 ++lData[j]; 484 } 485 } 486 487 opj_free(lCurrentData); 488 489 return OPJ_TRUE; 490 } 491 492 OPJ_BOOL opj_mct_decode_custom( 493 OPJ_BYTE * pDecodingData, 494 OPJ_UINT32 n, 495 OPJ_BYTE ** pData, 496 OPJ_UINT32 pNbComp, 497 OPJ_UINT32 isSigned) 498 { 499 OPJ_FLOAT32 * lMct; 500 OPJ_UINT32 i; 501 OPJ_UINT32 j; 502 OPJ_UINT32 k; 503 504 OPJ_FLOAT32 * lCurrentData = 00; 505 OPJ_FLOAT32 * lCurrentResult = 00; 506 OPJ_FLOAT32 ** lData = (OPJ_FLOAT32 **) pData; 507 508 OPJ_ARG_NOT_USED(isSigned); 509 510 lCurrentData = (OPJ_FLOAT32 *) opj_malloc (2 * pNbComp * sizeof(OPJ_FLOAT32)); 511 if (! lCurrentData) { 512 return OPJ_FALSE; 513 } 514 lCurrentResult = lCurrentData + pNbComp; 515 516 for (i = 0; i < n; ++i) { 517 lMct = (OPJ_FLOAT32 *) pDecodingData; 518 for (j=0;j<pNbComp;++j) { 519 lCurrentData[j] = (OPJ_FLOAT32) (*(lData[j])); 520 } 521 for (j=0;j<pNbComp;++j) { 522 lCurrentResult[j] = 0; 523 for (k=0;k<pNbComp;++k) { 524 lCurrentResult[j] += *(lMct++) * lCurrentData[k]; 525 } 526 *(lData[j]++) = (OPJ_FLOAT32) (lCurrentResult[j]); 527 } 528 } 529 opj_free(lCurrentData); 530 return OPJ_TRUE; 531 } 532 533 void opj_calculate_norms( OPJ_FLOAT64 * pNorms, 534 OPJ_UINT32 pNbComps, 535 OPJ_FLOAT32 * pMatrix) 536 { 537 OPJ_UINT32 i,j,lIndex; 538 OPJ_FLOAT32 lCurrentValue; 539 OPJ_FLOAT64 * lNorms = (OPJ_FLOAT64 *) pNorms; 540 OPJ_FLOAT32 * lMatrix = (OPJ_FLOAT32 *) pMatrix; 541 542 for (i=0;i<pNbComps;++i) { 543 lNorms[i] = 0; 544 lIndex = i; 545 546 for (j=0;j<pNbComps;++j) { 547 lCurrentValue = lMatrix[lIndex]; 548 lIndex += pNbComps; 549 lNorms[i] += lCurrentValue * lCurrentValue; 550 } 551 lNorms[i] = sqrt(lNorms[i]); 552 } 553 } 554
