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