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