MagickCore 7.1.1-43
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*/
63struct _MatrixInfo
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
89 SemaphoreInfo
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
211 MatrixInfo
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 a u s s J o r d a n E l i m i n a t i o n %
419% %
420% %
421% %
422%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
423%
424% GaussJordanElimination() returns a matrix in reduced row echelon form,
425% while simultaneously reducing and thus solving the augmented results
426% matrix.
427%
428% See also http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
429%
430% The format of the GaussJordanElimination method is:
431%
432% MagickBooleanType GaussJordanElimination(double **matrix,
433% double **vectors,const size_t rank,const size_t number_vectors)
434%
435% A description of each parameter follows:
436%
437% o matrix: the matrix to be reduced, as an 'array of row pointers'.
438%
439% o vectors: the additional matrix argumenting the matrix for row reduction.
440% Producing an 'array of column vectors'.
441%
442% o rank: The size of the matrix (both rows and columns).
443% Also represents the number terms that need to be solved.
444%
445% o number_vectors: Number of vectors columns, argumenting the above matrix.
446% Usually 1, but can be more for more complex equation solving.
447%
448% Note that the 'matrix' is given as a 'array of row pointers' of rank size.
449% That is values can be assigned as matrix[row][column] where 'row' is
450% typically the equation, and 'column' is the term of the equation.
451% That is the matrix is in the form of a 'row first array'.
452%
453% However 'vectors' is a 'array of column pointers' which can have any number
454% of columns, with each column array the same 'rank' size as 'matrix'.
455%
456% This allows for simpler handling of the results, especially is only one
457% column 'vector' is all that is required to produce the desired solution.
458%
459% For example, the 'vectors' can consist of a pointer to a simple array of
460% doubles. when only one set of simultaneous equations is to be solved from
461% the given set of coefficient weighted terms.
462%
463% double **matrix = AcquireMagickMatrix(8UL,8UL);
464% double coefficients[8];
465% ...
466% GaussJordanElimination(matrix, &coefficients, 8UL, 1UL);
467%
468% However by specifying more 'columns' (as an 'array of vector columns',
469% you can use this function to solve a set of 'separable' equations.
470%
471% For example a distortion function where u = U(x,y) v = V(x,y)
472% And the functions U() and V() have separate coefficients, but are being
473% generated from a common x,y->u,v data set.
474%
475% Another example is generation of a color gradient from a set of colors at
476% specific coordinates, such as a list x,y -> r,g,b,a.
477%
478% You can also use the 'vectors' to generate an inverse of the given 'matrix'
479% though as a 'column first array' rather than a 'row first array'. For
480% details see http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
481%
482*/
483MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix,
484 double **vectors,const size_t rank,const size_t number_vectors)
485{
486#define GaussJordanSwap(x,y) \
487{ \
488 if ((x) != (y)) \
489 { \
490 (x)+=(y); \
491 (y)=(x)-(y); \
492 (x)=(x)-(y); \
493 } \
494}
495
496 double
497 max,
498 scale;
499
500 ssize_t
501 i,
502 j,
503 k;
504
505 ssize_t
506 column,
507 *columns,
508 *pivots,
509 row,
510 *rows;
511
512 columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns));
513 rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows));
514 pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots));
515 if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) ||
516 (pivots == (ssize_t *) NULL))
517 {
518 if (pivots != (ssize_t *) NULL)
519 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
520 if (columns != (ssize_t *) NULL)
521 columns=(ssize_t *) RelinquishMagickMemory(columns);
522 if (rows != (ssize_t *) NULL)
523 rows=(ssize_t *) RelinquishMagickMemory(rows);
524 return(MagickFalse);
525 }
526 (void) memset(columns,0,rank*sizeof(*columns));
527 (void) memset(rows,0,rank*sizeof(*rows));
528 (void) memset(pivots,0,rank*sizeof(*pivots));
529 column=0;
530 row=0;
531 for (i=0; i < (ssize_t) rank; i++)
532 {
533 max=0.0;
534 for (j=0; j < (ssize_t) rank; j++)
535 if (pivots[j] != 1)
536 {
537 for (k=0; k < (ssize_t) rank; k++)
538 if (pivots[k] != 0)
539 {
540 if (pivots[k] > 1)
541 return(MagickFalse);
542 }
543 else
544 if (fabs(matrix[j][k]) >= max)
545 {
546 max=fabs(matrix[j][k]);
547 row=j;
548 column=k;
549 }
550 }
551 pivots[column]++;
552 if (row != column)
553 {
554 for (k=0; k < (ssize_t) rank; k++)
555 GaussJordanSwap(matrix[row][k],matrix[column][k]);
556 for (k=0; k < (ssize_t) number_vectors; k++)
557 GaussJordanSwap(vectors[k][row],vectors[k][column]);
558 }
559 rows[i]=row;
560 columns[i]=column;
561 if (matrix[column][column] == 0.0)
562 return(MagickFalse); /* singularity */
563 scale=PerceptibleReciprocal(matrix[column][column]);
564 matrix[column][column]=1.0;
565 for (j=0; j < (ssize_t) rank; j++)
566 matrix[column][j]*=scale;
567 for (j=0; j < (ssize_t) number_vectors; j++)
568 vectors[j][column]*=scale;
569 for (j=0; j < (ssize_t) rank; j++)
570 if (j != column)
571 {
572 scale=matrix[j][column];
573 matrix[j][column]=0.0;
574 for (k=0; k < (ssize_t) rank; k++)
575 matrix[j][k]-=scale*matrix[column][k];
576 for (k=0; k < (ssize_t) number_vectors; k++)
577 vectors[k][j]-=scale*vectors[k][column];
578 }
579 }
580 for (j=(ssize_t) rank-1; j >= 0; j--)
581 if (columns[j] != rows[j])
582 for (i=0; i < (ssize_t) rank; i++)
583 GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]);
584 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
585 rows=(ssize_t *) RelinquishMagickMemory(rows);
586 columns=(ssize_t *) RelinquishMagickMemory(columns);
587 return(MagickTrue);
588}
589␌
590/*
591%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
592% %
593% %
594% %
595% G e t M a t r i x C o l u m n s %
596% %
597% %
598% %
599%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
600%
601% GetMatrixColumns() returns the number of columns in the matrix.
602%
603% The format of the GetMatrixColumns method is:
604%
605% size_t GetMatrixColumns(const MatrixInfo *matrix_info)
606%
607% A description of each parameter follows:
608%
609% o matrix_info: the matrix.
610%
611*/
612MagickExport size_t GetMatrixColumns(const MatrixInfo *matrix_info)
613{
614 assert(matrix_info != (MatrixInfo *) NULL);
615 assert(matrix_info->signature == MagickCoreSignature);
616 return(matrix_info->columns);
617}
618␌
619/*
620%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
621% %
622% %
623% %
624% G e t M a t r i x E l e m e n t %
625% %
626% %
627% %
628%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
629%
630% GetMatrixElement() returns the specified element in the matrix.
631%
632% The format of the GetMatrixElement method is:
633%
634% MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
635% const ssize_t x,const ssize_t y,void *value)
636%
637% A description of each parameter follows:
638%
639% o matrix_info: the matrix columns.
640%
641% o x: the matrix x-offset.
642%
643% o y: the matrix y-offset.
644%
645% o value: return the matrix element in this buffer.
646%
647*/
648
649static inline ssize_t EdgeX(const ssize_t x,const size_t columns)
650{
651 if (x < 0L)
652 return(0L);
653 if (x >= (ssize_t) columns)
654 return((ssize_t) (columns-1));
655 return(x);
656}
657
658static inline ssize_t EdgeY(const ssize_t y,const size_t rows)
659{
660 if (y < 0L)
661 return(0L);
662 if (y >= (ssize_t) rows)
663 return((ssize_t) (rows-1));
664 return(y);
665}
666
667static inline MagickOffsetType ReadMatrixElements(
668 const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset,
669 const MagickSizeType length,unsigned char *magick_restrict buffer)
670{
671 MagickOffsetType
672 i;
673
674 ssize_t
675 count;
676
677#if !defined(MAGICKCORE_HAVE_PREAD)
678 LockSemaphoreInfo(matrix_info->semaphore);
679 if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
680 {
681 UnlockSemaphoreInfo(matrix_info->semaphore);
682 return((MagickOffsetType) -1);
683 }
684#endif
685 count=0;
686 for (i=0; i < (MagickOffsetType) length; i+=count)
687 {
688#if !defined(MAGICKCORE_HAVE_PREAD)
689 count=read(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
690 (MagickSizeType) MagickMaxBufferExtent));
691#else
692 count=pread(matrix_info->file,buffer+i,(size_t) MagickMin(length-
693 (MagickSizeType) i,(MagickSizeType) MagickMaxBufferExtent),offset+i);
694#endif
695 if (count <= 0)
696 {
697 count=0;
698 if (errno != EINTR)
699 break;
700 }
701 }
702#if !defined(MAGICKCORE_HAVE_PREAD)
703 UnlockSemaphoreInfo(matrix_info->semaphore);
704#endif
705 return(i);
706}
707
708MagickExport MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
709 const ssize_t x,const ssize_t y,void *value)
710{
711 MagickOffsetType
712 count,
713 i;
714
715 assert(matrix_info != (const MatrixInfo *) NULL);
716 assert(matrix_info->signature == MagickCoreSignature);
717 i=EdgeY(y,matrix_info->rows)*(MagickOffsetType) matrix_info->columns+
718 EdgeX(x,matrix_info->columns);
719 if (matrix_info->type != DiskCache)
720 {
721 (void) memcpy(value,(unsigned char *) matrix_info->elements+i*
722 (MagickOffsetType) matrix_info->stride,matrix_info->stride);
723 return(MagickTrue);
724 }
725 count=ReadMatrixElements(matrix_info,i*(MagickOffsetType) matrix_info->stride,
726 matrix_info->stride,(unsigned char *) value);
727 if (count != (MagickOffsetType) matrix_info->stride)
728 return(MagickFalse);
729 return(MagickTrue);
730}
731␌
732/*
733%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
734% %
735% %
736% %
737% G e t M a t r i x R o w s %
738% %
739% %
740% %
741%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
742%
743% GetMatrixRows() returns the number of rows in the matrix.
744%
745% The format of the GetMatrixRows method is:
746%
747% size_t GetMatrixRows(const MatrixInfo *matrix_info)
748%
749% A description of each parameter follows:
750%
751% o matrix_info: the matrix.
752%
753*/
754MagickExport size_t GetMatrixRows(const MatrixInfo *matrix_info)
755{
756 assert(matrix_info != (const MatrixInfo *) NULL);
757 assert(matrix_info->signature == MagickCoreSignature);
758 return(matrix_info->rows);
759}
760␌
761/*
762%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
763% %
764% %
765% %
766+ L e a s t S q u a r e s A d d T e r m s %
767% %
768% %
769% %
770%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
771%
772% LeastSquaresAddTerms() adds one set of terms and associate results to the
773% given matrix and vectors for solving using least-squares function fitting.
774%
775% The format of the AcquireMagickMatrix method is:
776%
777% void LeastSquaresAddTerms(double **matrix,double **vectors,
778% const double *terms,const double *results,const size_t rank,
779% const size_t number_vectors);
780%
781% A description of each parameter follows:
782%
783% o matrix: the square matrix to add given terms/results to.
784%
785% o vectors: the result vectors to add terms/results to.
786%
787% o terms: the pre-calculated terms (without the unknown coefficient
788% weights) that forms the equation being added.
789%
790% o results: the result(s) that should be generated from the given terms
791% weighted by the yet-to-be-solved coefficients.
792%
793% o rank: the rank or size of the dimensions of the square matrix.
794% Also the length of vectors, and number of terms being added.
795%
796% o number_vectors: Number of result vectors, and number or results being
797% added. Also represents the number of separable systems of equations
798% that is being solved.
799%
800% Example of use...
801%
802% 2 dimensional Affine Equations (which are separable)
803% c0*x + c2*y + c4*1 => u
804% c1*x + c3*y + c5*1 => v
805%
806% double **matrix = AcquireMagickMatrix(3UL,3UL);
807% double **vectors = AcquireMagickMatrix(2UL,3UL);
808% double terms[3], results[2];
809% ...
810% for each given x,y -> u,v
811% terms[0] = x;
812% terms[1] = y;
813% terms[2] = 1;
814% results[0] = u;
815% results[1] = v;
816% LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL);
817% ...
818% if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) {
819% c0 = vectors[0][0];
820% c2 = vectors[0][1];
821% c4 = vectors[0][2];
822% c1 = vectors[1][0];
823% c3 = vectors[1][1];
824% c5 = vectors[1][2];
825% }
826% else
827% printf("Matrix unsolvable\n");
828% RelinquishMagickMatrix(matrix,3UL);
829% RelinquishMagickMatrix(vectors,2UL);
830%
831*/
832MagickPrivate void LeastSquaresAddTerms(double **matrix,double **vectors,
833 const double *terms,const double *results,const size_t rank,
834 const size_t number_vectors)
835{
836 ssize_t
837 i,
838 j;
839
840 for (j=0; j < (ssize_t) rank; j++)
841 {
842 for (i=0; i < (ssize_t) rank; i++)
843 matrix[i][j]+=terms[i]*terms[j];
844 for (i=0; i < (ssize_t) number_vectors; i++)
845 vectors[i][j]+=results[i]*terms[j];
846 }
847}
848␌
849/*
850%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
851% %
852% %
853% %
854% M a t r i x T o I m a g e %
855% %
856% %
857% %
858%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
859%
860% MatrixToImage() returns a matrix as an image. The matrix elements must be
861% of type double otherwise nonsense is returned.
862%
863% The format of the MatrixToImage method is:
864%
865% Image *MatrixToImage(const MatrixInfo *matrix_info,
866% ExceptionInfo *exception)
867%
868% A description of each parameter follows:
869%
870% o matrix_info: the matrix.
871%
872% o exception: return any errors or warnings in this structure.
873%
874*/
875MagickExport Image *MatrixToImage(const MatrixInfo *matrix_info,
876 ExceptionInfo *exception)
877{
878 CacheView
879 *image_view;
880
881 double
882 max_value,
883 min_value,
884 scale_factor;
885
886 Image
887 *image;
888
889 MagickBooleanType
890 status;
891
892 ssize_t
893 y;
894
895 assert(matrix_info != (const MatrixInfo *) NULL);
896 assert(matrix_info->signature == MagickCoreSignature);
897 assert(exception != (ExceptionInfo *) NULL);
898 assert(exception->signature == MagickCoreSignature);
899 if (matrix_info->stride < sizeof(double))
900 return((Image *) NULL);
901 /*
902 Determine range of matrix.
903 */
904 (void) GetMatrixElement(matrix_info,0,0,&min_value);
905 max_value=min_value;
906 for (y=0; y < (ssize_t) matrix_info->rows; y++)
907 {
908 ssize_t
909 x;
910
911 for (x=0; x < (ssize_t) matrix_info->columns; x++)
912 {
913 double
914 value;
915
916 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
917 continue;
918 if (value < min_value)
919 min_value=value;
920 else
921 if (value > max_value)
922 max_value=value;
923 }
924 }
925 if ((min_value == 0.0) && (max_value == 0.0))
926 scale_factor=0;
927 else
928 if (min_value == max_value)
929 {
930 scale_factor=(double) QuantumRange/min_value;
931 min_value=0;
932 }
933 else
934 scale_factor=(double) QuantumRange/(max_value-min_value);
935 /*
936 Convert matrix to image.
937 */
938 image=AcquireImage((ImageInfo *) NULL,exception);
939 image->columns=matrix_info->columns;
940 image->rows=matrix_info->rows;
941 image->colorspace=GRAYColorspace;
942 status=MagickTrue;
943 image_view=AcquireAuthenticCacheView(image,exception);
944#if defined(MAGICKCORE_OPENMP_SUPPORT)
945 #pragma omp parallel for schedule(static) shared(status) \
946 magick_number_threads(image,image,image->rows,2)
947#endif
948 for (y=0; y < (ssize_t) image->rows; y++)
949 {
950 double
951 value;
952
953 Quantum
954 *q;
955
956 ssize_t
957 x;
958
959 if (status == MagickFalse)
960 continue;
961 q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
962 if (q == (Quantum *) NULL)
963 {
964 status=MagickFalse;
965 continue;
966 }
967 for (x=0; x < (ssize_t) image->columns; x++)
968 {
969 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
970 continue;
971 value=scale_factor*(value-min_value);
972 *q=ClampToQuantum(value);
973 q+=(ptrdiff_t) GetPixelChannels(image);
974 }
975 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
976 status=MagickFalse;
977 }
978 image_view=DestroyCacheView(image_view);
979 if (status == MagickFalse)
980 image=DestroyImage(image);
981 return(image);
982}
983␌
984/*
985%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
986% %
987% %
988% %
989% N u l l M a t r i x %
990% %
991% %
992% %
993%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
994%
995% NullMatrix() sets all elements of the matrix to zero.
996%
997% The format of the memset method is:
998%
999% MagickBooleanType *NullMatrix(MatrixInfo *matrix_info)
1000%
1001% A description of each parameter follows:
1002%
1003% o matrix_info: the matrix.
1004%
1005*/
1006MagickExport MagickBooleanType NullMatrix(MatrixInfo *matrix_info)
1007{
1008 ssize_t
1009 x;
1010
1011 ssize_t
1012 count,
1013 y;
1014
1015 unsigned char
1016 value;
1017
1018 assert(matrix_info != (const MatrixInfo *) NULL);
1019 assert(matrix_info->signature == MagickCoreSignature);
1020 if (matrix_info->type != DiskCache)
1021 {
1022 (void) memset(matrix_info->elements,0,(size_t)
1023 matrix_info->length);
1024 return(MagickTrue);
1025 }
1026 value=0;
1027 (void) lseek(matrix_info->file,0,SEEK_SET);
1028 for (y=0; y < (ssize_t) matrix_info->rows; y++)
1029 {
1030 for (x=0; x < (ssize_t) matrix_info->length; x++)
1031 {
1032 count=write(matrix_info->file,&value,sizeof(value));
1033 if (count != (ssize_t) sizeof(value))
1034 break;
1035 }
1036 if (x < (ssize_t) matrix_info->length)
1037 break;
1038 }
1039 return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue);
1040}
1041␌
1042/*
1043%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1044% %
1045% %
1046% %
1047% R e l i n q u i s h M a g i c k M a t r i x %
1048% %
1049% %
1050% %
1051%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1052%
1053% RelinquishMagickMatrix() frees the previously acquired matrix (array of
1054% pointers to arrays of doubles).
1055%
1056% The format of the RelinquishMagickMatrix method is:
1057%
1058% double **RelinquishMagickMatrix(double **matrix,
1059% const size_t number_rows)
1060%
1061% A description of each parameter follows:
1062%
1063% o matrix: the matrix to relinquish
1064%
1065% o number_rows: the first dimension of the acquired matrix (number of
1066% pointers)
1067%
1068*/
1069MagickExport double **RelinquishMagickMatrix(double **matrix,
1070 const size_t number_rows)
1071{
1072 ssize_t
1073 i;
1074
1075 if (matrix == (double **) NULL )
1076 return(matrix);
1077 for (i=0; i < (ssize_t) number_rows; i++)
1078 matrix[i]=(double *) RelinquishMagickMemory(matrix[i]);
1079 matrix=(double **) RelinquishMagickMemory(matrix);
1080 return(matrix);
1081}
1082␌
1083/*
1084%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1085% %
1086% %
1087% %
1088% S e t M a t r i x E l e m e n t %
1089% %
1090% %
1091% %
1092%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1093%
1094% SetMatrixElement() sets the specified element in the matrix.
1095%
1096% The format of the SetMatrixElement method is:
1097%
1098% MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
1099% const ssize_t x,const ssize_t y,void *value)
1100%
1101% A description of each parameter follows:
1102%
1103% o matrix_info: the matrix columns.
1104%
1105% o x: the matrix x-offset.
1106%
1107% o y: the matrix y-offset.
1108%
1109% o value: set the matrix element to this value.
1110%
1111*/
1112
1113MagickExport MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
1114 const ssize_t x,const ssize_t y,const void *value)
1115{
1116 MagickOffsetType
1117 count,
1118 i;
1119
1120 assert(matrix_info != (const MatrixInfo *) NULL);
1121 assert(matrix_info->signature == MagickCoreSignature);
1122 i=y*(MagickOffsetType) matrix_info->columns+x;
1123 if ((i < 0) ||
1124 (((MagickSizeType) i*matrix_info->stride) >= matrix_info->length))
1125 return(MagickFalse);
1126 if (matrix_info->type != DiskCache)
1127 {
1128 (void) memcpy((unsigned char *) matrix_info->elements+i*
1129 (MagickOffsetType) matrix_info->stride,value,matrix_info->stride);
1130 return(MagickTrue);
1131 }
1132 count=WriteMatrixElements(matrix_info,i*(MagickOffsetType)
1133 matrix_info->stride,matrix_info->stride,(unsigned char *) value);
1134 if (count != (MagickOffsetType) matrix_info->stride)
1135 return(MagickFalse);
1136 return(MagickTrue);
1137}
_CacheView
Definition cache-view.c:66
_ExceptionInfo
Definition exception.h:102
_ImageInfo
Definition image.h:359
_Image
Definition image.h:132
_MatrixInfo
Definition matrix.c:64