/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % M M AAA TTTTT RRRR IIIII X X % % MM MM A A T R R I X X % % M M M AAAAA T RRRR I X % % M M A A T R R I X X % % M M A A T R R IIIII X X % % % % % % ImageMagick Matrix Methods % % % % Software Design % % John Cristy % % August 2007 % % % % % % Copyright 1999-2007 ImageMagick Studio LLC, a non-profit organization % % dedicated to making software imaging solutions freely available. % % % % You may not use this file except in compliance with the License. You may % % obtain a copy of the License at % % % % http://www.imagemagick.org/script/license.php % % % % Unless required by applicable law or agreed to in writing, software % % distributed under the License is distributed on an "AS IS" BASIS, % % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. % % See the License for the specific language governing permissions and % % limitations under the License. % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % */ /* Include declarations. */ #include "magick/studio.h" #include "magick/matrix.h" #include "magick/memory_.h" #include "magick/utility.h" /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % A c q u i r e M a g i c k M a t r i x % % % % % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % AcquireMagickMatrix() allocates and returns a matrix in the form of an % array of pointers to an array of doubles, with all values pre-set to zero. % % This is used to generate the two dimentional matrix, and vectors required % for the GaussJordanElimination() method below, solving some system of % simultanious equations. % The format of the AcquireMagickMatrix method is: % % double **AcquireMagickMatrix(const unsigned long nptrs, % const unsigned long size) % % A description of each parameter follows: % % o nptrs: The number pointers for the array of pointers % (first dimension) % % o size: The size of the array of doubles each pointer points to. % (second dimension) % */ MagickExport double **AcquireMagickMatrix(const unsigned long nptrs, const unsigned long size) { double **matrix; register unsigned long i, j; matrix=(double **) AcquireQuantumMemory(nptrs,sizeof(*matrix)); if (matrix == (double **) NULL) return((double **)NULL); for (i=0; i < nptrs; i++) { matrix[i]=(double *) AcquireQuantumMemory(size,sizeof(*matrix[i])); if (matrix[i] == (double *) NULL) { for (j=0; j < i; j++) matrix[j]=(double *) RelinquishMagickMemory(matrix[j]); matrix=(double **) RelinquishMagickMemory(matrix); return((double **) NULL); } /*(void) ResetMagickMemory(matrix[i],0,size*sizeof(*matrix[i])); */ for (j=0; j < size; j++) matrix[i][j] = 0.0; } return(matrix); } /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % G a u s s J o r d a n E l i m i n a t i o n % % % % % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % GaussJordanElimination() returns a matrix in reduced row echelon form, % while simultaneously reducing and thus solving the augumented results % matrix. % % See also http://en.wikipedia.org/wiki/Gauss-Jordan_elimination % % The format of the GaussJordanElimination method is: % % MagickBooleanType GaussJordanElimination(double **matrix, % double **vectors, const unsigned long rank, const unsigned long nvecs) % % A description of each parameter follows: % % o matrix: The matrix to be reduced, as an 'array of row pointers'. % % o vectors: The additional matrix argumenting the matrix for row reduction. % Producing an 'array of column vectors'. % % o rank: The size of the matrix (both rows and columns). % Also represents the number terms that need to be solved. % % o nvecs: Number of vectors columns, argumenting the above matrix. % Usally 1, but can be more for more complex equation solving. % % Note that the 'matrix' is given as a 'array of row pointers' of rank size. % That is values can be assigned as matrix[row][column] where 'row' is % typically the equation, and 'column' is the term of the equation. % That is the matrix is in the form of a 'row first array'. % % However 'vectors' is a 'array of column pointers' which can have any number % of columns, with each column array the same 'rank' size as 'matrix'. % % This allows for simpler handling of the results, especially is only one % column 'vector' is all that is required to produce the desired solution. % % For example, the 'vectors' can consist of a pointer to a simple array of % doubles. when only one set of simultanious equations is to be solved from % the given set of coefficient weighted terms. % % double **matrix = AcquireMagickMatrix(8UL,8UL); % double coefficents[8]; % ... % GaussJordanElimination(matrix, &coefficents, 8UL, 1UL); % % However by specifing more 'columns' (as an 'array of vector columns', % you can use this function to solve a set of 'separable' equations. % % For example a distortion function where u = U(x,y) v = V(x,y) % And the functions U() and V() have separate coefficents, but are being % generated from a common x,y->u,v data set. % % Another example is generation of a color gradient from a set of colors % at specific coordients, such as a list x,y -> r,g,b,a % (Reference to be added - Anthony) % % You can also use the 'vectors' to generate an inverse of the given 'matrix' % though as a 'column first array' rather than a 'row first array'. For % details see http://en.wikipedia.org/wiki/Gauss-Jordan_elimination % */ MagickExport MagickBooleanType GaussJordanElimination(double **matrix, double **vectors, const unsigned long rank, const unsigned long nvecs) { #define GaussJordanSwap(x,y) \ { \ if ((x) != (y)) \ { \ (x)+=(y); \ (y)=(x)-(y); \ (x)=(x)-(y); \ } \ } double max, scale; long column, *columns, *pivots, row, *rows; register long i, j, k; columns=(long *) AcquireQuantumMemory(rank,sizeof(*columns)); rows=(long *) AcquireQuantumMemory(rank,sizeof(*rows)); pivots=(long *) AcquireQuantumMemory(rank,sizeof(*pivots)); if ((rows == (long *) NULL) || (columns == (long *) NULL) || (pivots == (long *) NULL)) { if (pivots != (long *) NULL) pivots=(long *) RelinquishMagickMemory(pivots); if (columns != (long *) NULL) columns=(long *) RelinquishMagickMemory(columns); if (rows != (long *) NULL) rows=(long *) RelinquishMagickMemory(rows); return(MagickFalse); } (void) ResetMagickMemory(columns,0,rank*sizeof(*columns)); (void) ResetMagickMemory(rows,0,rank*sizeof(*rows)); (void) ResetMagickMemory(pivots,0,rank*sizeof(*pivots)); column=0; row=0; for (i=0; i < (long) rank; i++) { max=0.0; for (j=0; j < (long) rank; j++) if (pivots[j] != 1) { for (k=0; k < (long) rank; k++) if (pivots[k] != 0) { if (pivots[k] > 1) return(MagickFalse); } else if (fabs(matrix[j][k]) >= max) { max=fabs(matrix[j][k]); row=j; column=k; } } pivots[column]++; if (row != column) { for (k=0; k < (long) rank; k++) GaussJordanSwap(matrix[row][k],matrix[column][k]); for (k=0; k < (long) nvecs; k++) GaussJordanSwap(vectors[k][row],vectors[k][column]); } rows[i]=row; columns[i]=column; if (matrix[column][column] == 0.0) return(MagickFalse); /* sigularity */ scale=1.0/matrix[column][column]; matrix[column][column]=1.0; for (j=0; j < (long) rank; j++) matrix[column][j]*=scale; for (j=0; j < (long) nvecs; j++) vectors[j][column]*=scale; for (j=0; j < (long) rank; j++) if (j != column) { scale=matrix[j][column]; matrix[j][column]=0.0; for (k=0; k < (long) rank; k++) matrix[j][k]-=scale*matrix[column][k]; for (k=0; k < (long) nvecs; k++) vectors[k][j]-=scale*vectors[k][column]; } } for (j=(long) rank-1; j >= 0; j--) if (columns[j] != rows[j]) for (i=0; i < (long) rank; i++) GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]); pivots=(long *) RelinquishMagickMemory(pivots); rows=(long *) RelinquishMagickMemory(rows); columns=(long *) RelinquishMagickMemory(columns); return(MagickTrue); } /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % L e a s t S q u a r e s A d d T e r m s % % % % % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % LeastSquaresAddTerms() adds one set of terms and associate results to the % given matrix and vectors for solving using least-squares function fitting. % % The format of the AcquireMagickMatrix method is: % % void LeastSquaresAddTerms(double **matrix,double **vectors, % const double *terms, const double *results, % const unsigned long rank, const unsigned long nvecs); % % A description of each parameter follows: % % o matrix: The square matrix to add given terms/results to. % % o vectors: The result vectors to add terms/results to. % % o terms: The pre-calculated terms (without the unknown coefficent % weights) that forms the equation being added. % % o results: The result(s) that should be generated from the given terms % weighted by the yet-to-be-solved coefficents. % % o rank: The rank or size of the dimentions of the square matrix. % Also the length of vectors, and number of terms being added. % % o nvecs: Number of result vectors, and number or results being added. % Also represents the number of separable systems of equations % that is being solved. % % Example of use... % % // 2 dimentional Affine Equations (which are separable) % // c0*x + c2*y + c4*1 => u % // c1*x + c3*y + c5*1 => v % % double **matrix = AcquireMagickMatrix(3UL,3UL); % double **vectors = AcquireMagickMatrix(2UL,3UL); % double terms[3], results[2]; % ... % //for each given x,y -> u,v % terms[0] = x; % terms[1] = y; % terms[2] = 1; % results[0] = u; % results[1] = v; % LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL); % ... % if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) { % c0 = vectors[0][0]; % c2 = vectors[0][1]; % c4 = vectors[0][2]; % c1 = vectors[1][0]; % c3 = vectors[1][1]; % c5 = vectors[1][2]; % } % else % printf("Matrix unsolvable\n); % RelinquishMagickMatrix(matrix,3UL); % RelinquishMagickMatrix(vectors,2UL); % */ MagickExport void LeastSquaresAddTerms(double **matrix,double **vectors, const double *terms, const double *results, const unsigned long rank, const unsigned long nvecs) { register unsigned long i, j; for(j=0; j