From: Mickael Savinaud Date: Wed, 24 Oct 2012 15:19:51 +0000 (+0000) Subject: [trunk] add the support of complex mct encoding when we setup the j2k encoder X-Git-Tag: version.2.0~90 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=76947f007466ea68cddabd915633d9670dabdb98;p=openjpeg [trunk] add the support of complex mct encoding when we setup the j2k encoder --- diff --git a/src/lib/openjp2/CMakeLists.txt b/src/lib/openjp2/CMakeLists.txt index 33281588..ed8f681c 100644 --- a/src/lib/openjp2/CMakeLists.txt +++ b/src/lib/openjp2/CMakeLists.txt @@ -15,6 +15,7 @@ set(OPENJPEG_SRCS ${CMAKE_CURRENT_SOURCE_DIR}/dwt.c ${CMAKE_CURRENT_SOURCE_DIR}/event.c ${CMAKE_CURRENT_SOURCE_DIR}/image.c + ${CMAKE_CURRENT_SOURCE_DIR}/invert.c ${CMAKE_CURRENT_SOURCE_DIR}/j2k.c ${CMAKE_CURRENT_SOURCE_DIR}/jp2.c ${CMAKE_CURRENT_SOURCE_DIR}/jpt.c diff --git a/src/lib/openjp2/invert.c b/src/lib/openjp2/invert.c new file mode 100644 index 00000000..31088952 --- /dev/null +++ b/src/lib/openjp2/invert.c @@ -0,0 +1,286 @@ +/* + * Copyright (c) 2008, Jerome Fimes, Communications & Systemes + * 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. + */ + +#include "opj_includes.h" + +/** + * LUP decomposition + */ +static opj_bool opj_lupDecompose(OPJ_FLOAT32 * matrix, + OPJ_UINT32 * permutations, + OPJ_FLOAT32 * p_swap_area, + OPJ_UINT32 n); +/** + * LUP solving + */ +static void opj_lupSolve(OPJ_FLOAT32 * pResult, + OPJ_FLOAT32* pMatrix, + OPJ_FLOAT32* pVector, + OPJ_UINT32* pPermutations, + OPJ_UINT32 n, + OPJ_FLOAT32 * p_intermediate_data); + +/** + *LUP inversion (call with the result of lupDecompose) + */ +static void opj_lupInvert ( OPJ_FLOAT32 * pSrcMatrix, + OPJ_FLOAT32 * pDestMatrix, + OPJ_UINT32 n, + OPJ_UINT32 * pPermutations, + OPJ_FLOAT32 * p_src_temp, + OPJ_FLOAT32 * p_dest_temp, + OPJ_FLOAT32 * p_swap_area); + +/* +========================================================== + Matric inversion interface +========================================================== +*/ +/** + * Matrix inversion. + */ +opj_bool opj_matrix_inversion_f(OPJ_FLOAT32 * pSrcMatrix, + OPJ_FLOAT32 * pDestMatrix, + OPJ_UINT32 n) +{ + OPJ_BYTE * l_data = 00; + OPJ_UINT32 l_permutation_size = n * sizeof(OPJ_UINT32); + OPJ_UINT32 l_swap_size = n * sizeof(OPJ_FLOAT32); + OPJ_UINT32 l_total_size = l_permutation_size + 3 * l_swap_size; + OPJ_UINT32 * lPermutations = 00; + OPJ_FLOAT32 * l_double_data = 00; + + l_data = (OPJ_BYTE *) opj_malloc(l_total_size); + if (l_data == 0) { + return OPJ_FALSE; + } + lPermutations = (OPJ_UINT32 *) l_data; + l_double_data = (OPJ_FLOAT32 *) (l_data + l_permutation_size); + memset(lPermutations,0,l_permutation_size); + + if(! opj_lupDecompose(pSrcMatrix,lPermutations,l_double_data,n)) { + opj_free(l_data); + return OPJ_FALSE; + } + + opj_lupInvert(pSrcMatrix,pDestMatrix,n,lPermutations,l_double_data,l_double_data + n,l_double_data + 2*n); + opj_free(l_data); + + return OPJ_TRUE; +} + + +/* +========================================================== + Local functions +========================================================== +*/ +opj_bool opj_lupDecompose(OPJ_FLOAT32 * matrix,OPJ_UINT32 * permutations, + OPJ_FLOAT32 * p_swap_area, + OPJ_UINT32 n) +{ + OPJ_UINT32 * tmpPermutations = permutations; + OPJ_UINT32 * dstPermutations; + OPJ_UINT32 k2=0,t; + OPJ_FLOAT32 temp; + OPJ_UINT32 i,j,k; + OPJ_FLOAT32 p; + OPJ_UINT32 lLastColum = n - 1; + OPJ_UINT32 lSwapSize = n * sizeof(OPJ_FLOAT32); + OPJ_FLOAT32 * lTmpMatrix = matrix; + OPJ_FLOAT32 * lColumnMatrix,* lDestMatrix; + OPJ_UINT32 offset = 1; + OPJ_UINT32 lStride = n-1; + + //initialize permutations + for (i = 0; i < n; ++i) + { + *tmpPermutations++ = i; + } + // now make a pivot with colum switch + tmpPermutations = permutations; + for (k = 0; k < lLastColum; ++k) { + p = 0.0; + + // take the middle element + lColumnMatrix = lTmpMatrix + k; + + // make permutation with the biggest value in the column + for (i = k; i < n; ++i) { + temp = ((*lColumnMatrix > 0) ? *lColumnMatrix : -(*lColumnMatrix)); + if (temp > p) { + p = temp; + k2 = i; + } + // next line + lColumnMatrix += n; + } + + // a whole rest of 0 -> non singular + if (p == 0.0) { + return OPJ_FALSE; + } + + // should we permute ? + if (k2 != k) { + //exchange of line + // k2 > k + dstPermutations = tmpPermutations + k2 - k; + // swap indices + t = *tmpPermutations; + *tmpPermutations = *dstPermutations; + *dstPermutations = t; + + // and swap entire line. + lColumnMatrix = lTmpMatrix + (k2 - k) * n; + memcpy(p_swap_area,lColumnMatrix,lSwapSize); + memcpy(lColumnMatrix,lTmpMatrix,lSwapSize); + memcpy(lTmpMatrix,p_swap_area,lSwapSize); + } + + // now update data in the rest of the line and line after + lDestMatrix = lTmpMatrix + k; + lColumnMatrix = lDestMatrix + n; + // take the middle element + temp = *(lDestMatrix++); + + // now compute up data (i.e. coeff up of the diagonal). + for (i = offset; i < n; ++i) { + //lColumnMatrix; + // divide the lower column elements by the diagonal value + + // matrix[i][k] /= matrix[k][k]; + // p = matrix[i][k] + p = *lColumnMatrix / temp; + *(lColumnMatrix++) = p; + + for (j = /* k + 1 */ offset; j < n; ++j) { + // matrix[i][j] -= matrix[i][k] * matrix[k][j]; + *(lColumnMatrix++) -= p * (*(lDestMatrix++)); + } + // come back to the k+1th element + lDestMatrix -= lStride; + // go to kth element of the next line + lColumnMatrix += k; + } + + // offset is now k+2 + ++offset; + // 1 element less for stride + --lStride; + // next line + lTmpMatrix+=n; + // next permutation element + ++tmpPermutations; + } + return OPJ_TRUE; +} + +void opj_lupSolve (OPJ_FLOAT32 * pResult, + OPJ_FLOAT32 * pMatrix, + OPJ_FLOAT32 * pVector, + OPJ_UINT32* pPermutations, + OPJ_UINT32 n,OPJ_FLOAT32 * p_intermediate_data) +{ + OPJ_UINT32 i,j; + OPJ_FLOAT32 sum; + OPJ_FLOAT32 u; + OPJ_UINT32 lStride = n+1; + OPJ_FLOAT32 * lCurrentPtr; + OPJ_FLOAT32 * lIntermediatePtr; + OPJ_FLOAT32 * lDestPtr; + OPJ_FLOAT32 * lTmpMatrix; + OPJ_FLOAT32 * lLineMatrix = pMatrix; + OPJ_FLOAT32 * lBeginPtr = pResult + n - 1; + OPJ_FLOAT32 * lGeneratedData; + OPJ_UINT32 * lCurrentPermutationPtr = pPermutations; + + + lIntermediatePtr = p_intermediate_data; + lGeneratedData = p_intermediate_data + n - 1; + + for (i = 0; i < n; ++i) { + sum = 0.0; + lCurrentPtr = p_intermediate_data; + lTmpMatrix = lLineMatrix; + for (j = 1; j <= i; ++j) + { + // sum += matrix[i][j-1] * y[j-1]; + sum += (*(lTmpMatrix++)) * (*(lCurrentPtr++)); + } + //y[i] = pVector[pPermutations[i]] - sum; + *(lIntermediatePtr++) = pVector[*(lCurrentPermutationPtr++)] - sum; + lLineMatrix += n; + } + + // we take the last point of the matrix + lLineMatrix = pMatrix + n*n - 1; + + // and we take after the last point of the destination vector + lDestPtr = pResult + n; + + for (i = n - 1; i != -1 ; --i) { + sum = 0.0; + lTmpMatrix = lLineMatrix; + u = *(lTmpMatrix++); + lCurrentPtr = lDestPtr--; + for (j = i + 1; j < n; ++j) { + // sum += matrix[i][j] * x[j] + sum += (*(lTmpMatrix++)) * (*(lCurrentPtr++)); + } + //x[i] = (y[i] - sum) / u; + *(lBeginPtr--) = (*(lGeneratedData--) - sum) / u; + lLineMatrix -= lStride; + } +} + + +void opj_lupInvert (OPJ_FLOAT32 * pSrcMatrix, + OPJ_FLOAT32 * pDestMatrix, + OPJ_UINT32 n, + OPJ_UINT32 * pPermutations, + OPJ_FLOAT32 * p_src_temp, + OPJ_FLOAT32 * p_dest_temp, + OPJ_FLOAT32 * p_swap_area ) +{ + OPJ_UINT32 j,i; + OPJ_FLOAT32 * lCurrentPtr; + OPJ_FLOAT32 * lLineMatrix = pDestMatrix; + OPJ_UINT32 lSwapSize = n * sizeof(OPJ_FLOAT32); + + for (j = 0; j < n; ++j) { + lCurrentPtr = lLineMatrix++; + memset(p_src_temp,0,lSwapSize); + p_src_temp[j] = 1.0; + opj_lupSolve(p_dest_temp,pSrcMatrix,p_src_temp, pPermutations, n , p_swap_area); + + for (i = 0; i < n; ++i) { + *(lCurrentPtr) = p_dest_temp[i]; + lCurrentPtr+=n; + } + } +} + diff --git a/src/lib/openjp2/invert.h b/src/lib/openjp2/invert.h new file mode 100644 index 00000000..669a68e3 --- /dev/null +++ b/src/lib/openjp2/invert.h @@ -0,0 +1,58 @@ +/* + * Copyright (c) 2008, Jerome Fimes, Communications & Systemes + * 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. + */ + +#ifndef __INVERT_H +#define __INVERT_H +/** +@file invert.h +@brief Implementation of the matrix inversion + +The function in INVERT.H compute a matrix inversion with a LUP method +*/ + +/** @defgroup INVERT INVERT - Implementation of a matrix inversion */ +/*@{*/ +/** @name Exported functions */ +/*@{*/ +/* ----------------------------------------------------------------------- */ + +/** + * Calculates a n x n double matrix inversion with a LUP method. Data is aligned, rows after rows (or columns after columns). + * The function does not take ownership of any memory block, data must be fred by the user. + * + * @param pSrcMatrix the matrix to invert. + * @param pDestMatrix data to store the inverted matrix. + * @return OPJ_TRUE if the inversion is successful, OPJ_FALSE if the matrix is singular. + */ +opj_bool opj_matrix_inversion_f(OPJ_FLOAT32 * pSrcMatrix, + OPJ_FLOAT32 * pDestMatrix, + OPJ_UINT32 n); +/* ----------------------------------------------------------------------- */ +/*@}*/ + +/*@}*/ + +#endif /* __INVERT_H */ \ No newline at end of file diff --git a/src/lib/openjp2/j2k.c b/src/lib/openjp2/j2k.c index 5481a6c6..82bc0627 100644 --- a/src/lib/openjp2/j2k.c +++ b/src/lib/openjp2/j2k.c @@ -5909,36 +5909,31 @@ void opj_j2k_setup_encoder( opj_j2k_v2_t *p_j2k, tcp->tccps = (opj_tccp_t*) opj_calloc(image->numcomps, sizeof(opj_tccp_t)); if (parameters->mct_data) { + + OPJ_UINT32 lMctSize = image->numcomps * image->numcomps * sizeof(OPJ_FLOAT32); + OPJ_FLOAT32 * lTmpBuf = (OPJ_FLOAT32*)opj_malloc(lMctSize); + OPJ_INT32 * l_dc_shift = (OPJ_INT32 *) ((OPJ_BYTE *) parameters->mct_data + lMctSize); - opj_event_msg_v2(p_manager, EVT_ERROR, "MCT not supported for now\n"); - return; + tcp->mct = 2; + tcp->m_mct_coding_matrix = (OPJ_FLOAT32*)opj_malloc(lMctSize); + memcpy(tcp->m_mct_coding_matrix,parameters->mct_data,lMctSize); + memcpy(lTmpBuf,parameters->mct_data,lMctSize); - /* TODO MSD : merge v2 add invert.c or used a external lib ? - OPJ_UINT32 lMctSize = image->numcomps * image->numcomps * sizeof(OPJ_FLOAT32); - OPJ_FLOAT32 * lTmpBuf = (OPJ_FLOAT32*)opj_malloc(lMctSize); - OPJ_INT32 * l_dc_shift = (OPJ_INT32 *) ((OPJ_BYTE *) parameters->mct_data + lMctSize); + tcp->m_mct_decoding_matrix = (OPJ_FLOAT32*)opj_malloc(lMctSize); + assert(opj_matrix_inversion_f(lTmpBuf,(tcp->m_mct_decoding_matrix),image->numcomps)); - tcp->mct = 2; - tcp->m_mct_coding_matrix = (OPJ_FLOAT32*)opj_malloc(lMctSize); - memcpy(tcp->m_mct_coding_matrix,parameters->mct_data,lMctSize); - memcpy(lTmpBuf,parameters->mct_data,lMctSize); + tcp->mct_norms = (OPJ_FLOAT64*) + opj_malloc(image->numcomps * sizeof(OPJ_FLOAT64)); - tcp->m_mct_decoding_matrix = (OPJ_FLOAT32*)opj_malloc(lMctSize); - assert(opj_matrix_inversion_f(lTmpBuf,(tcp->m_mct_decoding_matrix),image->numcomps)); + opj_calculate_norms(tcp->mct_norms,image->numcomps,tcp->m_mct_decoding_matrix); + opj_free(lTmpBuf); - tcp->mct_norms = (OPJ_FLOAT64*) - opj_malloc(image->numcomps * sizeof(OPJ_FLOAT64)); - - opj_calculate_norms(tcp->mct_norms,image->numcomps,tcp->m_mct_decoding_matrix); - opj_free(lTmpBuf); - - for (i = 0; i < image->numcomps; i++) { - opj_tccp_t *tccp = &tcp->tccps[i]; - tccp->m_dc_level_shift = l_dc_shift[i]; - } + for (i = 0; i < image->numcomps; i++) { + opj_tccp_t *tccp = &tcp->tccps[i]; + tccp->m_dc_level_shift = l_dc_shift[i]; + } - opj_j2k_setup_mct_encoding(tcp,image); - */ + opj_j2k_setup_mct_encoding(tcp,image); } else { for (i = 0; i < image->numcomps; i++) { diff --git a/src/lib/openjp2/opj_includes.h b/src/lib/openjp2/opj_includes.h index 427bf8b4..c31ae4cf 100644 --- a/src/lib/openjp2/opj_includes.h +++ b/src/lib/openjp2/opj_includes.h @@ -153,6 +153,7 @@ static INLINE long lrintf(float f){ #include "cio.h" #include "image.h" +#include "invert.h" #include "j2k.h" #include "jp2.h" #include "jpt.h"