diff --git a/MagickCore/matrix.c.orig b/MagickCore/matrix.c.orig deleted file mode 100644 index 9c1ea9d739..0000000000 --- a/MagickCore/matrix.c.orig +++ /dev/null @@ -1,1136 +0,0 @@ -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -% M M AAA TTTTT RRRR IIIII X X % -% MM MM A A T R R I X X % -% M M M AAAAA T RRRR I X % -% M M A A T R R I X X % -% M M A A T R R IIIII X X % -% % -% % -% MagickCore Matrix Methods % -% % -% Software Design % -% Cristy % -% August 2007 % -% % -% % -% Copyright @ 1999 ImageMagick Studio LLC, a non-profit organization % -% dedicated to making software imaging solutions freely available. % -% % -% You may not use this file except in compliance with the License. You may % -% obtain a copy of the License at % -% % -% https://imagemagick.org/script/license.php % -% % -% Unless required by applicable law or agreed to in writing, software % -% distributed under the License is distributed on an "AS IS" BASIS, % -% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. % -% See the License for the specific language governing permissions and % -% limitations under the License. % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% -*/ - -/* - Include declarations. -*/ -#include "MagickCore/studio.h" -#include "MagickCore/blob.h" -#include "MagickCore/blob-private.h" -#include "MagickCore/cache.h" -#include "MagickCore/exception.h" -#include "MagickCore/exception-private.h" -#include "MagickCore/image-private.h" -#include "MagickCore/matrix.h" -#include "MagickCore/matrix-private.h" -#include "MagickCore/memory_.h" -#include "MagickCore/nt-base-private.h" -#include "MagickCore/pixel-accessor.h" -#include "MagickCore/resource_.h" -#include "MagickCore/semaphore.h" -#include "MagickCore/thread-private.h" -#include "MagickCore/utility.h" - -/* - Typedef declaration. -*/ -struct _MatrixInfo -{ - CacheType - type; - - size_t - columns, - rows, - stride; - - MagickSizeType - length; - - MagickBooleanType - mapped, - synchronize; - - char - path[MagickPathExtent]; - - int - file; - - void - *elements; - - SemaphoreInfo - *semaphore; - - size_t - signature; -}; - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -% A c q u i r e M a t r i x I n f o % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% AcquireMatrixInfo() allocates the ImageInfo structure. -% -% The format of the AcquireMatrixInfo method is: -% -% MatrixInfo *AcquireMatrixInfo(const size_t columns,const size_t rows, -% const size_t stride,ExceptionInfo *exception) -% -% A description of each parameter follows: -% -% o columns: the matrix columns. -% -% o rows: the matrix rows. -% -% o stride: the matrix stride. -% -% o exception: return any errors or warnings in this structure. -% -*/ - -#if defined(SIGBUS) -static void MatrixSignalHandler(int magick_unused(status)) -{ - magick_unreferenced(status); - ThrowFatalException(CacheFatalError,"UnableToExtendMatrixCache"); -} -#endif - -static inline MagickOffsetType WriteMatrixElements( - const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset, - const MagickSizeType length,const unsigned char *magick_restrict buffer) -{ - MagickOffsetType - i; - - ssize_t - count; - -#if !defined(MAGICKCORE_HAVE_PWRITE) - LockSemaphoreInfo(matrix_info->semaphore); - if (lseek(matrix_info->file,offset,SEEK_SET) < 0) - { - UnlockSemaphoreInfo(matrix_info->semaphore); - return((MagickOffsetType) -1); - } -#endif - count=0; - for (i=0; i < (MagickOffsetType) length; i+=count) - { -#if !defined(MAGICKCORE_HAVE_PWRITE) - count=write(matrix_info->file,buffer+i,(size_t) MagickMin(length- - (MagickSizeType) i,(MagickSizeType) MagickMaxBufferExtent)); -#else - count=pwrite(matrix_info->file,buffer+i,(size_t) MagickMin(length- - (MagickSizeType) i,(MagickSizeType) MagickMaxBufferExtent),offset+i); -#endif - if (count <= 0) - { - count=0; - if (errno != EINTR) - break; - } - } -#if !defined(MAGICKCORE_HAVE_PWRITE) - UnlockSemaphoreInfo(matrix_info->semaphore); -#endif - return(i); -} - -static MagickBooleanType SetMatrixExtent( - MatrixInfo *magick_restrict matrix_info,MagickSizeType length) -{ - MagickOffsetType - count, - extent, - offset; - - if (length != (MagickSizeType) ((MagickOffsetType) length)) - return(MagickFalse); - offset=(MagickOffsetType) lseek(matrix_info->file,0,SEEK_END); - if (offset < 0) - return(MagickFalse); - if ((MagickSizeType) offset >= length) - return(MagickTrue); - extent=(MagickOffsetType) length-1; - count=WriteMatrixElements(matrix_info,extent,1,(const unsigned char *) ""); -#if defined(MAGICKCORE_HAVE_POSIX_FALLOCATE) - if (matrix_info->synchronize != MagickFalse) - (void) posix_fallocate(matrix_info->file,offset+1,extent-offset); -#endif -#if defined(SIGBUS) - (void) signal(SIGBUS,MatrixSignalHandler); -#endif - return(count != (MagickOffsetType) 1 ? MagickFalse : MagickTrue); -} - -MagickExport MatrixInfo *AcquireMatrixInfo(const size_t columns, - const size_t rows,const size_t stride,ExceptionInfo *exception) -{ - char - *synchronize; - - MagickBooleanType - status; - - MatrixInfo - *matrix_info; - - matrix_info=(MatrixInfo *) AcquireMagickMemory(sizeof(*matrix_info)); - if (matrix_info == (MatrixInfo *) NULL) - return((MatrixInfo *) NULL); - (void) memset(matrix_info,0,sizeof(*matrix_info)); - matrix_info->signature=MagickCoreSignature; - matrix_info->columns=columns; - matrix_info->rows=rows; - matrix_info->stride=stride; - matrix_info->semaphore=AcquireSemaphoreInfo(); - synchronize=GetEnvironmentValue("MAGICK_SYNCHRONIZE"); - if (synchronize != (const char *) NULL) - { - matrix_info->synchronize=IsStringTrue(synchronize); - synchronize=DestroyString(synchronize); - } - matrix_info->length=(MagickSizeType) columns*rows*stride; - if (matrix_info->columns != (size_t) (matrix_info->length/rows/stride)) - { - (void) ThrowMagickException(exception,GetMagickModule(),CacheError, - "CacheResourcesExhausted","`%s'","matrix cache"); - return(DestroyMatrixInfo(matrix_info)); - } - matrix_info->type=MemoryCache; - status=AcquireMagickResource(AreaResource,matrix_info->length); - if ((status != MagickFalse) && - (matrix_info->length == (MagickSizeType) ((size_t) matrix_info->length))) - { - status=AcquireMagickResource(MemoryResource,matrix_info->length); - if (status != MagickFalse) - { - matrix_info->mapped=MagickFalse; - matrix_info->elements=AcquireMagickMemory((size_t) - matrix_info->length); - if (matrix_info->elements == NULL) - { - matrix_info->mapped=MagickTrue; - matrix_info->elements=MapBlob(-1,IOMode,0,(size_t) - matrix_info->length); - } - if (matrix_info->elements == (unsigned short *) NULL) - RelinquishMagickResource(MemoryResource,matrix_info->length); - } - } - matrix_info->file=(-1); - if (matrix_info->elements == (unsigned short *) NULL) - { - status=AcquireMagickResource(DiskResource,matrix_info->length); - if (status == MagickFalse) - { - (void) ThrowMagickException(exception,GetMagickModule(),CacheError, - "CacheResourcesExhausted","`%s'","matrix cache"); - return(DestroyMatrixInfo(matrix_info)); - } - matrix_info->type=DiskCache; - matrix_info->file=AcquireUniqueFileResource(matrix_info->path); - if (matrix_info->file == -1) - return(DestroyMatrixInfo(matrix_info)); - status=AcquireMagickResource(MapResource,matrix_info->length); - if (status != MagickFalse) - { - status=SetMatrixExtent(matrix_info,matrix_info->length); - if (status != MagickFalse) - matrix_info->elements=(void *) MapBlob(matrix_info->file,IOMode,0, - (size_t) matrix_info->length); - if (matrix_info->elements != NULL) - matrix_info->type=MapCache; - else - RelinquishMagickResource(MapResource,matrix_info->length); - } - } - return(matrix_info); -} - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -% A c q u i r e M a g i c k M a t r i x % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% AcquireMagickMatrix() allocates and returns a matrix in the form of an -% array of pointers to an array of doubles, with all values pre-set to zero. -% -% This used to generate the two dimensional matrix, and vectors required -% for the GaussJordanElimination() method below, solving some system of -% simultaneous equations. -% -% The format of the AcquireMagickMatrix method is: -% -% double **AcquireMagickMatrix(const size_t number_rows, -% const size_t size) -% -% A description of each parameter follows: -% -% o number_rows: the number pointers for the array of pointers -% (first dimension). -% -% o size: the size of the array of doubles each pointer points to -% (second dimension). -% -*/ -MagickExport double **AcquireMagickMatrix(const size_t number_rows, - const size_t size) -{ - double - **matrix; - - ssize_t - i, - j; - - matrix=(double **) AcquireQuantumMemory(number_rows,sizeof(*matrix)); - if (matrix == (double **) NULL) - return((double **) NULL); - for (i=0; i < (ssize_t) number_rows; i++) - { - matrix[i]=(double *) AcquireQuantumMemory(size,sizeof(*matrix[i])); - if (matrix[i] == (double *) NULL) - { - for (j=0; j < i; j++) - matrix[j]=(double *) RelinquishMagickMemory(matrix[j]); - matrix=(double **) RelinquishMagickMemory(matrix); - return((double **) NULL); - } - for (j=0; j < (ssize_t) size; j++) - matrix[i][j]=0.0; - } - return(matrix); -} - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -% D e s t r o y M a t r i x I n f o % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% DestroyMatrixInfo() dereferences a matrix, deallocating memory associated -% with the matrix. -% -% The format of the DestroyImage method is: -% -% MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info) -% -% A description of each parameter follows: -% -% o matrix_info: the matrix. -% -*/ -MagickExport MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info) -{ - assert(matrix_info != (MatrixInfo *) NULL); - assert(matrix_info->signature == MagickCoreSignature); - LockSemaphoreInfo(matrix_info->semaphore); - switch (matrix_info->type) - { - case MemoryCache: - { - if (matrix_info->mapped == MagickFalse) - matrix_info->elements=RelinquishMagickMemory(matrix_info->elements); - else - { - (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length); - matrix_info->elements=(unsigned short *) NULL; - } - RelinquishMagickResource(MemoryResource,matrix_info->length); - break; - } - case MapCache: - { - (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length); - matrix_info->elements=NULL; - RelinquishMagickResource(MapResource,matrix_info->length); - magick_fallthrough; - } - case DiskCache: - { - if (matrix_info->file != -1) - (void) close(matrix_info->file); - (void) RelinquishUniqueFileResource(matrix_info->path); - RelinquishMagickResource(DiskResource,matrix_info->length); - break; - } - default: - break; - } - UnlockSemaphoreInfo(matrix_info->semaphore); - RelinquishSemaphoreInfo(&matrix_info->semaphore); - return((MatrixInfo *) RelinquishMagickMemory(matrix_info)); -} - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -+ G a u s s J o r d a n E l i m i n a t i o n % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% GaussJordanElimination() returns a matrix in reduced row echelon form, -% while simultaneously reducing and thus solving the augmented results -% matrix. -% -% See also http://en.wikipedia.org/wiki/Gauss-Jordan_elimination -% -% The format of the GaussJordanElimination method is: -% -% MagickBooleanType GaussJordanElimination(double **matrix, -% double **vectors,const size_t rank,const size_t number_vectors) -% -% A description of each parameter follows: -% -% o matrix: the matrix to be reduced, as an 'array of row pointers'. -% -% o vectors: the additional matrix argumenting the matrix for row reduction. -% Producing an 'array of column vectors'. -% -% o rank: The size of the matrix (both rows and columns). -% Also represents the number terms that need to be solved. -% -% o number_vectors: Number of vectors columns, argumenting the above matrix. -% Usually 1, but can be more for more complex equation solving. -% -% Note that the 'matrix' is given as a 'array of row pointers' of rank size. -% That is values can be assigned as matrix[row][column] where 'row' is -% typically the equation, and 'column' is the term of the equation. -% That is the matrix is in the form of a 'row first array'. -% -% However 'vectors' is a 'array of column pointers' which can have any number -% of columns, with each column array the same 'rank' size as 'matrix'. -% -% This allows for simpler handling of the results, especially is only one -% column 'vector' is all that is required to produce the desired solution. -% -% For example, the 'vectors' can consist of a pointer to a simple array of -% doubles. when only one set of simultaneous equations is to be solved from -% the given set of coefficient weighted terms. -% -% double **matrix = AcquireMagickMatrix(8UL,8UL); -% double coefficients[8]; -% ... -% GaussJordanElimination(matrix, &coefficients, 8UL, 1UL); -% -% However by specifying more 'columns' (as an 'array of vector columns', -% you can use this function to solve a set of 'separable' equations. -% -% For example a distortion function where u = U(x,y) v = V(x,y) -% And the functions U() and V() have separate coefficients, but are being -% generated from a common x,y->u,v data set. -% -% Another example is generation of a color gradient from a set of colors at -% specific coordinates, such as a list x,y -> r,g,b,a. -% -% You can also use the 'vectors' to generate an inverse of the given 'matrix' -% though as a 'column first array' rather than a 'row first array'. For -% details see http://en.wikipedia.org/wiki/Gauss-Jordan_elimination -% -*/ -MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix, - double **vectors,const size_t rank,const size_t number_vectors) -{ -#define GaussJordanSwap(x,y) \ -{ \ - if ((x) != (y)) \ - { \ - (x)+=(y); \ - (y)=(x)-(y); \ - (x)=(x)-(y); \ - } \ -} - - double - max, - scale; - - ssize_t - i, - j, - k; - - ssize_t - column, - *columns, - *pivots, - row, - *rows; - - columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns)); - rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows)); - pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots)); - if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) || - (pivots == (ssize_t *) NULL)) - { - if (pivots != (ssize_t *) NULL) - pivots=(ssize_t *) RelinquishMagickMemory(pivots); - if (columns != (ssize_t *) NULL) - columns=(ssize_t *) RelinquishMagickMemory(columns); - if (rows != (ssize_t *) NULL) - rows=(ssize_t *) RelinquishMagickMemory(rows); - return(MagickFalse); - } - (void) memset(columns,0,rank*sizeof(*columns)); - (void) memset(rows,0,rank*sizeof(*rows)); - (void) memset(pivots,0,rank*sizeof(*pivots)); - column=0; - row=0; - for (i=0; i < (ssize_t) rank; i++) - { - max=0.0; - for (j=0; j < (ssize_t) rank; j++) - if (pivots[j] != 1) - { - for (k=0; k < (ssize_t) rank; k++) - if (pivots[k] != 0) - { - if (pivots[k] > 1) - return(MagickFalse); - } - else - if (fabs(matrix[j][k]) >= max) - { - max=fabs(matrix[j][k]); - row=j; - column=k; - } - } - pivots[column]++; - if (row != column) - { - for (k=0; k < (ssize_t) rank; k++) - GaussJordanSwap(matrix[row][k],matrix[column][k]); - for (k=0; k < (ssize_t) number_vectors; k++) - GaussJordanSwap(vectors[k][row],vectors[k][column]); - } - rows[i]=row; - columns[i]=column; - if (matrix[column][column] == 0.0) - return(MagickFalse); /* singularity */ - scale=PerceptibleReciprocal(matrix[column][column]); - matrix[column][column]=1.0; - for (j=0; j < (ssize_t) rank; j++) - matrix[column][j]*=scale; - for (j=0; j < (ssize_t) number_vectors; j++) - vectors[j][column]*=scale; - for (j=0; j < (ssize_t) rank; j++) - if (j != column) - { - scale=matrix[j][column]; - matrix[j][column]=0.0; - for (k=0; k < (ssize_t) rank; k++) - matrix[j][k]-=scale*matrix[column][k]; - for (k=0; k < (ssize_t) number_vectors; k++) - vectors[k][j]-=scale*vectors[k][column]; - } - } - for (j=(ssize_t) rank-1; j >= 0; j--) - if (columns[j] != rows[j]) - for (i=0; i < (ssize_t) rank; i++) - GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]); - pivots=(ssize_t *) RelinquishMagickMemory(pivots); - rows=(ssize_t *) RelinquishMagickMemory(rows); - columns=(ssize_t *) RelinquishMagickMemory(columns); - return(MagickTrue); -} - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -% G e t M a t r i x C o l u m n s % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% GetMatrixColumns() returns the number of columns in the matrix. -% -% The format of the GetMatrixColumns method is: -% -% size_t GetMatrixColumns(const MatrixInfo *matrix_info) -% -% A description of each parameter follows: -% -% o matrix_info: the matrix. -% -*/ -MagickExport size_t GetMatrixColumns(const MatrixInfo *matrix_info) -{ - assert(matrix_info != (MatrixInfo *) NULL); - assert(matrix_info->signature == MagickCoreSignature); - return(matrix_info->columns); -} - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -% G e t M a t r i x E l e m e n t % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% GetMatrixElement() returns the specified element in the matrix. -% -% The format of the GetMatrixElement method is: -% -% MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info, -% const ssize_t x,const ssize_t y,void *value) -% -% A description of each parameter follows: -% -% o matrix_info: the matrix columns. -% -% o x: the matrix x-offset. -% -% o y: the matrix y-offset. -% -% o value: return the matrix element in this buffer. -% -*/ - -static inline ssize_t EdgeX(const ssize_t x,const size_t columns) -{ - if (x < 0L) - return(0L); - if (x >= (ssize_t) columns) - return((ssize_t) (columns-1)); - return(x); -} - -static inline ssize_t EdgeY(const ssize_t y,const size_t rows) -{ - if (y < 0L) - return(0L); - if (y >= (ssize_t) rows) - return((ssize_t) (rows-1)); - return(y); -} - -static inline MagickOffsetType ReadMatrixElements( - const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset, - const MagickSizeType length,unsigned char *magick_restrict buffer) -{ - MagickOffsetType - i; - - ssize_t - count; - -#if !defined(MAGICKCORE_HAVE_PREAD) - LockSemaphoreInfo(matrix_info->semaphore); - if (lseek(matrix_info->file,offset,SEEK_SET) < 0) - { - UnlockSemaphoreInfo(matrix_info->semaphore); - return((MagickOffsetType) -1); - } -#endif - count=0; - for (i=0; i < (MagickOffsetType) length; i+=count) - { -#if !defined(MAGICKCORE_HAVE_PREAD) - count=read(matrix_info->file,buffer+i,(size_t) MagickMin(length-i, - (MagickSizeType) MagickMaxBufferExtent)); -#else - count=pread(matrix_info->file,buffer+i,(size_t) MagickMin(length- - (MagickSizeType) i,(MagickSizeType) MagickMaxBufferExtent),offset+i); -#endif - if (count <= 0) - { - count=0; - if (errno != EINTR) - break; - } - } -#if !defined(MAGICKCORE_HAVE_PREAD) - UnlockSemaphoreInfo(matrix_info->semaphore); -#endif - return(i); -} - -MagickExport MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info, - const ssize_t x,const ssize_t y,void *value) -{ - MagickOffsetType - count, - i; - - assert(matrix_info != (const MatrixInfo *) NULL); - assert(matrix_info->signature == MagickCoreSignature); - i=EdgeY(y,matrix_info->rows)*(MagickOffsetType) matrix_info->columns+ - EdgeX(x,matrix_info->columns); - if (matrix_info->type != DiskCache) - { - (void) memcpy(value,(unsigned char *) matrix_info->elements+i* - (MagickOffsetType) matrix_info->stride,matrix_info->stride); - return(MagickTrue); - } - count=ReadMatrixElements(matrix_info,i*(MagickOffsetType) matrix_info->stride, - matrix_info->stride,(unsigned char *) value); - if (count != (MagickOffsetType) matrix_info->stride) - return(MagickFalse); - return(MagickTrue); -} - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -% G e t M a t r i x R o w s % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% GetMatrixRows() returns the number of rows in the matrix. -% -% The format of the GetMatrixRows method is: -% -% size_t GetMatrixRows(const MatrixInfo *matrix_info) -% -% A description of each parameter follows: -% -% o matrix_info: the matrix. -% -*/ -MagickExport size_t GetMatrixRows(const MatrixInfo *matrix_info) -{ - assert(matrix_info != (const MatrixInfo *) NULL); - assert(matrix_info->signature == MagickCoreSignature); - return(matrix_info->rows); -} - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -+ L e a s t S q u a r e s A d d T e r m s % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% LeastSquaresAddTerms() adds one set of terms and associate results to the -% given matrix and vectors for solving using least-squares function fitting. -% -% The format of the AcquireMagickMatrix method is: -% -% void LeastSquaresAddTerms(double **matrix,double **vectors, -% const double *terms,const double *results,const size_t rank, -% const size_t number_vectors); -% -% A description of each parameter follows: -% -% o matrix: the square matrix to add given terms/results to. -% -% o vectors: the result vectors to add terms/results to. -% -% o terms: the pre-calculated terms (without the unknown coefficient -% weights) that forms the equation being added. -% -% o results: the result(s) that should be generated from the given terms -% weighted by the yet-to-be-solved coefficients. -% -% o rank: the rank or size of the dimensions of the square matrix. -% Also the length of vectors, and number of terms being added. -% -% o number_vectors: Number of result vectors, and number or results being -% added. Also represents the number of separable systems of equations -% that is being solved. -% -% Example of use... -% -% 2 dimensional Affine Equations (which are separable) -% c0*x + c2*y + c4*1 => u -% c1*x + c3*y + c5*1 => v -% -% double **matrix = AcquireMagickMatrix(3UL,3UL); -% double **vectors = AcquireMagickMatrix(2UL,3UL); -% double terms[3], results[2]; -% ... -% for each given x,y -> u,v -% terms[0] = x; -% terms[1] = y; -% terms[2] = 1; -% results[0] = u; -% results[1] = v; -% LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL); -% ... -% if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) { -% c0 = vectors[0][0]; -% c2 = vectors[0][1]; -% c4 = vectors[0][2]; -% c1 = vectors[1][0]; -% c3 = vectors[1][1]; -% c5 = vectors[1][2]; -% } -% else -% printf("Matrix unsolvable\n"); -% RelinquishMagickMatrix(matrix,3UL); -% RelinquishMagickMatrix(vectors,2UL); -% -*/ -MagickPrivate void LeastSquaresAddTerms(double **matrix,double **vectors, - const double *terms,const double *results,const size_t rank, - const size_t number_vectors) -{ - ssize_t - i, - j; - - for (j=0; j < (ssize_t) rank; j++) - { - for (i=0; i < (ssize_t) rank; i++) - matrix[i][j]+=terms[i]*terms[j]; - for (i=0; i < (ssize_t) number_vectors; i++) - vectors[i][j]+=results[i]*terms[j]; - } -} - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -% M a t r i x T o I m a g e % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% MatrixToImage() returns a matrix as an image. The matrix elements must be -% of type double otherwise nonsense is returned. -% -% The format of the MatrixToImage method is: -% -% Image *MatrixToImage(const MatrixInfo *matrix_info, -% ExceptionInfo *exception) -% -% A description of each parameter follows: -% -% o matrix_info: the matrix. -% -% o exception: return any errors or warnings in this structure. -% -*/ -MagickExport Image *MatrixToImage(const MatrixInfo *matrix_info, - ExceptionInfo *exception) -{ - CacheView - *image_view; - - double - max_value, - min_value, - scale_factor; - - Image - *image; - - MagickBooleanType - status; - - ssize_t - y; - - assert(matrix_info != (const MatrixInfo *) NULL); - assert(matrix_info->signature == MagickCoreSignature); - assert(exception != (ExceptionInfo *) NULL); - assert(exception->signature == MagickCoreSignature); - if (matrix_info->stride < sizeof(double)) - return((Image *) NULL); - /* - Determine range of matrix. - */ - (void) GetMatrixElement(matrix_info,0,0,&min_value); - max_value=min_value; - for (y=0; y < (ssize_t) matrix_info->rows; y++) - { - ssize_t - x; - - for (x=0; x < (ssize_t) matrix_info->columns; x++) - { - double - value; - - if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse) - continue; - if (value < min_value) - min_value=value; - else - if (value > max_value) - max_value=value; - } - } - if ((min_value == 0.0) && (max_value == 0.0)) - scale_factor=0; - else - if (min_value == max_value) - { - scale_factor=(double) QuantumRange/min_value; - min_value=0; - } - else - scale_factor=(double) QuantumRange/(max_value-min_value); - /* - Convert matrix to image. - */ - image=AcquireImage((ImageInfo *) NULL,exception); - image->columns=matrix_info->columns; - image->rows=matrix_info->rows; - image->colorspace=GRAYColorspace; - status=MagickTrue; - image_view=AcquireAuthenticCacheView(image,exception); -#if defined(MAGICKCORE_OPENMP_SUPPORT) - #pragma omp parallel for schedule(static) shared(status) \ - magick_number_threads(image,image,image->rows,2) -#endif - for (y=0; y < (ssize_t) image->rows; y++) - { - double - value; - - Quantum - *q; - - ssize_t - x; - - if (status == MagickFalse) - continue; - q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception); - if (q == (Quantum *) NULL) - { - status=MagickFalse; - continue; - } - for (x=0; x < (ssize_t) image->columns; x++) - { - if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse) - continue; - value=scale_factor*(value-min_value); - *q=ClampToQuantum(value); - q+=(ptrdiff_t) GetPixelChannels(image); - } - if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse) - status=MagickFalse; - } - image_view=DestroyCacheView(image_view); - if (status == MagickFalse) - image=DestroyImage(image); - return(image); -} - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -% N u l l M a t r i x % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% NullMatrix() sets all elements of the matrix to zero. -% -% The format of the memset method is: -% -% MagickBooleanType *NullMatrix(MatrixInfo *matrix_info) -% -% A description of each parameter follows: -% -% o matrix_info: the matrix. -% -*/ -MagickExport MagickBooleanType NullMatrix(MatrixInfo *matrix_info) -{ - ssize_t - x; - - ssize_t - count, - y; - - unsigned char - value; - - assert(matrix_info != (const MatrixInfo *) NULL); - assert(matrix_info->signature == MagickCoreSignature); - if (matrix_info->type != DiskCache) - { - (void) memset(matrix_info->elements,0,(size_t) - matrix_info->length); - return(MagickTrue); - } - value=0; - (void) lseek(matrix_info->file,0,SEEK_SET); - for (y=0; y < (ssize_t) matrix_info->rows; y++) - { - for (x=0; x < (ssize_t) matrix_info->length; x++) - { - count=write(matrix_info->file,&value,sizeof(value)); - if (count != (ssize_t) sizeof(value)) - break; - } - if (x < (ssize_t) matrix_info->length) - break; - } - return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue); -} - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -% R e l i n q u i s h M a g i c k M a t r i x % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% RelinquishMagickMatrix() frees the previously acquired matrix (array of -% pointers to arrays of doubles). -% -% The format of the RelinquishMagickMatrix method is: -% -% double **RelinquishMagickMatrix(double **matrix, -% const size_t number_rows) -% -% A description of each parameter follows: -% -% o matrix: the matrix to relinquish -% -% o number_rows: the first dimension of the acquired matrix (number of -% pointers) -% -*/ -MagickExport double **RelinquishMagickMatrix(double **matrix, - const size_t number_rows) -{ - ssize_t - i; - - if (matrix == (double **) NULL ) - return(matrix); - for (i=0; i < (ssize_t) number_rows; i++) - matrix[i]=(double *) RelinquishMagickMemory(matrix[i]); - matrix=(double **) RelinquishMagickMemory(matrix); - return(matrix); -} - -/* -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% % -% % -% S e t M a t r i x E l e m e n t % -% % -% % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% SetMatrixElement() sets the specified element in the matrix. -% -% The format of the SetMatrixElement method is: -% -% MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info, -% const ssize_t x,const ssize_t y,void *value) -% -% A description of each parameter follows: -% -% o matrix_info: the matrix columns. -% -% o x: the matrix x-offset. -% -% o y: the matrix y-offset. -% -% o value: set the matrix element to this value. -% -*/ - -MagickExport MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info, - const ssize_t x,const ssize_t y,const void *value) -{ - MagickOffsetType - count, - i; - - assert(matrix_info != (const MatrixInfo *) NULL); - assert(matrix_info->signature == MagickCoreSignature); - i=y*(MagickOffsetType) matrix_info->columns+x; - if ((i < 0) || - (((MagickSizeType) i*matrix_info->stride) >= matrix_info->length)) - return(MagickFalse); - if (matrix_info->type != DiskCache) - { - (void) memcpy((unsigned char *) matrix_info->elements+i* - (MagickOffsetType) matrix_info->stride,value,matrix_info->stride); - return(MagickTrue); - } - count=WriteMatrixElements(matrix_info,i*(MagickOffsetType) - matrix_info->stride,matrix_info->stride,(unsigned char *) value); - if (count != (MagickOffsetType) matrix_info->stride) - return(MagickFalse); - return(MagickTrue); -}