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