src/lib/openjp2/mct.c

   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* OPJ_RESTRICT c0,
  78     OPJ_INT32* OPJ_RESTRICT c1,
  79     OPJ_INT32* OPJ_RESTRICT c2,
  80     OPJ_SIZE_T n)
  81 {
  82     OPJ_SIZE_T i;
  83     const OPJ_SIZE_T len = n;
  84     /* buffer are aligned on 16 bytes */
  85     assert(((size_t)c0 & 0xf) == 0);
  86     assert(((size_t)c1 & 0xf) == 0);
  87     assert(((size_t)c2 & 0xf) == 0);
  88
  89     for (i = 0; i < (len & ~3U); i += 4) {
  90         __m128i y, u, v;
  91         __m128i r = _mm_load_si128((const __m128i *) & (c0[i]));
  92         __m128i g = _mm_load_si128((const __m128i *) & (c1[i]));
  93         __m128i b = _mm_load_si128((const __m128i *) & (c2[i]));
  94         y = _mm_add_epi32(g, g);
  95         y = _mm_add_epi32(y, b);
  96         y = _mm_add_epi32(y, r);
  97         y = _mm_srai_epi32(y, 2);
  98         u = _mm_sub_epi32(b, g);
  99         v = _mm_sub_epi32(r, g);
 100         _mm_store_si128((__m128i *) & (c0[i]), y);
 101         _mm_store_si128((__m128i *) & (c1[i]), u);
 102         _mm_store_si128((__m128i *) & (c2[i]), v);
 103     }
 104
 105     for (; i < len; ++i) {
 106         OPJ_INT32 r = c0[i];
 107         OPJ_INT32 g = c1[i];
 108         OPJ_INT32 b = c2[i];
 109         OPJ_INT32 y = (r + (g * 2) + b) >> 2;
 110         OPJ_INT32 u = b - g;
 111         OPJ_INT32 v = r - g;
 112         c0[i] = y;
 113         c1[i] = u;
 114         c2[i] = v;
 115     }
 116 }
 117 #else
 118 void opj_mct_encode(
 119     OPJ_INT32* OPJ_RESTRICT c0,
 120     OPJ_INT32* OPJ_RESTRICT c1,
 121     OPJ_INT32* OPJ_RESTRICT c2,
 122     OPJ_SIZE_T n)
 123 {
 124     OPJ_SIZE_T i;
 125     const OPJ_SIZE_T len = n;
 126
 127     for (i = 0; i < len; ++i) {
 128         OPJ_INT32 r = c0[i];
 129         OPJ_INT32 g = c1[i];
 130         OPJ_INT32 b = c2[i];
 131         OPJ_INT32 y = (r + (g * 2) + b) >> 2;
 132         OPJ_INT32 u = b - g;
 133         OPJ_INT32 v = r - g;
 134         c0[i] = y;
 135         c1[i] = u;
 136         c2[i] = v;
 137     }
 138 }
 139 #endif
 140
 141 /* <summary> */
 142 /* Inverse reversible MCT. */
 143 /* </summary> */
 144 #ifdef __SSE2__
 145 void opj_mct_decode(
 146     OPJ_INT32* OPJ_RESTRICT c0,
 147     OPJ_INT32* OPJ_RESTRICT c1,
 148     OPJ_INT32* OPJ_RESTRICT c2,
 149     OPJ_SIZE_T n)
 150 {
 151     OPJ_SIZE_T i;
 152     const OPJ_SIZE_T len = n;
 153
 154     for (i = 0; i < (len & ~3U); i += 4) {
 155         __m128i r, g, b;
 156         __m128i y = _mm_load_si128((const __m128i *) & (c0[i]));
 157         __m128i u = _mm_load_si128((const __m128i *) & (c1[i]));
 158         __m128i v = _mm_load_si128((const __m128i *) & (c2[i]));
 159         g = y;
 160         g = _mm_sub_epi32(g, _mm_srai_epi32(_mm_add_epi32(u, v), 2));
 161         r = _mm_add_epi32(v, g);
 162         b = _mm_add_epi32(u, g);
 163         _mm_store_si128((__m128i *) & (c0[i]), r);
 164         _mm_store_si128((__m128i *) & (c1[i]), g);
 165         _mm_store_si128((__m128i *) & (c2[i]), b);
 166     }
 167     for (; i < len; ++i) {
 168         OPJ_INT32 y = c0[i];
 169         OPJ_INT32 u = c1[i];
 170         OPJ_INT32 v = c2[i];
 171         OPJ_INT32 g = y - ((u + v) >> 2);
 172         OPJ_INT32 r = v + g;
 173         OPJ_INT32 b = u + g;
 174         c0[i] = r;
 175         c1[i] = g;
 176         c2[i] = b;
 177     }
 178 }
 179 #else
 180 void opj_mct_decode(
 181     OPJ_INT32* OPJ_RESTRICT c0,
 182     OPJ_INT32* OPJ_RESTRICT c1,
 183     OPJ_INT32* OPJ_RESTRICT c2,
 184     OPJ_SIZE_T n)
 185 {
 186     OPJ_SIZE_T i;
 187     for (i = 0; i < n; ++i) {
 188         OPJ_INT32 y = c0[i];
 189         OPJ_INT32 u = c1[i];
 190         OPJ_INT32 v = c2[i];
 191         OPJ_INT32 g = y - ((u + v) >> 2);
 192         OPJ_INT32 r = v + g;
 193         OPJ_INT32 b = u + g;
 194         c0[i] = r;
 195         c1[i] = g;
 196         c2[i] = b;
 197     }
 198 }
 199 #endif
 200
 201 /* <summary> */
 202 /* Get norm of basis function of reversible MCT. */
 203 /* </summary> */
 204 OPJ_FLOAT64 opj_mct_getnorm(OPJ_UINT32 compno)
 205 {
 206     return opj_mct_norms[compno];
 207 }
 208
 209 /* <summary> */
 210 /* Forward irreversible MCT. */
 211 /* </summary> */
 212 void opj_mct_encode_real(
 213     OPJ_FLOAT32* OPJ_RESTRICT c0,
 214     OPJ_FLOAT32* OPJ_RESTRICT c1,
 215     OPJ_FLOAT32* OPJ_RESTRICT c2,
 216     OPJ_SIZE_T n)
 217 {
 218     OPJ_SIZE_T i;
 219     for (i = 0; i < n; ++i) {
 220         OPJ_FLOAT32 r = c0[i];
 221         OPJ_FLOAT32 g = c1[i];
 222         OPJ_FLOAT32 b = c2[i];
 223         OPJ_FLOAT32 y = 0.299f * r + 0.587f * g + 0.114f * b;
 224         OPJ_FLOAT32 u = -0.16875f * r - 0.331260f * g + 0.5f * b;
 225         OPJ_FLOAT32 v = 0.5f * r - 0.41869f * g - 0.08131f * b;
 226         c0[i] = y;
 227         c1[i] = u;
 228         c2[i] = v;
 229     }
 230 }
 231
 232 /* <summary> */
 233 /* Inverse irreversible MCT. */
 234 /* </summary> */
 235 void opj_mct_decode_real(
 236     OPJ_FLOAT32* OPJ_RESTRICT c0,
 237     OPJ_FLOAT32* OPJ_RESTRICT c1,
 238     OPJ_FLOAT32* OPJ_RESTRICT c2,
 239     OPJ_SIZE_T n)
 240 {
 241     OPJ_SIZE_T i;
 242 #ifdef __SSE__
 243     __m128 vrv, vgu, vgv, vbu;
 244     vrv = _mm_set1_ps(1.402f);
 245     vgu = _mm_set1_ps(0.34413f);
 246     vgv = _mm_set1_ps(0.71414f);
 247     vbu = _mm_set1_ps(1.772f);
 248     for (i = 0; i < (n >> 3); ++i) {
 249         __m128 vy, vu, vv;
 250         __m128 vr, vg, vb;
 251
 252         vy = _mm_load_ps(c0);
 253         vu = _mm_load_ps(c1);
 254         vv = _mm_load_ps(c2);
 255         vr = _mm_add_ps(vy, _mm_mul_ps(vv, vrv));
 256         vg = _mm_sub_ps(_mm_sub_ps(vy, _mm_mul_ps(vu, vgu)), _mm_mul_ps(vv, vgv));
 257         vb = _mm_add_ps(vy, _mm_mul_ps(vu, vbu));
 258         _mm_store_ps(c0, vr);
 259         _mm_store_ps(c1, vg);
 260         _mm_store_ps(c2, vb);
 261         c0 += 4;
 262         c1 += 4;
 263         c2 += 4;
 264
 265         vy = _mm_load_ps(c0);
 266         vu = _mm_load_ps(c1);
 267         vv = _mm_load_ps(c2);
 268         vr = _mm_add_ps(vy, _mm_mul_ps(vv, vrv));
 269         vg = _mm_sub_ps(_mm_sub_ps(vy, _mm_mul_ps(vu, vgu)), _mm_mul_ps(vv, vgv));
 270         vb = _mm_add_ps(vy, _mm_mul_ps(vu, vbu));
 271         _mm_store_ps(c0, vr);
 272         _mm_store_ps(c1, vg);
 273         _mm_store_ps(c2, vb);
 274         c0 += 4;
 275         c1 += 4;
 276         c2 += 4;
 277     }
 278     n &= 7;
 279 #endif
 280     for (i = 0; i < n; ++i) {
 281         OPJ_FLOAT32 y = c0[i];
 282         OPJ_FLOAT32 u = c1[i];
 283         OPJ_FLOAT32 v = c2[i];
 284         OPJ_FLOAT32 r = y + (v * 1.402f);
 285         OPJ_FLOAT32 g = y - (u * 0.34413f) - (v * (0.71414f));
 286         OPJ_FLOAT32 b = y + (u * 1.772f);
 287         c0[i] = r;
 288         c1[i] = g;
 289         c2[i] = b;
 290     }
 291 }
 292
 293 /* <summary> */
 294 /* Get norm of basis function of irreversible MCT. */
 295 /* </summary> */
 296 OPJ_FLOAT64 opj_mct_getnorm_real(OPJ_UINT32 compno)
 297 {
 298     return opj_mct_norms_real[compno];
 299 }
 300
 301
 302 OPJ_BOOL opj_mct_encode_custom(
 303     OPJ_BYTE * pCodingdata,
 304     OPJ_SIZE_T n,
 305     OPJ_BYTE ** pData,
 306     OPJ_UINT32 pNbComp,
 307     OPJ_UINT32 isSigned)
 308 {
 309     OPJ_FLOAT32 * lMct = (OPJ_FLOAT32 *) pCodingdata;
 310     OPJ_SIZE_T i;
 311     OPJ_UINT32 j;
 312     OPJ_UINT32 k;
 313     OPJ_UINT32 lNbMatCoeff = pNbComp * pNbComp;
 314     OPJ_INT32 * lCurrentData = 00;
 315     OPJ_INT32 * lCurrentMatrix = 00;
 316     OPJ_INT32 ** lData = (OPJ_INT32 **) pData;
 317     OPJ_UINT32 lMultiplicator = 1 << 13;
 318     OPJ_INT32 * lMctPtr;
 319
 320     OPJ_ARG_NOT_USED(isSigned);
 321
 322     lCurrentData = (OPJ_INT32 *) opj_malloc((pNbComp + lNbMatCoeff) * sizeof(
 323             OPJ_INT32));
 324     if (! lCurrentData) {
 325         return OPJ_FALSE;
 326     }
 327
 328     lCurrentMatrix = lCurrentData + pNbComp;
 329
 330     for (i = 0; i < lNbMatCoeff; ++i) {
 331         lCurrentMatrix[i] = (OPJ_INT32)(*(lMct++) * (OPJ_FLOAT32)lMultiplicator);
 332     }
 333
 334     for (i = 0; i < n; ++i)  {
 335         lMctPtr = lCurrentMatrix;
 336         for (j = 0; j < pNbComp; ++j) {
 337             lCurrentData[j] = (*(lData[j]));
 338         }
 339
 340         for (j = 0; j < pNbComp; ++j) {
 341             *(lData[j]) = 0;
 342             for (k = 0; k < pNbComp; ++k) {
 343                 *(lData[j]) += opj_int_fix_mul(*lMctPtr, lCurrentData[k]);
 344                 ++lMctPtr;
 345             }
 346
 347             ++lData[j];
 348         }
 349     }
 350
 351     opj_free(lCurrentData);
 352
 353     return OPJ_TRUE;
 354 }
 355
 356 OPJ_BOOL opj_mct_decode_custom(
 357     OPJ_BYTE * pDecodingData,
 358     OPJ_SIZE_T n,
 359     OPJ_BYTE ** pData,
 360     OPJ_UINT32 pNbComp,
 361     OPJ_UINT32 isSigned)
 362 {
 363     OPJ_FLOAT32 * lMct;
 364     OPJ_SIZE_T i;
 365     OPJ_UINT32 j;
 366     OPJ_UINT32 k;
 367
 368     OPJ_FLOAT32 * lCurrentData = 00;
 369     OPJ_FLOAT32 * lCurrentResult = 00;
 370     OPJ_FLOAT32 ** lData = (OPJ_FLOAT32 **) pData;
 371
 372     OPJ_ARG_NOT_USED(isSigned);
 373
 374     lCurrentData = (OPJ_FLOAT32 *) opj_malloc(2 * pNbComp * sizeof(OPJ_FLOAT32));
 375     if (! lCurrentData) {
 376         return OPJ_FALSE;
 377     }
 378     lCurrentResult = lCurrentData + pNbComp;
 379
 380     for (i = 0; i < n; ++i) {
 381         lMct = (OPJ_FLOAT32 *) pDecodingData;
 382         for (j = 0; j < pNbComp; ++j) {
 383             lCurrentData[j] = (OPJ_FLOAT32)(*(lData[j]));
 384         }
 385         for (j = 0; j < pNbComp; ++j) {
 386             lCurrentResult[j] = 0;
 387             for (k = 0; k < pNbComp; ++k) {
 388                 lCurrentResult[j] += *(lMct++) * lCurrentData[k];
 389             }
 390             *(lData[j]++) = (OPJ_FLOAT32)(lCurrentResult[j]);
 391         }
 392     }
 393     opj_free(lCurrentData);
 394     return OPJ_TRUE;
 395 }
 396
 397 void opj_calculate_norms(OPJ_FLOAT64 * pNorms,
 398                          OPJ_UINT32 pNbComps,
 399                          OPJ_FLOAT32 * pMatrix)
 400 {
 401     OPJ_UINT32 i, j, lIndex;
 402     OPJ_FLOAT32 lCurrentValue;
 403     OPJ_FLOAT64 * lNorms = (OPJ_FLOAT64 *) pNorms;
 404     OPJ_FLOAT32 * lMatrix = (OPJ_FLOAT32 *) pMatrix;
 405
 406     for (i = 0; i < pNbComps; ++i) {
 407         lNorms[i] = 0;
 408         lIndex = i;
 409
 410         for (j = 0; j < pNbComps; ++j) {
 411             lCurrentValue = lMatrix[lIndex];
 412             lIndex += pNbComps;
 413             lNorms[i] += lCurrentValue * lCurrentValue;
 414         }
 415         lNorms[i] = sqrt(lNorms[i]);
 416     }
 417 }