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