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