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