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