granicus.if.org Git - imagemagick/blob - MagickCore/matrix.c

   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                                  %
  16 %                              John Cristy                                    %
  17 %                              August 2007                                    %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
  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 \f
  39 /*
  40   Include declarations.
  41 */
  42 #include "MagickCore/studio.h"
  43 #include "MagickCore/matrix.h"
  44 #include "MagickCore/matrix-private.h"
  45 #include "MagickCore/pixel-private.h"
  46 #include "MagickCore/memory_.h"
  47 \f
  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 %
  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.
  67 %
  68 %  The format of the AcquireMagickMatrix method is:
  69 %
  70 %      double **AcquireMagickMatrix(const size_t number_rows,
  71 %        const size_t size)
  72 %
  73 %  A description of each parameter follows:
  74 %
  75 %    o number_rows: the number pointers for the array of pointers
  76 %      (first dimension).
  77 %
  78 %    o size: the size of the array of doubles each pointer points to
  79 %      (second dimension).
  80 %
  81 */
  82 MagickExport double **AcquireMagickMatrix(const size_t number_rows,
  83   const size_t size)
  84 {
  85   double
  86     **matrix;
  87
  88   register ssize_t
  89     i,
  90     j;
  91
  92   matrix=(double **) AcquireQuantumMemory(number_rows,sizeof(*matrix));
  93   if (matrix == (double **) NULL)
  94     return((double **) NULL);
  95   for (i=0; i < (ssize_t) number_rows; i++)
  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     }
 105     for (j=0; j < (ssize_t) size; j++)
 106       matrix[i][j]=0.0;
 107   }
 108   return(matrix);
 109 }
 110 \f
 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 %
 130 %      MagickBooleanType GaussJordanElimination(double **matrix,
 131 %        double **vectors,const size_t rank,const size_t number_vectors)
 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 %
 140 %    o rank:  The size of the square matrix (both rows and columns).
 141 %             Also represents the number terms that need to be solved.
 142 %
 143 %    o number_vectors: Number of vectors columns, argumenting the above matrix.
 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'.
 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.
 157 %
 158 %  This allows for simpler handling of the results, especially is only one
 159 %  column 'vector' is all that is required to produce the desired solution
 160 %  for that specific set of equations.
 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);
 170 %
 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.
 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
 180 %
 181 %  See LeastSquaresAddTerms() below for such an example.
 182 %
 183 %  You can also use the 'vectors' to generate an inverse of the given 'matrix'
 184 %  though as a 'column first array' rather than a 'row first array' (matrix
 185 %  is transposed).
 186 %
 187 %  For details of this process see...
 188 %     http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
 189 %
 190 */
 191 MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix,
 192   double **vectors,const size_t rank,const size_t number_vectors)
 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
 208   register ssize_t
 209     i,
 210     j,
 211     k;
 212
 213   ssize_t
 214     column,
 215     *columns,
 216     *pivots,
 217     row,
 218     *rows;
 219
 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))
 225     {
 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);
 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;
 239   for (i=0; i < (ssize_t) rank; i++)
 240   {
 241     max=0.0;
 242     for (j=0; j < (ssize_t) rank; j++)
 243       if (pivots[j] != 1)
 244         {
 245           for (k=0; k < (ssize_t) rank; k++)
 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       {
 262         for (k=0; k < (ssize_t) rank; k++)
 263           GaussJordanSwap(matrix[row][k],matrix[column][k]);
 264         for (k=0; k < (ssize_t) number_vectors; k++)
 265           GaussJordanSwap(vectors[k][row],vectors[k][column]);
 266       }
 267     rows[i]=row;
 268     columns[i]=column;
 269     if (matrix[column][column] == 0.0)
 270       return(MagickFalse);  /* singularity */
 271     scale=MagickEpsilonReciprocal(matrix[column][column]);
 272     matrix[column][column]=1.0;
 273     for (j=0; j < (ssize_t) rank; j++)
 274       matrix[column][j]*=scale;
 275     for (j=0; j < (ssize_t) number_vectors; j++)
 276       vectors[j][column]*=scale;
 277     for (j=0; j < (ssize_t) rank; j++)
 278       if (j != column)
 279         {
 280           scale=matrix[j][column];
 281           matrix[j][column]=0.0;
 282           for (k=0; k < (ssize_t) rank; k++)
 283             matrix[j][k]-=scale*matrix[column][k];
 284           for (k=0; k < (ssize_t) number_vectors; k++)
 285             vectors[k][j]-=scale*vectors[k][column];
 286         }
 287   }
 288   for (j=(ssize_t) rank-1; j >= 0; j--)
 289     if (columns[j] != rows[j])
 290       for (i=0; i < (ssize_t) rank; i++)
 291         GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]);
 292   pivots=(ssize_t *) RelinquishMagickMemory(pivots);
 293   rows=(ssize_t *) RelinquishMagickMemory(rows);
 294   columns=(ssize_t *) RelinquishMagickMemory(columns);
 295   return(MagickTrue);
 296 }
 297 \f
 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,
 315 %        const double *terms,const double *results,const size_t rank,
 316 %        const size_t number_vectors);
 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
 325 %      weights) that forms the equation being added.
 326 %
 327 %    o results: the result(s) that should be generated from the given terms
 328 %      weighted by the yet-to-be-solved coefficents.
 329 %
 330 %    o rank: the rank or size of the dimentions of the square matrix.
 331 %      Also the length of vectors, and number of terms being added.
 332 %
 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.
 336 %
 337 %  Example of use...
 338 %
 339 %     2 dimensional Affine Equations (which are separable)
 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];
 357 %       c2 = vectors[0][1];  %* weights to calculate u from any given x,y *%
 358 %       c4 = vectors[0][2];
 359 %       c1 = vectors[1][0];
 360 %       c3 = vectors[1][1];  %* weights for calculate v from any given x,y *%
 361 %       c5 = vectors[1][2];
 362 %     }
 363 %     else
 364 %       printf("Matrix unsolvable\n);
 365 %     RelinquishMagickMatrix(matrix,3UL);
 366 %     RelinquishMagickMatrix(vectors,2UL);
 367 %
 368 % More examples can be found in "distort.c"
 369 %
 370 */
 371 MagickPrivate void LeastSquaresAddTerms(double **matrix,double **vectors,
 372   const double *terms,const double *results,const size_t rank,
 373   const size_t number_vectors)
 374 {
 375   register ssize_t
 376     i,
 377     j;
 378
 379   for (j=0; j < (ssize_t) rank; j++)
 380   {
 381     for (i=0; i < (ssize_t) rank; i++)
 382       matrix[i][j]+=terms[i]*terms[j];
 383     for (i=0; i < (ssize_t) number_vectors; i++)
 384       vectors[i][j]+=results[i]*terms[j];
 385   }
 386 }
 387 \f
 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,
 405 %        const size_t number_rows)
 406 %
 407 %  A description of each parameter follows:
 408 %
 409 %    o matrix: the matrix to relinquish
 410 %
 411 %    o number_rows: the first dimension of the acquired matrix (number of
 412 %      pointers)
 413 %
 414 */
 415 MagickExport double **RelinquishMagickMatrix(double **matrix,
 416   const size_t number_rows)
 417 {
 418   register ssize_t
 419     i;
 420
 421   if (matrix == (double **) NULL )
 422     return(matrix);
 423   for (i=0; i < (ssize_t) number_rows; i++)
 424      matrix[i]=(double *) RelinquishMagickMemory(matrix[i]);
 425   matrix=(double **) RelinquishMagickMemory(matrix);
 426   return(matrix);
 427 }