cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 1 | /* |
| 2 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 3 | % % |
| 4 | % % |
| 5 | % % |
| 6 | % M M AAA TTTTT RRRR IIIII X X % |
| 7 | % MM MM A A T R R I X X % |
| 8 | % M M M AAAAA T RRRR I X % |
| 9 | % M M A A T R R I X X % |
| 10 | % M M A A T R R IIIII X X % |
| 11 | % % |
| 12 | % % |
| 13 | % MagickCore Matrix Methods % |
| 14 | % % |
| 15 | % Software Design % |
cristy | de984cd | 2013-12-01 14:49:27 +0000 | [diff] [blame^] | 16 | % Cristy % |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 17 | % August 2007 % |
| 18 | % % |
| 19 | % % |
cristy | fe676ee | 2013-11-18 13:03:38 +0000 | [diff] [blame] | 20 | % Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization % |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 21 | % dedicated to making software imaging solutions freely available. % |
| 22 | % % |
| 23 | % You may not use this file except in compliance with the License. You may % |
| 24 | % obtain a copy of the License at % |
| 25 | % % |
| 26 | % http://www.imagemagick.org/script/license.php % |
| 27 | % % |
| 28 | % Unless required by applicable law or agreed to in writing, software % |
| 29 | % distributed under the License is distributed on an "AS IS" BASIS, % |
| 30 | % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. % |
| 31 | % See the License for the specific language governing permissions and % |
| 32 | % limitations under the License. % |
| 33 | % % |
| 34 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 35 | % |
| 36 | % |
| 37 | */ |
| 38 | |
| 39 | /* |
| 40 | Include declarations. |
| 41 | */ |
cristy | 4c08aed | 2011-07-01 19:47:50 +0000 | [diff] [blame] | 42 | #include "MagickCore/studio.h" |
| 43 | #include "MagickCore/matrix.h" |
cristy | d1dd6e4 | 2011-09-04 01:46:08 +0000 | [diff] [blame] | 44 | #include "MagickCore/matrix-private.h" |
cristy | 35f1530 | 2012-06-07 14:59:02 +0000 | [diff] [blame] | 45 | #include "MagickCore/pixel-private.h" |
cristy | 4c08aed | 2011-07-01 19:47:50 +0000 | [diff] [blame] | 46 | #include "MagickCore/memory_.h" |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 47 | |
| 48 | /* |
| 49 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 50 | % % |
| 51 | % % |
| 52 | % % |
| 53 | % A c q u i r e M a g i c k M a t r i x % |
| 54 | % % |
| 55 | % % |
| 56 | % % |
| 57 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 58 | % |
| 59 | % AcquireMagickMatrix() allocates and returns a matrix in the form of an |
| 60 | % array of pointers to an array of doubles, with all values pre-set to zero. |
| 61 | % |
anthony | fd706f9 | 2012-01-19 04:22:02 +0000 | [diff] [blame] | 62 | % This used to generate the two dimensional matrix, that can be referenced |
| 63 | % using the simple C-code of the form "matrix[y][x]". |
| 64 | % |
| 65 | % This matrix is typically used for perform for the GaussJordanElimination() |
| 66 | % method below, solving some system of simultanious equations. |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 67 | % |
| 68 | % The format of the AcquireMagickMatrix method is: |
| 69 | % |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 70 | % double **AcquireMagickMatrix(const size_t number_rows, |
| 71 | % const size_t size) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 72 | % |
| 73 | % A description of each parameter follows: |
| 74 | % |
cristy | 74908cf | 2010-05-17 19:43:54 +0000 | [diff] [blame] | 75 | % o number_rows: the number pointers for the array of pointers |
cristy | 1ad491d | 2010-05-17 19:45:27 +0000 | [diff] [blame] | 76 | % (first dimension). |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 77 | % |
cristy | 1ad491d | 2010-05-17 19:45:27 +0000 | [diff] [blame] | 78 | % o size: the size of the array of doubles each pointer points to |
| 79 | % (second dimension). |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 80 | % |
| 81 | */ |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 82 | MagickExport double **AcquireMagickMatrix(const size_t number_rows, |
| 83 | const size_t size) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 84 | { |
| 85 | double |
cristy | 74908cf | 2010-05-17 19:43:54 +0000 | [diff] [blame] | 86 | **matrix; |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 87 | |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 88 | register ssize_t |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 89 | i, |
| 90 | j; |
| 91 | |
cristy | 74908cf | 2010-05-17 19:43:54 +0000 | [diff] [blame] | 92 | matrix=(double **) AcquireQuantumMemory(number_rows,sizeof(*matrix)); |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 93 | if (matrix == (double **) NULL) |
cristy | 2629532 | 2011-09-22 00:50:36 +0000 | [diff] [blame] | 94 | return((double **) NULL); |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 95 | for (i=0; i < (ssize_t) number_rows; i++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 96 | { |
| 97 | matrix[i]=(double *) AcquireQuantumMemory(size,sizeof(*matrix[i])); |
| 98 | if (matrix[i] == (double *) NULL) |
| 99 | { |
| 100 | for (j=0; j < i; j++) |
| 101 | matrix[j]=(double *) RelinquishMagickMemory(matrix[j]); |
| 102 | matrix=(double **) RelinquishMagickMemory(matrix); |
| 103 | return((double **) NULL); |
| 104 | } |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 105 | for (j=0; j < (ssize_t) size; j++) |
cristy | 74908cf | 2010-05-17 19:43:54 +0000 | [diff] [blame] | 106 | matrix[i][j]=0.0; |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 107 | } |
| 108 | return(matrix); |
| 109 | } |
| 110 | |
| 111 | /* |
| 112 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 113 | % % |
| 114 | % % |
| 115 | % % |
| 116 | % G a u s s J o r d a n E l i m i n a t i o n % |
| 117 | % % |
| 118 | % % |
| 119 | % % |
| 120 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 121 | % |
| 122 | % GaussJordanElimination() returns a matrix in reduced row echelon form, |
| 123 | % while simultaneously reducing and thus solving the augumented results |
| 124 | % matrix. |
| 125 | % |
| 126 | % See also http://en.wikipedia.org/wiki/Gauss-Jordan_elimination |
| 127 | % |
| 128 | % The format of the GaussJordanElimination method is: |
| 129 | % |
cristy | 2629532 | 2011-09-22 00:50:36 +0000 | [diff] [blame] | 130 | % MagickBooleanType GaussJordanElimination(double **matrix, |
| 131 | % double **vectors,const size_t rank,const size_t number_vectors) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 132 | % |
| 133 | % A description of each parameter follows: |
| 134 | % |
| 135 | % o matrix: the matrix to be reduced, as an 'array of row pointers'. |
| 136 | % |
| 137 | % o vectors: the additional matrix argumenting the matrix for row reduction. |
| 138 | % Producing an 'array of column vectors'. |
| 139 | % |
anthony | fd706f9 | 2012-01-19 04:22:02 +0000 | [diff] [blame] | 140 | % o rank: The size of the square matrix (both rows and columns). |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 141 | % Also represents the number terms that need to be solved. |
| 142 | % |
cristy | 74908cf | 2010-05-17 19:43:54 +0000 | [diff] [blame] | 143 | % o number_vectors: Number of vectors columns, argumenting the above matrix. |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 144 | % Usally 1, but can be more for more complex equation solving. |
| 145 | % |
| 146 | % Note that the 'matrix' is given as a 'array of row pointers' of rank size. |
| 147 | % That is values can be assigned as matrix[row][column] where 'row' is |
| 148 | % typically the equation, and 'column' is the term of the equation. |
| 149 | % That is the matrix is in the form of a 'row first array'. |
| 150 | % |
| 151 | % However 'vectors' is a 'array of column pointers' which can have any number |
| 152 | % of columns, with each column array the same 'rank' size as 'matrix'. |
anthony | fd706f9 | 2012-01-19 04:22:02 +0000 | [diff] [blame] | 153 | % It is assigned vector[column][row] where 'column' is the specific |
| 154 | % 'result' and 'row' is the 'values' for that answer. After processing |
| 155 | % the same vector array contains the 'weights' (answers) for each of the |
| 156 | % 'separatable' results. |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 157 | % |
| 158 | % This allows for simpler handling of the results, especially is only one |
anthony | fd706f9 | 2012-01-19 04:22:02 +0000 | [diff] [blame] | 159 | % column 'vector' is all that is required to produce the desired solution |
| 160 | % for that specific set of equations. |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 161 | % |
| 162 | % For example, the 'vectors' can consist of a pointer to a simple array of |
| 163 | % doubles. when only one set of simultanious equations is to be solved from |
| 164 | % the given set of coefficient weighted terms. |
| 165 | % |
| 166 | % double **matrix = AcquireMagickMatrix(8UL,8UL); |
| 167 | % double coefficents[8]; |
| 168 | % ... |
| 169 | % GaussJordanElimination(matrix, &coefficents, 8UL, 1UL); |
anthony | 34364f4 | 2011-10-21 05:31:53 +0000 | [diff] [blame] | 170 | % |
anthony | b16e143 | 2011-10-21 05:06:56 +0000 | [diff] [blame] | 171 | % However by specifing more 'columns' (as an 'array of vector columns'), |
| 172 | % you can use this function to solve multiple sets of 'separable' equations. |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 173 | % |
| 174 | % For example a distortion function where u = U(x,y) v = V(x,y) |
| 175 | % And the functions U() and V() have separate coefficents, but are being |
| 176 | % generated from a common x,y->u,v data set. |
| 177 | % |
| 178 | % Another example is generation of a color gradient from a set of colors |
| 179 | % at specific coordients, such as a list x,y -> r,g,b,a |
anthony | fd706f9 | 2012-01-19 04:22:02 +0000 | [diff] [blame] | 180 | % |
| 181 | % See LeastSquaresAddTerms() below for such an example. |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 182 | % |
| 183 | % You can also use the 'vectors' to generate an inverse of the given 'matrix' |
anthony | b16e143 | 2011-10-21 05:06:56 +0000 | [diff] [blame] | 184 | % though as a 'column first array' rather than a 'row first array' (matrix |
anthony | fd706f9 | 2012-01-19 04:22:02 +0000 | [diff] [blame] | 185 | % is transposed). |
| 186 | % |
| 187 | % For details of this process see... |
anthony | 34364f4 | 2011-10-21 05:31:53 +0000 | [diff] [blame] | 188 | % http://en.wikipedia.org/wiki/Gauss-Jordan_elimination |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 189 | % |
| 190 | */ |
cristy | d1dd6e4 | 2011-09-04 01:46:08 +0000 | [diff] [blame] | 191 | MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix, |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 192 | double **vectors,const size_t rank,const size_t number_vectors) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 193 | { |
| 194 | #define GaussJordanSwap(x,y) \ |
| 195 | { \ |
| 196 | if ((x) != (y)) \ |
| 197 | { \ |
| 198 | (x)+=(y); \ |
| 199 | (y)=(x)-(y); \ |
| 200 | (x)=(x)-(y); \ |
| 201 | } \ |
| 202 | } |
| 203 | |
| 204 | double |
| 205 | max, |
| 206 | scale; |
| 207 | |
cristy | 9d314ff | 2011-03-09 01:30:28 +0000 | [diff] [blame] | 208 | register ssize_t |
| 209 | i, |
| 210 | j, |
| 211 | k; |
| 212 | |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 213 | ssize_t |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 214 | column, |
| 215 | *columns, |
| 216 | *pivots, |
| 217 | row, |
| 218 | *rows; |
| 219 | |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 220 | columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns)); |
| 221 | rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows)); |
| 222 | pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots)); |
| 223 | if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) || |
| 224 | (pivots == (ssize_t *) NULL)) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 225 | { |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 226 | if (pivots != (ssize_t *) NULL) |
| 227 | pivots=(ssize_t *) RelinquishMagickMemory(pivots); |
| 228 | if (columns != (ssize_t *) NULL) |
| 229 | columns=(ssize_t *) RelinquishMagickMemory(columns); |
| 230 | if (rows != (ssize_t *) NULL) |
| 231 | rows=(ssize_t *) RelinquishMagickMemory(rows); |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 232 | return(MagickFalse); |
| 233 | } |
| 234 | (void) ResetMagickMemory(columns,0,rank*sizeof(*columns)); |
| 235 | (void) ResetMagickMemory(rows,0,rank*sizeof(*rows)); |
| 236 | (void) ResetMagickMemory(pivots,0,rank*sizeof(*pivots)); |
| 237 | column=0; |
| 238 | row=0; |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 239 | for (i=0; i < (ssize_t) rank; i++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 240 | { |
| 241 | max=0.0; |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 242 | for (j=0; j < (ssize_t) rank; j++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 243 | if (pivots[j] != 1) |
| 244 | { |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 245 | for (k=0; k < (ssize_t) rank; k++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 246 | if (pivots[k] != 0) |
| 247 | { |
| 248 | if (pivots[k] > 1) |
| 249 | return(MagickFalse); |
| 250 | } |
| 251 | else |
| 252 | if (fabs(matrix[j][k]) >= max) |
| 253 | { |
| 254 | max=fabs(matrix[j][k]); |
| 255 | row=j; |
| 256 | column=k; |
| 257 | } |
| 258 | } |
| 259 | pivots[column]++; |
| 260 | if (row != column) |
| 261 | { |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 262 | for (k=0; k < (ssize_t) rank; k++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 263 | GaussJordanSwap(matrix[row][k],matrix[column][k]); |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 264 | for (k=0; k < (ssize_t) number_vectors; k++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 265 | GaussJordanSwap(vectors[k][row],vectors[k][column]); |
| 266 | } |
| 267 | rows[i]=row; |
| 268 | columns[i]=column; |
| 269 | if (matrix[column][column] == 0.0) |
anthony | 34364f4 | 2011-10-21 05:31:53 +0000 | [diff] [blame] | 270 | return(MagickFalse); /* singularity */ |
cristy | 3e3ec3a | 2012-11-03 23:11:06 +0000 | [diff] [blame] | 271 | scale=PerceptibleReciprocal(matrix[column][column]); |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 272 | matrix[column][column]=1.0; |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 273 | for (j=0; j < (ssize_t) rank; j++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 274 | matrix[column][j]*=scale; |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 275 | for (j=0; j < (ssize_t) number_vectors; j++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 276 | vectors[j][column]*=scale; |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 277 | for (j=0; j < (ssize_t) rank; j++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 278 | if (j != column) |
| 279 | { |
| 280 | scale=matrix[j][column]; |
| 281 | matrix[j][column]=0.0; |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 282 | for (k=0; k < (ssize_t) rank; k++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 283 | matrix[j][k]-=scale*matrix[column][k]; |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 284 | for (k=0; k < (ssize_t) number_vectors; k++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 285 | vectors[k][j]-=scale*vectors[k][column]; |
| 286 | } |
| 287 | } |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 288 | for (j=(ssize_t) rank-1; j >= 0; j--) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 289 | if (columns[j] != rows[j]) |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 290 | for (i=0; i < (ssize_t) rank; i++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 291 | GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]); |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 292 | pivots=(ssize_t *) RelinquishMagickMemory(pivots); |
| 293 | rows=(ssize_t *) RelinquishMagickMemory(rows); |
| 294 | columns=(ssize_t *) RelinquishMagickMemory(columns); |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 295 | return(MagickTrue); |
| 296 | } |
| 297 | |
| 298 | /* |
| 299 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 300 | % % |
| 301 | % % |
| 302 | % % |
| 303 | % L e a s t S q u a r e s A d d T e r m s % |
| 304 | % % |
| 305 | % % |
| 306 | % % |
| 307 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 308 | % |
| 309 | % LeastSquaresAddTerms() adds one set of terms and associate results to the |
| 310 | % given matrix and vectors for solving using least-squares function fitting. |
| 311 | % |
| 312 | % The format of the AcquireMagickMatrix method is: |
| 313 | % |
| 314 | % void LeastSquaresAddTerms(double **matrix,double **vectors, |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 315 | % const double *terms,const double *results,const size_t rank, |
| 316 | % const size_t number_vectors); |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 317 | % |
| 318 | % A description of each parameter follows: |
| 319 | % |
| 320 | % o matrix: the square matrix to add given terms/results to. |
| 321 | % |
| 322 | % o vectors: the result vectors to add terms/results to. |
| 323 | % |
| 324 | % o terms: the pre-calculated terms (without the unknown coefficent |
cristy | 2629532 | 2011-09-22 00:50:36 +0000 | [diff] [blame] | 325 | % weights) that forms the equation being added. |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 326 | % |
| 327 | % o results: the result(s) that should be generated from the given terms |
cristy | 2629532 | 2011-09-22 00:50:36 +0000 | [diff] [blame] | 328 | % weighted by the yet-to-be-solved coefficents. |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 329 | % |
cristy | 4c08aed | 2011-07-01 19:47:50 +0000 | [diff] [blame] | 330 | % o rank: the rank or size of the dimentions of the square matrix. |
cristy | 2629532 | 2011-09-22 00:50:36 +0000 | [diff] [blame] | 331 | % Also the length of vectors, and number of terms being added. |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 332 | % |
cristy | 1ad491d | 2010-05-17 19:45:27 +0000 | [diff] [blame] | 333 | % o number_vectors: Number of result vectors, and number or results being |
| 334 | % added. Also represents the number of separable systems of equations |
| 335 | % that is being solved. |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 336 | % |
| 337 | % Example of use... |
| 338 | % |
glennrp | 2489f53 | 2011-06-25 03:02:43 +0000 | [diff] [blame] | 339 | % 2 dimensional Affine Equations (which are separable) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 340 | % c0*x + c2*y + c4*1 => u |
| 341 | % c1*x + c3*y + c5*1 => v |
| 342 | % |
| 343 | % double **matrix = AcquireMagickMatrix(3UL,3UL); |
| 344 | % double **vectors = AcquireMagickMatrix(2UL,3UL); |
| 345 | % double terms[3], results[2]; |
| 346 | % ... |
| 347 | % for each given x,y -> u,v |
| 348 | % terms[0] = x; |
| 349 | % terms[1] = y; |
| 350 | % terms[2] = 1; |
| 351 | % results[0] = u; |
| 352 | % results[1] = v; |
| 353 | % LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL); |
| 354 | % ... |
| 355 | % if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) { |
| 356 | % c0 = vectors[0][0]; |
anthony | fd706f9 | 2012-01-19 04:22:02 +0000 | [diff] [blame] | 357 | % c2 = vectors[0][1]; %* weights to calculate u from any given x,y *% |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 358 | % c4 = vectors[0][2]; |
| 359 | % c1 = vectors[1][0]; |
anthony | fd706f9 | 2012-01-19 04:22:02 +0000 | [diff] [blame] | 360 | % c3 = vectors[1][1]; %* weights for calculate v from any given x,y *% |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 361 | % c5 = vectors[1][2]; |
| 362 | % } |
| 363 | % else |
| 364 | % printf("Matrix unsolvable\n); |
| 365 | % RelinquishMagickMatrix(matrix,3UL); |
| 366 | % RelinquishMagickMatrix(vectors,2UL); |
| 367 | % |
anthony | fd706f9 | 2012-01-19 04:22:02 +0000 | [diff] [blame] | 368 | % More examples can be found in "distort.c" |
| 369 | % |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 370 | */ |
cristy | d1dd6e4 | 2011-09-04 01:46:08 +0000 | [diff] [blame] | 371 | MagickPrivate void LeastSquaresAddTerms(double **matrix,double **vectors, |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 372 | const double *terms,const double *results,const size_t rank, |
| 373 | const size_t number_vectors) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 374 | { |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 375 | register ssize_t |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 376 | i, |
| 377 | j; |
| 378 | |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 379 | for (j=0; j < (ssize_t) rank; j++) |
cristy | 74908cf | 2010-05-17 19:43:54 +0000 | [diff] [blame] | 380 | { |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 381 | for (i=0; i < (ssize_t) rank; i++) |
cristy | 74908cf | 2010-05-17 19:43:54 +0000 | [diff] [blame] | 382 | matrix[i][j]+=terms[i]*terms[j]; |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 383 | for (i=0; i < (ssize_t) number_vectors; i++) |
cristy | 74908cf | 2010-05-17 19:43:54 +0000 | [diff] [blame] | 384 | vectors[i][j]+=results[i]*terms[j]; |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 385 | } |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 386 | } |
| 387 | |
| 388 | /* |
| 389 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 390 | % % |
| 391 | % % |
| 392 | % % |
| 393 | % R e l i n q u i s h M a g i c k M a t r i x % |
| 394 | % % |
| 395 | % % |
| 396 | % % |
| 397 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 398 | % |
| 399 | % RelinquishMagickMatrix() frees the previously acquired matrix (array of |
| 400 | % pointers to arrays of doubles). |
| 401 | % |
| 402 | % The format of the RelinquishMagickMatrix method is: |
| 403 | % |
| 404 | % double **RelinquishMagickMatrix(double **matrix, |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 405 | % const size_t number_rows) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 406 | % |
| 407 | % A description of each parameter follows: |
| 408 | % |
| 409 | % o matrix: the matrix to relinquish |
| 410 | % |
cristy | 1ad491d | 2010-05-17 19:45:27 +0000 | [diff] [blame] | 411 | % o number_rows: the first dimension of the acquired matrix (number of |
| 412 | % pointers) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 413 | % |
| 414 | */ |
| 415 | MagickExport double **RelinquishMagickMatrix(double **matrix, |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 416 | const size_t number_rows) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 417 | { |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 418 | register ssize_t |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 419 | i; |
| 420 | |
| 421 | if (matrix == (double **) NULL ) |
| 422 | return(matrix); |
cristy | bb50337 | 2010-05-27 20:51:26 +0000 | [diff] [blame] | 423 | for (i=0; i < (ssize_t) number_rows; i++) |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 424 | matrix[i]=(double *) RelinquishMagickMemory(matrix[i]); |
| 425 | matrix=(double **) RelinquishMagickMemory(matrix); |
cristy | 3ed852e | 2009-09-05 21:47:34 +0000 | [diff] [blame] | 426 | return(matrix); |
| 427 | } |