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