2 * Copyright (c) 2008, Jerome Fimes, Communications & Systemes <jerome.fimes@c-s.fr>
\r
3 * All rights reserved.
\r
5 * Redistribution and use in source and binary forms, with or without
\r
6 * modification, are permitted provided that the following conditions
\r
8 * 1. Redistributions of source code must retain the above copyright
\r
9 * notice, this list of conditions and the following disclaimer.
\r
10 * 2. Redistributions in binary form must reproduce the above copyright
\r
11 * notice, this list of conditions and the following disclaimer in the
\r
12 * documentation and/or other materials provided with the distribution.
\r
14 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
\r
15 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
\r
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
\r
17 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
\r
18 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
\r
19 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
\r
20 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
\r
21 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
\r
22 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
\r
23 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
\r
24 * POSSIBILITY OF SUCH DAMAGE.
\r
28 #include "opj_malloc.h"
\r
31 bool opj_lupDecompose(OPJ_FLOAT32 * matrix,OPJ_UINT32 * permutations, OPJ_FLOAT32 * p_swap_area,OPJ_UINT32 n);
\r
32 void opj_lupSolve(OPJ_FLOAT32 * pResult, OPJ_FLOAT32* pMatrix, OPJ_FLOAT32* pVector, OPJ_UINT32* pPermutations, OPJ_UINT32 n,OPJ_FLOAT32 * p_intermediate_data);
\r
33 void opj_lupInvert (OPJ_FLOAT32 * pSrcMatrix,
\r
34 OPJ_FLOAT32 * pDestMatrix,
\r
36 OPJ_UINT32 * pPermutations,
\r
37 OPJ_FLOAT32 * p_src_temp,
\r
38 OPJ_FLOAT32 * p_dest_temp,
\r
39 OPJ_FLOAT32 * p_swap_area);
\r
44 bool opj_matrix_inversion_f(OPJ_FLOAT32 * pSrcMatrix,OPJ_FLOAT32 * pDestMatrix, OPJ_UINT32 n)
\r
46 OPJ_BYTE * l_data = 00;
\r
47 OPJ_UINT32 l_permutation_size = n * sizeof(OPJ_UINT32);
\r
48 OPJ_UINT32 l_swap_size = n * sizeof(OPJ_FLOAT32);
\r
49 OPJ_UINT32 l_total_size = l_permutation_size + 3 * l_swap_size;
\r
50 OPJ_UINT32 * lPermutations = 00;
\r
51 OPJ_FLOAT32 * l_double_data = 00;
\r
53 l_data = (OPJ_BYTE *) opj_malloc(l_total_size);
\r
59 lPermutations = (OPJ_UINT32 *) l_data;
\r
60 l_double_data = (OPJ_FLOAT32 *) (l_data + l_permutation_size);
\r
61 memset(lPermutations,0,l_permutation_size);
\r
64 (! opj_lupDecompose(pSrcMatrix,lPermutations,l_double_data,n))
\r
69 opj_lupInvert(pSrcMatrix,pDestMatrix,n,lPermutations,l_double_data,l_double_data + n,l_double_data + 2*n);
\r
78 bool opj_lupDecompose(OPJ_FLOAT32 * matrix,OPJ_UINT32 * permutations, OPJ_FLOAT32 * p_swap_area,OPJ_UINT32 n)
\r
80 OPJ_UINT32 * tmpPermutations = permutations;
\r
81 OPJ_UINT32 * dstPermutations;
\r
86 OPJ_UINT32 lLastColum = n - 1;
\r
87 OPJ_UINT32 lSwapSize = n * sizeof(OPJ_FLOAT32);
\r
88 OPJ_FLOAT32 * lTmpMatrix = matrix;
\r
89 OPJ_FLOAT32 * lColumnMatrix,* lDestMatrix;
\r
90 OPJ_UINT32 offset = 1;
\r
91 OPJ_UINT32 lStride = n-1;
\r
93 //initialize permutations
\r
95 (i = 0; i < n; ++i)
\r
97 *tmpPermutations++ = i;
\r
102 // now make a pivot with colum switch
\r
103 tmpPermutations = permutations;
\r
105 (k = 0; k < lLastColum; ++k)
\r
109 // take the middle element
\r
110 lColumnMatrix = lTmpMatrix + k;
\r
112 // make permutation with the biggest value in the column
\r
114 (i = k; i < n; ++i)
\r
116 temp = ((*lColumnMatrix > 0) ? *lColumnMatrix : -(*lColumnMatrix));
\r
124 lColumnMatrix += n;
\r
127 // a whole rest of 0 -> non singular
\r
134 // should we permute ?
\r
140 dstPermutations = tmpPermutations + k2 - k;
\r
142 t = *tmpPermutations;
\r
143 *tmpPermutations = *dstPermutations;
\r
144 *dstPermutations = t;
\r
146 // and swap entire line.
\r
147 lColumnMatrix = lTmpMatrix + (k2 - k) * n;
\r
148 memcpy(p_swap_area,lColumnMatrix,lSwapSize);
\r
149 memcpy(lColumnMatrix,lTmpMatrix,lSwapSize);
\r
150 memcpy(lTmpMatrix,p_swap_area,lSwapSize);
\r
153 // now update data in the rest of the line and line after
\r
154 lDestMatrix = lTmpMatrix + k;
\r
155 lColumnMatrix = lDestMatrix + n;
\r
156 // take the middle element
\r
157 temp = *(lDestMatrix++);
\r
159 // now compute up data (i.e. coeff up of the diagonal).
\r
160 for (i = offset; i < n; ++i)
\r
163 // divide the lower column elements by the diagonal value
\r
165 // matrix[i][k] /= matrix[k][k];
\r
166 // p = matrix[i][k]
\r
167 p = *lColumnMatrix / temp;
\r
168 *(lColumnMatrix++) = p;
\r
170 (j = /* k + 1 */ offset; j < n; ++j)
\r
172 // matrix[i][j] -= matrix[i][k] * matrix[k][j];
\r
173 *(lColumnMatrix++) -= p * (*(lDestMatrix++));
\r
175 // come back to the k+1th element
\r
176 lDestMatrix -= lStride;
\r
177 // go to kth element of the next line
\r
178 lColumnMatrix += k;
\r
180 // offset is now k+2
\r
182 // 1 element less for stride
\r
186 // next permutation element
\r
197 void opj_lupSolve (OPJ_FLOAT32 * pResult, OPJ_FLOAT32 * pMatrix, OPJ_FLOAT32 * pVector, OPJ_UINT32* pPermutations, OPJ_UINT32 n,OPJ_FLOAT32 * p_intermediate_data)
\r
202 OPJ_UINT32 lStride = n+1;
\r
203 OPJ_FLOAT32 * lCurrentPtr;
\r
204 OPJ_FLOAT32 * lIntermediatePtr;
\r
205 OPJ_FLOAT32 * lDestPtr;
\r
206 OPJ_FLOAT32 * lTmpMatrix;
\r
207 OPJ_FLOAT32 * lLineMatrix = pMatrix;
\r
208 OPJ_FLOAT32 * lBeginPtr = pResult + n - 1;
\r
209 OPJ_FLOAT32 * lGeneratedData;
\r
210 OPJ_UINT32 * lCurrentPermutationPtr = pPermutations;
\r
213 lIntermediatePtr = p_intermediate_data;
\r
214 lGeneratedData = p_intermediate_data + n - 1;
\r
217 (i = 0; i < n; ++i)
\r
220 lCurrentPtr = p_intermediate_data;
\r
221 lTmpMatrix = lLineMatrix;
\r
223 (j = 1; j <= i; ++j)
\r
225 // sum += matrix[i][j-1] * y[j-1];
\r
226 sum += (*(lTmpMatrix++)) * (*(lCurrentPtr++));
\r
228 //y[i] = pVector[pPermutations[i]] - sum;
\r
229 *(lIntermediatePtr++) = pVector[*(lCurrentPermutationPtr++)] - sum;
\r
233 // we take the last point of the matrix
\r
234 lLineMatrix = pMatrix + n*n - 1;
\r
236 // and we take after the last point of the destination vector
\r
237 lDestPtr = pResult + n;
\r
240 (i = n - 1; i != -1 ; --i)
\r
243 lTmpMatrix = lLineMatrix;
\r
244 u = *(lTmpMatrix++);
\r
245 lCurrentPtr = lDestPtr--;
\r
247 (j = i + 1; j < n; ++j)
\r
249 // sum += matrix[i][j] * x[j]
\r
250 sum += (*(lTmpMatrix++)) * (*(lCurrentPtr++));
\r
252 //x[i] = (y[i] - sum) / u;
\r
253 *(lBeginPtr--) = (*(lGeneratedData--) - sum) / u;
\r
254 lLineMatrix -= lStride;
\r
258 /** LUP inversion (call with the result of lupDecompose)
\r
260 void opj_lupInvert (
\r
261 OPJ_FLOAT32 * pSrcMatrix,
\r
262 OPJ_FLOAT32 * pDestMatrix,
\r
264 OPJ_UINT32 * pPermutations,
\r
265 OPJ_FLOAT32 * p_src_temp,
\r
266 OPJ_FLOAT32 * p_dest_temp,
\r
267 OPJ_FLOAT32 * p_swap_area
\r
271 OPJ_FLOAT32 * lCurrentPtr;
\r
272 OPJ_FLOAT32 * lLineMatrix = pDestMatrix;
\r
273 OPJ_UINT32 lSwapSize = n * sizeof(OPJ_FLOAT32);
\r
276 (j = 0; j < n; ++j)
\r
278 lCurrentPtr = lLineMatrix++;
\r
279 memset(p_src_temp,0,lSwapSize);
\r
280 p_src_temp[j] = 1.0;
\r
281 opj_lupSolve(p_dest_temp,pSrcMatrix,p_src_temp, pPermutations, n , p_swap_area);
\r
284 (i = 0; i < n; ++i)
\r
286 *(lCurrentPtr) = p_dest_temp[i];
\r