QRScanner/mobile/jni/thirdparty/openjpeg/libopenjpeg/mct.c

   1 /*
   2  * Copyright (c) 2002-2007, Communications and Remote Sensing Laboratory, Universite catholique de Louvain (UCL), Belgium
   3  * Copyright (c) 2002-2007, Professor Benoit Macq
   4  * Copyright (c) 2001-2003, David Janssens
   5  * Copyright (c) 2002-2003, Yannick Verschueren
   6  * Copyright (c) 2003-2007, Francois-Olivier Devaux and Antonin Descampe
   7  * Copyright (c) 2005, Herve Drolon, FreeImage Team
   8  * Copyright (c) 2008;2011-2012, Centre National d'Etudes Spatiales (CNES), France
   9  * Copyright (c) 2012, CS Systemes d'Information, France
  10  * All rights reserved.
  11  *
  12  * Redistribution and use in source and binary forms, with or without
  13  * modification, are permitted provided that the following conditions
  14  * are met:
  15  * 1. Redistributions of source code must retain the above copyright
  16  *    notice, this list of conditions and the following disclaimer.
  17  * 2. Redistributions in binary form must reproduce the above copyright
  18  *    notice, this list of conditions and the following disclaimer in the
  19  *    documentation and/or other materials provided with the distribution.
  20  *
  21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
  22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  31  * POSSIBILITY OF SUCH DAMAGE.
  32  */
  33
  34 #ifdef __SSE__
  35 #include <xmmintrin.h>
  36 #endif
  37
  38 #include "opj_includes.h"
  39
  40 /* <summary> */
  41 /* This table contains the norms of the basis function of the reversible MCT. */
  42 /* </summary> */
  43 static const OPJ_FLOAT64 opj_mct_norms[3] = { 1.732, .8292, .8292 };
  44
  45 /* <summary> */
  46 /* This table contains the norms of the basis function of the irreversible MCT. */
  47 /* </summary> */
  48 static const OPJ_FLOAT64 opj_mct_norms_real[3] = { 1.732, 1.805, 1.573 };
  49
  50 const OPJ_FLOAT64 * opj_mct_get_mct_norms ()
  51 {
  52         return opj_mct_norms;
  53 }
  54
  55 const OPJ_FLOAT64 * opj_mct_get_mct_norms_real ()
  56 {
  57         return opj_mct_norms_real;
  58 }
  59
  60 /* <summary> */
  61 /* Foward reversible MCT. */
  62 /* </summary> */
  63 void opj_mct_encode(
  64                 OPJ_INT32* restrict c0,
  65                 OPJ_INT32* restrict c1,
  66                 OPJ_INT32* restrict c2,
  67                 OPJ_UINT32 n)
  68 {
  69         OPJ_UINT32 i;
  70         for(i = 0; i < n; ++i) {
  71                 OPJ_INT32 r = c0[i];
  72                 OPJ_INT32 g = c1[i];
  73                 OPJ_INT32 b = c2[i];
  74                 OPJ_INT32 y = (r + (g * 2) + b) >> 2;
  75                 OPJ_INT32 u = b - g;
  76                 OPJ_INT32 v = r - g;
  77                 c0[i] = y;
  78                 c1[i] = u;
  79                 c2[i] = v;
  80         }
  81 }
  82
  83 /* <summary> */
  84 /* Inverse reversible MCT. */
  85 /* </summary> */
  86 void opj_mct_decode(
  87                 OPJ_INT32* restrict c0,
  88                 OPJ_INT32* restrict c1,
  89                 OPJ_INT32* restrict c2,
  90                 OPJ_UINT32 n)
  91 {
  92         OPJ_UINT32 i;
  93         for (i = 0; i < n; ++i) {
  94                 OPJ_INT32 y = c0[i];
  95                 OPJ_INT32 u = c1[i];
  96                 OPJ_INT32 v = c2[i];
  97                 OPJ_INT32 g = y - ((u + v) >> 2);
  98                 OPJ_INT32 r = v + g;
  99                 OPJ_INT32 b = u + g;
 100                 c0[i] = r;
 101                 c1[i] = g;
 102                 c2[i] = b;
 103         }
 104 }
 105
 106 /* <summary> */
 107 /* Get norm of basis function of reversible MCT. */
 108 /* </summary> */
 109 OPJ_FLOAT64 opj_mct_getnorm(OPJ_UINT32 compno) {
 110         return opj_mct_norms[compno];
 111 }
 112
 113 /* <summary> */
 114 /* Foward irreversible MCT. */
 115 /* </summary> */
 116 void opj_mct_encode_real(
 117                 OPJ_INT32* restrict c0,
 118                 OPJ_INT32* restrict c1,
 119                 OPJ_INT32* restrict c2,
 120                 OPJ_UINT32 n)
 121 {
 122         OPJ_UINT32 i;
 123         for(i = 0; i < n; ++i) {
 124                 OPJ_INT32 r = c0[i];
 125                 OPJ_INT32 g = c1[i];
 126                 OPJ_INT32 b = c2[i];
 127                 OPJ_INT32 y =  opj_int_fix_mul(r, 2449) + opj_int_fix_mul(g, 4809) + opj_int_fix_mul(b, 934);
 128                 OPJ_INT32 u = -opj_int_fix_mul(r, 1382) - opj_int_fix_mul(g, 2714) + opj_int_fix_mul(b, 4096);
 129                 OPJ_INT32 v =  opj_int_fix_mul(r, 4096) - opj_int_fix_mul(g, 3430) - opj_int_fix_mul(b, 666);
 130                 c0[i] = y;
 131                 c1[i] = u;
 132                 c2[i] = v;
 133         }
 134 }
 135
 136 /* <summary> */
 137 /* Inverse irreversible MCT. */
 138 /* </summary> */
 139 void opj_mct_decode_real(
 140                 OPJ_FLOAT32* restrict c0,
 141                 OPJ_FLOAT32* restrict c1,
 142                 OPJ_FLOAT32* restrict c2,
 143                 OPJ_UINT32 n)
 144 {
 145         OPJ_UINT32 i;
 146 #ifdef __SSE__
 147         __m128 vrv, vgu, vgv, vbu;
 148         vrv = _mm_set1_ps(1.402f);
 149         vgu = _mm_set1_ps(0.34413f);
 150         vgv = _mm_set1_ps(0.71414f);
 151         vbu = _mm_set1_ps(1.772f);
 152         for (i = 0; i < (n >> 3); ++i) {
 153                 __m128 vy, vu, vv;
 154                 __m128 vr, vg, vb;
 155
 156                 vy = _mm_load_ps(c0);
 157                 vu = _mm_load_ps(c1);
 158                 vv = _mm_load_ps(c2);
 159                 vr = _mm_add_ps(vy, _mm_mul_ps(vv, vrv));
 160                 vg = _mm_sub_ps(_mm_sub_ps(vy, _mm_mul_ps(vu, vgu)), _mm_mul_ps(vv, vgv));
 161                 vb = _mm_add_ps(vy, _mm_mul_ps(vu, vbu));
 162                 _mm_store_ps(c0, vr);
 163                 _mm_store_ps(c1, vg);
 164                 _mm_store_ps(c2, vb);
 165                 c0 += 4;
 166                 c1 += 4;
 167                 c2 += 4;
 168
 169                 vy = _mm_load_ps(c0);
 170                 vu = _mm_load_ps(c1);
 171                 vv = _mm_load_ps(c2);
 172                 vr = _mm_add_ps(vy, _mm_mul_ps(vv, vrv));
 173                 vg = _mm_sub_ps(_mm_sub_ps(vy, _mm_mul_ps(vu, vgu)), _mm_mul_ps(vv, vgv));
 174                 vb = _mm_add_ps(vy, _mm_mul_ps(vu, vbu));
 175                 _mm_store_ps(c0, vr);
 176                 _mm_store_ps(c1, vg);
 177                 _mm_store_ps(c2, vb);
 178                 c0 += 4;
 179                 c1 += 4;
 180                 c2 += 4;
 181         }
 182         n &= 7;
 183 #endif
 184         for(i = 0; i < n; ++i) {
 185                 OPJ_FLOAT32 y = c0[i];
 186                 OPJ_FLOAT32 u = c1[i];
 187                 OPJ_FLOAT32 v = c2[i];
 188                 OPJ_FLOAT32 r = y + (v * 1.402f);
 189                 OPJ_FLOAT32 g = y - (u * 0.34413f) - (v * (0.71414f));
 190                 OPJ_FLOAT32 b = y + (u * 1.772f);
 191                 c0[i] = r;
 192                 c1[i] = g;
 193                 c2[i] = b;
 194         }
 195 }
 196
 197 /* <summary> */
 198 /* Get norm of basis function of irreversible MCT. */
 199 /* </summary> */
 200 OPJ_FLOAT64 opj_mct_getnorm_real(OPJ_UINT32 compno) {
 201         return opj_mct_norms_real[compno];
 202 }
 203
 204
 205 OPJ_BOOL opj_mct_encode_custom(
 206                                            OPJ_BYTE * pCodingdata,
 207                                            OPJ_UINT32 n,
 208                                            OPJ_BYTE ** pData,
 209                                            OPJ_UINT32 pNbComp,
 210                                            OPJ_UINT32 isSigned)
 211 {
 212         OPJ_FLOAT32 * lMct = (OPJ_FLOAT32 *) pCodingdata;
 213         OPJ_UINT32 i;
 214         OPJ_UINT32 j;
 215         OPJ_UINT32 k;
 216         OPJ_UINT32 lNbMatCoeff = pNbComp * pNbComp;
 217         OPJ_INT32 * lCurrentData = 00;
 218         OPJ_INT32 * lCurrentMatrix = 00;
 219         OPJ_INT32 ** lData = (OPJ_INT32 **) pData;
 220         OPJ_UINT32 lMultiplicator = 1 << 13;
 221         OPJ_INT32 * lMctPtr;
 222
 223     OPJ_ARG_NOT_USED(isSigned);
 224
 225         lCurrentData = (OPJ_INT32 *) opj_malloc((pNbComp + lNbMatCoeff) * sizeof(OPJ_INT32));
 226         if (! lCurrentData) {
 227                 return OPJ_FALSE;
 228         }
 229
 230         lCurrentMatrix = lCurrentData + pNbComp;
 231
 232         for (i =0;i<lNbMatCoeff;++i) {
 233                 lCurrentMatrix[i] = (OPJ_INT32) (*(lMct++) * lMultiplicator);
 234         }
 235
 236         for (i = 0; i < n; ++i)  {
 237                 lMctPtr = lCurrentMatrix;
 238                 for (j=0;j<pNbComp;++j) {
 239                         lCurrentData[j] = (*(lData[j]));
 240                 }
 241
 242                 for (j=0;j<pNbComp;++j) {
 243                         *(lData[j]) = 0;
 244                         for (k=0;k<pNbComp;++k) {
 245                                 *(lData[j]) += opj_int_fix_mul(*lMctPtr, lCurrentData[k]);
 246                                 ++lMctPtr;
 247                         }
 248
 249                         ++lData[j];
 250                 }
 251         }
 252
 253         opj_free(lCurrentData);
 254
 255         return OPJ_TRUE;
 256 }
 257
 258 OPJ_BOOL opj_mct_decode_custom(
 259                                            OPJ_BYTE * pDecodingData,
 260                                            OPJ_UINT32 n,
 261                                            OPJ_BYTE ** pData,
 262                                            OPJ_UINT32 pNbComp,
 263                                            OPJ_UINT32 isSigned)
 264 {
 265         OPJ_FLOAT32 * lMct;
 266         OPJ_UINT32 i;
 267         OPJ_UINT32 j;
 268         OPJ_UINT32 k;
 269
 270         OPJ_FLOAT32 * lCurrentData = 00;
 271         OPJ_FLOAT32 * lCurrentResult = 00;
 272         OPJ_FLOAT32 ** lData = (OPJ_FLOAT32 **) pData;
 273
 274     OPJ_ARG_NOT_USED(isSigned);
 275
 276         lCurrentData = (OPJ_FLOAT32 *) opj_malloc (2 * pNbComp * sizeof(OPJ_FLOAT32));
 277         if (! lCurrentData) {
 278                 return OPJ_FALSE;
 279         }
 280         lCurrentResult = lCurrentData + pNbComp;
 281
 282         for (i = 0; i < n; ++i) {
 283                 lMct = (OPJ_FLOAT32 *) pDecodingData;
 284                 for (j=0;j<pNbComp;++j) {
 285                         lCurrentData[j] = (OPJ_FLOAT32) (*(lData[j]));
 286                 }
 287                 for (j=0;j<pNbComp;++j) {
 288                         lCurrentResult[j] = 0;
 289                         for     (k=0;k<pNbComp;++k)     {
 290                                 lCurrentResult[j] += *(lMct++) * lCurrentData[k];
 291                         }
 292                         *(lData[j]++) = (OPJ_FLOAT32) (lCurrentResult[j]);
 293                 }
 294         }
 295         opj_free(lCurrentData);
 296         return OPJ_TRUE;
 297 }
 298
 299 void opj_calculate_norms(       OPJ_FLOAT64 * pNorms,
 300                                                         OPJ_UINT32 pNbComps,
 301                                                         OPJ_FLOAT32 * pMatrix)
 302 {
 303         OPJ_UINT32 i,j,lIndex;
 304         OPJ_FLOAT32 lCurrentValue;
 305         OPJ_FLOAT64 * lNorms = (OPJ_FLOAT64 *) pNorms;
 306         OPJ_FLOAT32 * lMatrix = (OPJ_FLOAT32 *) pMatrix;
 307
 308         for     (i=0;i<pNbComps;++i) {
 309                 lNorms[i] = 0;
 310                 lIndex = i;
 311
 312                 for     (j=0;j<pNbComps;++j) {
 313                         lCurrentValue = lMatrix[lIndex];
 314                         lIndex += pNbComps;
 315                         lNorms[i] += lCurrentValue * lCurrentValue;
 316                 }
 317                 lNorms[i] = sqrt(lNorms[i]);
 318         }
 319 }