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 %                                 Cristy                                      %
  17 %                              August 2007                                    %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2017 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 %    https://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/blob.h"
  44 #include "MagickCore/blob-private.h"
  45 #include "MagickCore/cache.h"
  46 #include "MagickCore/exception.h"
  47 #include "MagickCore/exception-private.h"
  48 #include "MagickCore/image-private.h"
  49 #include "MagickCore/matrix.h"
  50 #include "MagickCore/matrix-private.h"
  51 #include "MagickCore/memory_.h"
  52 #include "MagickCore/pixel-accessor.h"
  53 #include "MagickCore/pixel-private.h"
  54 #include "MagickCore/resource_.h"
  55 #include "MagickCore/semaphore.h"
  56 #include "MagickCore/thread-private.h"
  57 #include "MagickCore/utility.h"
  58 \f
  59 /*
  60   Typedef declaration.
  61 */
  62 struct _MatrixInfo
  63 {
  64   CacheType
  65     type;
  66
  67   size_t
  68     columns,
  69     rows,
  70     stride;
  71
  72   MagickSizeType
  73     length;
  74
  75   MagickBooleanType
  76     mapped,
  77     synchronize;
  78
  79   char
  80     path[MagickPathExtent];
  81
  82   int
  83     file;
  84
  85   void
  86     *elements;
  87
  88   SemaphoreInfo
  89     *semaphore;
  90
  91   size_t
  92     signature;
  93 };
  94 \f
  95 /*
  96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  97 %                                                                             %
  98 %                                                                             %
  99 %                                                                             %
 100 %   A c q u i r e M a t r i x I n f o                                         %
 101 %                                                                             %
 102 %                                                                             %
 103 %                                                                             %
 104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 105 %
 106 %  AcquireMatrixInfo() allocates the ImageInfo structure.
 107 %
 108 %  The format of the AcquireMatrixInfo method is:
 109 %
 110 %      MatrixInfo *AcquireMatrixInfo(const size_t columns,const size_t rows,
 111 %        const size_t stride,ExceptionInfo *exception)
 112 %
 113 %  A description of each parameter follows:
 114 %
 115 %    o columns: the matrix columns.
 116 %
 117 %    o rows: the matrix rows.
 118 %
 119 %    o stride: the matrix stride.
 120 %
 121 %    o exception: return any errors or warnings in this structure.
 122 %
 123 */
 124
 125 #if defined(SIGBUS)
 126 static void MatrixSignalHandler(int status)
 127 {
 128   ThrowFatalException(CacheFatalError,"UnableToExtendMatrixCache");
 129 }
 130 #endif
 131
 132 static inline MagickOffsetType WriteMatrixElements(
 133   const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset,
 134   const MagickSizeType length,const unsigned char *magick_restrict buffer)
 135 {
 136   register MagickOffsetType
 137     i;
 138
 139   ssize_t
 140     count;
 141
 142 #if !defined(MAGICKCORE_HAVE_PWRITE)
 143   LockSemaphoreInfo(matrix_info->semaphore);
 144   if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
 145     {
 146       UnlockSemaphoreInfo(matrix_info->semaphore);
 147       return((MagickOffsetType) -1);
 148     }
 149 #endif
 150   count=0;
 151   for (i=0; i < (MagickOffsetType) length; i+=count)
 152   {
 153 #if !defined(MAGICKCORE_HAVE_PWRITE)
 154     count=write(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
 155       (MagickSizeType) SSIZE_MAX));
 156 #else
 157     count=pwrite(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
 158       (MagickSizeType) SSIZE_MAX),(off_t) (offset+i));
 159 #endif
 160     if (count <= 0)
 161       {
 162         count=0;
 163         if (errno != EINTR)
 164           break;
 165       }
 166   }
 167 #if !defined(MAGICKCORE_HAVE_PWRITE)
 168   UnlockSemaphoreInfo(matrix_info->semaphore);
 169 #endif
 170   return(i);
 171 }
 172
 173 static MagickBooleanType SetMatrixExtent(
 174   MatrixInfo *magick_restrict matrix_info,MagickSizeType length)
 175 {
 176   MagickOffsetType
 177     count,
 178     extent,
 179     offset;
 180
 181   if (length != (MagickSizeType) ((MagickOffsetType) length))
 182     return(MagickFalse);
 183   offset=(MagickOffsetType) lseek(matrix_info->file,0,SEEK_END);
 184   if (offset < 0)
 185     return(MagickFalse);
 186   if ((MagickSizeType) offset >= length)
 187     return(MagickTrue);
 188   extent=(MagickOffsetType) length-1;
 189   count=WriteMatrixElements(matrix_info,extent,1,(const unsigned char *) "");
 190 #if defined(MAGICKCORE_HAVE_POSIX_FALLOCATE)
 191   if (matrix_info->synchronize != MagickFalse)
 192     (void) posix_fallocate(matrix_info->file,offset+1,extent-offset);
 193 #endif
 194 #if defined(SIGBUS)
 195   (void) signal(SIGBUS,MatrixSignalHandler);
 196 #endif
 197   return(count != (MagickOffsetType) 1 ? MagickFalse : MagickTrue);
 198 }
 199
 200 MagickExport MatrixInfo *AcquireMatrixInfo(const size_t columns,
 201   const size_t rows,const size_t stride,ExceptionInfo *exception)
 202 {
 203   char
 204     *synchronize;
 205
 206   MagickBooleanType
 207     status;
 208
 209   MatrixInfo
 210     *matrix_info;
 211
 212   matrix_info=(MatrixInfo *) AcquireMagickMemory(sizeof(*matrix_info));
 213   if (matrix_info == (MatrixInfo *) NULL)
 214     return((MatrixInfo *) NULL);
 215   (void) ResetMagickMemory(matrix_info,0,sizeof(*matrix_info));
 216   matrix_info->signature=MagickCoreSignature;
 217   matrix_info->columns=columns;
 218   matrix_info->rows=rows;
 219   matrix_info->stride=stride;
 220   matrix_info->semaphore=AcquireSemaphoreInfo();
 221   synchronize=GetEnvironmentValue("MAGICK_SYNCHRONIZE");
 222   if (synchronize != (const char *) NULL)
 223     {
 224       matrix_info->synchronize=IsStringTrue(synchronize);
 225       synchronize=DestroyString(synchronize);
 226     }
 227   matrix_info->length=(MagickSizeType) columns*rows*stride;
 228   if (matrix_info->columns != (size_t) (matrix_info->length/rows/stride))
 229     {
 230       (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
 231         "CacheResourcesExhausted","`%s'","matrix cache");
 232       return(DestroyMatrixInfo(matrix_info));
 233     }
 234   matrix_info->type=MemoryCache;
 235   status=AcquireMagickResource(AreaResource,matrix_info->length);
 236   if ((status != MagickFalse) &&
 237       (matrix_info->length == (MagickSizeType) ((size_t) matrix_info->length)))
 238     {
 239       status=AcquireMagickResource(MemoryResource,matrix_info->length);
 240       if (status != MagickFalse)
 241         {
 242           matrix_info->mapped=MagickFalse;
 243           matrix_info->elements=AcquireMagickMemory((size_t)
 244             matrix_info->length);
 245           if (matrix_info->elements == NULL)
 246             {
 247               matrix_info->mapped=MagickTrue;
 248               matrix_info->elements=MapBlob(-1,IOMode,0,(size_t)
 249                 matrix_info->length);
 250             }
 251           if (matrix_info->elements == (unsigned short *) NULL)
 252             RelinquishMagickResource(MemoryResource,matrix_info->length);
 253         }
 254     }
 255   matrix_info->file=(-1);
 256   if (matrix_info->elements == (unsigned short *) NULL)
 257     {
 258       status=AcquireMagickResource(DiskResource,matrix_info->length);
 259       if (status == MagickFalse)
 260         {
 261           (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
 262             "CacheResourcesExhausted","`%s'","matrix cache");
 263           return(DestroyMatrixInfo(matrix_info));
 264         }
 265       matrix_info->type=DiskCache;
 266       matrix_info->file=AcquireUniqueFileResource(matrix_info->path);
 267       if (matrix_info->file == -1)
 268         return(DestroyMatrixInfo(matrix_info));
 269       status=AcquireMagickResource(MapResource,matrix_info->length);
 270       if (status != MagickFalse)
 271         {
 272           status=SetMatrixExtent(matrix_info,matrix_info->length);
 273           if (status != MagickFalse)
 274             matrix_info->elements=(void *) MapBlob(matrix_info->file,IOMode,0,
 275               (size_t) matrix_info->length);
 276           if (matrix_info->elements != NULL)
 277             matrix_info->type=MapCache;
 278           else
 279             RelinquishMagickResource(MapResource,matrix_info->length);
 280         }
 281     }
 282   return(matrix_info);
 283 }
 284 \f
 285 /*
 286 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 287 %                                                                             %
 288 %                                                                             %
 289 %                                                                             %
 290 %   A c q u i r e M a g i c k M a t r i x                                     %
 291 %                                                                             %
 292 %                                                                             %
 293 %                                                                             %
 294 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 295 %
 296 %  AcquireMagickMatrix() allocates and returns a matrix in the form of an
 297 %  array of pointers to an array of doubles, with all values pre-set to zero.
 298 %
 299 %  This used to generate the two dimensional matrix, and vectors required
 300 %  for the GaussJordanElimination() method below, solving some system of
 301 %  simultanious equations.
 302 %
 303 %  The format of the AcquireMagickMatrix method is:
 304 %
 305 %      double **AcquireMagickMatrix(const size_t number_rows,
 306 %        const size_t size)
 307 %
 308 %  A description of each parameter follows:
 309 %
 310 %    o number_rows: the number pointers for the array of pointers
 311 %      (first dimension).
 312 %
 313 %    o size: the size of the array of doubles each pointer points to
 314 %      (second dimension).
 315 %
 316 */
 317 MagickExport double **AcquireMagickMatrix(const size_t number_rows,
 318   const size_t size)
 319 {
 320   double
 321     **matrix;
 322
 323   register ssize_t
 324     i,
 325     j;
 326
 327   matrix=(double **) AcquireQuantumMemory(number_rows,sizeof(*matrix));
 328   if (matrix == (double **) NULL)
 329     return((double **) NULL);
 330   for (i=0; i < (ssize_t) number_rows; i++)
 331   {
 332     matrix[i]=(double *) AcquireQuantumMemory(size,sizeof(*matrix[i]));
 333     if (matrix[i] == (double *) NULL)
 334       {
 335         for (j=0; j < i; j++)
 336           matrix[j]=(double *) RelinquishMagickMemory(matrix[j]);
 337         matrix=(double **) RelinquishMagickMemory(matrix);
 338         return((double **) NULL);
 339       }
 340     for (j=0; j < (ssize_t) size; j++)
 341       matrix[i][j]=0.0;
 342   }
 343   return(matrix);
 344 }
 345 \f
 346 /*
 347 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 348 %                                                                             %
 349 %                                                                             %
 350 %                                                                             %
 351 %   D e s t r o y M a t r i x I n f o                                         %
 352 %                                                                             %
 353 %                                                                             %
 354 %                                                                             %
 355 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 356 %
 357 %  DestroyMatrixInfo() dereferences a matrix, deallocating memory associated
 358 %  with the matrix.
 359 %
 360 %  The format of the DestroyImage method is:
 361 %
 362 %      MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info)
 363 %
 364 %  A description of each parameter follows:
 365 %
 366 %    o matrix_info: the matrix.
 367 %
 368 */
 369 MagickExport MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info)
 370 {
 371   assert(matrix_info != (MatrixInfo *) NULL);
 372   assert(matrix_info->signature == MagickCoreSignature);
 373   LockSemaphoreInfo(matrix_info->semaphore);
 374   switch (matrix_info->type)
 375   {
 376     case MemoryCache:
 377     {
 378       if (matrix_info->mapped == MagickFalse)
 379         matrix_info->elements=RelinquishMagickMemory(matrix_info->elements);
 380       else
 381         {
 382           (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length);
 383           matrix_info->elements=(unsigned short *) NULL;
 384         }
 385       RelinquishMagickResource(MemoryResource,matrix_info->length);
 386       break;
 387     }
 388     case MapCache:
 389     {
 390       (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length);
 391       matrix_info->elements=NULL;
 392       RelinquishMagickResource(MapResource,matrix_info->length);
 393     }
 394     case DiskCache:
 395     {
 396       if (matrix_info->file != -1)
 397         (void) close(matrix_info->file);
 398       (void) RelinquishUniqueFileResource(matrix_info->path);
 399       RelinquishMagickResource(DiskResource,matrix_info->length);
 400       break;
 401     }
 402     default:
 403       break;
 404   }
 405   UnlockSemaphoreInfo(matrix_info->semaphore);
 406   RelinquishSemaphoreInfo(&matrix_info->semaphore);
 407   return((MatrixInfo *) RelinquishMagickMemory(matrix_info));
 408 }
 409 \f
 410 /*
 411 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 412 %                                                                             %
 413 %                                                                             %
 414 %                                                                             %
 415 +   G a u s s J o r d a n E l i m i n a t i o n                               %
 416 %                                                                             %
 417 %                                                                             %
 418 %                                                                             %
 419 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 420 %
 421 %  GaussJordanElimination() returns a matrix in reduced row echelon form,
 422 %  while simultaneously reducing and thus solving the augumented results
 423 %  matrix.
 424 %
 425 %  See also  http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
 426 %
 427 %  The format of the GaussJordanElimination method is:
 428 %
 429 %      MagickBooleanType GaussJordanElimination(double **matrix,
 430 %        double **vectors,const size_t rank,const size_t number_vectors)
 431 %
 432 %  A description of each parameter follows:
 433 %
 434 %    o matrix: the matrix to be reduced, as an 'array of row pointers'.
 435 %
 436 %    o vectors: the additional matrix argumenting the matrix for row reduction.
 437 %             Producing an 'array of column vectors'.
 438 %
 439 %    o rank:  The size of the matrix (both rows and columns).
 440 %             Also represents the number terms that need to be solved.
 441 %
 442 %    o number_vectors: Number of vectors columns, argumenting the above matrix.
 443 %             Usally 1, but can be more for more complex equation solving.
 444 %
 445 %  Note that the 'matrix' is given as a 'array of row pointers' of rank size.
 446 %  That is values can be assigned as   matrix[row][column]   where 'row' is
 447 %  typically the equation, and 'column' is the term of the equation.
 448 %  That is the matrix is in the form of a 'row first array'.
 449 %
 450 %  However 'vectors' is a 'array of column pointers' which can have any number
 451 %  of columns, with each column array the same 'rank' size as 'matrix'.
 452 %
 453 %  This allows for simpler handling of the results, especially is only one
 454 %  column 'vector' is all that is required to produce the desired solution.
 455 %
 456 %  For example, the 'vectors' can consist of a pointer to a simple array of
 457 %  doubles.  when only one set of simultanious equations is to be solved from
 458 %  the given set of coefficient weighted terms.
 459 %
 460 %     double **matrix = AcquireMagickMatrix(8UL,8UL);
 461 %     double coefficents[8];
 462 %     ...
 463 %     GaussJordanElimination(matrix, &coefficents, 8UL, 1UL);
 464 %
 465 %  However by specifing more 'columns' (as an 'array of vector columns',
 466 %  you can use this function to solve a set of 'separable' equations.
 467 %
 468 %  For example a distortion function where    u = U(x,y)   v = V(x,y)
 469 %  And the functions U() and V() have separate coefficents, but are being
 470 %  generated from a common x,y->u,v  data set.
 471 %
 472 %  Another example is generation of a color gradient from a set of colors at
 473 %  specific coordients, such as a list x,y -> r,g,b,a.
 474 %
 475 %  You can also use the 'vectors' to generate an inverse of the given 'matrix'
 476 %  though as a 'column first array' rather than a 'row first array'. For
 477 %  details see http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
 478 %
 479 */
 480 MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix,
 481   double **vectors,const size_t rank,const size_t number_vectors)
 482 {
 483 #define GaussJordanSwap(x,y) \
 484 { \
 485   if ((x) != (y)) \
 486     { \
 487       (x)+=(y); \
 488       (y)=(x)-(y); \
 489       (x)=(x)-(y); \
 490     } \
 491 }
 492
 493   double
 494     max,
 495     scale;
 496
 497   register ssize_t
 498     i,
 499     j,
 500     k;
 501
 502   ssize_t
 503     column,
 504     *columns,
 505     *pivots,
 506     row,
 507     *rows;
 508
 509   columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns));
 510   rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows));
 511   pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots));
 512   if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) ||
 513       (pivots == (ssize_t *) NULL))
 514     {
 515       if (pivots != (ssize_t *) NULL)
 516         pivots=(ssize_t *) RelinquishMagickMemory(pivots);
 517       if (columns != (ssize_t *) NULL)
 518         columns=(ssize_t *) RelinquishMagickMemory(columns);
 519       if (rows != (ssize_t *) NULL)
 520         rows=(ssize_t *) RelinquishMagickMemory(rows);
 521       return(MagickFalse);
 522     }
 523   (void) ResetMagickMemory(columns,0,rank*sizeof(*columns));
 524   (void) ResetMagickMemory(rows,0,rank*sizeof(*rows));
 525   (void) ResetMagickMemory(pivots,0,rank*sizeof(*pivots));
 526   column=0;
 527   row=0;
 528   for (i=0; i < (ssize_t) rank; i++)
 529   {
 530     max=0.0;
 531     for (j=0; j < (ssize_t) rank; j++)
 532       if (pivots[j] != 1)
 533         {
 534           for (k=0; k < (ssize_t) rank; k++)
 535             if (pivots[k] != 0)
 536               {
 537                 if (pivots[k] > 1)
 538                   return(MagickFalse);
 539               }
 540             else
 541               if (fabs(matrix[j][k]) >= max)
 542                 {
 543                   max=fabs(matrix[j][k]);
 544                   row=j;
 545                   column=k;
 546                 }
 547         }
 548     pivots[column]++;
 549     if (row != column)
 550       {
 551         for (k=0; k < (ssize_t) rank; k++)
 552           GaussJordanSwap(matrix[row][k],matrix[column][k]);
 553         for (k=0; k < (ssize_t) number_vectors; k++)
 554           GaussJordanSwap(vectors[k][row],vectors[k][column]);
 555       }
 556     rows[i]=row;
 557     columns[i]=column;
 558     if (matrix[column][column] == 0.0)
 559       return(MagickFalse);  /* sigularity */
 560     scale=PerceptibleReciprocal(matrix[column][column]);
 561     matrix[column][column]=1.0;
 562     for (j=0; j < (ssize_t) rank; j++)
 563       matrix[column][j]*=scale;
 564     for (j=0; j < (ssize_t) number_vectors; j++)
 565       vectors[j][column]*=scale;
 566     for (j=0; j < (ssize_t) rank; j++)
 567       if (j != column)
 568         {
 569           scale=matrix[j][column];
 570           matrix[j][column]=0.0;
 571           for (k=0; k < (ssize_t) rank; k++)
 572             matrix[j][k]-=scale*matrix[column][k];
 573           for (k=0; k < (ssize_t) number_vectors; k++)
 574             vectors[k][j]-=scale*vectors[k][column];
 575         }
 576   }
 577   for (j=(ssize_t) rank-1; j >= 0; j--)
 578     if (columns[j] != rows[j])
 579       for (i=0; i < (ssize_t) rank; i++)
 580         GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]);
 581   pivots=(ssize_t *) RelinquishMagickMemory(pivots);
 582   rows=(ssize_t *) RelinquishMagickMemory(rows);
 583   columns=(ssize_t *) RelinquishMagickMemory(columns);
 584   return(MagickTrue);
 585 }
 586 \f
 587 /*
 588 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 589 %                                                                             %
 590 %                                                                             %
 591 %                                                                             %
 592 %   G e t M a t r i x C o l u m n s                                           %
 593 %                                                                             %
 594 %                                                                             %
 595 %                                                                             %
 596 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 597 %
 598 %  GetMatrixColumns() returns the number of columns in the matrix.
 599 %
 600 %  The format of the GetMatrixColumns method is:
 601 %
 602 %      size_t GetMatrixColumns(const MatrixInfo *matrix_info)
 603 %
 604 %  A description of each parameter follows:
 605 %
 606 %    o matrix_info: the matrix.
 607 %
 608 */
 609 MagickExport size_t GetMatrixColumns(const MatrixInfo *matrix_info)
 610 {
 611   assert(matrix_info != (MatrixInfo *) NULL);
 612   assert(matrix_info->signature == MagickCoreSignature);
 613   return(matrix_info->columns);
 614 }
 615 \f
 616 /*
 617 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 618 %                                                                             %
 619 %                                                                             %
 620 %                                                                             %
 621 %   G e t M a t r i x E l e m e n t                                           %
 622 %                                                                             %
 623 %                                                                             %
 624 %                                                                             %
 625 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 626 %
 627 %  GetMatrixElement() returns the specifed element in the matrix.
 628 %
 629 %  The format of the GetMatrixElement method is:
 630 %
 631 %      MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
 632 %        const ssize_t x,const ssize_t y,void *value)
 633 %
 634 %  A description of each parameter follows:
 635 %
 636 %    o matrix_info: the matrix columns.
 637 %
 638 %    o x: the matrix x-offset.
 639 %
 640 %    o y: the matrix y-offset.
 641 %
 642 %    o value: return the matrix element in this buffer.
 643 %
 644 */
 645
 646 static inline ssize_t EdgeX(const ssize_t x,const size_t columns)
 647 {
 648   if (x < 0L)
 649     return(0L);
 650   if (x >= (ssize_t) columns)
 651     return((ssize_t) (columns-1));
 652   return(x);
 653 }
 654
 655 static inline ssize_t EdgeY(const ssize_t y,const size_t rows)
 656 {
 657   if (y < 0L)
 658     return(0L);
 659   if (y >= (ssize_t) rows)
 660     return((ssize_t) (rows-1));
 661   return(y);
 662 }
 663
 664 static inline MagickOffsetType ReadMatrixElements(
 665   const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset,
 666   const MagickSizeType length,unsigned char *magick_restrict buffer)
 667 {
 668   register MagickOffsetType
 669     i;
 670
 671   ssize_t
 672     count;
 673
 674 #if !defined(MAGICKCORE_HAVE_PREAD)
 675   LockSemaphoreInfo(matrix_info->semaphore);
 676   if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
 677     {
 678       UnlockSemaphoreInfo(matrix_info->semaphore);
 679       return((MagickOffsetType) -1);
 680     }
 681 #endif
 682   count=0;
 683   for (i=0; i < (MagickOffsetType) length; i+=count)
 684   {
 685 #if !defined(MAGICKCORE_HAVE_PREAD)
 686     count=read(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
 687       (MagickSizeType) SSIZE_MAX));
 688 #else
 689     count=pread(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
 690       (MagickSizeType) SSIZE_MAX),(off_t) (offset+i));
 691 #endif
 692     if (count <= 0)
 693       {
 694         count=0;
 695         if (errno != EINTR)
 696           break;
 697       }
 698   }
 699 #if !defined(MAGICKCORE_HAVE_PREAD)
 700   UnlockSemaphoreInfo(matrix_info->semaphore);
 701 #endif
 702   return(i);
 703 }
 704
 705 MagickExport MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
 706   const ssize_t x,const ssize_t y,void *value)
 707 {
 708   MagickOffsetType
 709     count,
 710     i;
 711
 712   assert(matrix_info != (const MatrixInfo *) NULL);
 713   assert(matrix_info->signature == MagickCoreSignature);
 714   i=(MagickOffsetType) EdgeY(y,matrix_info->rows)*matrix_info->columns+
 715     EdgeX(x,matrix_info->columns);
 716   if (matrix_info->type != DiskCache)
 717     {
 718       (void) memcpy(value,(unsigned char *) matrix_info->elements+i*
 719         matrix_info->stride,matrix_info->stride);
 720       return(MagickTrue);
 721     }
 722   count=ReadMatrixElements(matrix_info,i*matrix_info->stride,
 723     matrix_info->stride,(unsigned char *) value);
 724   if (count != (MagickOffsetType) matrix_info->stride)
 725     return(MagickFalse);
 726   return(MagickTrue);
 727 }
 728 \f
 729 /*
 730 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 731 %                                                                             %
 732 %                                                                             %
 733 %                                                                             %
 734 %   G e t M a t r i x R o w s                                                 %
 735 %                                                                             %
 736 %                                                                             %
 737 %                                                                             %
 738 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 739 %
 740 %  GetMatrixRows() returns the number of rows in the matrix.
 741 %
 742 %  The format of the GetMatrixRows method is:
 743 %
 744 %      size_t GetMatrixRows(const MatrixInfo *matrix_info)
 745 %
 746 %  A description of each parameter follows:
 747 %
 748 %    o matrix_info: the matrix.
 749 %
 750 */
 751 MagickExport size_t GetMatrixRows(const MatrixInfo *matrix_info)
 752 {
 753   assert(matrix_info != (const MatrixInfo *) NULL);
 754   assert(matrix_info->signature == MagickCoreSignature);
 755   return(matrix_info->rows);
 756 }
 757 \f
 758 /*
 759 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 760 %                                                                             %
 761 %                                                                             %
 762 %                                                                             %
 763 +   L e a s t S q u a r e s A d d T e r m s                                   %
 764 %                                                                             %
 765 %                                                                             %
 766 %                                                                             %
 767 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 768 %
 769 %  LeastSquaresAddTerms() adds one set of terms and associate results to the
 770 %  given matrix and vectors for solving using least-squares function fitting.
 771 %
 772 %  The format of the AcquireMagickMatrix method is:
 773 %
 774 %      void LeastSquaresAddTerms(double **matrix,double **vectors,
 775 %        const double *terms,const double *results,const size_t rank,
 776 %        const size_t number_vectors);
 777 %
 778 %  A description of each parameter follows:
 779 %
 780 %    o matrix: the square matrix to add given terms/results to.
 781 %
 782 %    o vectors: the result vectors to add terms/results to.
 783 %
 784 %    o terms: the pre-calculated terms (without the unknown coefficent
 785 %             weights) that forms the equation being added.
 786 %
 787 %    o results: the result(s) that should be generated from the given terms
 788 %               weighted by the yet-to-be-solved coefficents.
 789 %
 790 %    o rank: the rank or size of the dimensions of the square matrix.
 791 %            Also the length of vectors, and number of terms being added.
 792 %
 793 %    o number_vectors: Number of result vectors, and number or results being
 794 %      added.  Also represents the number of separable systems of equations
 795 %      that is being solved.
 796 %
 797 %  Example of use...
 798 %
 799 %     2 dimensional Affine Equations (which are separable)
 800 %         c0*x + c2*y + c4*1 => u
 801 %         c1*x + c3*y + c5*1 => v
 802 %
 803 %     double **matrix = AcquireMagickMatrix(3UL,3UL);
 804 %     double **vectors = AcquireMagickMatrix(2UL,3UL);
 805 %     double terms[3], results[2];
 806 %     ...
 807 %     for each given x,y -> u,v
 808 %        terms[0] = x;
 809 %        terms[1] = y;
 810 %        terms[2] = 1;
 811 %        results[0] = u;
 812 %        results[1] = v;
 813 %        LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL);
 814 %     ...
 815 %     if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) {
 816 %       c0 = vectors[0][0];
 817 %       c2 = vectors[0][1];
 818 %       c4 = vectors[0][2];
 819 %       c1 = vectors[1][0];
 820 %       c3 = vectors[1][1];
 821 %       c5 = vectors[1][2];
 822 %     }
 823 %     else
 824 %       printf("Matrix unsolvable\n);
 825 %     RelinquishMagickMatrix(matrix,3UL);
 826 %     RelinquishMagickMatrix(vectors,2UL);
 827 %
 828 */
 829 MagickPrivate void LeastSquaresAddTerms(double **matrix,double **vectors,
 830   const double *terms,const double *results,const size_t rank,
 831   const size_t number_vectors)
 832 {
 833   register ssize_t
 834     i,
 835     j;
 836
 837   for (j=0; j < (ssize_t) rank; j++)
 838   {
 839     for (i=0; i < (ssize_t) rank; i++)
 840       matrix[i][j]+=terms[i]*terms[j];
 841     for (i=0; i < (ssize_t) number_vectors; i++)
 842       vectors[i][j]+=results[i]*terms[j];
 843   }
 844 }
 845 \f
 846 /*
 847 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 848 %                                                                             %
 849 %                                                                             %
 850 %                                                                             %
 851 %   M a t r i x T o I m a g e                                                 %
 852 %                                                                             %
 853 %                                                                             %
 854 %                                                                             %
 855 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 856 %
 857 %  MatrixToImage() returns a matrix as an image.  The matrix elements must be
 858 %  of type double otherwise nonsense is returned.
 859 %
 860 %  The format of the MatrixToImage method is:
 861 %
 862 %      Image *MatrixToImage(const MatrixInfo *matrix_info,
 863 %        ExceptionInfo *exception)
 864 %
 865 %  A description of each parameter follows:
 866 %
 867 %    o matrix_info: the matrix.
 868 %
 869 %    o exception: return any errors or warnings in this structure.
 870 %
 871 */
 872 MagickExport Image *MatrixToImage(const MatrixInfo *matrix_info,
 873   ExceptionInfo *exception)
 874 {
 875   CacheView
 876     *image_view;
 877
 878   double
 879     max_value,
 880     min_value,
 881     scale_factor,
 882     value;
 883
 884   Image
 885     *image;
 886
 887   MagickBooleanType
 888     status;
 889
 890   ssize_t
 891     y;
 892
 893   assert(matrix_info != (const MatrixInfo *) NULL);
 894   assert(matrix_info->signature == MagickCoreSignature);
 895   assert(exception != (ExceptionInfo *) NULL);
 896   assert(exception->signature == MagickCoreSignature);
 897   if (matrix_info->stride < sizeof(double))
 898     return((Image *) NULL);
 899   /*
 900     Determine range of matrix.
 901   */
 902   (void) GetMatrixElement(matrix_info,0,0,&value);
 903   min_value=value;
 904   max_value=value;
 905   for (y=0; y < (ssize_t) matrix_info->rows; y++)
 906   {
 907     register ssize_t
 908       x;
 909
 910     for (x=0; x < (ssize_t) matrix_info->columns; x++)
 911     {
 912       if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
 913         continue;
 914       if (value < min_value)
 915         min_value=value;
 916       else
 917         if (value > max_value)
 918           max_value=value;
 919     }
 920   }
 921   if ((min_value == 0.0) && (max_value == 0.0))
 922     scale_factor=0;
 923   else
 924     if (min_value == max_value)
 925       {
 926         scale_factor=(double) QuantumRange/min_value;
 927         min_value=0;
 928       }
 929     else
 930       scale_factor=(double) QuantumRange/(max_value-min_value);
 931   /*
 932     Convert matrix to image.
 933   */
 934   image=AcquireImage((ImageInfo *) NULL,exception);
 935   image->columns=matrix_info->columns;
 936   image->rows=matrix_info->rows;
 937   image->colorspace=GRAYColorspace;
 938   status=MagickTrue;
 939   image_view=AcquireAuthenticCacheView(image,exception);
 940 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 941   #pragma omp parallel for schedule(static,4) shared(status) \
 942     magick_number_threads(image,image,image->rows,1)
 943 #endif
 944   for (y=0; y < (ssize_t) image->rows; y++)
 945   {
 946     double
 947       value;
 948
 949     register Quantum
 950       *q;
 951
 952     register ssize_t
 953       x;
 954
 955     if (status == MagickFalse)
 956       continue;
 957     q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 958     if (q == (Quantum *) NULL)
 959       {
 960         status=MagickFalse;
 961         continue;
 962       }
 963     for (x=0; x < (ssize_t) image->columns; x++)
 964     {
 965       if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
 966         continue;
 967       value=scale_factor*(value-min_value);
 968       *q=ClampToQuantum(value);
 969       q+=GetPixelChannels(image);
 970     }
 971     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 972       status=MagickFalse;
 973   }
 974   image_view=DestroyCacheView(image_view);
 975   if (status == MagickFalse)
 976     image=DestroyImage(image);
 977   return(image);
 978 }
 979 \f
 980 /*
 981 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 982 %                                                                             %
 983 %                                                                             %
 984 %                                                                             %
 985 %   N u l l M a t r i x                                                       %
 986 %                                                                             %
 987 %                                                                             %
 988 %                                                                             %
 989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 990 %
 991 %  NullMatrix() sets all elements of the matrix to zero.
 992 %
 993 %  The format of the ResetMagickMemory method is:
 994 %
 995 %      MagickBooleanType *NullMatrix(MatrixInfo *matrix_info)
 996 %
 997 %  A description of each parameter follows:
 998 %
 999 %    o matrix_info: the matrix.
1000 %
1001 */
1002 MagickExport MagickBooleanType NullMatrix(MatrixInfo *matrix_info)
1003 {
1004   register ssize_t
1005     x;
1006
1007   ssize_t
1008     count,
1009     y;
1010
1011   unsigned char
1012     value;
1013
1014   assert(matrix_info != (const MatrixInfo *) NULL);
1015   assert(matrix_info->signature == MagickCoreSignature);
1016   if (matrix_info->type != DiskCache)
1017     {
1018       (void) ResetMagickMemory(matrix_info->elements,0,(size_t)
1019         matrix_info->length);
1020       return(MagickTrue);
1021     }
1022   value=0;
1023   (void) lseek(matrix_info->file,0,SEEK_SET);
1024   for (y=0; y < (ssize_t) matrix_info->rows; y++)
1025   {
1026     for (x=0; x < (ssize_t) matrix_info->length; x++)
1027     {
1028       count=write(matrix_info->file,&value,sizeof(value));
1029       if (count != (ssize_t) sizeof(value))
1030         break;
1031     }
1032     if (x < (ssize_t) matrix_info->length)
1033       break;
1034   }
1035   return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue);
1036 }
1037 \f
1038 /*
1039 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1040 %                                                                             %
1041 %                                                                             %
1042 %                                                                             %
1043 %   R e l i n q u i s h M a g i c k M a t r i x                               %
1044 %                                                                             %
1045 %                                                                             %
1046 %                                                                             %
1047 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1048 %
1049 %  RelinquishMagickMatrix() frees the previously acquired matrix (array of
1050 %  pointers to arrays of doubles).
1051 %
1052 %  The format of the RelinquishMagickMatrix method is:
1053 %
1054 %      double **RelinquishMagickMatrix(double **matrix,
1055 %        const size_t number_rows)
1056 %
1057 %  A description of each parameter follows:
1058 %
1059 %    o matrix: the matrix to relinquish
1060 %
1061 %    o number_rows: the first dimension of the acquired matrix (number of
1062 %      pointers)
1063 %
1064 */
1065 MagickExport double **RelinquishMagickMatrix(double **matrix,
1066   const size_t number_rows)
1067 {
1068   register ssize_t
1069     i;
1070
1071   if (matrix == (double **) NULL )
1072     return(matrix);
1073   for (i=0; i < (ssize_t) number_rows; i++)
1074     matrix[i]=(double *) RelinquishMagickMemory(matrix[i]);
1075   matrix=(double **) RelinquishMagickMemory(matrix);
1076   return(matrix);
1077 }
1078 \f
1079 /*
1080 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1081 %                                                                             %
1082 %                                                                             %
1083 %                                                                             %
1084 %   S e t M a t r i x E l e m e n t                                           %
1085 %                                                                             %
1086 %                                                                             %
1087 %                                                                             %
1088 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1089 %
1090 %  SetMatrixElement() sets the specifed element in the matrix.
1091 %
1092 %  The format of the SetMatrixElement method is:
1093 %
1094 %      MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
1095 %        const ssize_t x,const ssize_t y,void *value)
1096 %
1097 %  A description of each parameter follows:
1098 %
1099 %    o matrix_info: the matrix columns.
1100 %
1101 %    o x: the matrix x-offset.
1102 %
1103 %    o y: the matrix y-offset.
1104 %
1105 %    o value: set the matrix element to this value.
1106 %
1107 */
1108
1109 MagickExport MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
1110   const ssize_t x,const ssize_t y,const void *value)
1111 {
1112   MagickOffsetType
1113     count,
1114     i;
1115
1116   assert(matrix_info != (const MatrixInfo *) NULL);
1117   assert(matrix_info->signature == MagickCoreSignature);
1118   i=(MagickOffsetType) y*matrix_info->columns+x;
1119   if ((i < 0) ||
1120       ((MagickSizeType) (i*matrix_info->stride) >= matrix_info->length))
1121     return(MagickFalse);
1122   if (matrix_info->type != DiskCache)
1123     {
1124       (void) memcpy((unsigned char *) matrix_info->elements+i*
1125         matrix_info->stride,value,matrix_info->stride);
1126       return(MagickTrue);
1127     }
1128   count=WriteMatrixElements(matrix_info,i*matrix_info->stride,
1129     matrix_info->stride,(unsigned char *) value);
1130   if (count != (MagickOffsetType) matrix_info->stride)
1131     return(MagickFalse);
1132   return(MagickTrue);
1133 }