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