MagickCore  7.0.7
Convert, Edit, Or Compose Bitmap Images
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-2018 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
22 % %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
25 % %
26 % https://www.imagemagick.org/script/license.php %
27 % %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
33 % %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 */
38 
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  ThrowFatalException(CacheFatalError,"UnableToExtendMatrixCache");
129 }
130 #endif
131 
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 
175 {
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 
201  const size_t rows,const size_t stride,ExceptionInfo *exception)
202 {
203  char
204  *synchronize;
205 
207  status;
208 
209  MatrixInfo
210  *matrix_info;
211 
212  matrix_info=(MatrixInfo *) AcquireMagickMemory(sizeof(*matrix_info));
213  if (matrix_info == (MatrixInfo *) NULL)
214  return((MatrixInfo *) NULL);
215  (void) ResetMagickMemory(matrix_info,0,sizeof(*matrix_info));
216  matrix_info->signature=MagickCoreSignature;
217  matrix_info->columns=columns;
218  matrix_info->rows=rows;
219  matrix_info->stride=stride;
220  matrix_info->semaphore=AcquireSemaphoreInfo();
221  synchronize=GetEnvironmentValue("MAGICK_SYNCHRONIZE");
222  if (synchronize != (const char *) NULL)
223  {
224  matrix_info->synchronize=IsStringTrue(synchronize);
225  synchronize=DestroyString(synchronize);
226  }
227  matrix_info->length=(MagickSizeType) columns*rows*stride;
228  if (matrix_info->columns != (size_t) (matrix_info->length/rows/stride))
229  {
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)
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  {
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
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 */
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 */
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  }
386  break;
387  }
388  case MapCache:
389  {
390  (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length);
391  matrix_info->elements=NULL;
393  }
394  case DiskCache:
395  {
396  if (matrix_info->file != -1)
397  (void) close(matrix_info->file);
398  (void) RelinquishUniqueFileResource(matrix_info->path);
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 */
481  double **vectors,const size_t rank,const size_t number_vectors)
482 {
483 #define GaussJordanSwap(x,y) \
484 { \
485  if ((x) != (y)) \
486  { \
487  (x)+=(y); \
488  (y)=(x)-(y); \
489  (x)=(x)-(y); \
490  } \
491 }
492 
493  double
494  max,
495  scale;
496 
497  register ssize_t
498  i,
499  j,
500  k;
501 
502  ssize_t
503  column,
504  *columns,
505  *pivots,
506  row,
507  *rows;
508 
509  columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns));
510  rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows));
511  pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots));
512  if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) ||
513  (pivots == (ssize_t *) NULL))
514  {
515  if (pivots != (ssize_t *) NULL)
516  pivots=(ssize_t *) RelinquishMagickMemory(pivots);
517  if (columns != (ssize_t *) NULL)
518  columns=(ssize_t *) RelinquishMagickMemory(columns);
519  if (rows != (ssize_t *) NULL)
520  rows=(ssize_t *) RelinquishMagickMemory(rows);
521  return(MagickFalse);
522  }
523  (void) ResetMagickMemory(columns,0,rank*sizeof(*columns));
524  (void) ResetMagickMemory(rows,0,rank*sizeof(*rows));
525  (void) ResetMagickMemory(pivots,0,rank*sizeof(*pivots));
526  column=0;
527  row=0;
528  for (i=0; i < (ssize_t) rank; i++)
529  {
530  max=0.0;
531  for (j=0; j < (ssize_t) rank; j++)
532  if (pivots[j] != 1)
533  {
534  for (k=0; k < (ssize_t) rank; k++)
535  if (pivots[k] != 0)
536  {
537  if (pivots[k] > 1)
538  return(MagickFalse);
539  }
540  else
541  if (fabs(matrix[j][k]) >= max)
542  {
543  max=fabs(matrix[j][k]);
544  row=j;
545  column=k;
546  }
547  }
548  pivots[column]++;
549  if (row != column)
550  {
551  for (k=0; k < (ssize_t) rank; k++)
552  GaussJordanSwap(matrix[row][k],matrix[column][k]);
553  for (k=0; k < (ssize_t) number_vectors; k++)
554  GaussJordanSwap(vectors[k][row],vectors[k][column]);
555  }
556  rows[i]=row;
557  columns[i]=column;
558  if (matrix[column][column] == 0.0)
559  return(MagickFalse); /* sigularity */
560  scale=PerceptibleReciprocal(matrix[column][column]);
561  matrix[column][column]=1.0;
562  for (j=0; j < (ssize_t) rank; j++)
563  matrix[column][j]*=scale;
564  for (j=0; j < (ssize_t) number_vectors; j++)
565  vectors[j][column]*=scale;
566  for (j=0; j < (ssize_t) rank; j++)
567  if (j != column)
568  {
569  scale=matrix[j][column];
570  matrix[j][column]=0.0;
571  for (k=0; k < (ssize_t) rank; k++)
572  matrix[j][k]-=scale*matrix[column][k];
573  for (k=0; k < (ssize_t) number_vectors; k++)
574  vectors[k][j]-=scale*vectors[k][column];
575  }
576  }
577  for (j=(ssize_t) rank-1; j >= 0; j--)
578  if (columns[j] != rows[j])
579  for (i=0; i < (ssize_t) rank; i++)
580  GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]);
581  pivots=(ssize_t *) RelinquishMagickMemory(pivots);
582  rows=(ssize_t *) RelinquishMagickMemory(rows);
583  columns=(ssize_t *) RelinquishMagickMemory(columns);
584  return(MagickTrue);
585 }
586 
587 /*
588 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
589 % %
590 % %
591 % %
592 % G e t M a t r i x C o l u m n s %
593 % %
594 % %
595 % %
596 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
597 %
598 % GetMatrixColumns() returns the number of columns in the matrix.
599 %
600 % The format of the GetMatrixColumns method is:
601 %
602 % size_t GetMatrixColumns(const MatrixInfo *matrix_info)
603 %
604 % A description of each parameter follows:
605 %
606 % o matrix_info: the matrix.
607 %
608 */
609 MagickExport size_t GetMatrixColumns(const MatrixInfo *matrix_info)
610 {
611  assert(matrix_info != (MatrixInfo *) NULL);
612  assert(matrix_info->signature == MagickCoreSignature);
613  return(matrix_info->columns);
614 }
615 
616 /*
617 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
618 % %
619 % %
620 % %
621 % G e t M a t r i x E l e m e n t %
622 % %
623 % %
624 % %
625 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
626 %
627 % GetMatrixElement() returns the specifed element in the matrix.
628 %
629 % The format of the GetMatrixElement method is:
630 %
631 % MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
632 % const ssize_t x,const ssize_t y,void *value)
633 %
634 % A description of each parameter follows:
635 %
636 % o matrix_info: the matrix columns.
637 %
638 % o x: the matrix x-offset.
639 %
640 % o y: the matrix y-offset.
641 %
642 % o value: return the matrix element in this buffer.
643 %
644 */
645 
646 static inline ssize_t EdgeX(const ssize_t x,const size_t columns)
647 {
648  if (x < 0L)
649  return(0L);
650  if (x >= (ssize_t) columns)
651  return((ssize_t) (columns-1));
652  return(x);
653 }
654 
655 static inline ssize_t EdgeY(const ssize_t y,const size_t rows)
656 {
657  if (y < 0L)
658  return(0L);
659  if (y >= (ssize_t) rows)
660  return((ssize_t) (rows-1));
661  return(y);
662 }
663 
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 
706  const ssize_t x,const ssize_t y,void *value)
707 {
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 */
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 */
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 */
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 
888  status;
889 
890  ssize_t
891  y;
892 
893  assert(matrix_info != (const MatrixInfo *) NULL);
894  assert(matrix_info->signature == MagickCoreSignature);
895  assert(exception != (ExceptionInfo *) NULL);
896  assert(exception->signature == MagickCoreSignature);
897  if (matrix_info->stride < sizeof(double))
898  return((Image *) NULL);
899  /*
900  Determine range of matrix.
901  */
902  (void) GetMatrixElement(matrix_info,0,0,&value);
903  min_value=value;
904  max_value=value;
905  for (y=0; y < (ssize_t) matrix_info->rows; y++)
906  {
907  register ssize_t
908  x;
909 
910  for (x=0; x < (ssize_t) matrix_info->columns; x++)
911  {
912  if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
913  continue;
914  if (value < min_value)
915  min_value=value;
916  else
917  if (value > max_value)
918  max_value=value;
919  }
920  }
921  if ((min_value == 0.0) && (max_value == 0.0))
922  scale_factor=0;
923  else
924  if (min_value == max_value)
925  {
926  scale_factor=(double) QuantumRange/min_value;
927  min_value=0;
928  }
929  else
930  scale_factor=(double) QuantumRange/(max_value-min_value);
931  /*
932  Convert matrix to image.
933  */
934  image=AcquireImage((ImageInfo *) NULL,exception);
935  image->columns=matrix_info->columns;
936  image->rows=matrix_info->rows;
937  image->colorspace=GRAYColorspace;
938  status=MagickTrue;
939  image_view=AcquireAuthenticCacheView(image,exception);
940 #if defined(MAGICKCORE_OPENMP_SUPPORT)
941  #pragma omp parallel for schedule(static,4) shared(status) \
942  magick_number_threads(image,image,image->rows,1)
943 #endif
944  for (y=0; y < (ssize_t) image->rows; y++)
945  {
946  double
947  value;
948 
949  register Quantum
950  *q;
951 
952  register ssize_t
953  x;
954 
955  if (status == MagickFalse)
956  continue;
957  q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
958  if (q == (Quantum *) NULL)
959  {
960  status=MagickFalse;
961  continue;
962  }
963  for (x=0; x < (ssize_t) image->columns; x++)
964  {
965  if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
966  continue;
967  value=scale_factor*(value-min_value);
968  *q=ClampToQuantum(value);
969  q+=GetPixelChannels(image);
970  }
971  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
972  status=MagickFalse;
973  }
974  image_view=DestroyCacheView(image_view);
975  if (status == MagickFalse)
976  image=DestroyImage(image);
977  return(image);
978 }
979 
980 /*
981 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
982 % %
983 % %
984 % %
985 % N u l l M a t r i x %
986 % %
987 % %
988 % %
989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
990 %
991 % NullMatrix() sets all elements of the matrix to zero.
992 %
993 % The format of the ResetMagickMemory method is:
994 %
995 % MagickBooleanType *NullMatrix(MatrixInfo *matrix_info)
996 %
997 % A description of each parameter follows:
998 %
999 % o matrix_info: the matrix.
1000 %
1001 */
1003 {
1004  register ssize_t
1005  x;
1006 
1007  ssize_t
1008  count,
1009  y;
1010 
1011  unsigned char
1012  value;
1013 
1014  assert(matrix_info != (const MatrixInfo *) NULL);
1015  assert(matrix_info->signature == MagickCoreSignature);
1016  if (matrix_info->type != DiskCache)
1017  {
1018  (void) ResetMagickMemory(matrix_info->elements,0,(size_t)
1019  matrix_info->length);
1020  return(MagickTrue);
1021  }
1022  value=0;
1023  (void) lseek(matrix_info->file,0,SEEK_SET);
1024  for (y=0; y < (ssize_t) matrix_info->rows; y++)
1025  {
1026  for (x=0; x < (ssize_t) matrix_info->length; x++)
1027  {
1028  count=write(matrix_info->file,&value,sizeof(value));
1029  if (count != (ssize_t) sizeof(value))
1030  break;
1031  }
1032  if (x < (ssize_t) matrix_info->length)
1033  break;
1034  }
1035  return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue);
1036 }
1037 
1038 /*
1039 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1040 % %
1041 % %
1042 % %
1043 % R e l i n q u i s h M a g i c k M a t r i x %
1044 % %
1045 % %
1046 % %
1047 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1048 %
1049 % RelinquishMagickMatrix() frees the previously acquired matrix (array of
1050 % pointers to arrays of doubles).
1051 %
1052 % The format of the RelinquishMagickMatrix method is:
1053 %
1054 % double **RelinquishMagickMatrix(double **matrix,
1055 % const size_t number_rows)
1056 %
1057 % A description of each parameter follows:
1058 %
1059 % o matrix: the matrix to relinquish
1060 %
1061 % o number_rows: the first dimension of the acquired matrix (number of
1062 % pointers)
1063 %
1064 */
1065 MagickExport double **RelinquishMagickMatrix(double **matrix,
1066  const size_t number_rows)
1067 {
1068  register ssize_t
1069  i;
1070 
1071  if (matrix == (double **) NULL )
1072  return(matrix);
1073  for (i=0; i < (ssize_t) number_rows; i++)
1074  matrix[i]=(double *) RelinquishMagickMemory(matrix[i]);
1075  matrix=(double **) RelinquishMagickMemory(matrix);
1076  return(matrix);
1077 }
1078 
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 
1110  const ssize_t x,const ssize_t y,const void *value)
1111 {
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 }
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:829
MagickExport void UnlockSemaphoreInfo(SemaphoreInfo *semaphore_info)
Definition: semaphore.c:450
#define ThrowFatalException(severity, tag)
MagickExport MagickBooleanType NullMatrix(MatrixInfo *matrix_info)
Definition: matrix.c:1002
size_t signature
Definition: exception.h:123
MagickExport SemaphoreInfo * AcquireSemaphoreInfo(void)
Definition: semaphore.c:192
Definition: blob.h:31
MagickBooleanType synchronize
Definition: matrix.c:76
MagickExport size_t GetMatrixColumns(const MatrixInfo *matrix_info)
Definition: matrix.c:609
MagickExport void RelinquishMagickResource(const ResourceType type, const MagickSizeType size)
Definition: resource.c:932
MagickExport MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info, const ssize_t x, const ssize_t y, void *value)
Definition: matrix.c:705
MagickExport MagickBooleanType AcquireMagickResource(const ResourceType type, const MagickSizeType size)
Definition: resource.c:169
ssize_t MagickOffsetType
Definition: magick-type.h:127
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:1109
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:156
int file
Definition: matrix.c:83
static double PerceptibleReciprocal(const double x)
MagickExport int AcquireUniqueFileResource(char *path)
Definition: resource.c:549
MagickExport void * ResetMagickMemory(void *memory, int byte, const size_t size)
Definition: memory.c:1164
#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:529
MagickExport MagickBooleanType RelinquishUniqueFileResource(const char *path)
Definition: resource.c:1131
size_t MagickSizeType
Definition: magick-type.h:128
#define MagickPathExtent
MagickExport Image * MatrixToImage(const MatrixInfo *matrix_info, ExceptionInfo *exception)
Definition: matrix.c:872
MagickExport MagickBooleanType IsStringTrue(const char *value)
Definition: string.c:1464
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:200
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:1058
static MagickOffsetType ReadMatrixElements(const MatrixInfo *magick_restrict matrix_info, const MagickOffsetType offset, const MagickSizeType length, unsigned char *magick_restrict buffer)
Definition: matrix.c:664
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:480
MagickExport Image * AcquireImage(const ImageInfo *image_info, ExceptionInfo *exception)
Definition: image.c:154
MagickExport char * GetEnvironmentValue(const char *name)
Definition: string.c:1250
size_t columns
Definition: matrix.c:68
MagickExport double ** RelinquishMagickMatrix(double **matrix, const size_t number_rows)
Definition: matrix.c:1065
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:132
static Quantum ClampToQuantum(const MagickRealType value)
Definition: quantum.h:84
unsigned short Quantum
Definition: magick-type.h:82
static MagickBooleanType SetMatrixExtent(MatrixInfo *magick_restrict matrix_info, MagickSizeType length)
Definition: matrix.c:173
MagickExport double ** AcquireMagickMatrix(const size_t number_rows, const size_t size)
Definition: matrix.c:317
MagickExport char * DestroyString(char *string)
Definition: string.c:810
MagickExport void * AcquireMagickMemory(const size_t size)
Definition: memory.c:458
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:27
static ssize_t EdgeY(const ssize_t y, const size_t rows)
Definition: matrix.c:655
static ssize_t EdgeX(const ssize_t x, const size_t columns)
Definition: matrix.c:646
MagickExport void * RelinquishMagickMemory(void *memory)
Definition: memory.c:1038
MagickExport MagickBooleanType UnmapBlob(void *, const size_t)
Definition: blob.c:5220
#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:1182
ColorspaceType colorspace
Definition: image.h:157
CacheType type
Definition: matrix.c:65
#define QuantumRange
Definition: magick-type.h:83
MagickExport size_t GetMatrixRows(const MatrixInfo *matrix_info)
Definition: matrix.c:751
MagickExport MatrixInfo * DestroyMatrixInfo(MatrixInfo *matrix_info)
Definition: matrix.c:369
MagickExport void * MapBlob(int, const MapMode, const MagickOffsetType, const size_t)