From: cristy Date: Sun, 9 Mar 2014 01:06:05 +0000 (+0000) Subject: (no commit message) X-Git-Tag: 7.0.1-0~2602 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=106efa18c0522058b3a89b2d3ce035245c452bce;p=imagemagick --- diff --git a/MagickCore/matrix.c b/MagickCore/matrix.c index 7b0dcc9fd..855f751a3 100644 --- a/MagickCore/matrix.c +++ b/MagickCore/matrix.c @@ -40,10 +40,155 @@ 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/matrix.h" -#include "MagickCore/matrix-private.h" -#include "MagickCore/pixel-private.h" #include "MagickCore/memory_.h" +#include "MagickCore/pixel-private.h" +#include "MagickCore/resource_.h" +#include "MagickCore/utility.h" + +/* + Typedef declaration. +*/ +struct _MatrixInfo +{ + CacheType + type; + + size_t + columns, + rows, + stride; + + MagickSizeType + length; + + MagickBooleanType + mapped; + + char + path[MaxTextExtent]; + + int + file; + + void + *elements; + + 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. +% +*/ +MagickExport MatrixInfo *AcquireMatrixInfo(const size_t columns, + const size_t rows,const size_t stride,ExceptionInfo *exception) +{ + MagickBooleanType + status; + + MatrixInfo + *matrix_info; + + matrix_info=(MatrixInfo *) AcquireMagickMemory(sizeof(*matrix_info)); + if (matrix_info == (MatrixInfo *) NULL) + return((MatrixInfo *) NULL); + (void) ResetMagickMemory(matrix_info,0,sizeof(*matrix_info)); + matrix_info->signature=MagickSignature; + matrix_info->columns=columns; + matrix_info->rows=rows; + matrix_info->stride=stride; + 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; + (void) AcquireMagickResource(MemoryResource,matrix_info->length); + 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=ResetMatrixInfo(matrix_info); + 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); +} /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -59,11 +204,9 @@ % 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, that can be referenced -% using the simple C-code of the form "matrix[y][x]". -% -% This matrix is typically used for perform for the GaussJordanElimination() -% method below, solving some system of simultanious equations. +% This used to generate the two dimensional matrix, and vectors required +% for the GaussJordanElimination() method below, solving some system of +% simultanious equations. % % The format of the AcquireMagickMatrix method is: % @@ -91,7 +234,7 @@ MagickExport double **AcquireMagickMatrix(const size_t number_rows, matrix=(double **) AcquireQuantumMemory(number_rows,sizeof(*matrix)); if (matrix == (double **) NULL) - return((double **) NULL); + return((double **)NULL); for (i=0; i < (ssize_t) number_rows; i++) { matrix[i]=(double *) AcquireQuantumMemory(size,sizeof(*matrix[i])); @@ -113,6 +256,67 @@ MagickExport double **AcquireMagickMatrix(const size_t number_rows, % % % % % % +% 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 == MagickSignature); + 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); + } + 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; + } + 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 % % % % % @@ -127,8 +331,8 @@ MagickExport double **AcquireMagickMatrix(const size_t number_rows, % % The format of the GaussJordanElimination method is: % -% MagickBooleanType GaussJordanElimination(double **matrix, -% double **vectors,const size_t rank,const size_t number_vectors) +% MagickBooleanType GaussJordanElimination(double **matrix,double **vectors, +% const size_t rank,const size_t number_vectors) % % A description of each parameter follows: % @@ -137,7 +341,7 @@ MagickExport double **AcquireMagickMatrix(const size_t number_rows, % o vectors: the additional matrix argumenting the matrix for row reduction. % Producing an 'array of column vectors'. % -% o rank: The size of the square matrix (both rows and columns). +% 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. @@ -150,14 +354,9 @@ MagickExport double **AcquireMagickMatrix(const size_t number_rows, % % 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'. -% It is assigned vector[column][row] where 'column' is the specific -% 'result' and 'row' is the 'values' for that answer. After processing -% the same vector array contains the 'weights' (answers) for each of the -% 'separatable' results. % % 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 that specific set of equations. +% 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 simultanious equations is to be solved from @@ -168,8 +367,8 @@ MagickExport double **AcquireMagickMatrix(const size_t number_rows, % ... % GaussJordanElimination(matrix, &coefficents, 8UL, 1UL); % -% However by specifing more 'columns' (as an 'array of vector columns'), -% you can use this function to solve multiple sets of 'separable' equations. +% However by specifing 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 coefficents, but are being @@ -177,18 +376,14 @@ MagickExport double **AcquireMagickMatrix(const size_t number_rows, % % Another example is generation of a color gradient from a set of colors % at specific coordients, such as a list x,y -> r,g,b,a -% -% See LeastSquaresAddTerms() below for such an example. +% (Reference to be added - Anthony) % % 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' (matrix -% is transposed). -% -% For details of this process see... -% http://en.wikipedia.org/wiki/Gauss-Jordan_elimination +% 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, +MagickExport MagickBooleanType GaussJordanElimination(double **matrix, double **vectors,const size_t rank,const size_t number_vectors) { #define GaussJordanSwap(x,y) \ @@ -267,7 +462,7 @@ MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix, rows[i]=row; columns[i]=column; if (matrix[column][column] == 0.0) - return(MagickFalse); /* singularity */ + return(MagickFalse); /* sigularity */ scale=PerceptibleReciprocal(matrix[column][column]); matrix[column][column]=1.0; for (j=0; j < (ssize_t) rank; j++) @@ -300,6 +495,173 @@ MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix, % % % % % % +% 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 == MagickSignature); + return(matrix_info->columns); +} + +/* +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% % +% % +% G e t M a t r i x E l e m e n t % +% % +% % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% GetMatrixElement() returns the specifed 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 size_t MagickMin(const size_t x,const size_t y) +{ + if (x < y) + return(x); + return(y); +} + +static inline ssize_t ReadMatrixInfo(const MatrixInfo *matrix_info, + const MagickOffsetType offset,const size_t length,unsigned char *buffer) +{ + register ssize_t + i; + + ssize_t + count; + +#if !defined(MAGICKCORE_HAVE_PPREAD) +#if defined(MAGICKCORE_OPENMP_SUPPORT) + #pragma omp critical (MagickCore_ReadMatrixInfo) +#endif + { + i=(-1); + if (lseek(matrix_info->file,offset,SEEK_SET) >= 0) + { +#endif + count=0; + for (i=0; i < (ssize_t) length; i+=count) + { +#if !defined(MAGICKCORE_HAVE_PPREAD) + count=read(matrix_info->file,buffer+i,MagickMin(length-i,(size_t) + SSIZE_MAX)); +#else + count=pread(matrix_info->file,buffer+i,MagickMin(length-i,(size_t) + SSIZE_MAX),offset+i); +#endif + if (count > 0) + continue; + count=0; + if (errno != EINTR) + { + i=(-1); + break; + } + } +#if !defined(MAGICKCORE_HAVE_PPREAD) + } + } +#endif + return(i); +} + +MagickExport MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info, + const ssize_t x,const ssize_t y,void *value) +{ + MagickOffsetType + i; + + ssize_t + count; + + assert(matrix_info != (const MatrixInfo *) NULL); + assert(matrix_info->signature == MagickSignature); + i=(MagickOffsetType) matrix_info->rows*x+y; + if ((i < 0) || + ((MagickSizeType) (i*matrix_info->stride) >= matrix_info->length)) + return(MagickFalse); + if (matrix_info->type != DiskCache) + { + (void) memcpy(value,matrix_info->elements+i*matrix_info->stride, + matrix_info->stride); + return(MagickTrue); + } + count=ReadMatrixInfo(matrix_info,i*matrix_info->stride,matrix_info->stride, + value); + if (count != 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 == MagickSignature); + return(matrix_info->rows); +} + +/* +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% % +% % % L e a s t S q u a r e s A d d T e r m s % % % % % @@ -322,13 +684,13 @@ MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix, % o vectors: the result vectors to add terms/results to. % % o terms: the pre-calculated terms (without the unknown coefficent -% weights) that forms the equation being added. +% 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 coefficents. +% weighted by the yet-to-be-solved coefficents. % -% o rank: the rank or size of the dimentions of the square matrix. -% Also the length of vectors, and number of terms being added. +% 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 @@ -354,10 +716,10 @@ MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix, % ... % if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) { % c0 = vectors[0][0]; -% c2 = vectors[0][1]; %* weights to calculate u from any given x,y *% +% c2 = vectors[0][1]; % c4 = vectors[0][2]; % c1 = vectors[1][0]; -% c3 = vectors[1][1]; %* weights for calculate v from any given x,y *% +% c3 = vectors[1][1]; % c5 = vectors[1][2]; % } % else @@ -365,10 +727,8 @@ MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix, % RelinquishMagickMatrix(matrix,3UL); % RelinquishMagickMatrix(vectors,2UL); % -% More examples can be found in "distort.c" -% */ -MagickPrivate void LeastSquaresAddTerms(double **matrix,double **vectors, +MagickExport void LeastSquaresAddTerms(double **matrix,double **vectors, const double *terms,const double *results,const size_t rank, const size_t number_vectors) { @@ -383,6 +743,7 @@ MagickPrivate void LeastSquaresAddTerms(double **matrix,double **vectors, for (i=0; i < (ssize_t) number_vectors; i++) vectors[i][j]+=results[i]*terms[j]; } + return; } /* @@ -425,3 +786,163 @@ MagickExport double **RelinquishMagickMatrix(double **matrix, matrix=(double **) RelinquishMagickMemory(matrix); return(matrix); } + +/* +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% % +% % +% R e s e t M a g i c k M e m o r y % +% % +% % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% ResetMagickMemory() sets all elements of the matrix to zero. +% +% The format of the ResetMagickMemory method is: +% +% MagickBoolenType *ResetMatrixInfo(MatrixInfo *matrix_info) +% +% A description of each parameter follows: +% +% o matrix_info: the matrix. +% +*/ +MagickExport MagickBooleanType ResetMatrixInfo(MatrixInfo *matrix_info) +{ + register ssize_t + x; + + ssize_t + count, + y; + + unsigned char + value; + + assert(matrix_info != (const MatrixInfo *) NULL); + assert(matrix_info->signature == MagickSignature); + if (matrix_info->type != DiskCache) + { + (void) ResetMagickMemory(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); +} + +/* +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% % +% % +% S e t M a t r i x E l e m e n t % +% % +% % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% SetMatrixElement() sets the specifed 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. +% +*/ + +static inline ssize_t WriteMatrixInfo(const MatrixInfo *matrix_info, + const MagickOffsetType offset,const size_t length,const unsigned char *buffer) +{ + register ssize_t + i; + + ssize_t + count; + + i=0; +#if !defined(MAGICKCORE_HAVE_PWRITE) +#if defined(MAGICKCORE_OPENMP_SUPPORT) + #pragma omp critical (MagickCore_WriteMatrixInfo) +#endif + { + if (lseek(matrix_info->file,offset,SEEK_SET) >= 0) + { +#endif + count=0; + for (i=0; i < (ssize_t) length; i+=count) + { +#if !defined(MAGICKCORE_HAVE_PWRITE) + count=write(matrix_info->file,buffer+i,MagickMin(length-i,(size_t) + SSIZE_MAX)); +#else + count=pwrite(matrix_info->file,buffer+i,MagickMin(length-i,(size_t) + SSIZE_MAX),offset+i); +#endif + if (count > 0) + continue; + count=0; + if (errno != EINTR) + { + i=(-1); + break; + } + } +#if !defined(MAGICKCORE_HAVE_PWRITE) + } + } +#endif + return(i); +} + +MagickExport MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info, + const ssize_t x,const ssize_t y,const void *value) +{ + MagickOffsetType + i; + + ssize_t + count; + + assert(matrix_info != (const MatrixInfo *) NULL); + assert(matrix_info->signature == MagickSignature); + i=(MagickOffsetType) matrix_info->rows*x+y; + if ((i < 0) || + ((MagickSizeType) (i*matrix_info->stride) >= matrix_info->length)) + return(MagickFalse); + if (matrix_info->type != DiskCache) + { + (void) memcpy(matrix_info->elements+i*matrix_info->stride,value, + matrix_info->stride); + return(MagickTrue); + } + count=WriteMatrixInfo(matrix_info,i*matrix_info->stride,matrix_info->stride, + value); + if (count != (ssize_t) matrix_info->stride) + return(MagickFalse); + return(MagickTrue); +} diff --git a/MagickCore/matrix.h b/MagickCore/matrix.h index 61d99c530..6da31b81c 100644 --- a/MagickCore/matrix.h +++ b/MagickCore/matrix.h @@ -13,7 +13,7 @@ See the License for the specific language governing permissions and limitations under the License. - MagickCore graphic resample methods. + MagickCore matrix methods. */ #ifndef _MAGICKCORE_MATRIX_H #define _MAGICKCORE_MATRIX_H @@ -22,10 +22,31 @@ extern "C" { #endif +typedef struct _MatrixInfo + MatrixInfo; + extern MagickExport double **AcquireMagickMatrix(const size_t,const size_t), **RelinquishMagickMatrix(double **,const size_t); +extern MagickExport MagickBooleanType + GaussJordanElimination(double **,double **,const size_t,const size_t), + GetMatrixElement(const MatrixInfo *,const ssize_t,const ssize_t,void *), + ResetMatrixInfo(MatrixInfo *), + SetMatrixElement(const MatrixInfo *,const ssize_t,const ssize_t,const void *); + +MagickExport MatrixInfo + *AcquireMatrixInfo(const size_t,const size_t,const size_t,ExceptionInfo *), + *DestroyMatrixInfo(MatrixInfo *); + +MagickExport size_t + GetMatrixColumns(const MatrixInfo *), + GetMatrixRows(const MatrixInfo *); + +extern MagickExport void + LeastSquaresAddTerms(double **,double **,const double *,const double *, + const size_t,const size_t); + #if defined(__cplusplus) || defined(c_plusplus) } #endif diff --git a/MagickCore/shear.c b/MagickCore/shear.c index 62ab2354e..59125fc2d 100644 --- a/MagickCore/shear.c +++ b/MagickCore/shear.c @@ -63,6 +63,7 @@ #include "MagickCore/geometry.h" #include "MagickCore/image.h" #include "MagickCore/image-private.h" +#include "MagickCore/matrix.h" #include "MagickCore/memory_.h" #include "MagickCore/list.h" #include "MagickCore/monitor.h" @@ -217,335 +218,27 @@ static MagickBooleanType CropToFitImage(Image **image, % */ -typedef struct _RadonInfo +static void RadonProjection(const Image *image,MatrixInfo *source_cells, + MatrixInfo *destination_cells,const ssize_t sign,size_t *projection) { - CacheType - type; - - size_t - width, - height; - - MagickSizeType - length; - - MagickBooleanType - mapped; - - char - path[MaxTextExtent]; - - int - file; - - unsigned short - *cells; -} RadonInfo; - -static RadonInfo *DestroyRadonInfo(RadonInfo *radon_info) -{ - assert(radon_info != (RadonInfo *) NULL); - switch (radon_info->type) - { - case MemoryCache: - { - if (radon_info->mapped == MagickFalse) - radon_info->cells=(unsigned short *) RelinquishMagickMemory( - radon_info->cells); - else - { - (void) UnmapBlob(radon_info->cells,(size_t) radon_info->length); - radon_info->cells=(unsigned short *) NULL; - } - RelinquishMagickResource(MemoryResource,radon_info->length); - break; - } - case MapCache: - { - (void) UnmapBlob(radon_info->cells,(size_t) radon_info->length); - radon_info->cells=(unsigned short *) NULL; - RelinquishMagickResource(MapResource,radon_info->length); - } - case DiskCache: - { - if (radon_info->file != -1) - (void) close(radon_info->file); - (void) RelinquishUniqueFileResource(radon_info->path); - RelinquishMagickResource(DiskResource,radon_info->length); - break; - } - default: - break; - } - return((RadonInfo *) RelinquishMagickMemory(radon_info)); -} - -static MagickBooleanType ResetRadonCells(RadonInfo *radon_info) -{ - register ssize_t - x; - - ssize_t - count, - y; - - unsigned short - value; - - if (radon_info->type != DiskCache) - { - (void) ResetMagickMemory(radon_info->cells,0,(size_t) radon_info->length); - return(MagickTrue); - } - value=0; - (void) lseek(radon_info->file,0,SEEK_SET); - for (y=0; y < (ssize_t) radon_info->height; y++) - { - for (x=0; x < (ssize_t) radon_info->width; x++) - { - count=write(radon_info->file,&value,sizeof(*radon_info->cells)); - if (count != (ssize_t) sizeof(*radon_info->cells)) - break; - } - if (x < (ssize_t) radon_info->width) - break; - } - return(y < (ssize_t) radon_info->height ? MagickFalse : MagickTrue); -} - -static RadonInfo *AcquireRadonInfo(const Image *image,const size_t width, - const size_t height,ExceptionInfo *exception) -{ - MagickBooleanType - status; - - RadonInfo - *radon_info; - - radon_info=(RadonInfo *) AcquireMagickMemory(sizeof(*radon_info)); - if (radon_info == (RadonInfo *) NULL) - return((RadonInfo *) NULL); - (void) ResetMagickMemory(radon_info,0,sizeof(*radon_info)); - radon_info->width=width; - radon_info->height=height; - radon_info->length=(MagickSizeType) width*height*sizeof(*radon_info->cells); - radon_info->type=MemoryCache; - status=AcquireMagickResource(AreaResource,radon_info->length); - if ((status != MagickFalse) && - (radon_info->length == (MagickSizeType) ((size_t) radon_info->length))) - { - status=AcquireMagickResource(MemoryResource,radon_info->length); - if (status != MagickFalse) - { - radon_info->mapped=MagickFalse; - radon_info->cells=(unsigned short *) AcquireMagickMemory((size_t) - radon_info->length); - if (radon_info->cells == (unsigned short *) NULL) - { - radon_info->mapped=MagickTrue; - radon_info->cells=(unsigned short *) MapBlob(-1,IOMode,0,(size_t) - radon_info->length); - } - if (radon_info->cells == (unsigned short *) NULL) - RelinquishMagickResource(MemoryResource,radon_info->length); - } - } - radon_info->file=(-1); - if (radon_info->cells == (unsigned short *) NULL) - { - status=AcquireMagickResource(DiskResource,radon_info->length); - if (status == MagickFalse) - { - (void) ThrowMagickException(exception,GetMagickModule(),CacheError, - "CacheResourcesExhausted","`%s'",image->filename); - return(DestroyRadonInfo(radon_info)); - } - radon_info->type=DiskCache; - (void) AcquireMagickResource(MemoryResource,radon_info->length); - radon_info->file=AcquireUniqueFileResource(radon_info->path); - if (radon_info->file == -1) - return(DestroyRadonInfo(radon_info)); - status=AcquireMagickResource(MapResource,radon_info->length); - if (status != MagickFalse) - { - status=ResetRadonCells(radon_info); - if (status != MagickFalse) - { - radon_info->cells=(unsigned short *) MapBlob(radon_info->file, - IOMode,0,(size_t) radon_info->length); - if (radon_info->cells != (unsigned short *) NULL) - radon_info->type=MapCache; - else - RelinquishMagickResource(MapResource,radon_info->length); - } - } - } - return(radon_info); -} - -static inline size_t MagickMin(const size_t x,const size_t y) -{ - if (x < y) - return(x); - return(y); -} - -static inline ssize_t ReadRadonCell(const RadonInfo *radon_info, - const MagickOffsetType offset,const size_t length,unsigned char *buffer) -{ - register ssize_t - i; - - ssize_t - count; - -#if !defined(MAGICKCORE_HAVE_PPREAD) -#if defined(MAGICKCORE_OPENMP_SUPPORT) - #pragma omp critical (MagickCore_ReadRadonCell) -#endif - { - i=(-1); - if (lseek(radon_info->file,offset,SEEK_SET) >= 0) - { -#endif - count=0; - for (i=0; i < (ssize_t) length; i+=count) - { -#if !defined(MAGICKCORE_HAVE_PPREAD) - count=read(radon_info->file,buffer+i,MagickMin(length-i,(size_t) - SSIZE_MAX)); -#else - count=pread(radon_info->file,buffer+i,MagickMin(length-i,(size_t) - SSIZE_MAX),offset+i); -#endif - if (count > 0) - continue; - count=0; - if (errno != EINTR) - { - i=(-1); - break; - } - } -#if !defined(MAGICKCORE_HAVE_PPREAD) - } - } -#endif - return(i); -} - -static inline ssize_t WriteRadonCell(const RadonInfo *radon_info, - const MagickOffsetType offset,const size_t length,const unsigned char *buffer) -{ - register ssize_t - i; - - ssize_t - count; - - i=0; -#if !defined(MAGICKCORE_HAVE_PWRITE) -#if defined(MAGICKCORE_OPENMP_SUPPORT) - #pragma omp critical (MagickCore_WriteRadonCell) -#endif - { - if (lseek(radon_info->file,offset,SEEK_SET) >= 0) - { -#endif - count=0; - for (i=0; i < (ssize_t) length; i+=count) - { -#if !defined(MAGICKCORE_HAVE_PWRITE) - count=write(radon_info->file,buffer+i,MagickMin(length-i,(size_t) - SSIZE_MAX)); -#else - count=pwrite(radon_info->file,buffer+i,MagickMin(length-i,(size_t) - SSIZE_MAX),offset+i); -#endif - if (count > 0) - continue; - count=0; - if (errno != EINTR) - { - i=(-1); - break; - } - } -#if !defined(MAGICKCORE_HAVE_PWRITE) - } - } -#endif - return(i); -} - -static inline unsigned short GetRadonCell(const RadonInfo *radon_info, - const ssize_t x,const ssize_t y) -{ - MagickOffsetType - i; - - unsigned short - value; - - i=(MagickOffsetType) radon_info->height*x+y; - if ((i < 0) || - ((MagickSizeType) (i*sizeof(*radon_info->cells)) >= radon_info->length)) - return(0); - if (radon_info->type != DiskCache) - return(radon_info->cells[i]); - value=0; - (void) ReadRadonCell(radon_info,i*sizeof(*radon_info->cells), - sizeof(*radon_info->cells),(unsigned char *) &value); - return(value); -} - -static inline MagickBooleanType SetRadonCell(const RadonInfo *radon_info, - const ssize_t x,const ssize_t y,const unsigned short value) -{ - MagickOffsetType - i; - - ssize_t - count; - - i=(MagickOffsetType) radon_info->height*x+y; - if ((i < 0) || - ((MagickSizeType) (i*sizeof(*radon_info->cells)) >= radon_info->length)) - return(MagickFalse); - if (radon_info->type != DiskCache) - { - radon_info->cells[i]=value; - return(MagickTrue); - } - count=WriteRadonCell(radon_info,i*sizeof(*radon_info->cells), - sizeof(*radon_info->cells),(const unsigned char *) &value); - if (count != (ssize_t) sizeof(*radon_info->cells)) - return(MagickFalse); - return(MagickTrue); -} - -static void RadonProjection(const Image *image,RadonInfo *source_cells, - RadonInfo *destination_cells,const ssize_t sign,size_t *projection) -{ - RadonInfo + MatrixInfo *swap; - register ssize_t - x; - - register RadonInfo + register MatrixInfo *p, *q; + register ssize_t + x; + size_t step; - assert(image != (Image *) NULL); p=source_cells; q=destination_cells; - for (step=1; step < p->width; step*=2) + for (step=1; step < GetMatrixColumns(p); step*=2) { - for (x=0; x < (ssize_t) p->width; x+=2*(ssize_t) step) + for (x=0; x < (ssize_t) GetMatrixColumns(p); x+=2*(ssize_t) step) { register ssize_t i; @@ -554,30 +247,46 @@ static void RadonProjection(const Image *image,RadonInfo *source_cells, y; unsigned short - cell; + cell, + neighbor; for (i=0; i < (ssize_t) step; i++) { - for (y=0; y < (ssize_t) (p->height-i-1); y++) + for (y=0; y < (ssize_t) (GetMatrixRows(p)-i-1); y++) { - cell=GetRadonCell(p,x+i,y); - (void) SetRadonCell(q,x+2*i,y,cell+GetRadonCell(p,x+i+(ssize_t) - step,y+i)); - (void) SetRadonCell(q,x+2*i+1,y,cell+GetRadonCell(p,x+i+(ssize_t) - step,y+i+1)); + if (GetMatrixElement(p,x+i,y,&cell) == MagickFalse) + continue; + if (GetMatrixElement(p,x+i+step,y+i,&neighbor) == MagickFalse) + continue; + neighbor+=cell; + if (SetMatrixElement(q,x+2*i,y,&neighbor) == MagickFalse) + continue; + if (GetMatrixElement(p,x+i+step,y+i+1,&neighbor) == MagickFalse) + continue; + neighbor+=cell; + if (SetMatrixElement(q,x+2*i+1,y,&neighbor) == MagickFalse) + continue; } - for ( ; y < (ssize_t) (p->height-i); y++) + for ( ; y < (ssize_t) (GetMatrixRows(p)-i); y++) { - cell=GetRadonCell(p,x+i,y); - (void) SetRadonCell(q,x+2*i,y,cell+GetRadonCell(p,x+i+(ssize_t) step, - y+i)); - (void) SetRadonCell(q,x+2*i+1,y,cell); + if (GetMatrixElement(p,x+i,y,&cell) == MagickFalse) + continue; + if (GetMatrixElement(p,x+i+step,y+i,&neighbor) == MagickFalse) + continue; + neighbor+=cell; + if (SetMatrixElement(q,x+2*i,y,&neighbor) == MagickFalse) + continue; + if (SetMatrixElement(q,x+2*i+1,y,&cell) == MagickFalse) + continue; } - for ( ; y < (ssize_t) p->height; y++) + for ( ; y < (ssize_t) GetMatrixRows(p); y++) { - cell=GetRadonCell(p,x+i,y); - (void) SetRadonCell(q,x+2*i,y,cell); - (void) SetRadonCell(q,x+2*i+1,y,cell); + if (GetMatrixElement(p,x+i,y,&cell) == MagickFalse) + continue; + if (SetMatrixElement(q,x+2*i,y,&cell) == MagickFalse) + continue; + if (SetMatrixElement(q,x+2*i+1,y,&cell) == MagickFalse) + continue; } } } @@ -589,7 +298,7 @@ static void RadonProjection(const Image *image,RadonInfo *source_cells, #pragma omp parallel for schedule(static,4) \ magick_threads(image,image,1,1) #endif - for (x=0; x < (ssize_t) p->width; x++) + for (x=0; x < (ssize_t) GetMatrixColumns(p); x++) { register ssize_t y; @@ -598,15 +307,23 @@ static void RadonProjection(const Image *image,RadonInfo *source_cells, sum; sum=0; - for (y=0; y < (ssize_t) (p->height-1); y++) + for (y=0; y < (ssize_t) (GetMatrixRows(p)-1); y++) { ssize_t delta; - delta=GetRadonCell(p,x,y)-(ssize_t) GetRadonCell(p,x,y+1); + unsigned short + cell, + neighbor; + + if (GetMatrixElement(p,x,y,&cell) == MagickFalse) + continue; + if (GetMatrixElement(p,x,y+1,&neighbor) == MagickFalse) + continue; + delta=(ssize_t) cell-(ssize_t) neighbor; sum+=delta*delta; } - projection[p->width+sign*x-1]=sum; + projection[GetMatrixColumns(p)+sign*x-1]=sum; } } @@ -616,13 +333,13 @@ static MagickBooleanType RadonTransform(const Image *image, CacheView *image_view; - MagickBooleanType - status; - - RadonInfo + MatrixInfo *destination_cells, *source_cells; + MagickBooleanType + status; + register ssize_t i; @@ -640,21 +357,23 @@ static MagickBooleanType RadonTransform(const Image *image, bits[256]; for (width=1; width < ((image->columns+7)/8); width<<=1) ; - source_cells=AcquireRadonInfo(image,width,image->rows,exception); - destination_cells=AcquireRadonInfo(image,width,image->rows,exception); - if ((source_cells == (RadonInfo *) NULL) || - (destination_cells == (RadonInfo *) NULL)) + source_cells=AcquireMatrixInfo(width,image->rows,sizeof(unsigned short), + exception); + destination_cells=AcquireMatrixInfo(width,image->rows,sizeof(unsigned short), + exception); + if ((source_cells == (MatrixInfo *) NULL) || + (destination_cells == (MatrixInfo *) NULL)) { - if (destination_cells != (RadonInfo *) NULL) - destination_cells=DestroyRadonInfo(destination_cells); - if (source_cells != (RadonInfo *) NULL) - source_cells=DestroyRadonInfo(source_cells); + if (destination_cells != (MatrixInfo *) NULL) + destination_cells=DestroyMatrixInfo(destination_cells); + if (source_cells != (MatrixInfo *) NULL) + source_cells=DestroyMatrixInfo(source_cells); return(MagickFalse); } - if (ResetRadonCells(source_cells) == MagickFalse) + if (ResetMatrixInfo(source_cells) == MagickFalse) { - destination_cells=DestroyRadonInfo(destination_cells); - source_cells=DestroyRadonInfo(source_cells); + destination_cells=DestroyMatrixInfo(destination_cells); + source_cells=DestroyMatrixInfo(source_cells); return(MagickFalse); } for (i=0; i < 256; i++) @@ -668,7 +387,7 @@ static MagickBooleanType RadonTransform(const Image *image, image_view=AcquireVirtualCacheView(image,exception); #if defined(MAGICKCORE_OPENMP_SUPPORT) #pragma omp parallel for schedule(static,4) shared(status) \ - magick_threads(image,image,image->rows,1) + magick_threads(image,image,1,1) #endif for (y=0; y < (ssize_t) image->rows; y++) { @@ -683,6 +402,9 @@ static MagickBooleanType RadonTransform(const Image *image, bit, byte; + unsigned short + value; + if (status == MagickFalse) continue; p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); @@ -697,12 +419,15 @@ static MagickBooleanType RadonTransform(const Image *image, for (x=0; x < (ssize_t) image->columns; x++) { byte<<=1; - if (GetPixelIntensity(image,p) < threshold) + if (((MagickRealType) GetPixelRed(image,p) < threshold) || + ((MagickRealType) GetPixelGreen(image,p) < threshold) || + ((MagickRealType) GetPixelBlue(image,p) < threshold)) byte|=0x01; bit++; if (bit == 8) { - (void) SetRadonCell(source_cells,--i,y,bits[byte]); + value=bits[byte]; + (void) SetMatrixElement(source_cells,--i,y,&value); bit=0; byte=0; } @@ -711,14 +436,15 @@ static MagickBooleanType RadonTransform(const Image *image, if (bit != 0) { byte<<=(8-bit); - (void) SetRadonCell(source_cells,--i,y,bits[byte]); + value=bits[byte]; + (void) SetMatrixElement(source_cells,--i,y,&value); } } RadonProjection(image,source_cells,destination_cells,-1,projection); - (void) ResetRadonCells(source_cells); + (void) ResetMatrixInfo(source_cells); #if defined(MAGICKCORE_OPENMP_SUPPORT) #pragma omp parallel for schedule(static,4) shared(status) \ - magick_threads(image,image,1,1) + magick_threads(image,image,image->rows,1) #endif for (y=0; y < (ssize_t) image->rows; y++) { @@ -733,6 +459,9 @@ static MagickBooleanType RadonTransform(const Image *image, bit, byte; + unsigned short + value; + if (status == MagickFalse) continue; p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); @@ -747,12 +476,15 @@ static MagickBooleanType RadonTransform(const Image *image, for (x=0; x < (ssize_t) image->columns; x++) { byte<<=1; - if (GetPixelIntensity(image,p) < threshold) + if (((MagickRealType) GetPixelRed(image,p) < threshold) || + ((MagickRealType) GetPixelGreen(image,p) < threshold) || + ((MagickRealType) GetPixelBlue(image,p) < threshold)) byte|=0x01; bit++; if (bit == 8) { - (void) SetRadonCell(source_cells,i++,y,bits[byte]); + value=bits[byte]; + (void) SetMatrixElement(source_cells,i++,y,&value); bit=0; byte=0; } @@ -761,13 +493,14 @@ static MagickBooleanType RadonTransform(const Image *image, if (bit != 0) { byte<<=(8-bit); - (void) SetRadonCell(source_cells,i++,y,bits[byte]); + value=bits[byte]; + (void) SetMatrixElement(source_cells,i++,y,&value); } } RadonProjection(image,source_cells,destination_cells,1,projection); image_view=DestroyCacheView(image_view); - destination_cells=DestroyRadonInfo(destination_cells); - source_cells=DestroyRadonInfo(source_cells); + destination_cells=DestroyMatrixInfo(destination_cells); + source_cells=DestroyMatrixInfo(source_cells); return(MagickTrue); }