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