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