MagickCore  6.9.13-51
Convert, Edit, Or Compose Bitmap Images
matrix.c
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % M M AAA TTTTT RRRR IIIII X X %
7 % MM MM A A T R R I X X %
8 % M M M AAAAA T RRRR I X %
9 % M M A A T R R I X X %
10 % M M A A T R R IIIII X X %
11 % %
12 % %
13 % MagickCore Matrix Methods %
14 % %
15 % Software Design %
16 % Cristy %
17 % August 2007 %
18 % %
19 % %
20 % Copyright 1999 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/license/ %
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 "magick/studio.h"
43 #include "magick/blob.h"
44 #include "magick/blob-private.h"
45 #include "magick/exception.h"
46 #include "magick/exception-private.h"
47 #include "magick/image-private.h"
48 #include "magick/matrix.h"
49 #include "magick/memory_.h"
50 #include "magick/memory-private.h"
51 #include "magick/pixel-private.h"
52 #include "magick/resource_.h"
53 #include "magick/semaphore.h"
54 #include "magick/thread-private.h"
55 #include "magick/utility.h"
56 ␌
57 /*
58  Typedef declaration.
59 */
60 struct _MatrixInfo
61 {
62  CacheType
63  type;
64 
65  size_t
66  columns,
67  rows,
68  stride;
69 
70  MagickSizeType
71  length;
72 
73  MagickBooleanType
74  mapped,
75  synchronize;
76 
77  char
78  path[MaxTextExtent];
79 
80  int
81  file;
82 
83  void
84  *elements;
85 
86  SemaphoreInfo
87  *semaphore;
88 
89  size_t
90  signature;
91 };
92 ␌
93 /*
94 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
95 % %
96 % %
97 % %
98 % A c q u i r e M a t r i x I n f o %
99 % %
100 % %
101 % %
102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
103 %
104 % AcquireMatrixInfo() allocates the ImageInfo structure.
105 %
106 % The format of the AcquireMatrixInfo method is:
107 %
108 % MatrixInfo *AcquireMatrixInfo(const size_t columns,const size_t rows,
109 % const size_t stride,ExceptionInfo *exception)
110 %
111 % A description of each parameter follows:
112 %
113 % o columns: the matrix columns.
114 %
115 % o rows: the matrix rows.
116 %
117 % o stride: the matrix stride.
118 %
119 % o exception: return any errors or warnings in this structure.
120 %
121 */
122 
123 #if defined(SIGBUS)
124 static void MatrixSignalHandler(int magick_unused(status))
125 {
126  magick_unreferenced(status);
127  ThrowFatalException(CacheFatalError,"UnableToExtendMatrixCache");
128 }
129 #endif
130 
131 static inline MagickOffsetType WriteMatrixElements(
132  const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset,
133  const MagickSizeType length,const unsigned char *magick_restrict buffer)
134 {
135  MagickOffsetType
136  i;
137 
138  ssize_t
139  count;
140 
141 #if !defined(MAGICKCORE_HAVE_PWRITE)
142  LockSemaphoreInfo(matrix_info->semaphore);
143  if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
144  {
145  UnlockSemaphoreInfo(matrix_info->semaphore);
146  return((MagickOffsetType) -1);
147  }
148 #endif
149  count=0;
150  for (i=0; i < (MagickOffsetType) length; i+=count)
151  {
152 #if !defined(MAGICKCORE_HAVE_PWRITE)
153  count=write(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
154  (MagickSizeType) MagickMaxBufferExtent));
155 #else
156  count=pwrite(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
157  (MagickSizeType) MagickMaxBufferExtent),(off_t) (offset+i));
158 #endif
159  if (count <= 0)
160  {
161  count=0;
162  if (errno != EINTR)
163  break;
164  }
165  }
166 #if !defined(MAGICKCORE_HAVE_PWRITE)
167  UnlockSemaphoreInfo(matrix_info->semaphore);
168 #endif
169  return(i);
170 }
171 
172 static MagickBooleanType SetMatrixExtent(
173  MatrixInfo *magick_restrict matrix_info,MagickSizeType length)
174 {
175  MagickOffsetType
176  count,
177  extent,
178  offset;
179 
180  if (length != (MagickSizeType) ((MagickOffsetType) length))
181  return(MagickFalse);
182  offset=(MagickOffsetType) lseek(matrix_info->file,0,SEEK_END);
183  if (offset < 0)
184  return(MagickFalse);
185  if ((MagickSizeType) offset >= length)
186  return(MagickTrue);
187  extent=(MagickOffsetType) length-1;
188  count=WriteMatrixElements(matrix_info,extent,1,(const unsigned char *) "");
189 #if defined(MAGICKCORE_HAVE_POSIX_FALLOCATE)
190  if (matrix_info->synchronize != MagickFalse)
191  (void) posix_fallocate(matrix_info->file,offset+1,extent-offset);
192 #endif
193 #if defined(SIGBUS)
194  (void) signal(SIGBUS,MatrixSignalHandler);
195 #endif
196  return(count != (MagickOffsetType) 1 ? MagickFalse : MagickTrue);
197 }
198 
199 MagickExport MatrixInfo *AcquireMatrixInfo(const size_t columns,
200  const size_t rows,const size_t stride,ExceptionInfo *exception)
201 {
202  char
203  *synchronize;
204 
205  MagickBooleanType
206  status;
207 
208  MatrixInfo
209  *matrix_info;
210 
211  matrix_info=(MatrixInfo *) AcquireMagickMemory(sizeof(*matrix_info));
212  if (matrix_info == (MatrixInfo *) NULL)
213  return((MatrixInfo *) NULL);
214  (void) memset(matrix_info,0,sizeof(*matrix_info));
215  matrix_info->signature=MagickCoreSignature;
216  matrix_info->columns=columns;
217  matrix_info->rows=rows;
218  matrix_info->stride=stride;
219  matrix_info->semaphore=AllocateSemaphoreInfo();
220  synchronize=GetEnvironmentValue("MAGICK_SYNCHRONIZE");
221  if (synchronize != (const char *) NULL)
222  {
223  matrix_info->synchronize=IsStringTrue(synchronize);
224  synchronize=DestroyString(synchronize);
225  }
226  matrix_info->length=(MagickSizeType) columns*rows*stride;
227  if (matrix_info->columns != (size_t) (matrix_info->length/rows/stride))
228  {
229  (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
230  "CacheResourcesExhausted","`%s'","matrix cache");
231  return(DestroyMatrixInfo(matrix_info));
232  }
233  matrix_info->type=MemoryCache;
234  status=AcquireMagickResource(AreaResource,matrix_info->length);
235  if ((status != MagickFalse) &&
236  (matrix_info->length == (MagickSizeType) ((size_t) matrix_info->length)) &&
237  ((size_t) matrix_info->length <= GetMaxMemoryRequest()))
238  {
239  status=AcquireMagickResource(MemoryResource,matrix_info->length);
240  if (status != MagickFalse)
241  {
242  matrix_info->mapped=MagickFalse;
243  matrix_info->elements=MagickAssumeAligned(AcquireAlignedMemory(1,
244  (size_t) 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 % simultaneous 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  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 */
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=RelinquishAlignedMemory(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  magick_fallthrough;
394  }
395  case DiskCache:
396  {
397  if (matrix_info->file != -1)
398  (void) close(matrix_info->file);
399  (void) RelinquishUniqueFileResource(matrix_info->path);
400  RelinquishMagickResource(DiskResource,matrix_info->length);
401  break;
402  }
403  default:
404  break;
405  }
406  UnlockSemaphoreInfo(matrix_info->semaphore);
407  DestroySemaphoreInfo(&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 augmented 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). Also represents
441 % the number terms that need to be solved.
442 %
443 % o number_vectors: Number of vectors columns, argumenting the above matrix.
444 % Usually 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 simultaneous 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', you
467 % 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 coordinates, 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 */
481 MagickExport MagickBooleanType GaussJordanElimination(double **matrix,
482  double **vectors,const size_t rank,const size_t number_vectors)
483 {
484 #define GaussJordanSwap(x,y) \
485 { \
486  double temp = (x); \
487  (x)=(y); \
488  (y)=temp; \
489 }
490 #define GaussJordanSwapLD(x,y) \
491 { \
492  long double temp = (x); \
493  (x)=(y); \
494  (y)=temp; \
495 }
496 #define ThrowGaussJordanException() \
497 { \
498  for (i=0; i < (ssize_t) rank; i++) \
499  hp_matrix[i]=(long double *) RelinquishMagickMemory(hp_matrix[i]); \
500  hp_matrix=(long double **) RelinquishMagickMemory(hp_matrix); \
501  if (pivots != (ssize_t *) NULL) \
502  pivots=(ssize_t *) RelinquishMagickMemory(pivots); \
503  if (rows != (ssize_t *) NULL) \
504  rows=(ssize_t *) RelinquishMagickMemory(rows); \
505  if (columns != (ssize_t *) NULL) \
506  columns=(ssize_t *) RelinquishMagickMemory(columns); \
507  return(MagickFalse); \
508 }
509 
510  long double
511  **hp_matrix = (long double **) NULL,
512  scale;
513 
514  ssize_t
515  column,
516  *columns = (ssize_t *) NULL,
517  i,
518  j,
519  *pivots = (ssize_t *) NULL,
520  row,
521  *rows = (ssize_t *) NULL;
522 
523  /*
524  Allocate high precision matrix.
525  */
526  hp_matrix=(long double **) AcquireQuantumMemory(rank,sizeof(*hp_matrix));
527  if (hp_matrix == (long double **) NULL)
528  return(MagickFalse);
529  for (i=0; i < (ssize_t) rank; i++)
530  {
531  hp_matrix[i]=(long double *) AcquireQuantumMemory(rank,
532  sizeof(*hp_matrix[i]));
533  if (hp_matrix[i] == (long double *) NULL)
534  ThrowGaussJordanException();
535  for (j=0; j < (ssize_t) rank; j++)
536  hp_matrix[i][j]=(long double)matrix[i][j];
537  }
538  columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns));
539  rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows));
540  pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots));
541  if ((columns == (ssize_t *) NULL) || (rows == (ssize_t *) NULL) ||
542  (pivots == (ssize_t *) NULL))
543  ThrowGaussJordanException();
544  (void) memset(columns,0,rank*sizeof(*columns));
545  (void) memset(rows,0,rank*sizeof(*rows));
546  (void) memset(pivots,0,rank*sizeof(*pivots));
547  for (i=0; i < (ssize_t) rank; i++)
548  {
549  long double
550  max = 0.0;
551 
552  ssize_t
553  k;
554 
555  /*
556  Partial pivoting: find the largest absolute value in the unreduced
557  submatrix.
558  */
559  column=(-1);
560  row=(-1);
561  for (j=0; j < (ssize_t) rank; j++)
562  if (pivots[j] != 1)
563  for (k=0; k < (ssize_t) rank; k++)
564  if ((pivots[k] == 0) && (fabsl(hp_matrix[j][k]) > max))
565  {
566  max=fabsl(hp_matrix[j][k]);
567  row=j;
568  column=k;
569  }
570  if ((column == -1) || (row == -1) || (fabsl(max) < LDBL_MIN))
571  ThrowGaussJordanException(); /* Singular or nearly singular matrix */
572  pivots[column]++;
573  if (row != column)
574  {
575  for (k=0; k < (ssize_t) rank; k++)
576  GaussJordanSwapLD(hp_matrix[row][k],hp_matrix[column][k]);
577  for (k=0; k < (ssize_t) number_vectors; k++)
578  GaussJordanSwap(vectors[k][row],vectors[k][column]);
579  }
580  rows[i]=row;
581  columns[i]=column;
582  if (fabsl(hp_matrix[column][column]) < LDBL_MIN)
583  ThrowGaussJordanException(); /* Singular matrix */
584  scale=1.0L/hp_matrix[column][column];
585  hp_matrix[column][column]=1.0;
586  for (j=0; j < (ssize_t) rank; j++)
587  hp_matrix[column][j]*=scale;
588  for (j=0; j < (ssize_t) number_vectors; j++)
589  vectors[j][column]*=(double) scale;
590  for (j=0; j < (ssize_t) rank; j++)
591  if (j != column)
592  {
593  scale=hp_matrix[j][column];
594  hp_matrix[j][column]=0.0;
595  for (k=0; k < (ssize_t) rank; k++)
596  hp_matrix[j][k]-=scale*hp_matrix[column][k];
597  for (k=0; k < (ssize_t) number_vectors; k++)
598  vectors[k][j]-=(double)(scale*(long double) vectors[k][column]);
599  }
600  }
601  for (j=(ssize_t) rank-1; j >= 0; j--)
602  if (columns[j] != rows[j])
603  for (i=0; i < (ssize_t) rank; i++)
604  GaussJordanSwapLD(hp_matrix[i][columns[j]],hp_matrix[i][rows[j]]);
605  /*
606  Copy back the result to the original matrix.
607  */
608  for (i=0; i < (ssize_t) rank; i++)
609  for (j=0; j < (ssize_t) rank; j++)
610  matrix[i][j]=(double)hp_matrix[i][j];
611  /*
612  Free resources.
613  */
614  for (i=0; i < (ssize_t) rank; i++)
615  hp_matrix[i]=(long double *) RelinquishMagickMemory(hp_matrix[i]);
616  hp_matrix=(long double **) RelinquishMagickMemory(hp_matrix);
617  pivots=(ssize_t *) RelinquishMagickMemory(pivots);
618  rows=(ssize_t *) RelinquishMagickMemory(rows);
619  columns=(ssize_t *) RelinquishMagickMemory(columns);
620  return(MagickTrue);
621 }
622 ␌
623 /*
624 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
625 % %
626 % %
627 % %
628 % G e t M a t r i x C o l u m n s %
629 % %
630 % %
631 % %
632 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
633 %
634 % GetMatrixColumns() returns the number of columns in the matrix.
635 %
636 % The format of the GetMatrixColumns method is:
637 %
638 % size_t GetMatrixColumns(const MatrixInfo *matrix_info)
639 %
640 % A description of each parameter follows:
641 %
642 % o matrix_info: the matrix.
643 %
644 */
645 MagickExport size_t GetMatrixColumns(const MatrixInfo *matrix_info)
646 {
647  assert(matrix_info != (MatrixInfo *) NULL);
648  assert(matrix_info->signature == MagickCoreSignature);
649  return(matrix_info->columns);
650 }
651 ␌
652 /*
653 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
654 % %
655 % %
656 % %
657 % G e t M a t r i x E l e m e n t %
658 % %
659 % %
660 % %
661 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
662 %
663 % GetMatrixElement() returns the specified element in the matrix.
664 %
665 % The format of the GetMatrixElement method is:
666 %
667 % MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
668 % const ssize_t x,const ssize_t y,void *value)
669 %
670 % A description of each parameter follows:
671 %
672 % o matrix_info: the matrix columns.
673 %
674 % o x: the matrix x-offset.
675 %
676 % o y: the matrix y-offset.
677 %
678 % o value: return the matrix element in this buffer.
679 %
680 */
681 
682 static inline ssize_t EdgeX(const ssize_t x,const size_t columns)
683 {
684  if (x < 0L)
685  return(0L);
686  if (x >= (ssize_t) columns)
687  return((ssize_t) (columns-1));
688  return(x);
689 }
690 
691 static inline ssize_t EdgeY(const ssize_t y,const size_t rows)
692 {
693  if (y < 0L)
694  return(0L);
695  if (y >= (ssize_t) rows)
696  return((ssize_t) (rows-1));
697  return(y);
698 }
699 
700 static inline MagickOffsetType ReadMatrixElements(
701  const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset,
702  const MagickSizeType length,unsigned char *magick_restrict buffer)
703 {
704  MagickOffsetType
705  i;
706 
707  ssize_t
708  count;
709 
710 #if !defined(MAGICKCORE_HAVE_PREAD)
711  LockSemaphoreInfo(matrix_info->semaphore);
712  if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
713  {
714  UnlockSemaphoreInfo(matrix_info->semaphore);
715  return((MagickOffsetType) -1);
716  }
717 #endif
718  count=0;
719  for (i=0; i < (MagickOffsetType) length; i+=count)
720  {
721 #if !defined(MAGICKCORE_HAVE_PREAD)
722  count=read(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
723  (MagickSizeType) MagickMaxBufferExtent));
724 #else
725  count=pread(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
726  (MagickSizeType) MagickMaxBufferExtent),(off_t) (offset+i));
727 #endif
728  if (count <= 0)
729  {
730  count=0;
731  if (errno != EINTR)
732  break;
733  }
734  }
735 #if !defined(MAGICKCORE_HAVE_PREAD)
736  UnlockSemaphoreInfo(matrix_info->semaphore);
737 #endif
738  return(i);
739 }
740 
741 MagickExport MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
742  const ssize_t x,const ssize_t y,void *value)
743 {
744  MagickOffsetType
745  count,
746  i;
747 
748  assert(matrix_info != (const MatrixInfo *) NULL);
749  assert(matrix_info->signature == MagickCoreSignature);
750  i=(MagickOffsetType) EdgeY(y,matrix_info->rows)*matrix_info->columns+
751  EdgeX(x,matrix_info->columns);
752  if (matrix_info->type != DiskCache)
753  {
754  (void) memcpy(value,(unsigned char *) matrix_info->elements+i*
755  matrix_info->stride,matrix_info->stride);
756  return(MagickTrue);
757  }
758  count=ReadMatrixElements(matrix_info,i*matrix_info->stride,
759  matrix_info->stride,(unsigned char *) value);
760  if (count != (MagickOffsetType) matrix_info->stride)
761  return(MagickFalse);
762  return(MagickTrue);
763 }
764 ␌
765 /*
766 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
767 % %
768 % %
769 % %
770 % G e t M a t r i x R o w s %
771 % %
772 % %
773 % %
774 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
775 %
776 % GetMatrixRows() returns the number of rows in the matrix.
777 %
778 % The format of the GetMatrixRows method is:
779 %
780 % size_t GetMatrixRows(const MatrixInfo *matrix_info)
781 %
782 % A description of each parameter follows:
783 %
784 % o matrix_info: the matrix.
785 %
786 */
787 MagickExport size_t GetMatrixRows(const MatrixInfo *matrix_info)
788 {
789  assert(matrix_info != (const MatrixInfo *) NULL);
790  assert(matrix_info->signature == MagickCoreSignature);
791  return(matrix_info->rows);
792 }
793 ␌
794 /*
795 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
796 % %
797 % %
798 % %
799 % L e a s t S q u a r e s A d d T e r m s %
800 % %
801 % %
802 % %
803 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
804 %
805 % LeastSquaresAddTerms() adds one set of terms and associate results to the
806 % given matrix and vectors for solving using least-squares function fitting.
807 %
808 % The format of the AcquireMagickMatrix method is:
809 %
810 % void LeastSquaresAddTerms(double **matrix,double **vectors,
811 % const double *terms,const double *results,const size_t rank,
812 % const size_t number_vectors);
813 %
814 % A description of each parameter follows:
815 %
816 % o matrix: the square matrix to add given terms/results to.
817 %
818 % o vectors: the result vectors to add terms/results to.
819 %
820 % o terms: the pre-calculated terms (without the unknown coefficient
821 % weights) that forms the equation being added.
822 %
823 % o results: the result(s) that should be generated from the given terms
824 % weighted by the yet-to-be-solved coefficients.
825 %
826 % o rank: the rank or size of the dimensions of the square matrix.
827 % Also the length of vectors, and number of terms being added.
828 %
829 % o number_vectors: Number of result vectors, and number or results being
830 % added. Also represents the number of separable systems of equations
831 % that is being solved.
832 %
833 % Example of use...
834 %
835 % 2 dimensional Affine Equations (which are separable)
836 % c0*x + c2*y + c4*1 => u
837 % c1*x + c3*y + c5*1 => v
838 %
839 % double **matrix = AcquireMagickMatrix(3UL,3UL);
840 % double **vectors = AcquireMagickMatrix(2UL,3UL);
841 % double terms[3], results[2];
842 % ...
843 % for each given x,y -> u,v
844 % terms[0] = x;
845 % terms[1] = y;
846 % terms[2] = 1;
847 % results[0] = u;
848 % results[1] = v;
849 % LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL);
850 % ...
851 % if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) {
852 % c0 = vectors[0][0];
853 % c2 = vectors[0][1];
854 % c4 = vectors[0][2];
855 % c1 = vectors[1][0];
856 % c3 = vectors[1][1];
857 % c5 = vectors[1][2];
858 % }
859 % else
860 % printf("Matrix unsolvable\n);
861 % RelinquishMagickMatrix(matrix,3UL);
862 % RelinquishMagickMatrix(vectors,2UL);
863 %
864 */
865 MagickExport void LeastSquaresAddTerms(double **matrix,double **vectors,
866  const double *terms,const double *results,const size_t rank,
867  const size_t number_vectors)
868 {
869  ssize_t
870  i,
871  j;
872 
873  for (j=0; j < (ssize_t) rank; j++)
874  {
875  for (i=0; i < (ssize_t) rank; i++)
876  matrix[i][j]+=terms[i]*terms[j];
877  for (i=0; i < (ssize_t) number_vectors; i++)
878  vectors[i][j]+=results[i]*terms[j];
879  }
880 }
881 ␌
882 /*
883 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
884 % %
885 % %
886 % %
887 % M a t r i x T o I m a g e %
888 % %
889 % %
890 % %
891 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
892 %
893 % MatrixToImage() returns a matrix as an image. The matrix elements must be
894 % of type double otherwise nonsense is returned.
895 %
896 % The format of the MatrixToImage method is:
897 %
898 % Image *MatrixToImage(const MatrixInfo *matrix_info,
899 % ExceptionInfo *exception)
900 %
901 % A description of each parameter follows:
902 %
903 % o matrix_info: the matrix.
904 %
905 % o exception: return any errors or warnings in this structure.
906 %
907 */
908 MagickExport Image *MatrixToImage(const MatrixInfo *matrix_info,
909  ExceptionInfo *exception)
910 {
911  CacheView
912  *image_view;
913 
914  double
915  max_value,
916  min_value,
917  scale_factor,
918  value;
919 
920  Image
921  *image;
922 
923  MagickBooleanType
924  status;
925 
926  ssize_t
927  y;
928 
929  assert(matrix_info != (const MatrixInfo *) NULL);
930  assert(matrix_info->signature == MagickCoreSignature);
931  assert(exception != (ExceptionInfo *) NULL);
932  assert(exception->signature == MagickCoreSignature);
933  if (matrix_info->stride < sizeof(double))
934  return((Image *) NULL);
935  /*
936  Determine range of matrix.
937  */
938  (void) GetMatrixElement(matrix_info,0,0,&value);
939  min_value=value;
940  max_value=value;
941  for (y=0; y < (ssize_t) matrix_info->rows; y++)
942  {
943  ssize_t
944  x;
945 
946  for (x=0; x < (ssize_t) matrix_info->columns; x++)
947  {
948  if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
949  continue;
950  if (value < min_value)
951  min_value=value;
952  else
953  if (value > max_value)
954  max_value=value;
955  }
956  }
957  if ((min_value == 0.0) && (max_value == 0.0))
958  scale_factor=0;
959  else
960  if (min_value == max_value)
961  {
962  scale_factor=(double) QuantumRange/min_value;
963  min_value=0;
964  }
965  else
966  scale_factor=(double) QuantumRange/(max_value-min_value);
967  /*
968  Convert matrix to image.
969  */
970  image=AcquireImage((ImageInfo *) NULL);
971  image->columns=matrix_info->columns;
972  image->rows=matrix_info->rows;
973  image->colorspace=GRAYColorspace;
974  status=MagickTrue;
975  image_view=AcquireAuthenticCacheView(image,exception);
976 #if defined(MAGICKCORE_OPENMP_SUPPORT)
977  #pragma omp parallel for schedule(static) shared(status) \
978  magick_number_threads(image,image,image->rows,2)
979 #endif
980  for (y=0; y < (ssize_t) image->rows; y++)
981  {
982  double
983  value;
984 
985  PixelPacket
986  *q;
987 
988  ssize_t
989  x;
990 
991  if (status == MagickFalse)
992  continue;
993  q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
994  if (q == (PixelPacket *) NULL)
995  {
996  status=MagickFalse;
997  continue;
998  }
999  for (x=0; x < (ssize_t) image->columns; x++)
1000  {
1001  if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
1002  continue;
1003  value=scale_factor*(value-min_value);
1004  q->red=ClampToQuantum(value);
1005  q->green=q->red;
1006  q->blue=q->red;
1007  q++;
1008  }
1009  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1010  status=MagickFalse;
1011  }
1012  image_view=DestroyCacheView(image_view);
1013  if (status == MagickFalse)
1014  image=DestroyImage(image);
1015  return(image);
1016 }
1017 ␌
1018 /*
1019 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1020 % %
1021 % %
1022 % %
1023 % N u l l M a t r i x %
1024 % %
1025 % %
1026 % %
1027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1028 %
1029 % NullMatrix() sets all elements of the matrix to zero.
1030 %
1031 % The format of the memset method is:
1032 %
1033 % MagickBooleanType *NullMatrix(MatrixInfo *matrix_info)
1034 %
1035 % A description of each parameter follows:
1036 %
1037 % o matrix_info: the matrix.
1038 %
1039 */
1040 MagickExport MagickBooleanType NullMatrix(MatrixInfo *matrix_info)
1041 {
1042  ssize_t
1043  x;
1044 
1045  ssize_t
1046  count,
1047  y;
1048 
1049  unsigned char
1050  value;
1051 
1052  assert(matrix_info != (const MatrixInfo *) NULL);
1053  assert(matrix_info->signature == MagickCoreSignature);
1054  if (matrix_info->type != DiskCache)
1055  {
1056  (void) memset(matrix_info->elements,0,(size_t)
1057  matrix_info->length);
1058  return(MagickTrue);
1059  }
1060  value=0;
1061  (void) lseek(matrix_info->file,0,SEEK_SET);
1062  for (y=0; y < (ssize_t) matrix_info->rows; y++)
1063  {
1064  for (x=0; x < (ssize_t) matrix_info->length; x++)
1065  {
1066  count=write(matrix_info->file,&value,sizeof(value));
1067  if (count != (ssize_t) sizeof(value))
1068  break;
1069  }
1070  if (x < (ssize_t) matrix_info->length)
1071  break;
1072  }
1073  return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue);
1074 }
1075 ␌
1076 /*
1077 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1078 % %
1079 % %
1080 % %
1081 % R e l i n q u i s h M a g i c k M a t r i x %
1082 % %
1083 % %
1084 % %
1085 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1086 %
1087 % RelinquishMagickMatrix() frees the previously acquired matrix (array of
1088 % pointers to arrays of doubles).
1089 %
1090 % The format of the RelinquishMagickMatrix method is:
1091 %
1092 % double **RelinquishMagickMatrix(double **matrix,
1093 % const size_t number_rows)
1094 %
1095 % A description of each parameter follows:
1096 %
1097 % o matrix: the matrix to relinquish
1098 %
1099 % o number_rows: the first dimension of the acquired matrix (number of
1100 % pointers)
1101 %
1102 */
1103 MagickExport double **RelinquishMagickMatrix(double **matrix,
1104  const size_t number_rows)
1105 {
1106  ssize_t
1107  i;
1108 
1109  if (matrix == (double **) NULL )
1110  return(matrix);
1111  for (i=0; i < (ssize_t) number_rows; i++)
1112  matrix[i]=(double *) RelinquishMagickMemory(matrix[i]);
1113  matrix=(double **) RelinquishMagickMemory(matrix);
1114  return(matrix);
1115 }
1116 ␌
1117 /*
1118 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1119 % %
1120 % %
1121 % %
1122 % S e t M a t r i x E l e m e n t %
1123 % %
1124 % %
1125 % %
1126 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1127 %
1128 % SetMatrixElement() sets the specified element in the matrix.
1129 %
1130 % The format of the SetMatrixElement method is:
1131 %
1132 % MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
1133 % const ssize_t x,const ssize_t y,void *value)
1134 %
1135 % A description of each parameter follows:
1136 %
1137 % o matrix_info: the matrix columns.
1138 %
1139 % o x: the matrix x-offset.
1140 %
1141 % o y: the matrix y-offset.
1142 %
1143 % o value: set the matrix element to this value.
1144 %
1145 */
1146 
1147 MagickExport MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
1148  const ssize_t x,const ssize_t y,const void *value)
1149 {
1150  MagickOffsetType
1151  count,
1152  i;
1153 
1154  assert(matrix_info != (const MatrixInfo *) NULL);
1155  assert(matrix_info->signature == MagickCoreSignature);
1156  i=(MagickOffsetType) y*matrix_info->columns+x;
1157  if ((i < 0) ||
1158  ((MagickSizeType) (i*matrix_info->stride) >= matrix_info->length))
1159  return(MagickFalse);
1160  if (matrix_info->type != DiskCache)
1161  {
1162  (void) memcpy((unsigned char *) matrix_info->elements+i*
1163  matrix_info->stride,value,matrix_info->stride);
1164  return(MagickTrue);
1165  }
1166  count=WriteMatrixElements(matrix_info,i*matrix_info->stride,
1167  matrix_info->stride,(unsigned char *) value);
1168  if (count != (MagickOffsetType) matrix_info->stride)
1169  return(MagickFalse);
1170  return(MagickTrue);
1171 }
_CacheView
Definition: cache-view.c:66
_ExceptionInfo
Definition: exception.h:103
_ImageInfo
Definition: image.h:346
_Image
Definition: image.h:134
_MatrixInfo
Definition: matrix.c:61
_PixelPacket
Definition: pixel.h:132