granicus.if.org Git - imagemagick/blob - magick/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-2010 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 "magick/studio.h"
  43 #include "magick/matrix.h"
  44 #include "magick/memory_.h"
  45 #include "magick/utility.h"
  46 \f
  47 /*
  48 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  49 %                                                                             %
  50 %                                                                             %
  51 %                                                                             %
  52 %   A c q u i r e M a g i c k M a t r i x                                     %
  53 %                                                                             %
  54 %                                                                             %
  55 %                                                                             %
  56 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  57 %
  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.
  60 %
  61 %  This used to generate the two dimentional matrix, and vectors required
  62 %  for the GaussJordanElimination() method below, solving some system of
  63 %  simultanious equations.
  64 %
  65 %  The format of the AcquireMagickMatrix method is:
  66 %
  67 %      double **AcquireMagickMatrix(const unsigned long nptrs,
  68 %           const unsigned long size)
  69 %
  70 %  A description of each parameter follows:
  71 %
  72 %    o nptrs: the number pointers for the array of pointers
  73 %             (first dimension)
  74 %
  75 %    o size: the size of the array of doubles each pointer points to.
  76 %            (second dimension)
  77 %
  78 */
  79 MagickExport double **AcquireMagickMatrix(const unsigned long nptrs,
  80      const unsigned long size)
  81 {
  82   double
  83    **matrix;
  84
  85   register unsigned long
  86     i,
  87     j;
  88
  89   matrix=(double **) AcquireQuantumMemory(nptrs,sizeof(*matrix));
  90   if (matrix == (double **) NULL)
  91     return((double **)NULL);
  92
  93   for (i=0; i < nptrs; i++)
  94   {
  95     matrix[i]=(double *) AcquireQuantumMemory(size,sizeof(*matrix[i]));
  96     if (matrix[i] == (double *) NULL)
  97     {
  98       for (j=0; j < i; j++)
  99         matrix[j]=(double *) RelinquishMagickMemory(matrix[j]);
 100       matrix=(double **) RelinquishMagickMemory(matrix);
 101       return((double **) NULL);
 102     }
 103     /*(void) ResetMagickMemory(matrix[i],0,size*sizeof(*matrix[i])); */
 104     for (j=0; j < size; j++)
 105       matrix[i][j] = 0.0;
 106   }
 107   return(matrix);
 108 }
 109 \f
 110 /*
 111 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 112 %                                                                             %
 113 %                                                                             %
 114 %                                                                             %
 115 %   G a u s s J o r d a n E l i m i n a t i o n                               %
 116 %                                                                             %
 117 %                                                                             %
 118 %                                                                             %
 119 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 120 %
 121 %  GaussJordanElimination() returns a matrix in reduced row echelon form,
 122 %  while simultaneously reducing and thus solving the augumented results
 123 %  matrix.
 124 %
 125 %  See also  http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
 126 %
 127 %  The format of the GaussJordanElimination method is:
 128 %
 129 %      MagickBooleanType GaussJordanElimination(double **matrix,
 130 %        double **vectors, const unsigned long rank, const unsigned long nvecs)
 131 %
 132 %  A description of each parameter follows:
 133 %
 134 %    o matrix: the matrix to be reduced, as an 'array of row pointers'.
 135 %
 136 %    o vectors: the additional matrix argumenting the matrix for row reduction.
 137 %             Producing an 'array of column vectors'.
 138 %
 139 %    o rank:  The size of the matrix (both rows and columns).
 140 %             Also represents the number terms that need to be solved.
 141 %
 142 %    o nvecs: Number of vectors columns, argumenting the above matrix.
 143 %             Usally 1, but can be more for more complex equation solving.
 144 %
 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'.
 149 %
 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 %
 153 %  This allows for simpler handling of the results, especially is only one
 154 %  column 'vector' is all that is required to produce the desired solution.
 155 %
 156 %  For example, the 'vectors' can consist of a pointer to a simple array of
 157 %  doubles.  when only one set of simultanious equations is to be solved from
 158 %  the given set of coefficient weighted terms.
 159 %
 160 %     double **matrix = AcquireMagickMatrix(8UL,8UL);
 161 %     double coefficents[8];
 162 %     ...
 163 %     GaussJordanElimination(matrix, &coefficents, 8UL, 1UL);
 164 %
 165 %  However by specifing more 'columns' (as an 'array of vector columns',
 166 %  you can use this function to solve a set of 'separable' equations.
 167 %
 168 %  For example a distortion function where    u = U(x,y)   v = V(x,y)
 169 %  And the functions U() and V() have separate coefficents, but are being
 170 %  generated from a common x,y->u,v  data set.
 171 %
 172 %  Another example is generation of a color gradient from a set of colors
 173 %  at specific coordients, such as a list    x,y -> r,g,b,a
 174 %  (Reference to be added - Anthony)
 175 %
 176 %  You can also use the 'vectors' to generate an inverse of the given 'matrix'
 177 %  though as a 'column first array' rather than a 'row first array'. For
 178 %  details see    http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
 179 %
 180 */
 181 MagickExport MagickBooleanType GaussJordanElimination(double **matrix,
 182   double **vectors, const unsigned long rank, const unsigned long nvecs)
 183 {
 184 #define GaussJordanSwap(x,y) \
 185 { \
 186   if ((x) != (y)) \
 187     { \
 188       (x)+=(y); \
 189       (y)=(x)-(y); \
 190       (x)=(x)-(y); \
 191     } \
 192 }
 193
 194   double
 195     max,
 196     scale;
 197
 198   long
 199     column,
 200     *columns,
 201     *pivots,
 202     row,
 203     *rows;
 204
 205   register long
 206     i,
 207     j,
 208     k;
 209
 210   columns=(long *) AcquireQuantumMemory(rank,sizeof(*columns));
 211   rows=(long *) AcquireQuantumMemory(rank,sizeof(*rows));
 212   pivots=(long *) AcquireQuantumMemory(rank,sizeof(*pivots));
 213   if ((rows == (long *) NULL) || (columns == (long *) NULL) ||
 214       (pivots == (long *) NULL))
 215     {
 216       if (pivots != (long *) NULL)
 217         pivots=(long *) RelinquishMagickMemory(pivots);
 218       if (columns != (long *) NULL)
 219         columns=(long *) RelinquishMagickMemory(columns);
 220       if (rows != (long *) NULL)
 221         rows=(long *) RelinquishMagickMemory(rows);
 222       return(MagickFalse);
 223     }
 224   (void) ResetMagickMemory(columns,0,rank*sizeof(*columns));
 225   (void) ResetMagickMemory(rows,0,rank*sizeof(*rows));
 226   (void) ResetMagickMemory(pivots,0,rank*sizeof(*pivots));
 227   column=0;
 228   row=0;
 229   for (i=0; i < (long) rank; i++)
 230   {
 231     max=0.0;
 232     for (j=0; j < (long) rank; j++)
 233       if (pivots[j] != 1)
 234         {
 235           for (k=0; k < (long) rank; k++)
 236             if (pivots[k] != 0)
 237               {
 238                 if (pivots[k] > 1)
 239                   return(MagickFalse);
 240               }
 241             else
 242               if (fabs(matrix[j][k]) >= max)
 243                 {
 244                   max=fabs(matrix[j][k]);
 245                   row=j;
 246                   column=k;
 247                 }
 248         }
 249     pivots[column]++;
 250     if (row != column)
 251       {
 252         for (k=0; k < (long) rank; k++)
 253           GaussJordanSwap(matrix[row][k],matrix[column][k]);
 254         for (k=0; k < (long) nvecs; k++)
 255           GaussJordanSwap(vectors[k][row],vectors[k][column]);
 256       }
 257     rows[i]=row;
 258     columns[i]=column;
 259     if (matrix[column][column] == 0.0)
 260       return(MagickFalse);  /* sigularity */
 261     scale=1.0/matrix[column][column];
 262     matrix[column][column]=1.0;
 263     for (j=0; j < (long) rank; j++)
 264       matrix[column][j]*=scale;
 265     for (j=0; j < (long) nvecs; j++)
 266       vectors[j][column]*=scale;
 267     for (j=0; j < (long) rank; j++)
 268       if (j != column)
 269         {
 270           scale=matrix[j][column];
 271           matrix[j][column]=0.0;
 272           for (k=0; k < (long) rank; k++)
 273             matrix[j][k]-=scale*matrix[column][k];
 274           for (k=0; k < (long) nvecs; k++)
 275             vectors[k][j]-=scale*vectors[k][column];
 276         }
 277   }
 278   for (j=(long) rank-1; j >= 0; j--)
 279     if (columns[j] != rows[j])
 280       for (i=0; i < (long) rank; i++)
 281         GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]);
 282   pivots=(long *) RelinquishMagickMemory(pivots);
 283   rows=(long *) RelinquishMagickMemory(rows);
 284   columns=(long *) RelinquishMagickMemory(columns);
 285   return(MagickTrue);
 286 }
 287 \f
 288 /*
 289 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 290 %                                                                             %
 291 %                                                                             %
 292 %                                                                             %
 293 %   L e a s t S q u a r e s A d d T e r m s                                   %
 294 %                                                                             %
 295 %                                                                             %
 296 %                                                                             %
 297 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 298 %
 299 %  LeastSquaresAddTerms() adds one set of terms and associate results to the
 300 %  given matrix and vectors for solving using least-squares function fitting.
 301 %
 302 %  The format of the AcquireMagickMatrix method is:
 303 %
 304 %      void LeastSquaresAddTerms(double **matrix,double **vectors,
 305 %             const double *terms, const double *results,
 306 %             const unsigned long rank, const unsigned long nvecs);
 307 %
 308 %  A description of each parameter follows:
 309 %
 310 %    o matrix: the square matrix to add given terms/results to.
 311 %
 312 %    o vectors: the result vectors to add terms/results to.
 313 %
 314 %    o terms: the pre-calculated terms (without the unknown coefficent
 315 %             weights) that forms the equation being added.
 316 %
 317 %    o results: the result(s) that should be generated from the given terms
 318 %               weighted by the yet-to-be-solved coefficents.
 319 %
 320 %    o rank: the rank or size of the dimentions of the square matrix.
 321 %            Also the length of vectors, and number of terms being added.
 322 %
 323 %    o nvecs: Number of result vectors, and number or results being added.
 324 %             Also represents the number of separable systems of equations
 325 %             that is being solved.
 326 %
 327 %  Example of use...
 328 %
 329 %     2 dimentional Affine Equations (which are separable)
 330 %         c0*x + c2*y + c4*1 => u
 331 %         c1*x + c3*y + c5*1 => v
 332 %
 333 %     double **matrix = AcquireMagickMatrix(3UL,3UL);
 334 %     double **vectors = AcquireMagickMatrix(2UL,3UL);
 335 %     double terms[3], results[2];
 336 %     ...
 337 %     for each given x,y -> u,v
 338 %        terms[0] = x;
 339 %        terms[1] = y;
 340 %        terms[2] = 1;
 341 %        results[0] = u;
 342 %        results[1] = v;
 343 %        LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL);
 344 %     ...
 345 %     if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) {
 346 %       c0 = vectors[0][0];
 347 %       c2 = vectors[0][1];
 348 %       c4 = vectors[0][2];
 349 %       c1 = vectors[1][0];
 350 %       c3 = vectors[1][1];
 351 %       c5 = vectors[1][2];
 352 %     }
 353 %     else
 354 %       printf("Matrix unsolvable\n);
 355 %     RelinquishMagickMatrix(matrix,3UL);
 356 %     RelinquishMagickMatrix(vectors,2UL);
 357 %
 358 */
 359 MagickExport void LeastSquaresAddTerms(double **matrix,double **vectors,
 360      const double *terms, const double *results, const unsigned long rank,
 361      const unsigned long nvecs)
 362 {
 363   register unsigned long
 364     i,
 365     j;
 366
 367   for(j=0; j<rank; j++) {
 368     for(i=0; i<rank; i++)
 369       matrix[i][j] += terms[i] * terms[j];
 370     for(i=0; i<nvecs; i++)
 371       vectors[i][j] += results[i] * terms[j];
 372   }
 373
 374   return;
 375 }
 376 \f
 377 /*
 378 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 379 %                                                                             %
 380 %                                                                             %
 381 %                                                                             %
 382 %   R e l i n q u i s h M a g i c k M a t r i x                               %
 383 %                                                                             %
 384 %                                                                             %
 385 %                                                                             %
 386 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 387 %
 388 %  RelinquishMagickMatrix() frees the previously acquired matrix (array of
 389 %  pointers to arrays of doubles).
 390 %
 391 %  The format of the RelinquishMagickMatrix method is:
 392 %
 393 %      double **RelinquishMagickMatrix(double **matrix,
 394 %         const unsigned long nptrs)
 395 %
 396 %  A description of each parameter follows:
 397 %
 398 %    o matrix: the matrix to relinquish
 399 %
 400 %    o nptrs: the first dimention of the acquired matrix (number of pointers)
 401 %
 402 */
 403 MagickExport double **RelinquishMagickMatrix(double **matrix,
 404      const unsigned long nptrs)
 405 {
 406   register unsigned long
 407     i;
 408
 409   if (matrix == (double **) NULL )
 410     return(matrix);
 411
 412   for (i=0; i < nptrs; i++)
 413      matrix[i]=(double *) RelinquishMagickMemory(matrix[i]);
 414   matrix=(double **) RelinquishMagickMemory(matrix);
 415
 416   return(matrix);
 417 }
 418