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