2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 %
13 % MagickCore Matrix Methods %
20 % Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
26 % http://www.imagemagick.org/script/license.php %
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. %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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/matrix.h"
49 #include "MagickCore/memory_.h"
50 #include "MagickCore/pixel-accessor.h"
51 #include "MagickCore/pixel-private.h"
52 #include "MagickCore/resource_.h"
53 #include "MagickCore/semaphore.h"
54 #include "MagickCore/thread-private.h"
55 #include "MagickCore/utility.h"
94 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
98 % A c q u i r e M a t r i x I n f o %
102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
104 % AcquireMatrixInfo() allocates the ImageInfo structure.
106 % The format of the AcquireMatrixInfo method is:
108 % MatrixInfo *AcquireMatrixInfo(const size_t columns,const size_t rows,
109 % const size_t stride,ExceptionInfo *exception)
111 % A description of each parameter follows:
113 % o columns: the matrix columns.
115 % o rows: the matrix rows.
117 % o stride: the matrix stride.
119 % o exception: return any errors or warnings in this structure.
123 static inline MagickSizeType MagickMin(const MagickSizeType x,
124 const MagickSizeType y)
132 static void MatrixSignalHandler(int status)
134 ThrowFatalException(CacheFatalError,"UnableToExtendMatrixCache");
138 static inline MagickOffsetType WriteMatrixElements(
139 const MatrixInfo *restrict matrix_info,const MagickOffsetType offset,
140 const MagickSizeType length,const unsigned char *restrict buffer)
142 register MagickOffsetType
148 #if !defined(MAGICKCORE_HAVE_PWRITE)
149 LockSemaphoreInfo(matrix_info->semaphore);
150 if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
152 UnlockSemaphoreInfo(matrix_info->semaphore);
153 return((MagickOffsetType) -1);
157 for (i=0; i < (MagickOffsetType) length; i+=count)
159 #if !defined(MAGICKCORE_HAVE_PWRITE)
160 count=write(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
161 (MagickSizeType) SSIZE_MAX));
163 count=pwrite(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
164 (MagickSizeType) SSIZE_MAX),(off_t) (offset+i));
173 #if !defined(MAGICKCORE_HAVE_PWRITE)
174 UnlockSemaphoreInfo(matrix_info->semaphore);
179 static MagickBooleanType SetMatrixExtent(MatrixInfo *restrict matrix_info,
180 MagickSizeType length)
187 if (length != (MagickSizeType) ((MagickOffsetType) length))
189 offset=(MagickOffsetType) lseek(matrix_info->file,0,SEEK_END);
192 if ((MagickSizeType) offset >= length)
194 extent=(MagickOffsetType) length-1;
195 count=WriteMatrixElements(matrix_info,extent,1,(const unsigned char *) "");
196 #if defined(MAGICKCORE_HAVE_POSIX_FALLOCATE)
197 if (matrix_info->synchronize != MagickFalse)
202 status=posix_fallocate(matrix_info->file,offset+1,extent-offset);
208 (void) signal(SIGBUS,MatrixSignalHandler);
210 return(count != (MagickOffsetType) 1 ? MagickFalse : MagickTrue);
213 MagickExport MatrixInfo *AcquireMatrixInfo(const size_t columns,
214 const size_t rows,const size_t stride,ExceptionInfo *exception)
225 matrix_info=(MatrixInfo *) AcquireMagickMemory(sizeof(*matrix_info));
226 if (matrix_info == (MatrixInfo *) NULL)
227 return((MatrixInfo *) NULL);
228 (void) ResetMagickMemory(matrix_info,0,sizeof(*matrix_info));
229 matrix_info->signature=MagickSignature;
230 matrix_info->columns=columns;
231 matrix_info->rows=rows;
232 matrix_info->stride=stride;
233 matrix_info->semaphore=AcquireSemaphoreInfo();
234 synchronize=GetEnvironmentValue("MAGICK_SYNCHRONIZE");
235 if (synchronize != (const char *) NULL)
237 matrix_info->synchronize=IsStringTrue(synchronize);
238 synchronize=DestroyString(synchronize);
240 matrix_info->length=(MagickSizeType) columns*rows*stride;
241 if (matrix_info->columns != (size_t) (matrix_info->length/rows/stride))
243 (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
244 "CacheResourcesExhausted","`%s'","matrix cache");
245 return(DestroyMatrixInfo(matrix_info));
247 matrix_info->type=MemoryCache;
248 status=AcquireMagickResource(AreaResource,matrix_info->length);
249 if ((status != MagickFalse) &&
250 (matrix_info->length == (MagickSizeType) ((size_t) matrix_info->length)))
252 status=AcquireMagickResource(MemoryResource,matrix_info->length);
253 if (status != MagickFalse)
255 matrix_info->mapped=MagickFalse;
256 matrix_info->elements=AcquireMagickMemory((size_t)
257 matrix_info->length);
258 if (matrix_info->elements == NULL)
260 matrix_info->mapped=MagickTrue;
261 matrix_info->elements=MapBlob(-1,IOMode,0,(size_t)
262 matrix_info->length);
264 if (matrix_info->elements == (unsigned short *) NULL)
265 RelinquishMagickResource(MemoryResource,matrix_info->length);
268 matrix_info->file=(-1);
269 if (matrix_info->elements == (unsigned short *) NULL)
271 status=AcquireMagickResource(DiskResource,matrix_info->length);
272 if (status == MagickFalse)
274 (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
275 "CacheResourcesExhausted","`%s'","matrix cache");
276 return(DestroyMatrixInfo(matrix_info));
278 matrix_info->type=DiskCache;
279 (void) AcquireMagickResource(MemoryResource,matrix_info->length);
280 matrix_info->file=AcquireUniqueFileResource(matrix_info->path);
281 if (matrix_info->file == -1)
282 return(DestroyMatrixInfo(matrix_info));
283 status=AcquireMagickResource(MapResource,matrix_info->length);
284 if (status != MagickFalse)
286 status=SetMatrixExtent(matrix_info,matrix_info->length);
287 if (status != MagickFalse)
289 matrix_info->elements=(void *) MapBlob(matrix_info->file,IOMode,0,
290 (size_t) matrix_info->length);
291 if (matrix_info->elements != NULL)
292 matrix_info->type=MapCache;
294 RelinquishMagickResource(MapResource,matrix_info->length);
302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
306 % A c q u i r e M a g i c k M a t r i x %
310 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
312 % AcquireMagickMatrix() allocates and returns a matrix in the form of an
313 % array of pointers to an array of doubles, with all values pre-set to zero.
315 % This used to generate the two dimensional matrix, and vectors required
316 % for the GaussJordanElimination() method below, solving some system of
317 % simultanious equations.
319 % The format of the AcquireMagickMatrix method is:
321 % double **AcquireMagickMatrix(const size_t number_rows,
324 % A description of each parameter follows:
326 % o number_rows: the number pointers for the array of pointers
329 % o size: the size of the array of doubles each pointer points to
330 % (second dimension).
333 MagickExport double **AcquireMagickMatrix(const size_t number_rows,
343 matrix=(double **) AcquireQuantumMemory(number_rows,sizeof(*matrix));
344 if (matrix == (double **) NULL)
345 return((double **)NULL);
346 for (i=0; i < (ssize_t) number_rows; i++)
348 matrix[i]=(double *) AcquireQuantumMemory(size,sizeof(*matrix[i]));
349 if (matrix[i] == (double *) NULL)
351 for (j=0; j < i; j++)
352 matrix[j]=(double *) RelinquishMagickMemory(matrix[j]);
353 matrix=(double **) RelinquishMagickMemory(matrix);
354 return((double **) NULL);
356 for (j=0; j < (ssize_t) size; j++)
363 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
367 % D e s t r o y M a t r i x I n f o %
371 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
373 % DestroyMatrixInfo() dereferences a matrix, deallocating memory associated
376 % The format of the DestroyImage method is:
378 % MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info)
380 % A description of each parameter follows:
382 % o matrix_info: the matrix.
385 MagickExport MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info)
387 assert(matrix_info != (MatrixInfo *) NULL);
388 assert(matrix_info->signature == MagickSignature);
389 LockSemaphoreInfo(matrix_info->semaphore);
390 switch (matrix_info->type)
394 if (matrix_info->mapped == MagickFalse)
395 matrix_info->elements=RelinquishMagickMemory(matrix_info->elements);
398 (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length);
399 matrix_info->elements=(unsigned short *) NULL;
401 RelinquishMagickResource(MemoryResource,matrix_info->length);
406 (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length);
407 matrix_info->elements=NULL;
408 RelinquishMagickResource(MapResource,matrix_info->length);
412 if (matrix_info->file != -1)
413 (void) close(matrix_info->file);
414 (void) RelinquishUniqueFileResource(matrix_info->path);
415 RelinquishMagickResource(DiskResource,matrix_info->length);
421 UnlockSemaphoreInfo(matrix_info->semaphore);
422 RelinquishSemaphoreInfo(&matrix_info->semaphore);
423 return((MatrixInfo *) RelinquishMagickMemory(matrix_info));
427 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
431 % G a u s s J o r d a n E l i m i n a t i o n %
435 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
437 % GaussJordanElimination() returns a matrix in reduced row echelon form,
438 % while simultaneously reducing and thus solving the augumented results
441 % See also http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
443 % The format of the GaussJordanElimination method is:
445 % MagickBooleanType GaussJordanElimination(double **matrix,double **vectors,
446 % const size_t rank,const size_t number_vectors)
448 % A description of each parameter follows:
450 % o matrix: the matrix to be reduced, as an 'array of row pointers'.
452 % o vectors: the additional matrix argumenting the matrix for row reduction.
453 % Producing an 'array of column vectors'.
455 % o rank: The size of the matrix (both rows and columns).
456 % Also represents the number terms that need to be solved.
458 % o number_vectors: Number of vectors columns, argumenting the above matrix.
459 % Usally 1, but can be more for more complex equation solving.
461 % Note that the 'matrix' is given as a 'array of row pointers' of rank size.
462 % That is values can be assigned as matrix[row][column] where 'row' is
463 % typically the equation, and 'column' is the term of the equation.
464 % That is the matrix is in the form of a 'row first array'.
466 % However 'vectors' is a 'array of column pointers' which can have any number
467 % of columns, with each column array the same 'rank' size as 'matrix'.
469 % This allows for simpler handling of the results, especially is only one
470 % column 'vector' is all that is required to produce the desired solution.
472 % For example, the 'vectors' can consist of a pointer to a simple array of
473 % doubles. when only one set of simultanious equations is to be solved from
474 % the given set of coefficient weighted terms.
476 % double **matrix = AcquireMagickMatrix(8UL,8UL);
477 % double coefficents[8];
479 % GaussJordanElimination(matrix, &coefficents, 8UL, 1UL);
481 % However by specifing more 'columns' (as an 'array of vector columns',
482 % you can use this function to solve a set of 'separable' equations.
484 % For example a distortion function where u = U(x,y) v = V(x,y)
485 % And the functions U() and V() have separate coefficents, but are being
486 % generated from a common x,y->u,v data set.
488 % Another example is generation of a color gradient from a set of colors
489 % at specific coordients, such as a list x,y -> r,g,b,a
490 % (Reference to be added - Anthony)
492 % You can also use the 'vectors' to generate an inverse of the given 'matrix'
493 % though as a 'column first array' rather than a 'row first array'. For
494 % details see http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
497 MagickExport MagickBooleanType GaussJordanElimination(double **matrix,
498 double **vectors,const size_t rank,const size_t number_vectors)
500 #define GaussJordanSwap(x,y) \
526 columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns));
527 rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows));
528 pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots));
529 if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) ||
530 (pivots == (ssize_t *) NULL))
532 if (pivots != (ssize_t *) NULL)
533 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
534 if (columns != (ssize_t *) NULL)
535 columns=(ssize_t *) RelinquishMagickMemory(columns);
536 if (rows != (ssize_t *) NULL)
537 rows=(ssize_t *) RelinquishMagickMemory(rows);
540 (void) ResetMagickMemory(columns,0,rank*sizeof(*columns));
541 (void) ResetMagickMemory(rows,0,rank*sizeof(*rows));
542 (void) ResetMagickMemory(pivots,0,rank*sizeof(*pivots));
545 for (i=0; i < (ssize_t) rank; i++)
548 for (j=0; j < (ssize_t) rank; j++)
551 for (k=0; k < (ssize_t) rank; k++)
558 if (fabs(matrix[j][k]) >= max)
560 max=fabs(matrix[j][k]);
568 for (k=0; k < (ssize_t) rank; k++)
569 GaussJordanSwap(matrix[row][k],matrix[column][k]);
570 for (k=0; k < (ssize_t) number_vectors; k++)
571 GaussJordanSwap(vectors[k][row],vectors[k][column]);
575 if (matrix[column][column] == 0.0)
576 return(MagickFalse); /* sigularity */
577 scale=PerceptibleReciprocal(matrix[column][column]);
578 matrix[column][column]=1.0;
579 for (j=0; j < (ssize_t) rank; j++)
580 matrix[column][j]*=scale;
581 for (j=0; j < (ssize_t) number_vectors; j++)
582 vectors[j][column]*=scale;
583 for (j=0; j < (ssize_t) rank; j++)
586 scale=matrix[j][column];
587 matrix[j][column]=0.0;
588 for (k=0; k < (ssize_t) rank; k++)
589 matrix[j][k]-=scale*matrix[column][k];
590 for (k=0; k < (ssize_t) number_vectors; k++)
591 vectors[k][j]-=scale*vectors[k][column];
594 for (j=(ssize_t) rank-1; j >= 0; j--)
595 if (columns[j] != rows[j])
596 for (i=0; i < (ssize_t) rank; i++)
597 GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]);
598 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
599 rows=(ssize_t *) RelinquishMagickMemory(rows);
600 columns=(ssize_t *) RelinquishMagickMemory(columns);
605 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
609 % G e t M a t r i x C o l u m n s %
613 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
615 % GetMatrixColumns() returns the number of columns in the matrix.
617 % The format of the GetMatrixColumns method is:
619 % size_t GetMatrixColumns(const MatrixInfo *matrix_info)
621 % A description of each parameter follows:
623 % o matrix_info: the matrix.
626 MagickExport size_t GetMatrixColumns(const MatrixInfo *matrix_info)
628 assert(matrix_info != (MatrixInfo *) NULL);
629 assert(matrix_info->signature == MagickSignature);
630 return(matrix_info->columns);
634 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
638 % G e t M a t r i x E l e m e n t %
642 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
644 % GetMatrixElement() returns the specifed element in the matrix.
646 % The format of the GetMatrixElement method is:
648 % MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
649 % const ssize_t x,const ssize_t y,void *value)
651 % A description of each parameter follows:
653 % o matrix_info: the matrix columns.
655 % o x: the matrix x-offset.
657 % o y: the matrix y-offset.
659 % o value: return the matrix element in this buffer.
663 static inline ssize_t EdgeX(const ssize_t x,const size_t columns)
667 if (x >= (ssize_t) columns)
668 return((ssize_t) (columns-1));
672 static inline ssize_t EdgeY(const ssize_t y,const size_t rows)
676 if (y >= (ssize_t) rows)
677 return((ssize_t) (rows-1));
681 static inline MagickOffsetType ReadMatrixElements(
682 const MatrixInfo *restrict matrix_info,const MagickOffsetType offset,
683 const MagickSizeType length,unsigned char *restrict buffer)
685 register MagickOffsetType
691 #if !defined(MAGICKCORE_HAVE_PREAD)
692 LockSemaphoreInfo(matrix_info->semaphore);
693 if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
695 UnlockSemaphoreInfo(matrix_info->semaphore);
696 return((MagickOffsetType) -1);
700 for (i=0; i < (MagickOffsetType) length; i+=count)
702 #if !defined(MAGICKCORE_HAVE_PREAD)
703 count=read(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
704 (MagickSizeType) SSIZE_MAX));
706 count=pread(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
707 (MagickSizeType) SSIZE_MAX),(off_t) (offset+i));
716 #if !defined(MAGICKCORE_HAVE_PREAD)
717 UnlockSemaphoreInfo(matrix_info->semaphore);
722 MagickExport MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
723 const ssize_t x,const ssize_t y,void *value)
729 assert(matrix_info != (const MatrixInfo *) NULL);
730 assert(matrix_info->signature == MagickSignature);
731 i=(MagickOffsetType) EdgeY(y,matrix_info->rows)*matrix_info->columns+
732 EdgeX(x,matrix_info->columns);
733 if (matrix_info->type != DiskCache)
735 (void) memcpy(value,(unsigned char *) matrix_info->elements+i*
736 matrix_info->stride,matrix_info->stride);
739 count=ReadMatrixElements(matrix_info,i*matrix_info->stride,
740 matrix_info->stride,value);
741 if (count != (MagickOffsetType) matrix_info->stride)
747 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
751 % G e t M a t r i x R o w s %
755 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
757 % GetMatrixRows() returns the number of rows in the matrix.
759 % The format of the GetMatrixRows method is:
761 % size_t GetMatrixRows(const MatrixInfo *matrix_info)
763 % A description of each parameter follows:
765 % o matrix_info: the matrix.
768 MagickExport size_t GetMatrixRows(const MatrixInfo *matrix_info)
770 assert(matrix_info != (const MatrixInfo *) NULL);
771 assert(matrix_info->signature == MagickSignature);
772 return(matrix_info->rows);
776 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
780 % L e a s t S q u a r e s A d d T e r m s %
784 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
786 % LeastSquaresAddTerms() adds one set of terms and associate results to the
787 % given matrix and vectors for solving using least-squares function fitting.
789 % The format of the AcquireMagickMatrix method is:
791 % void LeastSquaresAddTerms(double **matrix,double **vectors,
792 % const double *terms,const double *results,const size_t rank,
793 % const size_t number_vectors);
795 % A description of each parameter follows:
797 % o matrix: the square matrix to add given terms/results to.
799 % o vectors: the result vectors to add terms/results to.
801 % o terms: the pre-calculated terms (without the unknown coefficent
802 % weights) that forms the equation being added.
804 % o results: the result(s) that should be generated from the given terms
805 % weighted by the yet-to-be-solved coefficents.
807 % o rank: the rank or size of the dimensions of the square matrix.
808 % Also the length of vectors, and number of terms being added.
810 % o number_vectors: Number of result vectors, and number or results being
811 % added. Also represents the number of separable systems of equations
812 % that is being solved.
816 % 2 dimensional Affine Equations (which are separable)
817 % c0*x + c2*y + c4*1 => u
818 % c1*x + c3*y + c5*1 => v
820 % double **matrix = AcquireMagickMatrix(3UL,3UL);
821 % double **vectors = AcquireMagickMatrix(2UL,3UL);
822 % double terms[3], results[2];
824 % for each given x,y -> u,v
830 % LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL);
832 % if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) {
833 % c0 = vectors[0][0];
834 % c2 = vectors[0][1];
835 % c4 = vectors[0][2];
836 % c1 = vectors[1][0];
837 % c3 = vectors[1][1];
838 % c5 = vectors[1][2];
841 % printf("Matrix unsolvable\n);
842 % RelinquishMagickMatrix(matrix,3UL);
843 % RelinquishMagickMatrix(vectors,2UL);
846 MagickExport void LeastSquaresAddTerms(double **matrix,double **vectors,
847 const double *terms,const double *results,const size_t rank,
848 const size_t number_vectors)
854 for (j=0; j < (ssize_t) rank; j++)
856 for (i=0; i < (ssize_t) rank; i++)
857 matrix[i][j]+=terms[i]*terms[j];
858 for (i=0; i < (ssize_t) number_vectors; i++)
859 vectors[i][j]+=results[i]*terms[j];
865 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
869 % M a t r i x T o I m a g e %
873 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
875 % MatrixToImage() returns a matrix as an image. The matrix elements must be
876 % of type double otherwise nonsense is returned.
878 % The format of the MatrixToImage method is:
880 % Image *MatrixToImage(const MatrixInfo *matrix_info,
881 % ExceptionInfo *exception)
883 % A description of each parameter follows:
885 % o matrix_info: the matrix.
887 % o exception: return any errors or warnings in this structure.
890 MagickExport Image *MatrixToImage(const MatrixInfo *matrix_info,
891 ExceptionInfo *exception)
911 assert(matrix_info != (const MatrixInfo *) NULL);
912 assert(matrix_info->signature == MagickSignature);
913 assert(exception != (ExceptionInfo *) NULL);
914 assert(exception->signature == MagickSignature);
915 if (matrix_info->stride < sizeof(double))
916 return((Image *) NULL);
918 Determine range of matrix.
920 (void) GetMatrixElement(matrix_info,0,0,&value);
923 for (y=0; y < (ssize_t) matrix_info->rows; y++)
928 for (x=0; x < (ssize_t) matrix_info->columns; x++)
930 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
932 if (value < min_value)
935 if (value > max_value)
939 if ((min_value == 0.0) && (max_value == 0.0))
942 if (min_value == max_value)
944 scale_factor=(double) QuantumRange/min_value;
948 scale_factor=(double) QuantumRange/(max_value-min_value);
950 Convert matrix to image.
952 image=AcquireImage((ImageInfo *) NULL,exception);
953 image->columns=matrix_info->columns;
954 image->rows=matrix_info->rows;
955 image->colorspace=GRAYColorspace;
957 image_view=AcquireAuthenticCacheView(image,exception);
958 #if defined(MAGICKCORE_OPENMP_SUPPORT)
959 #pragma omp parallel for schedule(static,4) shared(status) \
960 magick_threads(image,image,image->rows,1)
962 for (y=0; y < (ssize_t) image->rows; y++)
973 if (status == MagickFalse)
975 q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
976 if (q == (Quantum *) NULL)
981 for (x=0; x < (ssize_t) image->columns; x++)
983 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
985 value=scale_factor*(value-min_value);
986 *q=ClampToQuantum(value);
987 q+=GetPixelChannels(image);
989 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
992 image_view=DestroyCacheView(image_view);
993 if (status == MagickFalse)
994 image=DestroyImage(image);
999 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1003 % N u l l M a t r i x %
1007 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1009 % NullMatrix() sets all elements of the matrix to zero.
1011 % The format of the ResetMagickMemory method is:
1013 % MagickBooleanType *NullMatrix(MatrixInfo *matrix_info)
1015 % A description of each parameter follows:
1017 % o matrix_info: the matrix.
1020 MagickExport MagickBooleanType NullMatrix(MatrixInfo *matrix_info)
1032 assert(matrix_info != (const MatrixInfo *) NULL);
1033 assert(matrix_info->signature == MagickSignature);
1034 if (matrix_info->type != DiskCache)
1036 (void) ResetMagickMemory(matrix_info->elements,0,(size_t)
1037 matrix_info->length);
1041 (void) lseek(matrix_info->file,0,SEEK_SET);
1042 for (y=0; y < (ssize_t) matrix_info->rows; y++)
1044 for (x=0; x < (ssize_t) matrix_info->length; x++)
1046 count=write(matrix_info->file,&value,sizeof(value));
1047 if (count != (ssize_t) sizeof(value))
1050 if (x < (ssize_t) matrix_info->length)
1053 return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue);
1057 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1061 % R e l i n q u i s h M a g i c k M a t r i x %
1065 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1067 % RelinquishMagickMatrix() frees the previously acquired matrix (array of
1068 % pointers to arrays of doubles).
1070 % The format of the RelinquishMagickMatrix method is:
1072 % double **RelinquishMagickMatrix(double **matrix,
1073 % const size_t number_rows)
1075 % A description of each parameter follows:
1077 % o matrix: the matrix to relinquish
1079 % o number_rows: the first dimension of the acquired matrix (number of
1083 MagickExport double **RelinquishMagickMatrix(double **matrix,
1084 const size_t number_rows)
1089 if (matrix == (double **) NULL )
1091 for (i=0; i < (ssize_t) number_rows; i++)
1092 matrix[i]=(double *) RelinquishMagickMemory(matrix[i]);
1093 matrix=(double **) RelinquishMagickMemory(matrix);
1098 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1102 % S e t M a t r i x E l e m e n t %
1106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1108 % SetMatrixElement() sets the specifed element in the matrix.
1110 % The format of the SetMatrixElement method is:
1112 % MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
1113 % const ssize_t x,const ssize_t y,void *value)
1115 % A description of each parameter follows:
1117 % o matrix_info: the matrix columns.
1119 % o x: the matrix x-offset.
1121 % o y: the matrix y-offset.
1123 % o value: set the matrix element to this value.
1127 MagickExport MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
1128 const ssize_t x,const ssize_t y,const void *value)
1134 assert(matrix_info != (const MatrixInfo *) NULL);
1135 assert(matrix_info->signature == MagickSignature);
1136 i=(MagickOffsetType) y*matrix_info->columns+x;
1138 ((MagickSizeType) (i*matrix_info->stride) >= matrix_info->length))
1139 return(MagickFalse);
1140 if (matrix_info->type != DiskCache)
1142 (void) memcpy((unsigned char *) matrix_info->elements+i*
1143 matrix_info->stride,value,matrix_info->stride);
1146 count=WriteMatrixElements(matrix_info,i*matrix_info->stride,
1147 matrix_info->stride,value);
1148 if (count != (MagickOffsetType) matrix_info->stride)
1149 return(MagickFalse);