diff options
| author | Mathieu Malaterre <mathieu.malaterre@gmail.com> | 2012-09-28 08:11:41 +0000 |
|---|---|---|
| committer | Mathieu Malaterre <mathieu.malaterre@gmail.com> | 2012-09-28 08:11:41 +0000 |
| commit | d518970039a19a2a9b6d2bdd592cc88a43897bbb (patch) | |
| tree | 57bac2cf7e63e9352228231062763baac627563c /src/lib/openjp2/mct.c | |
| parent | 8363a6ab1e031bb4b2e40a92e56efd40fdab1aa1 (diff) | |
[trunk] Start FolderReorgProposal task
Update issue 177
Diffstat (limited to 'src/lib/openjp2/mct.c')
| -rw-r--r-- | src/lib/openjp2/mct.c | 333 |
1 files changed, 333 insertions, 0 deletions
diff --git a/src/lib/openjp2/mct.c b/src/lib/openjp2/mct.c new file mode 100644 index 00000000..fcc105e6 --- /dev/null +++ b/src/lib/openjp2/mct.c @@ -0,0 +1,333 @@ +/* + * Copyright (c) 2002-2007, Communications and Remote Sensing Laboratory, Universite catholique de Louvain (UCL), Belgium + * Copyright (c) 2002-2007, Professor Benoit Macq + * Copyright (c) 2001-2003, David Janssens + * Copyright (c) 2002-2003, Yannick Verschueren + * Copyright (c) 2003-2007, Francois-Olivier Devaux and Antonin Descampe + * Copyright (c) 2005, Herve Drolon, FreeImage Team + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS' + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#ifdef __SSE__ +#include <xmmintrin.h> +#endif + +#include "opj_includes.h" + +/* <summary> */ +/* This table contains the norms of the basis function of the reversible MCT. */ +/* </summary> */ +static const double mct_norms[3] = { 1.732, .8292, .8292 }; + +/* <summary> */ +/* This table contains the norms of the basis function of the irreversible MCT. */ +/* </summary> */ +static const double mct_norms_real[3] = { 1.732, 1.805, 1.573 }; + +const OPJ_FLOAT64 * get_mct_norms () +{ + return mct_norms; +} + +const OPJ_FLOAT64 * get_mct_norms_real () +{ + return mct_norms_real; +} + +/* <summary> */ +/* Foward reversible MCT. */ +/* </summary> */ +void mct_encode( + int* restrict c0, + int* restrict c1, + int* restrict c2, + int n) +{ + int i; + for(i = 0; i < n; ++i) { + int r = c0[i]; + int g = c1[i]; + int b = c2[i]; + int y = (r + (g * 2) + b) >> 2; + int u = b - g; + int v = r - g; + c0[i] = y; + c1[i] = u; + c2[i] = v; + } +} + +/* <summary> */ +/* Inverse reversible MCT. */ +/* </summary> */ +void mct_decode( + int* restrict c0, + int* restrict c1, + int* restrict c2, + int n) +{ + int i; + for (i = 0; i < n; ++i) { + int y = c0[i]; + int u = c1[i]; + int v = c2[i]; + int g = y - ((u + v) >> 2); + int r = v + g; + int b = u + g; + c0[i] = r; + c1[i] = g; + c2[i] = b; + } +} + +/* <summary> */ +/* Get norm of basis function of reversible MCT. */ +/* </summary> */ +double mct_getnorm(int compno) { + return mct_norms[compno]; +} + +/* <summary> */ +/* Foward irreversible MCT. */ +/* </summary> */ +void mct_encode_real( + int* restrict c0, + int* restrict c1, + int* restrict c2, + int n) +{ + int i; + for(i = 0; i < n; ++i) { + int r = c0[i]; + int g = c1[i]; + int b = c2[i]; + int y = fix_mul(r, 2449) + fix_mul(g, 4809) + fix_mul(b, 934); + int u = -fix_mul(r, 1382) - fix_mul(g, 2714) + fix_mul(b, 4096); + int v = fix_mul(r, 4096) - fix_mul(g, 3430) - fix_mul(b, 666); + c0[i] = y; + c1[i] = u; + c2[i] = v; + } +} + +/* <summary> */ +/* Inverse irreversible MCT. */ +/* </summary> */ +void mct_decode_real( + float* restrict c0, + float* restrict c1, + float* restrict c2, + int n) +{ + int i; +#ifdef __SSE__ + __m128 vrv, vgu, vgv, vbu; + vrv = _mm_set1_ps(1.402f); + vgu = _mm_set1_ps(0.34413f); + vgv = _mm_set1_ps(0.71414f); + vbu = _mm_set1_ps(1.772f); + for (i = 0; i < (n >> 3); ++i) { + __m128 vy, vu, vv; + __m128 vr, vg, vb; + + vy = _mm_load_ps(c0); + vu = _mm_load_ps(c1); + vv = _mm_load_ps(c2); + vr = _mm_add_ps(vy, _mm_mul_ps(vv, vrv)); + vg = _mm_sub_ps(_mm_sub_ps(vy, _mm_mul_ps(vu, vgu)), _mm_mul_ps(vv, vgv)); + vb = _mm_add_ps(vy, _mm_mul_ps(vu, vbu)); + _mm_store_ps(c0, vr); + _mm_store_ps(c1, vg); + _mm_store_ps(c2, vb); + c0 += 4; + c1 += 4; + c2 += 4; + + vy = _mm_load_ps(c0); + vu = _mm_load_ps(c1); + vv = _mm_load_ps(c2); + vr = _mm_add_ps(vy, _mm_mul_ps(vv, vrv)); + vg = _mm_sub_ps(_mm_sub_ps(vy, _mm_mul_ps(vu, vgu)), _mm_mul_ps(vv, vgv)); + vb = _mm_add_ps(vy, _mm_mul_ps(vu, vbu)); + _mm_store_ps(c0, vr); + _mm_store_ps(c1, vg); + _mm_store_ps(c2, vb); + c0 += 4; + c1 += 4; + c2 += 4; + } + n &= 7; +#endif + for(i = 0; i < n; ++i) { + float y = c0[i]; + float u = c1[i]; + float v = c2[i]; + float r = y + (v * 1.402f); + float g = y - (u * 0.34413f) - (v * (0.71414f)); + float b = y + (u * 1.772f); + c0[i] = r; + c1[i] = g; + c2[i] = b; + } +} + +/* <summary> */ +/* Get norm of basis function of irreversible MCT. */ +/* </summary> */ +double mct_getnorm_real(int compno) { + return mct_norms_real[compno]; +} + + +opj_bool mct_encode_custom( + /* MCT data */ + OPJ_BYTE * pCodingdata, + /* size of components */ + OPJ_UINT32 n, + /* components */ + OPJ_BYTE ** pData, + /* nb of components (i.e. size of pData) */ + OPJ_UINT32 pNbComp, + /* tells if the data is signed */ + OPJ_UINT32 isSigned) +{ + OPJ_FLOAT32 * lMct = (OPJ_FLOAT32 *) pCodingdata; + OPJ_UINT32 i; + OPJ_UINT32 j; + OPJ_UINT32 k; + OPJ_UINT32 lNbMatCoeff = pNbComp * pNbComp; + OPJ_INT32 * lCurrentData = 00; + OPJ_INT32 * lCurrentMatrix = 00; + OPJ_INT32 ** lData = (OPJ_INT32 **) pData; + OPJ_UINT32 lMultiplicator = 1 << 13; + OPJ_INT32 * lMctPtr; + + lCurrentData = (OPJ_INT32 *) opj_malloc((pNbComp + lNbMatCoeff) * sizeof(OPJ_INT32)); + if (! lCurrentData) { + return OPJ_FALSE; + } + + lCurrentMatrix = lCurrentData + pNbComp; + + for (i =0;i<lNbMatCoeff;++i) { + lCurrentMatrix[i] = (OPJ_INT32) (*(lMct++) * lMultiplicator); + } + + for (i = 0; i < n; ++i) { + lMctPtr = lCurrentMatrix; + for (j=0;j<pNbComp;++j) { + lCurrentData[j] = (*(lData[j])); + } + + for (j=0;j<pNbComp;++j) { + *(lData[j]) = 0; + for (k=0;k<pNbComp;++k) { + *(lData[j]) += fix_mul(*lMctPtr, lCurrentData[k]); + ++lMctPtr; + } + + ++lData[j]; + } + } + + opj_free(lCurrentData); + + return OPJ_TRUE; +} + +opj_bool mct_decode_custom( + /* MCT data */ + OPJ_BYTE * pDecodingData, + /* size of components */ + OPJ_UINT32 n, + /* components */ + OPJ_BYTE ** pData, + /* nb of components (i.e. size of pData) */ + OPJ_UINT32 pNbComp, + /* tells if the data is signed */ + OPJ_UINT32 isSigned) +{ + OPJ_FLOAT32 * lMct; + OPJ_UINT32 i; + OPJ_UINT32 j; + OPJ_UINT32 k; + + OPJ_FLOAT32 * lCurrentData = 00; + OPJ_FLOAT32 * lCurrentResult = 00; + OPJ_FLOAT32 ** lData = (OPJ_FLOAT32 **) pData; + + lCurrentData = (OPJ_FLOAT32 *) opj_malloc (2 * pNbComp * sizeof(OPJ_FLOAT32)); + if + (! lCurrentData) + { + return OPJ_FALSE; + } + lCurrentResult = lCurrentData + pNbComp; + + for + (i = 0; i < n; ++i) + { + lMct = (OPJ_FLOAT32 *) pDecodingData; + for + (j=0;j<pNbComp;++j) + { + lCurrentData[j] = (OPJ_FLOAT32) (*(lData[j])); + } + for + (j=0;j<pNbComp;++j) + { + lCurrentResult[j] = 0; + for + (k=0;k<pNbComp;++k) + { + lCurrentResult[j] += *(lMct++) * lCurrentData[k]; + } + *(lData[j]++) = (OPJ_FLOAT32) (lCurrentResult[j]); + } + } + opj_free(lCurrentData); + return OPJ_TRUE; +} + +void opj_calculate_norms( OPJ_FLOAT64 * pNorms, + OPJ_UINT32 pNbComps, + OPJ_FLOAT32 * pMatrix) +{ + OPJ_UINT32 i,j,lIndex; + OPJ_FLOAT32 lCurrentValue; + OPJ_FLOAT64 * lNorms = (OPJ_FLOAT64 *) pNorms; + OPJ_FLOAT32 * lMatrix = (OPJ_FLOAT32 *) pMatrix; + + for (i=0;i<pNbComps;++i) { + lNorms[i] = 0; + lIndex = i; + + for (j=0;j<pNbComps;++j) { + lCurrentValue = lMatrix[lIndex]; + lIndex += pNbComps; + lNorms[i] += lCurrentValue * lCurrentValue; + } + lNorms[i] = sqrt(lNorms[i]); + } +} |
