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