42#include "MagickCore/studio.h"
43#include "MagickCore/blob.h"
44#include "MagickCore/blob-private.h"
45#include "MagickCore/cache.h"
46#include "MagickCore/exception.h"
47#include "MagickCore/exception-private.h"
48#include "MagickCore/image-private.h"
49#include "MagickCore/matrix.h"
50#include "MagickCore/matrix-private.h"
51#include "MagickCore/memory_.h"
52#include "MagickCore/nt-base-private.h"
53#include "MagickCore/pixel-accessor.h"
54#include "MagickCore/resource_.h"
55#include "MagickCore/semaphore.h"
56#include "MagickCore/thread-private.h"
57#include "MagickCore/utility.h"
58#include "MagickCore/utility-private.h"
81 path[MagickPathExtent];
127static void MatrixSignalHandler(
int magick_unused(status))
129 magick_unreferenced(status);
130 ThrowFatalException(CacheFatalError,
"UnableToExtendMatrixCache");
134static inline MagickOffsetType WriteMatrixElements(
135 const MatrixInfo *magick_restrict matrix_info,
const MagickOffsetType offset,
136 const MagickSizeType length,
const unsigned char *magick_restrict buffer)
144#if !defined(MAGICKCORE_HAVE_PWRITE)
145 LockSemaphoreInfo(matrix_info->semaphore);
146 if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
148 UnlockSemaphoreInfo(matrix_info->semaphore);
149 return((MagickOffsetType) -1);
153 for (i=0; i < (MagickOffsetType) length; i+=count)
155#if !defined(MAGICKCORE_HAVE_PWRITE)
156 count=write(matrix_info->file,buffer+i,(
size_t) MagickMin(length-
157 (MagickSizeType) i,(MagickSizeType) MagickMaxBufferExtent));
159 count=pwrite(matrix_info->file,buffer+i,(
size_t) MagickMin(length-
160 (MagickSizeType) i,(MagickSizeType) MagickMaxBufferExtent),offset+i);
169#if !defined(MAGICKCORE_HAVE_PWRITE)
170 UnlockSemaphoreInfo(matrix_info->semaphore);
175static MagickBooleanType SetMatrixExtent(
176 MatrixInfo *magick_restrict matrix_info,MagickSizeType length)
183 if (length != (MagickSizeType) ((MagickOffsetType) length))
185 offset=(MagickOffsetType) lseek(matrix_info->file,0,SEEK_END);
188 if ((MagickSizeType) offset >= length)
190 extent=(MagickOffsetType) length-1;
191 count=WriteMatrixElements(matrix_info,extent,1,(
const unsigned char *)
"");
192#if defined(MAGICKCORE_HAVE_POSIX_FALLOCATE)
193 if (matrix_info->synchronize != MagickFalse)
194 (void) posix_fallocate(matrix_info->file,offset+1,extent-offset);
197 (void) signal(SIGBUS,MatrixSignalHandler);
199 return(count != (MagickOffsetType) 1 ? MagickFalse : MagickTrue);
202MagickExport
MatrixInfo *AcquireMatrixInfo(
const size_t columns,
203 const size_t rows,
const size_t stride,
ExceptionInfo *exception)
214 matrix_info=(
MatrixInfo *) AcquireMagickMemory(
sizeof(*matrix_info));
217 (void) memset(matrix_info,0,
sizeof(*matrix_info));
218 matrix_info->signature=MagickCoreSignature;
219 matrix_info->columns=columns;
220 matrix_info->rows=rows;
221 matrix_info->stride=stride;
222 matrix_info->semaphore=AcquireSemaphoreInfo();
223 synchronize=GetEnvironmentValue(
"MAGICK_SYNCHRONIZE");
224 if (synchronize != (
const char *) NULL)
226 matrix_info->synchronize=IsStringTrue(synchronize);
227 synchronize=DestroyString(synchronize);
229 matrix_info->length=(MagickSizeType) columns*rows*stride;
230 if (matrix_info->columns != (
size_t) (matrix_info->length/rows/stride))
232 (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
233 "CacheResourcesExhausted",
"`%s'",
"matrix cache");
234 return(DestroyMatrixInfo(matrix_info));
236 matrix_info->type=MemoryCache;
237 status=AcquireMagickResource(AreaResource,matrix_info->length);
238 if ((status != MagickFalse) &&
239 (matrix_info->length == (MagickSizeType) ((
size_t) matrix_info->length)))
241 status=AcquireMagickResource(MemoryResource,matrix_info->length);
242 if (status != MagickFalse)
244 matrix_info->mapped=MagickFalse;
245 matrix_info->elements=AcquireMagickMemory((
size_t)
246 matrix_info->length);
247 if (matrix_info->elements == NULL)
249 matrix_info->mapped=MagickTrue;
250 matrix_info->elements=MapBlob(-1,IOMode,0,(
size_t)
251 matrix_info->length);
253 if (matrix_info->elements == (
unsigned short *) NULL)
254 RelinquishMagickResource(MemoryResource,matrix_info->length);
257 matrix_info->file=(-1);
258 if (matrix_info->elements == (
unsigned short *) NULL)
260 status=AcquireMagickResource(DiskResource,matrix_info->length);
261 if (status == MagickFalse)
263 (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
264 "CacheResourcesExhausted",
"`%s'",
"matrix cache");
265 return(DestroyMatrixInfo(matrix_info));
267 matrix_info->type=DiskCache;
268 matrix_info->file=AcquireUniqueFileResource(matrix_info->path);
269 if (matrix_info->file == -1)
270 return(DestroyMatrixInfo(matrix_info));
271 status=AcquireMagickResource(MapResource,matrix_info->length);
272 if (status != MagickFalse)
274 status=SetMatrixExtent(matrix_info,matrix_info->length);
275 if (status != MagickFalse)
276 matrix_info->elements=(
void *) MapBlob(matrix_info->file,IOMode,0,
277 (
size_t) matrix_info->length);
278 if (matrix_info->elements != NULL)
279 matrix_info->type=MapCache;
281 RelinquishMagickResource(MapResource,matrix_info->length);
319MagickExport
double **AcquireMagickMatrix(
const size_t number_rows,
329 matrix=(
double **) AcquireQuantumMemory(number_rows,
sizeof(*matrix));
330 if (matrix == (
double **) NULL)
331 return((
double **) NULL);
332 for (i=0; i < (ssize_t) number_rows; i++)
334 matrix[i]=(
double *) AcquireQuantumMemory(size,
sizeof(*matrix[i]));
335 if (matrix[i] == (
double *) NULL)
337 for (j=0; j < i; j++)
338 matrix[j]=(
double *) RelinquishMagickMemory(matrix[j]);
339 matrix=(
double **) RelinquishMagickMemory(matrix);
340 return((
double **) NULL);
342 for (j=0; j < (ssize_t) size; j++)
374 assert(matrix_info->signature == MagickCoreSignature);
375 LockSemaphoreInfo(matrix_info->semaphore);
376 switch (matrix_info->type)
380 if (matrix_info->mapped == MagickFalse)
381 matrix_info->elements=RelinquishMagickMemory(matrix_info->elements);
384 (void) UnmapBlob(matrix_info->elements,(
size_t) matrix_info->length);
385 matrix_info->elements=(
unsigned short *) NULL;
387 RelinquishMagickResource(MemoryResource,matrix_info->length);
392 (void) UnmapBlob(matrix_info->elements,(
size_t) matrix_info->length);
393 matrix_info->elements=NULL;
394 RelinquishMagickResource(MapResource,matrix_info->length);
399 if (matrix_info->file != -1)
400 (void) close_utf8(matrix_info->file);
401 (void) RelinquishUniqueFileResource(matrix_info->path);
402 RelinquishMagickResource(DiskResource,matrix_info->length);
408 UnlockSemaphoreInfo(matrix_info->semaphore);
409 RelinquishSemaphoreInfo(&matrix_info->semaphore);
410 return((
MatrixInfo *) RelinquishMagickMemory(matrix_info));
435MagickExport
size_t GetMatrixColumns(
const MatrixInfo *matrix_info)
438 assert(matrix_info->signature == MagickCoreSignature);
439 return(matrix_info->columns);
472static inline ssize_t EdgeX(
const ssize_t x,
const size_t columns)
476 if (x >= (ssize_t) columns)
477 return((ssize_t) (columns-1));
481static inline ssize_t EdgeY(
const ssize_t y,
const size_t rows)
485 if (y >= (ssize_t) rows)
486 return((ssize_t) (rows-1));
490static inline MagickOffsetType ReadMatrixElements(
491 const MatrixInfo *magick_restrict matrix_info,
const MagickOffsetType offset,
492 const MagickSizeType length,
unsigned char *magick_restrict buffer)
500#if !defined(MAGICKCORE_HAVE_PREAD)
501 LockSemaphoreInfo(matrix_info->semaphore);
502 if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
504 UnlockSemaphoreInfo(matrix_info->semaphore);
505 return((MagickOffsetType) -1);
509 for (i=0; i < (MagickOffsetType) length; i+=count)
511#if !defined(MAGICKCORE_HAVE_PREAD)
512 count=read(matrix_info->file,buffer+i,(
size_t) MagickMin(length-i,
513 (MagickSizeType) MagickMaxBufferExtent));
515 count=pread(matrix_info->file,buffer+i,(
size_t) MagickMin(length-
516 (MagickSizeType) i,(MagickSizeType) MagickMaxBufferExtent),offset+i);
525#if !defined(MAGICKCORE_HAVE_PREAD)
526 UnlockSemaphoreInfo(matrix_info->semaphore);
531MagickExport MagickBooleanType GetMatrixElement(
const MatrixInfo *matrix_info,
532 const ssize_t x,
const ssize_t y,
void *value)
538 assert(matrix_info != (
const MatrixInfo *) NULL);
539 assert(matrix_info->signature == MagickCoreSignature);
540 i=EdgeY(y,matrix_info->rows)*(MagickOffsetType) matrix_info->columns+
541 EdgeX(x,matrix_info->columns);
542 if (matrix_info->type != DiskCache)
544 (void) memcpy(value,(
unsigned char *) matrix_info->elements+i*
545 (MagickOffsetType) matrix_info->stride,matrix_info->stride);
548 count=ReadMatrixElements(matrix_info,i*(MagickOffsetType) matrix_info->stride,
549 matrix_info->stride,(
unsigned char *) value);
550 if (count != (MagickOffsetType) matrix_info->stride)
577MagickExport
size_t GetMatrixRows(
const MatrixInfo *matrix_info)
579 assert(matrix_info != (
const MatrixInfo *) NULL);
580 assert(matrix_info->signature == MagickCoreSignature);
581 return(matrix_info->rows);
655MagickPrivate
void LeastSquaresAddTerms(
double **matrix,
double **vectors,
656 const double *terms,
const double *results,
const size_t rank,
657 const size_t number_vectors)
663 for (j=0; j < (ssize_t) rank; j++)
665 for (i=0; i < (ssize_t) rank; i++)
666 matrix[i][j]+=terms[i]*terms[j];
667 for (i=0; i < (ssize_t) number_vectors; i++)
668 vectors[i][j]+=results[i]*terms[j];
740MagickPrivate MagickBooleanType GaussJordanElimination(
double **matrix,
741 double **vectors,
const size_t rank,
const size_t number_vectors)
743#define GaussJordanSwap(x,y) \
749#define GaussJordanSwapLD(x,y) \
751 long double temp = (x); \
755#define ThrowGaussJordanException() \
757 for (i=0; i < (ssize_t) rank; i++) \
758 hp_matrix[i]=(long double *) RelinquishMagickMemory(hp_matrix[i]); \
759 hp_matrix=(long double **) RelinquishMagickMemory(hp_matrix); \
760 if (pivots != (ssize_t *) NULL) \
761 pivots=(ssize_t *) RelinquishMagickMemory(pivots); \
762 if (rows != (ssize_t *) NULL) \
763 rows=(ssize_t *) RelinquishMagickMemory(rows); \
764 if (columns != (ssize_t *) NULL) \
765 columns=(ssize_t *) RelinquishMagickMemory(columns); \
766 return(MagickFalse); \
770 **hp_matrix = (
long double **) NULL,
775 *columns = (ssize_t *) NULL,
778 *pivots = (ssize_t *) NULL,
780 *rows = (ssize_t *) NULL;
785 hp_matrix=(
long double **) AcquireQuantumMemory(rank,
sizeof(*hp_matrix));
786 if (hp_matrix == (
long double **) NULL)
788 for (i=0; i < (ssize_t) rank; i++)
790 hp_matrix[i]=(
long double *) AcquireQuantumMemory(rank,
791 sizeof(*hp_matrix[i]));
792 if (hp_matrix[i] == (
long double *) NULL)
793 ThrowGaussJordanException();
794 for (j=0; j < (ssize_t) rank; j++)
795 hp_matrix[i][j]=(
long double)matrix[i][j];
797 columns=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*columns));
798 rows=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*rows));
799 pivots=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*pivots));
800 if ((columns == (ssize_t *) NULL) || (rows == (ssize_t *) NULL) ||
801 (pivots == (ssize_t *) NULL))
802 ThrowGaussJordanException();
803 (void) memset(columns,0,rank*
sizeof(*columns));
804 (void) memset(rows,0,rank*
sizeof(*rows));
805 (void) memset(pivots,0,rank*
sizeof(*pivots));
806 for (i=0; i < (ssize_t) rank; i++)
820 for (j=0; j < (ssize_t) rank; j++)
822 for (k=0; k < (ssize_t) rank; k++)
823 if ((pivots[k] == 0) && (fabsl(hp_matrix[j][k]) > max))
825 max=fabsl(hp_matrix[j][k]);
829 if ((column == -1) || (row == -1) || (fabsl(max) < LDBL_MIN))
830 ThrowGaussJordanException();
834 for (k=0; k < (ssize_t) rank; k++)
835 GaussJordanSwapLD(hp_matrix[row][k],hp_matrix[column][k]);
836 for (k=0; k < (ssize_t) number_vectors; k++)
837 GaussJordanSwap(vectors[k][row],vectors[k][column]);
841 if (fabsl(hp_matrix[column][column]) < LDBL_MIN)
842 ThrowGaussJordanException();
843 scale=1.0L/hp_matrix[column][column];
844 hp_matrix[column][column]=1.0;
845 for (j=0; j < (ssize_t) rank; j++)
846 hp_matrix[column][j]*=scale;
847 for (j=0; j < (ssize_t) number_vectors; j++)
848 vectors[j][column]*=(
double) scale;
849 for (j=0; j < (ssize_t) rank; j++)
852 scale=hp_matrix[j][column];
853 hp_matrix[j][column]=0.0;
854 for (k=0; k < (ssize_t) rank; k++)
855 hp_matrix[j][k]-=scale*hp_matrix[column][k];
856 for (k=0; k < (ssize_t) number_vectors; k++)
857 vectors[k][j]-=(
double)(scale*(
long double) vectors[k][column]);
860 for (j=(ssize_t) rank-1; j >= 0; j--)
861 if (columns[j] != rows[j])
862 for (i=0; i < (ssize_t) rank; i++)
863 GaussJordanSwapLD(hp_matrix[i][columns[j]],hp_matrix[i][rows[j]]);
867 for (i=0; i < (ssize_t) rank; i++)
868 for (j=0; j < (ssize_t) rank; j++)
869 matrix[i][j]=(
double)hp_matrix[i][j];
873 for (i=0; i < (ssize_t) rank; i++)
874 hp_matrix[i]=(
long double *) RelinquishMagickMemory(hp_matrix[i]);
875 hp_matrix=(
long double **) RelinquishMagickMemory(hp_matrix);
876 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
877 rows=(ssize_t *) RelinquishMagickMemory(rows);
878 columns=(ssize_t *) RelinquishMagickMemory(columns);
928 assert(matrix_info != (
const MatrixInfo *) NULL);
929 assert(matrix_info->signature == MagickCoreSignature);
931 assert(exception->signature == MagickCoreSignature);
932 if (matrix_info->stride <
sizeof(
double))
933 return((
Image *) NULL);
937 (void) GetMatrixElement(matrix_info,0,0,&min_value);
939 for (y=0; y < (ssize_t) matrix_info->rows; y++)
944 for (x=0; x < (ssize_t) matrix_info->columns; x++)
949 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
951 if (value < min_value)
954 if (value > max_value)
958 if ((min_value == 0.0) && (max_value == 0.0))
961 if (min_value == max_value)
963 scale_factor=(double) QuantumRange/min_value;
967 scale_factor=(double) QuantumRange/(max_value-min_value);
971 image=AcquireImage((
ImageInfo *) NULL,exception);
972 image->columns=matrix_info->columns;
973 image->rows=matrix_info->rows;
974 image->colorspace=GRAYColorspace;
976 image_view=AcquireAuthenticCacheView(image,exception);
977#if defined(MAGICKCORE_OPENMP_SUPPORT)
978 #pragma omp parallel for schedule(static) shared(status) \
979 magick_number_threads(image,image,image->rows,2)
981 for (y=0; y < (ssize_t) image->rows; y++)
992 if (status == MagickFalse)
994 q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
995 if (q == (Quantum *) NULL)
1000 for (x=0; x < (ssize_t) image->columns; x++)
1002 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
1004 value=scale_factor*(value-min_value);
1005 *q=ClampToQuantum(value);
1006 q+=(ptrdiff_t) GetPixelChannels(image);
1008 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1011 image_view=DestroyCacheView(image_view);
1012 if (status == MagickFalse)
1013 image=DestroyImage(image);
1039MagickExport MagickBooleanType NullMatrix(
MatrixInfo *matrix_info)
1051 assert(matrix_info != (
const MatrixInfo *) NULL);
1052 assert(matrix_info->signature == MagickCoreSignature);
1053 if (matrix_info->type != DiskCache)
1055 (void) memset(matrix_info->elements,0,(
size_t)
1056 matrix_info->length);
1060 (void) lseek(matrix_info->file,0,SEEK_SET);
1061 for (y=0; y < (ssize_t) matrix_info->rows; y++)
1063 for (x=0; x < (ssize_t) matrix_info->length; x++)
1065 count=write(matrix_info->file,&value,
sizeof(value));
1066 if (count != (ssize_t)
sizeof(value))
1069 if (x < (ssize_t) matrix_info->length)
1072 return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue);
1102MagickExport
double **RelinquishMagickMatrix(
double **matrix,
1103 const size_t number_rows)
1108 if (matrix == (
double **) NULL )
1110 for (i=0; i < (ssize_t) number_rows; i++)
1111 matrix[i]=(
double *) RelinquishMagickMemory(matrix[i]);
1112 matrix=(
double **) RelinquishMagickMemory(matrix);
1146MagickExport MagickBooleanType SetMatrixElement(
const MatrixInfo *matrix_info,
1147 const ssize_t x,
const ssize_t y,
const void *value)
1153 assert(matrix_info != (
const MatrixInfo *) NULL);
1154 assert(matrix_info->signature == MagickCoreSignature);
1155 i=y*(MagickOffsetType) matrix_info->columns+x;
1157 (((MagickSizeType) i*matrix_info->stride) >= matrix_info->length))
1158 return(MagickFalse);
1159 if (matrix_info->type != DiskCache)
1161 (void) memcpy((
unsigned char *) matrix_info->elements+i*
1162 (MagickOffsetType) matrix_info->stride,value,matrix_info->stride);
1165 count=WriteMatrixElements(matrix_info,i*(MagickOffsetType)
1166 matrix_info->stride,matrix_info->stride,(
unsigned char *) value);
1167 if (count != (MagickOffsetType) matrix_info->stride)
1168 return(MagickFalse);