• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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-2019 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://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 
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 
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 
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)
MatrixSignalHandler(int status)126 static void MatrixSignalHandler(int status)
127 {
128   ThrowFatalException(CacheFatalError,"UnableToExtendMatrixCache");
129 }
130 #endif
131 
WriteMatrixElements(const MatrixInfo * magick_restrict matrix_info,const MagickOffsetType offset,const MagickSizeType length,const unsigned char * magick_restrict buffer)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 
SetMatrixExtent(MatrixInfo * magick_restrict matrix_info,MagickSizeType length)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 
AcquireMatrixInfo(const size_t columns,const size_t rows,const size_t stride,ExceptionInfo * exception)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) memset(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 
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 */
AcquireMagickMatrix(const size_t number_rows,const size_t size)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 
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 */
DestroyMatrixInfo(MatrixInfo * matrix_info)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 
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 */
GaussJordanElimination(double ** matrix,double ** vectors,const size_t rank,const size_t number_vectors)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) memset(columns,0,rank*sizeof(*columns));
524   (void) memset(rows,0,rank*sizeof(*rows));
525   (void) memset(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 
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 */
GetMatrixColumns(const MatrixInfo * matrix_info)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 
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 
EdgeX(const ssize_t x,const size_t columns)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 
EdgeY(const ssize_t y,const size_t rows)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 
ReadMatrixElements(const MatrixInfo * magick_restrict matrix_info,const MagickOffsetType offset,const MagickSizeType length,unsigned char * magick_restrict buffer)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 
GetMatrixElement(const MatrixInfo * matrix_info,const ssize_t x,const ssize_t y,void * value)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 
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 */
GetMatrixRows(const MatrixInfo * matrix_info)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 
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 */
LeastSquaresAddTerms(double ** matrix,double ** vectors,const double * terms,const double * results,const size_t rank,const size_t number_vectors)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 
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 */
MatrixToImage(const MatrixInfo * matrix_info,ExceptionInfo * exception)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) 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 
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 memset 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 */
NullMatrix(MatrixInfo * matrix_info)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) memset(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 
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 */
RelinquishMagickMatrix(double ** matrix,const size_t number_rows)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 
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 
SetMatrixElement(const MatrixInfo * matrix_info,const ssize_t x,const ssize_t y,const void * value)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 }
1134