42 #include "magick/studio.h"
43 #include "magick/blob.h"
44 #include "magick/blob-private.h"
45 #include "magick/exception.h"
46 #include "magick/exception-private.h"
47 #include "magick/image-private.h"
48 #include "magick/matrix.h"
49 #include "magick/memory_.h"
50 #include "magick/pixel-private.h"
51 #include "magick/resource_.h"
52 #include "magick/semaphore.h"
53 #include "magick/thread-private.h"
54 #include "magick/utility.h"
123 static void MatrixSignalHandler(
int status)
125 ThrowFatalException(CacheFatalError,
"UnableToExtendMatrixCache");
129 static inline MagickOffsetType WriteMatrixElements(
130 const MatrixInfo *magick_restrict matrix_info,
const MagickOffsetType offset,
131 const MagickSizeType length,
const unsigned char *magick_restrict buffer)
139 #if !defined(MAGICKCORE_HAVE_PWRITE)
140 LockSemaphoreInfo(matrix_info->semaphore);
141 if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
143 UnlockSemaphoreInfo(matrix_info->semaphore);
144 return((MagickOffsetType) -1);
148 for (i=0; i < (MagickOffsetType) length; i+=count)
150 #if !defined(MAGICKCORE_HAVE_PWRITE)
151 count=write(matrix_info->file,buffer+i,(
size_t) MagickMin(length-i,
152 (MagickSizeType) MAGICK_SSIZE_MAX));
154 count=pwrite(matrix_info->file,buffer+i,(
size_t) MagickMin(length-i,
155 (MagickSizeType) MAGICK_SSIZE_MAX),(off_t) (offset+i));
164 #if !defined(MAGICKCORE_HAVE_PWRITE)
165 UnlockSemaphoreInfo(matrix_info->semaphore);
170 static MagickBooleanType SetMatrixExtent(
171 MatrixInfo *magick_restrict matrix_info,MagickSizeType length)
178 if (length != (MagickSizeType) ((MagickOffsetType) length))
180 offset=(MagickOffsetType) lseek(matrix_info->file,0,SEEK_END);
183 if ((MagickSizeType) offset >= length)
185 extent=(MagickOffsetType) length-1;
186 count=WriteMatrixElements(matrix_info,extent,1,(
const unsigned char *)
"");
187 #if defined(MAGICKCORE_HAVE_POSIX_FALLOCATE)
188 if (matrix_info->synchronize != MagickFalse)
189 (void) posix_fallocate(matrix_info->file,offset+1,extent-offset);
192 (void) signal(SIGBUS,MatrixSignalHandler);
194 return(count != (MagickOffsetType) 1 ? MagickFalse : MagickTrue);
197 MagickExport
MatrixInfo *AcquireMatrixInfo(
const size_t columns,
198 const size_t rows,
const size_t stride,
ExceptionInfo *exception)
209 matrix_info=(
MatrixInfo *) AcquireMagickMemory(
sizeof(*matrix_info));
212 (void) memset(matrix_info,0,
sizeof(*matrix_info));
213 matrix_info->signature=MagickCoreSignature;
214 matrix_info->columns=columns;
215 matrix_info->rows=rows;
216 matrix_info->stride=stride;
217 matrix_info->semaphore=AllocateSemaphoreInfo();
218 synchronize=GetEnvironmentValue(
"MAGICK_SYNCHRONIZE");
219 if (synchronize != (
const char *) NULL)
221 matrix_info->synchronize=IsStringTrue(synchronize);
222 synchronize=DestroyString(synchronize);
224 matrix_info->length=(MagickSizeType) columns*rows*stride;
225 if (matrix_info->columns != (
size_t) (matrix_info->length/rows/stride))
227 (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
228 "CacheResourcesExhausted",
"`%s'",
"matrix cache");
229 return(DestroyMatrixInfo(matrix_info));
231 matrix_info->type=MemoryCache;
232 status=AcquireMagickResource(AreaResource,matrix_info->length);
233 if ((status != MagickFalse) &&
234 (matrix_info->length == (MagickSizeType) ((
size_t) matrix_info->length)))
236 status=AcquireMagickResource(MemoryResource,matrix_info->length);
237 if (status != MagickFalse)
239 matrix_info->mapped=MagickFalse;
240 matrix_info->elements=AcquireMagickMemory((
size_t)
241 matrix_info->length);
242 if (matrix_info->elements == NULL)
244 matrix_info->mapped=MagickTrue;
245 matrix_info->elements=MapBlob(-1,IOMode,0,(
size_t)
246 matrix_info->length);
248 if (matrix_info->elements == (
unsigned short *) NULL)
249 RelinquishMagickResource(MemoryResource,matrix_info->length);
252 matrix_info->file=(-1);
253 if (matrix_info->elements == (
unsigned short *) NULL)
255 status=AcquireMagickResource(DiskResource,matrix_info->length);
256 if (status == MagickFalse)
258 (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
259 "CacheResourcesExhausted",
"`%s'",
"matrix cache");
260 return(DestroyMatrixInfo(matrix_info));
262 matrix_info->type=DiskCache;
263 matrix_info->file=AcquireUniqueFileResource(matrix_info->path);
264 if (matrix_info->file == -1)
265 return(DestroyMatrixInfo(matrix_info));
266 status=AcquireMagickResource(MapResource,matrix_info->length);
267 if (status != MagickFalse)
269 status=SetMatrixExtent(matrix_info,matrix_info->length);
270 if (status != MagickFalse)
271 matrix_info->elements=(
void *) MapBlob(matrix_info->file,IOMode,0,
272 (
size_t) matrix_info->length);
273 if (matrix_info->elements != NULL)
274 matrix_info->type=MapCache;
276 RelinquishMagickResource(MapResource,matrix_info->length);
314 MagickExport
double **AcquireMagickMatrix(
const size_t number_rows,
324 matrix=(
double **) AcquireQuantumMemory(number_rows,
sizeof(*matrix));
325 if (matrix == (
double **) NULL)
326 return((
double **) NULL);
327 for (i=0; i < (ssize_t) number_rows; i++)
329 matrix[i]=(
double *) AcquireQuantumMemory(size,
sizeof(*matrix[i]));
330 if (matrix[i] == (
double *) NULL)
332 for (j=0; j < i; j++)
333 matrix[j]=(
double *) RelinquishMagickMemory(matrix[j]);
334 matrix=(
double **) RelinquishMagickMemory(matrix);
335 return((
double **) NULL);
337 for (j=0; j < (ssize_t) size; j++)
369 assert(matrix_info->signature == MagickCoreSignature);
370 LockSemaphoreInfo(matrix_info->semaphore);
371 switch (matrix_info->type)
375 if (matrix_info->mapped == MagickFalse)
376 matrix_info->elements=RelinquishMagickMemory(matrix_info->elements);
379 (void) UnmapBlob(matrix_info->elements,(
size_t) matrix_info->length);
380 matrix_info->elements=(
unsigned short *) NULL;
382 RelinquishMagickResource(MemoryResource,matrix_info->length);
387 (void) UnmapBlob(matrix_info->elements,(
size_t) matrix_info->length);
388 matrix_info->elements=NULL;
389 RelinquishMagickResource(MapResource,matrix_info->length);
393 if (matrix_info->file != -1)
394 (void) close(matrix_info->file);
395 (void) RelinquishUniqueFileResource(matrix_info->path);
396 RelinquishMagickResource(DiskResource,matrix_info->length);
402 UnlockSemaphoreInfo(matrix_info->semaphore);
403 DestroySemaphoreInfo(&matrix_info->semaphore);
404 return((
MatrixInfo *) RelinquishMagickMemory(matrix_info));
477 MagickExport MagickBooleanType GaussJordanElimination(
double **matrix,
478 double **vectors,
const size_t rank,
const size_t number_vectors)
480 #define GaussJordanSwap(x,y) \
506 columns=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*columns));
507 rows=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*rows));
508 pivots=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*pivots));
509 if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) ||
510 (pivots == (ssize_t *) NULL))
512 if (pivots != (ssize_t *) NULL)
513 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
514 if (columns != (ssize_t *) NULL)
515 columns=(ssize_t *) RelinquishMagickMemory(columns);
516 if (rows != (ssize_t *) NULL)
517 rows=(ssize_t *) RelinquishMagickMemory(rows);
520 (void) memset(columns,0,rank*
sizeof(*columns));
521 (void) memset(rows,0,rank*
sizeof(*rows));
522 (void) memset(pivots,0,rank*
sizeof(*pivots));
525 for (i=0; i < (ssize_t) rank; i++)
528 for (j=0; j < (ssize_t) rank; j++)
531 for (k=0; k < (ssize_t) rank; k++)
538 if (fabs(matrix[j][k]) >= max)
540 max=fabs(matrix[j][k]);
548 for (k=0; k < (ssize_t) rank; k++)
549 GaussJordanSwap(matrix[row][k],matrix[column][k]);
550 for (k=0; k < (ssize_t) number_vectors; k++)
551 GaussJordanSwap(vectors[k][row],vectors[k][column]);
555 if (matrix[column][column] == 0.0)
557 scale=PerceptibleReciprocal(matrix[column][column]);
558 matrix[column][column]=1.0;
559 for (j=0; j < (ssize_t) rank; j++)
560 matrix[column][j]*=scale;
561 for (j=0; j < (ssize_t) number_vectors; j++)
562 vectors[j][column]*=scale;
563 for (j=0; j < (ssize_t) rank; j++)
566 scale=matrix[j][column];
567 matrix[j][column]=0.0;
568 for (k=0; k < (ssize_t) rank; k++)
569 matrix[j][k]-=scale*matrix[column][k];
570 for (k=0; k < (ssize_t) number_vectors; k++)
571 vectors[k][j]-=scale*vectors[k][column];
574 for (j=(ssize_t) rank-1; j >= 0; j--)
575 if (columns[j] != rows[j])
576 for (i=0; i < (ssize_t) rank; i++)
577 GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]);
578 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
579 rows=(ssize_t *) RelinquishMagickMemory(rows);
580 columns=(ssize_t *) RelinquishMagickMemory(columns);
606 MagickExport
size_t GetMatrixColumns(
const MatrixInfo *matrix_info)
609 assert(matrix_info->signature == MagickCoreSignature);
610 return(matrix_info->columns);
643 static inline ssize_t EdgeX(
const ssize_t x,
const size_t columns)
647 if (x >= (ssize_t) columns)
648 return((ssize_t) (columns-1));
652 static inline ssize_t EdgeY(
const ssize_t y,
const size_t rows)
656 if (y >= (ssize_t) rows)
657 return((ssize_t) (rows-1));
661 static inline MagickOffsetType ReadMatrixElements(
662 const MatrixInfo *magick_restrict matrix_info,
const MagickOffsetType offset,
663 const MagickSizeType length,
unsigned char *magick_restrict buffer)
671 #if !defined(MAGICKCORE_HAVE_PREAD)
672 LockSemaphoreInfo(matrix_info->semaphore);
673 if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
675 UnlockSemaphoreInfo(matrix_info->semaphore);
676 return((MagickOffsetType) -1);
680 for (i=0; i < (MagickOffsetType) length; i+=count)
682 #if !defined(MAGICKCORE_HAVE_PREAD)
683 count=read(matrix_info->file,buffer+i,(
size_t) MagickMin(length-i,
684 (MagickSizeType) MAGICK_SSIZE_MAX));
686 count=pread(matrix_info->file,buffer+i,(
size_t) MagickMin(length-i,
687 (MagickSizeType) MAGICK_SSIZE_MAX),(off_t) (offset+i));
696 #if !defined(MAGICKCORE_HAVE_PREAD)
697 UnlockSemaphoreInfo(matrix_info->semaphore);
702 MagickExport MagickBooleanType GetMatrixElement(
const MatrixInfo *matrix_info,
703 const ssize_t x,
const ssize_t y,
void *value)
709 assert(matrix_info != (
const MatrixInfo *) NULL);
710 assert(matrix_info->signature == MagickCoreSignature);
711 i=(MagickOffsetType) EdgeY(y,matrix_info->rows)*matrix_info->columns+
712 EdgeX(x,matrix_info->columns);
713 if (matrix_info->type != DiskCache)
715 (void) memcpy(value,(
unsigned char *) matrix_info->elements+i*
716 matrix_info->stride,matrix_info->stride);
719 count=ReadMatrixElements(matrix_info,i*matrix_info->stride,
720 matrix_info->stride,(
unsigned char *) value);
721 if (count != (MagickOffsetType) matrix_info->stride)
748 MagickExport
size_t GetMatrixRows(
const MatrixInfo *matrix_info)
750 assert(matrix_info != (
const MatrixInfo *) NULL);
751 assert(matrix_info->signature == MagickCoreSignature);
752 return(matrix_info->rows);
826 MagickExport
void LeastSquaresAddTerms(
double **matrix,
double **vectors,
827 const double *terms,
const double *results,
const size_t rank,
828 const size_t number_vectors)
834 for (j=0; j < (ssize_t) rank; j++)
836 for (i=0; i < (ssize_t) rank; i++)
837 matrix[i][j]+=terms[i]*terms[j];
838 for (i=0; i < (ssize_t) number_vectors; i++)
839 vectors[i][j]+=results[i]*terms[j];
890 assert(matrix_info != (
const MatrixInfo *) NULL);
891 assert(matrix_info->signature == MagickCoreSignature);
893 assert(exception->signature == MagickCoreSignature);
894 if (matrix_info->stride <
sizeof(
double))
895 return((
Image *) NULL);
899 (void) GetMatrixElement(matrix_info,0,0,&value);
902 for (y=0; y < (ssize_t) matrix_info->rows; y++)
907 for (x=0; x < (ssize_t) matrix_info->columns; x++)
909 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
911 if (value < min_value)
914 if (value > max_value)
918 if ((min_value == 0.0) && (max_value == 0.0))
921 if (min_value == max_value)
923 scale_factor=(double) QuantumRange/min_value;
927 scale_factor=(double) QuantumRange/(max_value-min_value);
932 image->columns=matrix_info->columns;
933 image->rows=matrix_info->rows;
934 image->colorspace=GRAYColorspace;
936 image_view=AcquireAuthenticCacheView(image,exception);
937 #if defined(MAGICKCORE_OPENMP_SUPPORT)
938 #pragma omp parallel for schedule(static) shared(status) \
939 magick_number_threads(image,image,image->rows,1)
941 for (y=0; y < (ssize_t) image->rows; y++)
952 if (status == MagickFalse)
954 q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
960 for (x=0; x < (ssize_t) image->columns; x++)
962 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
964 value=scale_factor*(value-min_value);
965 q->red=ClampToQuantum(value);
970 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
973 image_view=DestroyCacheView(image_view);
974 if (status == MagickFalse)
975 image=DestroyImage(image);
1001 MagickExport MagickBooleanType NullMatrix(
MatrixInfo *matrix_info)
1013 assert(matrix_info != (
const MatrixInfo *) NULL);
1014 assert(matrix_info->signature == MagickCoreSignature);
1015 if (matrix_info->type != DiskCache)
1017 (void) memset(matrix_info->elements,0,(
size_t)
1018 matrix_info->length);
1022 (void) lseek(matrix_info->file,0,SEEK_SET);
1023 for (y=0; y < (ssize_t) matrix_info->rows; y++)
1025 for (x=0; x < (ssize_t) matrix_info->length; x++)
1027 count=write(matrix_info->file,&value,
sizeof(value));
1028 if (count != (ssize_t)
sizeof(value))
1031 if (x < (ssize_t) matrix_info->length)
1034 return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue);
1064 MagickExport
double **RelinquishMagickMatrix(
double **matrix,
1065 const size_t number_rows)
1070 if (matrix == (
double **) NULL )
1072 for (i=0; i < (ssize_t) number_rows; i++)
1073 matrix[i]=(
double *) RelinquishMagickMemory(matrix[i]);
1074 matrix=(
double **) RelinquishMagickMemory(matrix);
1108 MagickExport MagickBooleanType SetMatrixElement(
const MatrixInfo *matrix_info,
1109 const ssize_t x,
const ssize_t y,
const void *value)
1115 assert(matrix_info != (
const MatrixInfo *) NULL);
1116 assert(matrix_info->signature == MagickCoreSignature);
1117 i=(MagickOffsetType) y*matrix_info->columns+x;
1119 ((MagickSizeType) (i*matrix_info->stride) >= matrix_info->length))
1120 return(MagickFalse);
1121 if (matrix_info->type != DiskCache)
1123 (void) memcpy((
unsigned char *) matrix_info->elements+i*
1124 matrix_info->stride,value,matrix_info->stride);
1127 count=WriteMatrixElements(matrix_info,i*matrix_info->stride,
1128 matrix_info->stride,(
unsigned char *) value);
1129 if (count != (MagickOffsetType) matrix_info->stride)
1130 return(MagickFalse);