diff options
Diffstat (limited to 'src/lib/openjp2/invert.c')
| -rw-r--r-- | src/lib/openjp2/invert.c | 81 |
1 files changed, 42 insertions, 39 deletions
diff --git a/src/lib/openjp2/invert.c b/src/lib/openjp2/invert.c index 30651ba3..b05fabd4 100644 --- a/src/lib/openjp2/invert.c +++ b/src/lib/openjp2/invert.c @@ -32,7 +32,7 @@ static opj_bool opj_lupDecompose(OPJ_FLOAT32 * matrix, OPJ_UINT32 * permutations, OPJ_FLOAT32 * p_swap_area, - OPJ_UINT32 n); + OPJ_UINT32 nb_compo); /** * LUP solving */ @@ -40,7 +40,7 @@ static void opj_lupSolve(OPJ_FLOAT32 * pResult, OPJ_FLOAT32* pMatrix, OPJ_FLOAT32* pVector, OPJ_UINT32* pPermutations, - OPJ_UINT32 n, + OPJ_UINT32 nb_compo, OPJ_FLOAT32 * p_intermediate_data); /** @@ -48,7 +48,7 @@ static void opj_lupSolve(OPJ_FLOAT32 * pResult, */ static void opj_lupInvert ( OPJ_FLOAT32 * pSrcMatrix, OPJ_FLOAT32 * pDestMatrix, - OPJ_UINT32 n, + OPJ_UINT32 nb_compo, OPJ_UINT32 * pPermutations, OPJ_FLOAT32 * p_src_temp, OPJ_FLOAT32 * p_dest_temp, @@ -64,11 +64,11 @@ static void opj_lupInvert ( OPJ_FLOAT32 * pSrcMatrix, */ opj_bool opj_matrix_inversion_f(OPJ_FLOAT32 * pSrcMatrix, OPJ_FLOAT32 * pDestMatrix, - OPJ_UINT32 n) + OPJ_UINT32 nb_compo) { 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_permutation_size = nb_compo * sizeof(OPJ_UINT32); + OPJ_UINT32 l_swap_size = nb_compo * 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; @@ -81,12 +81,12 @@ opj_bool opj_matrix_inversion_f(OPJ_FLOAT32 * pSrcMatrix, 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)) { + if(! opj_lupDecompose(pSrcMatrix,lPermutations,l_double_data,nb_compo)) { 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_lupInvert(pSrcMatrix,pDestMatrix,nb_compo,lPermutations,l_double_data,l_double_data + nb_compo,l_double_data + 2*nb_compo); opj_free(l_data); return OPJ_TRUE; @@ -100,7 +100,7 @@ opj_bool opj_matrix_inversion_f(OPJ_FLOAT32 * pSrcMatrix, */ opj_bool opj_lupDecompose(OPJ_FLOAT32 * matrix,OPJ_UINT32 * permutations, OPJ_FLOAT32 * p_swap_area, - OPJ_UINT32 n) + OPJ_UINT32 nb_compo) { OPJ_UINT32 * tmpPermutations = permutations; OPJ_UINT32 * dstPermutations; @@ -108,15 +108,15 @@ opj_bool opj_lupDecompose(OPJ_FLOAT32 * matrix,OPJ_UINT32 * permutations, OPJ_FLOAT32 temp; OPJ_UINT32 i,j,k; OPJ_FLOAT32 p; - OPJ_UINT32 lLastColum = n - 1; - OPJ_UINT32 lSwapSize = n * sizeof(OPJ_FLOAT32); + OPJ_UINT32 lLastColum = nb_compo - 1; + OPJ_UINT32 lSwapSize = nb_compo * sizeof(OPJ_FLOAT32); OPJ_FLOAT32 * lTmpMatrix = matrix; OPJ_FLOAT32 * lColumnMatrix,* lDestMatrix; OPJ_UINT32 offset = 1; - OPJ_UINT32 lStride = n-1; + OPJ_UINT32 lStride = nb_compo-1; /*initialize permutations */ - for (i = 0; i < n; ++i) + for (i = 0; i < nb_compo; ++i) { *tmpPermutations++ = i; } @@ -129,14 +129,14 @@ opj_bool opj_lupDecompose(OPJ_FLOAT32 * matrix,OPJ_UINT32 * permutations, lColumnMatrix = lTmpMatrix + k; /* make permutation with the biggest value in the column */ - for (i = k; i < n; ++i) { + for (i = k; i < nb_compo; ++i) { temp = ((*lColumnMatrix > 0) ? *lColumnMatrix : -(*lColumnMatrix)); if (temp > p) { p = temp; k2 = i; } /* next line */ - lColumnMatrix += n; + lColumnMatrix += nb_compo; } /* a whole rest of 0 -> non singular */ @@ -155,7 +155,7 @@ opj_bool opj_lupDecompose(OPJ_FLOAT32 * matrix,OPJ_UINT32 * permutations, *dstPermutations = t; /* and swap entire line. */ - lColumnMatrix = lTmpMatrix + (k2 - k) * n; + lColumnMatrix = lTmpMatrix + (k2 - k) * nb_compo; memcpy(p_swap_area,lColumnMatrix,lSwapSize); memcpy(lColumnMatrix,lTmpMatrix,lSwapSize); memcpy(lTmpMatrix,p_swap_area,lSwapSize); @@ -163,12 +163,12 @@ opj_bool opj_lupDecompose(OPJ_FLOAT32 * matrix,OPJ_UINT32 * permutations, /* now update data in the rest of the line and line after */ lDestMatrix = lTmpMatrix + k; - lColumnMatrix = lDestMatrix + n; + lColumnMatrix = lDestMatrix + nb_compo; /* take the middle element */ temp = *(lDestMatrix++); /* now compute up data (i.e. coeff up of the diagonal). */ - for (i = offset; i < n; ++i) { + for (i = offset; i < nb_compo; ++i) { /*lColumnMatrix; */ /* divide the lower column elements by the diagonal value */ @@ -177,7 +177,7 @@ opj_bool opj_lupDecompose(OPJ_FLOAT32 * matrix,OPJ_UINT32 * permutations, p = *lColumnMatrix / temp; *(lColumnMatrix++) = p; - for (j = /* k + 1 */ offset; j < n; ++j) { + for (j = /* k + 1 */ offset; j < nb_compo; ++j) { /* matrix[i][j] -= matrix[i][k] * matrix[k][j]; */ *(lColumnMatrix++) -= p * (*(lDestMatrix++)); } @@ -192,7 +192,7 @@ opj_bool opj_lupDecompose(OPJ_FLOAT32 * matrix,OPJ_UINT32 * permutations, /* 1 element less for stride */ --lStride; /* next line */ - lTmpMatrix+=n; + lTmpMatrix+=nb_compo; /* next permutation element */ ++tmpPermutations; } @@ -203,26 +203,27 @@ 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 nb_compo,OPJ_FLOAT32 * p_intermediate_data) { - OPJ_UINT32 i,j; + OPJ_INT32 k; + OPJ_UINT32 i,j; OPJ_FLOAT32 sum; OPJ_FLOAT32 u; - OPJ_UINT32 lStride = n+1; + OPJ_UINT32 lStride = nb_compo+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 * lBeginPtr = pResult + nb_compo - 1; OPJ_FLOAT32 * lGeneratedData; OPJ_UINT32 * lCurrentPermutationPtr = pPermutations; lIntermediatePtr = p_intermediate_data; - lGeneratedData = p_intermediate_data + n - 1; + lGeneratedData = p_intermediate_data + nb_compo - 1; - for (i = 0; i < n; ++i) { + for (i = 0; i < nb_compo; ++i) { sum = 0.0; lCurrentPtr = p_intermediate_data; lTmpMatrix = lLineMatrix; @@ -233,25 +234,27 @@ void opj_lupSolve (OPJ_FLOAT32 * pResult, } /*y[i] = pVector[pPermutations[i]] - sum; */ *(lIntermediatePtr++) = pVector[*(lCurrentPermutationPtr++)] - sum; - lLineMatrix += n; + lLineMatrix += nb_compo; } /* we take the last point of the matrix */ - lLineMatrix = pMatrix + n*n - 1; + lLineMatrix = pMatrix + nb_compo*nb_compo - 1; /* and we take after the last point of the destination vector */ - lDestPtr = pResult + n; + lDestPtr = pResult + nb_compo; - for (i = n - 1; i != -1 ; --i) { + + assert(nb_compo != 0); + for (k = (OPJ_INT32)nb_compo - 1; k != -1 ; --k) { sum = 0.0; lTmpMatrix = lLineMatrix; u = *(lTmpMatrix++); lCurrentPtr = lDestPtr--; - for (j = i + 1; j < n; ++j) { - /* sum += matrix[i][j] * x[j] */ + for (j = k + 1; j < nb_compo; ++j) { + /* sum += matrix[k][j] * x[j] */ sum += (*(lTmpMatrix++)) * (*(lCurrentPtr++)); } - /*x[i] = (y[i] - sum) / u; */ + /*x[k] = (y[k] - sum) / u; */ *(lBeginPtr--) = (*(lGeneratedData--) - sum) / u; lLineMatrix -= lStride; } @@ -260,7 +263,7 @@ void opj_lupSolve (OPJ_FLOAT32 * pResult, void opj_lupInvert (OPJ_FLOAT32 * pSrcMatrix, OPJ_FLOAT32 * pDestMatrix, - OPJ_UINT32 n, + OPJ_UINT32 nb_compo, OPJ_UINT32 * pPermutations, OPJ_FLOAT32 * p_src_temp, OPJ_FLOAT32 * p_dest_temp, @@ -269,17 +272,17 @@ void opj_lupInvert (OPJ_FLOAT32 * pSrcMatrix, OPJ_UINT32 j,i; OPJ_FLOAT32 * lCurrentPtr; OPJ_FLOAT32 * lLineMatrix = pDestMatrix; - OPJ_UINT32 lSwapSize = n * sizeof(OPJ_FLOAT32); + OPJ_UINT32 lSwapSize = nb_compo * sizeof(OPJ_FLOAT32); - for (j = 0; j < n; ++j) { + for (j = 0; j < nb_compo; ++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); + opj_lupSolve(p_dest_temp,pSrcMatrix,p_src_temp, pPermutations, nb_compo , p_swap_area); - for (i = 0; i < n; ++i) { + for (i = 0; i < nb_compo; ++i) { *(lCurrentPtr) = p_dest_temp[i]; - lCurrentPtr+=n; + lCurrentPtr+=nb_compo; } } } |
