2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 %
13 % MagickCore Matrix Methods %
20 % Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
26 % http://www.imagemagick.org/script/license.php %
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. %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42 #include "MagickCore/studio.h"
43 #include "MagickCore/matrix.h"
44 #include "MagickCore/matrix-private.h"
45 #include "MagickCore/memory_.h"
48 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52 % A c q u i r e M a g i c k M a t r i x %
56 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58 % AcquireMagickMatrix() allocates and returns a matrix in the form of an
59 % array of pointers to an array of doubles, with all values pre-set to zero.
61 % This used to generate the two dimensional matrix, that can be referenced
62 % using the simple C-code of the form "matrix[y][x]".
64 % This matrix is typically used for perform for the GaussJordanElimination()
65 % method below, solving some system of simultanious equations.
67 % The format of the AcquireMagickMatrix method is:
69 % double **AcquireMagickMatrix(const size_t number_rows,
72 % A description of each parameter follows:
74 % o number_rows: the number pointers for the array of pointers
77 % o size: the size of the array of doubles each pointer points to
81 MagickExport double **AcquireMagickMatrix(const size_t number_rows,
91 matrix=(double **) AcquireQuantumMemory(number_rows,sizeof(*matrix));
92 if (matrix == (double **) NULL)
93 return((double **) NULL);
94 for (i=0; i < (ssize_t) number_rows; i++)
96 matrix[i]=(double *) AcquireQuantumMemory(size,sizeof(*matrix[i]));
97 if (matrix[i] == (double *) NULL)
100 matrix[j]=(double *) RelinquishMagickMemory(matrix[j]);
101 matrix=(double **) RelinquishMagickMemory(matrix);
102 return((double **) NULL);
104 for (j=0; j < (ssize_t) size; j++)
111 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
115 % G a u s s J o r d a n E l i m i n a t i o n %
119 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121 % GaussJordanElimination() returns a matrix in reduced row echelon form,
122 % while simultaneously reducing and thus solving the augumented results
125 % See also http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
127 % The format of the GaussJordanElimination method is:
129 % MagickBooleanType GaussJordanElimination(double **matrix,
130 % double **vectors,const size_t rank,const size_t number_vectors)
132 % A description of each parameter follows:
134 % o matrix: the matrix to be reduced, as an 'array of row pointers'.
136 % o vectors: the additional matrix argumenting the matrix for row reduction.
137 % Producing an 'array of column vectors'.
139 % o rank: The size of the square matrix (both rows and columns).
140 % Also represents the number terms that need to be solved.
142 % o number_vectors: Number of vectors columns, argumenting the above matrix.
143 % Usally 1, but can be more for more complex equation solving.
145 % Note that the 'matrix' is given as a 'array of row pointers' of rank size.
146 % That is values can be assigned as matrix[row][column] where 'row' is
147 % typically the equation, and 'column' is the term of the equation.
148 % That is the matrix is in the form of a 'row first array'.
150 % However 'vectors' is a 'array of column pointers' which can have any number
151 % of columns, with each column array the same 'rank' size as 'matrix'.
152 % It is assigned vector[column][row] where 'column' is the specific
153 % 'result' and 'row' is the 'values' for that answer. After processing
154 % the same vector array contains the 'weights' (answers) for each of the
155 % 'separatable' results.
157 % This allows for simpler handling of the results, especially is only one
158 % column 'vector' is all that is required to produce the desired solution
159 % for that specific set of equations.
161 % For example, the 'vectors' can consist of a pointer to a simple array of
162 % doubles. when only one set of simultanious equations is to be solved from
163 % the given set of coefficient weighted terms.
165 % double **matrix = AcquireMagickMatrix(8UL,8UL);
166 % double coefficents[8];
168 % GaussJordanElimination(matrix, &coefficents, 8UL, 1UL);
170 % However by specifing more 'columns' (as an 'array of vector columns'),
171 % you can use this function to solve multiple sets of 'separable' equations.
173 % For example a distortion function where u = U(x,y) v = V(x,y)
174 % And the functions U() and V() have separate coefficents, but are being
175 % generated from a common x,y->u,v data set.
177 % Another example is generation of a color gradient from a set of colors
178 % at specific coordients, such as a list x,y -> r,g,b,a
180 % See LeastSquaresAddTerms() below for such an example.
182 % You can also use the 'vectors' to generate an inverse of the given 'matrix'
183 % though as a 'column first array' rather than a 'row first array' (matrix
186 % For details of this process see...
187 % http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
190 MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix,
191 double **vectors,const size_t rank,const size_t number_vectors)
193 #define GaussJordanSwap(x,y) \
219 columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns));
220 rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows));
221 pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots));
222 if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) ||
223 (pivots == (ssize_t *) NULL))
225 if (pivots != (ssize_t *) NULL)
226 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
227 if (columns != (ssize_t *) NULL)
228 columns=(ssize_t *) RelinquishMagickMemory(columns);
229 if (rows != (ssize_t *) NULL)
230 rows=(ssize_t *) RelinquishMagickMemory(rows);
233 (void) ResetMagickMemory(columns,0,rank*sizeof(*columns));
234 (void) ResetMagickMemory(rows,0,rank*sizeof(*rows));
235 (void) ResetMagickMemory(pivots,0,rank*sizeof(*pivots));
238 for (i=0; i < (ssize_t) rank; i++)
241 for (j=0; j < (ssize_t) rank; j++)
244 for (k=0; k < (ssize_t) rank; k++)
251 if (fabs(matrix[j][k]) >= max)
253 max=fabs(matrix[j][k]);
261 for (k=0; k < (ssize_t) rank; k++)
262 GaussJordanSwap(matrix[row][k],matrix[column][k]);
263 for (k=0; k < (ssize_t) number_vectors; k++)
264 GaussJordanSwap(vectors[k][row],vectors[k][column]);
268 if (matrix[column][column] == 0.0)
269 return(MagickFalse); /* singularity */
270 scale=1.0/matrix[column][column];
271 matrix[column][column]=1.0;
272 for (j=0; j < (ssize_t) rank; j++)
273 matrix[column][j]*=scale;
274 for (j=0; j < (ssize_t) number_vectors; j++)
275 vectors[j][column]*=scale;
276 for (j=0; j < (ssize_t) rank; j++)
279 scale=matrix[j][column];
280 matrix[j][column]=0.0;
281 for (k=0; k < (ssize_t) rank; k++)
282 matrix[j][k]-=scale*matrix[column][k];
283 for (k=0; k < (ssize_t) number_vectors; k++)
284 vectors[k][j]-=scale*vectors[k][column];
287 for (j=(ssize_t) rank-1; j >= 0; j--)
288 if (columns[j] != rows[j])
289 for (i=0; i < (ssize_t) rank; i++)
290 GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]);
291 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
292 rows=(ssize_t *) RelinquishMagickMemory(rows);
293 columns=(ssize_t *) RelinquishMagickMemory(columns);
298 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
302 % L e a s t S q u a r e s A d d T e r m s %
306 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
308 % LeastSquaresAddTerms() adds one set of terms and associate results to the
309 % given matrix and vectors for solving using least-squares function fitting.
311 % The format of the AcquireMagickMatrix method is:
313 % void LeastSquaresAddTerms(double **matrix,double **vectors,
314 % const double *terms,const double *results,const size_t rank,
315 % const size_t number_vectors);
317 % A description of each parameter follows:
319 % o matrix: the square matrix to add given terms/results to.
321 % o vectors: the result vectors to add terms/results to.
323 % o terms: the pre-calculated terms (without the unknown coefficent
324 % weights) that forms the equation being added.
326 % o results: the result(s) that should be generated from the given terms
327 % weighted by the yet-to-be-solved coefficents.
329 % o rank: the rank or size of the dimentions of the square matrix.
330 % Also the length of vectors, and number of terms being added.
332 % o number_vectors: Number of result vectors, and number or results being
333 % added. Also represents the number of separable systems of equations
334 % that is being solved.
338 % 2 dimensional Affine Equations (which are separable)
339 % c0*x + c2*y + c4*1 => u
340 % c1*x + c3*y + c5*1 => v
342 % double **matrix = AcquireMagickMatrix(3UL,3UL);
343 % double **vectors = AcquireMagickMatrix(2UL,3UL);
344 % double terms[3], results[2];
346 % for each given x,y -> u,v
352 % LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL);
354 % if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) {
355 % c0 = vectors[0][0];
356 % c2 = vectors[0][1]; %* weights to calculate u from any given x,y *%
357 % c4 = vectors[0][2];
358 % c1 = vectors[1][0];
359 % c3 = vectors[1][1]; %* weights for calculate v from any given x,y *%
360 % c5 = vectors[1][2];
363 % printf("Matrix unsolvable\n);
364 % RelinquishMagickMatrix(matrix,3UL);
365 % RelinquishMagickMatrix(vectors,2UL);
367 % More examples can be found in "distort.c"
370 MagickPrivate void LeastSquaresAddTerms(double **matrix,double **vectors,
371 const double *terms,const double *results,const size_t rank,
372 const size_t number_vectors)
378 for (j=0; j < (ssize_t) rank; j++)
380 for (i=0; i < (ssize_t) rank; i++)
381 matrix[i][j]+=terms[i]*terms[j];
382 for (i=0; i < (ssize_t) number_vectors; i++)
383 vectors[i][j]+=results[i]*terms[j];
388 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
392 % R e l i n q u i s h M a g i c k M a t r i x %
396 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
398 % RelinquishMagickMatrix() frees the previously acquired matrix (array of
399 % pointers to arrays of doubles).
401 % The format of the RelinquishMagickMatrix method is:
403 % double **RelinquishMagickMatrix(double **matrix,
404 % const size_t number_rows)
406 % A description of each parameter follows:
408 % o matrix: the matrix to relinquish
410 % o number_rows: the first dimension of the acquired matrix (number of
414 MagickExport double **RelinquishMagickMatrix(double **matrix,
415 const size_t number_rows)
420 if (matrix == (double **) NULL )
422 for (i=0; i < (ssize_t) number_rows; i++)
423 matrix[i]=(double *) RelinquishMagickMemory(matrix[i]);
424 matrix=(double **) RelinquishMagickMemory(matrix);