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