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 44 #include "opj_includes.h" 45 46 /* <summary> */ 47 /* This table contains the norms of the basis function of the reversible MCT. */ 48 /* </summary> */ 49 static const OPJ_FLOAT64 opj_mct_norms[3] = { 1.732, .8292, .8292 }; 50 51 /* <summary> */ 52 /* This table contains the norms of the basis function of the irreversible MCT. */ 53 /* </summary> */ 54 static const OPJ_FLOAT64 opj_mct_norms_real[3] = { 1.732, 1.805, 1.573 }; 55 56 const OPJ_FLOAT64 * opj_mct_get_mct_norms () 57 { 58 return opj_mct_norms; 59 } 60 61 const OPJ_FLOAT64 * opj_mct_get_mct_norms_real () 62 { 63 return opj_mct_norms_real; 64 } 65 66 /* <summary> */ 67 /* Foward reversible MCT. */ 68 /* </summary> */ 69 void opj_mct_encode( 70 OPJ_INT32* restrict c0, 71 OPJ_INT32* restrict c1, 72 OPJ_INT32* restrict c2, 73 OPJ_UINT32 n) 74 { 75 OPJ_UINT32 i; 76 for(i = 0; i < n; ++i) { 77 OPJ_INT32 r = c0[i]; 78 OPJ_INT32 g = c1[i]; 79 OPJ_INT32 b = c2[i]; 80 OPJ_INT32 y = (r + (g * 2) + b) >> 2; 81 OPJ_INT32 u = b - g; 82 OPJ_INT32 v = r - g; 83 c0[i] = y; 84 c1[i] = u; 85 c2[i] = v; 86 } 87 } 88 89 /* <summary> */ 90 /* Inverse reversible MCT. */ 91 /* </summary> */ 92 void opj_mct_decode( 93 OPJ_INT32* restrict c0, 94 OPJ_INT32* restrict c1, 95 OPJ_INT32* restrict c2, 96 OPJ_UINT32 n) 97 { 98 OPJ_UINT32 i; 99 for (i = 0; i < n; ++i) { 100 OPJ_INT32 y = c0[i]; 101 OPJ_INT32 u = c1[i]; 102 OPJ_INT32 v = c2[i]; 103 OPJ_INT32 g = y - ((u + v) >> 2); 104 OPJ_INT32 r = v + g; 105 OPJ_INT32 b = u + g; 106 c0[i] = r; 107 c1[i] = g; 108 c2[i] = b; 109 } 110 } 111 112 /* <summary> */ 113 /* Get norm of basis function of reversible MCT. */ 114 /* </summary> */ 115 OPJ_FLOAT64 opj_mct_getnorm(OPJ_UINT32 compno) { 116 return opj_mct_norms[compno]; 117 } 118 119 /* <summary> */ 120 /* Foward irreversible MCT. */ 121 /* </summary> */ 122 void opj_mct_encode_real( 123 OPJ_INT32* restrict c0, 124 OPJ_INT32* restrict c1, 125 OPJ_INT32* restrict c2, 126 OPJ_UINT32 n) 127 { 128 OPJ_UINT32 i; 129 for(i = 0; i < n; ++i) { 130 OPJ_INT32 r = c0[i]; 131 OPJ_INT32 g = c1[i]; 132 OPJ_INT32 b = c2[i]; 133 OPJ_INT32 y = opj_int_fix_mul(r, 2449) + opj_int_fix_mul(g, 4809) + opj_int_fix_mul(b, 934); 134 OPJ_INT32 u = -opj_int_fix_mul(r, 1382) - opj_int_fix_mul(g, 2714) + opj_int_fix_mul(b, 4096); 135 OPJ_INT32 v = opj_int_fix_mul(r, 4096) - opj_int_fix_mul(g, 3430) - opj_int_fix_mul(b, 666); 136 c0[i] = y; 137 c1[i] = u; 138 c2[i] = v; 139 } 140 } 141 142 /* <summary> */ 143 /* Inverse irreversible MCT. */ 144 /* </summary> */ 145 void opj_mct_decode_real( 146 OPJ_FLOAT32* restrict c0, 147 OPJ_FLOAT32* restrict c1, 148 OPJ_FLOAT32* restrict c2, 149 OPJ_UINT32 n) 150 { 151 OPJ_UINT32 i; 152 #ifdef __SSE__ 153 // Mantis BUGID: 0056291. The address must be 16-byte aligned. 154 // TestFile: fuzz-signal_sigsegv_6e9e7f_5076_5265.pdf 155 if ((uintptr_t)c0 % 16 == 0 && (uintptr_t)c1 % 16 == 0 && (uintptr_t)c2 % 16 == 0){ 156 __m128 vrv, vgu, vgv, vbu; 157 vrv = _mm_set1_ps(1.402f); 158 vgu = _mm_set1_ps(0.34413f); 159 vgv = _mm_set1_ps(0.71414f); 160 vbu = _mm_set1_ps(1.772f); 161 for (i = 0; i < (n >> 3); ++i) { 162 __m128 vy, vu, vv; 163 __m128 vr, vg, vb; 164 165 vy = _mm_load_ps(c0); 166 vu = _mm_load_ps(c1); 167 vv = _mm_load_ps(c2); 168 vr = _mm_add_ps(vy, _mm_mul_ps(vv, vrv)); 169 vg = _mm_sub_ps(_mm_sub_ps(vy, _mm_mul_ps(vu, vgu)), _mm_mul_ps(vv, vgv)); 170 vb = _mm_add_ps(vy, _mm_mul_ps(vu, vbu)); 171 _mm_store_ps(c0, vr); 172 _mm_store_ps(c1, vg); 173 _mm_store_ps(c2, vb); 174 c0 += 4; 175 c1 += 4; 176 c2 += 4; 177 178 vy = _mm_load_ps(c0); 179 vu = _mm_load_ps(c1); 180 vv = _mm_load_ps(c2); 181 vr = _mm_add_ps(vy, _mm_mul_ps(vv, vrv)); 182 vg = _mm_sub_ps(_mm_sub_ps(vy, _mm_mul_ps(vu, vgu)), _mm_mul_ps(vv, vgv)); 183 vb = _mm_add_ps(vy, _mm_mul_ps(vu, vbu)); 184 _mm_store_ps(c0, vr); 185 _mm_store_ps(c1, vg); 186 _mm_store_ps(c2, vb); 187 c0 += 4; 188 c1 += 4; 189 c2 += 4; 190 } 191 n &= 7; 192 } else { 193 for(i = 0; i < n; ++i) { 194 OPJ_FLOAT32 y = c0[i]; 195 OPJ_FLOAT32 u = c1[i]; 196 OPJ_FLOAT32 v = c2[i]; 197 OPJ_FLOAT32 r = y + (v * 1.402f); 198 OPJ_FLOAT32 g = y - (u * 0.34413f) - (v * (0.71414f)); 199 OPJ_FLOAT32 b = y + (u * 1.772f); 200 c0[i] = r; 201 c1[i] = g; 202 c2[i] = b; 203 } 204 } 205 206 #endif 207 for(i = 0; i < n; ++i) { 208 OPJ_FLOAT32 y = c0[i]; 209 OPJ_FLOAT32 u = c1[i]; 210 OPJ_FLOAT32 v = c2[i]; 211 OPJ_FLOAT32 r = y + (v * 1.402f); 212 OPJ_FLOAT32 g = y - (u * 0.34413f) - (v * (0.71414f)); 213 OPJ_FLOAT32 b = y + (u * 1.772f); 214 c0[i] = r; 215 c1[i] = g; 216 c2[i] = b; 217 } 218 } 219 220 /* <summary> */ 221 /* Get norm of basis function of irreversible MCT. */ 222 /* </summary> */ 223 OPJ_FLOAT64 opj_mct_getnorm_real(OPJ_UINT32 compno) { 224 return opj_mct_norms_real[compno]; 225 } 226 227 228 OPJ_BOOL opj_mct_encode_custom( 229 OPJ_BYTE * pCodingdata, 230 OPJ_UINT32 n, 231 OPJ_BYTE ** pData, 232 OPJ_UINT32 pNbComp, 233 OPJ_UINT32 isSigned) 234 { 235 OPJ_FLOAT32 * lMct = (OPJ_FLOAT32 *) pCodingdata; 236 OPJ_UINT32 i; 237 OPJ_UINT32 j; 238 OPJ_UINT32 k; 239 OPJ_UINT32 lNbMatCoeff = pNbComp * pNbComp; 240 OPJ_INT32 * lCurrentData = 00; 241 OPJ_INT32 * lCurrentMatrix = 00; 242 OPJ_INT32 ** lData = (OPJ_INT32 **) pData; 243 OPJ_UINT32 lMultiplicator = 1 << 13; 244 OPJ_INT32 * lMctPtr; 245 246 OPJ_ARG_NOT_USED(isSigned); 247 248 lCurrentData = (OPJ_INT32 *) opj_malloc((pNbComp + lNbMatCoeff) * sizeof(OPJ_INT32)); 249 if (! lCurrentData) { 250 return OPJ_FALSE; 251 } 252 253 lCurrentMatrix = lCurrentData + pNbComp; 254 255 for (i =0;i<lNbMatCoeff;++i) { 256 lCurrentMatrix[i] = (OPJ_INT32) (*(lMct++) * (OPJ_FLOAT32)lMultiplicator); 257 } 258 259 for (i = 0; i < n; ++i) { 260 lMctPtr = lCurrentMatrix; 261 for (j=0;j<pNbComp;++j) { 262 lCurrentData[j] = (*(lData[j])); 263 } 264 265 for (j=0;j<pNbComp;++j) { 266 *(lData[j]) = 0; 267 for (k=0;k<pNbComp;++k) { 268 *(lData[j]) += opj_int_fix_mul(*lMctPtr, lCurrentData[k]); 269 ++lMctPtr; 270 } 271 272 ++lData[j]; 273 } 274 } 275 276 opj_free(lCurrentData); 277 278 return OPJ_TRUE; 279 } 280 281 OPJ_BOOL opj_mct_decode_custom( 282 OPJ_BYTE * pDecodingData, 283 OPJ_UINT32 n, 284 OPJ_BYTE ** pData, 285 OPJ_UINT32 pNbComp, 286 OPJ_UINT32 isSigned) 287 { 288 OPJ_FLOAT32 * lMct; 289 OPJ_UINT32 i; 290 OPJ_UINT32 j; 291 OPJ_UINT32 k; 292 293 OPJ_FLOAT32 * lCurrentData = 00; 294 OPJ_FLOAT32 * lCurrentResult = 00; 295 OPJ_FLOAT32 ** lData = (OPJ_FLOAT32 **) pData; 296 297 OPJ_ARG_NOT_USED(isSigned); 298 299 lCurrentData = (OPJ_FLOAT32 *) opj_malloc (2 * pNbComp * sizeof(OPJ_FLOAT32)); 300 if (! lCurrentData) { 301 return OPJ_FALSE; 302 } 303 lCurrentResult = lCurrentData + pNbComp; 304 305 for (i = 0; i < n; ++i) { 306 lMct = (OPJ_FLOAT32 *) pDecodingData; 307 for (j=0;j<pNbComp;++j) { 308 lCurrentData[j] = (OPJ_FLOAT32) (*(lData[j])); 309 } 310 for (j=0;j<pNbComp;++j) { 311 lCurrentResult[j] = 0; 312 for (k=0;k<pNbComp;++k) { 313 lCurrentResult[j] += *(lMct++) * lCurrentData[k]; 314 } 315 *(lData[j]++) = (OPJ_FLOAT32) (lCurrentResult[j]); 316 } 317 } 318 opj_free(lCurrentData); 319 return OPJ_TRUE; 320 } 321 322 void opj_calculate_norms( OPJ_FLOAT64 * pNorms, 323 OPJ_UINT32 pNbComps, 324 OPJ_FLOAT32 * pMatrix) 325 { 326 OPJ_UINT32 i,j,lIndex; 327 OPJ_FLOAT32 lCurrentValue; 328 OPJ_FLOAT64 * lNorms = (OPJ_FLOAT64 *) pNorms; 329 OPJ_FLOAT32 * lMatrix = (OPJ_FLOAT32 *) pMatrix; 330 331 for (i=0;i<pNbComps;++i) { 332 lNorms[i] = 0; 333 lIndex = i; 334 335 for (j=0;j<pNbComps;++j) { 336 lCurrentValue = lMatrix[lIndex]; 337 lIndex += pNbComps; 338 lNorms[i] += lCurrentValue * lCurrentValue; 339 } 340 lNorms[i] = sqrt(lNorms[i]); 341 } 342 } 343
