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/pixel-accessor.h"
53#include "MagickCore/resource_.h"
54#include "MagickCore/semaphore.h"
55#include "MagickCore/thread-private.h"
56#include "MagickCore/utility.h"
79 path[MagickPathExtent];
125static void MatrixSignalHandler(
int magick_unused(status))
127 magick_unreferenced(status);
128 ThrowFatalException(CacheFatalError,
"UnableToExtendMatrixCache");
132static inline MagickOffsetType WriteMatrixElements(
133 const MatrixInfo *magick_restrict matrix_info,
const MagickOffsetType offset,
134 const MagickSizeType length,
const unsigned char *magick_restrict buffer)
142#if !defined(MAGICKCORE_HAVE_PWRITE)
143 LockSemaphoreInfo(matrix_info->semaphore);
144 if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
146 UnlockSemaphoreInfo(matrix_info->semaphore);
147 return((MagickOffsetType) -1);
151 for (i=0; i < (MagickOffsetType) length; i+=count)
153#if !defined(MAGICKCORE_HAVE_PWRITE)
154 count=write(matrix_info->file,buffer+i,(
size_t) MagickMin(length-
155 (MagickSizeType) i,(MagickSizeType) MAGICK_SSIZE_MAX));
157 count=pwrite(matrix_info->file,buffer+i,(
size_t) MagickMin(length-
158 (MagickSizeType) i,(MagickSizeType) MAGICK_SSIZE_MAX),offset+i);
167#if !defined(MAGICKCORE_HAVE_PWRITE)
168 UnlockSemaphoreInfo(matrix_info->semaphore);
173static MagickBooleanType SetMatrixExtent(
174 MatrixInfo *magick_restrict matrix_info,MagickSizeType length)
181 if (length != (MagickSizeType) ((MagickOffsetType) length))
183 offset=(MagickOffsetType) lseek(matrix_info->file,0,SEEK_END);
186 if ((MagickSizeType) offset >= length)
188 extent=(MagickOffsetType) length-1;
189 count=WriteMatrixElements(matrix_info,extent,1,(
const unsigned char *)
"");
190#if defined(MAGICKCORE_HAVE_POSIX_FALLOCATE)
191 if (matrix_info->synchronize != MagickFalse)
192 (void) posix_fallocate(matrix_info->file,offset+1,extent-offset);
195 (void) signal(SIGBUS,MatrixSignalHandler);
197 return(count != (MagickOffsetType) 1 ? MagickFalse : MagickTrue);
200MagickExport
MatrixInfo *AcquireMatrixInfo(
const size_t columns,
201 const size_t rows,
const size_t stride,
ExceptionInfo *exception)
212 matrix_info=(
MatrixInfo *) AcquireMagickMemory(
sizeof(*matrix_info));
215 (void) memset(matrix_info,0,
sizeof(*matrix_info));
216 matrix_info->signature=MagickCoreSignature;
217 matrix_info->columns=columns;
218 matrix_info->rows=rows;
219 matrix_info->stride=stride;
220 matrix_info->semaphore=AcquireSemaphoreInfo();
221 synchronize=GetEnvironmentValue(
"MAGICK_SYNCHRONIZE");
222 if (synchronize != (
const char *) NULL)
224 matrix_info->synchronize=IsStringTrue(synchronize);
225 synchronize=DestroyString(synchronize);
227 matrix_info->length=(MagickSizeType) columns*rows*stride;
228 if (matrix_info->columns != (
size_t) (matrix_info->length/rows/stride))
230 (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
231 "CacheResourcesExhausted",
"`%s'",
"matrix cache");
232 return(DestroyMatrixInfo(matrix_info));
234 matrix_info->type=MemoryCache;
235 status=AcquireMagickResource(AreaResource,matrix_info->length);
236 if ((status != MagickFalse) &&
237 (matrix_info->length == (MagickSizeType) ((
size_t) matrix_info->length)))
239 status=AcquireMagickResource(MemoryResource,matrix_info->length);
240 if (status != MagickFalse)
242 matrix_info->mapped=MagickFalse;
243 matrix_info->elements=AcquireMagickMemory((
size_t)
244 matrix_info->length);
245 if (matrix_info->elements == NULL)
247 matrix_info->mapped=MagickTrue;
248 matrix_info->elements=MapBlob(-1,IOMode,0,(
size_t)
249 matrix_info->length);
251 if (matrix_info->elements == (
unsigned short *) NULL)
252 RelinquishMagickResource(MemoryResource,matrix_info->length);
255 matrix_info->file=(-1);
256 if (matrix_info->elements == (
unsigned short *) NULL)
258 status=AcquireMagickResource(DiskResource,matrix_info->length);
259 if (status == MagickFalse)
261 (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
262 "CacheResourcesExhausted",
"`%s'",
"matrix cache");
263 return(DestroyMatrixInfo(matrix_info));
265 matrix_info->type=DiskCache;
266 matrix_info->file=AcquireUniqueFileResource(matrix_info->path);
267 if (matrix_info->file == -1)
268 return(DestroyMatrixInfo(matrix_info));
269 status=AcquireMagickResource(MapResource,matrix_info->length);
270 if (status != MagickFalse)
272 status=SetMatrixExtent(matrix_info,matrix_info->length);
273 if (status != MagickFalse)
274 matrix_info->elements=(
void *) MapBlob(matrix_info->file,IOMode,0,
275 (
size_t) matrix_info->length);
276 if (matrix_info->elements != NULL)
277 matrix_info->type=MapCache;
279 RelinquishMagickResource(MapResource,matrix_info->length);
317MagickExport
double **AcquireMagickMatrix(
const size_t number_rows,
327 matrix=(
double **) AcquireQuantumMemory(number_rows,
sizeof(*matrix));
328 if (matrix == (
double **) NULL)
329 return((
double **) NULL);
330 for (i=0; i < (ssize_t) number_rows; i++)
332 matrix[i]=(
double *) AcquireQuantumMemory(size,
sizeof(*matrix[i]));
333 if (matrix[i] == (
double *) NULL)
335 for (j=0; j < i; j++)
336 matrix[j]=(
double *) RelinquishMagickMemory(matrix[j]);
337 matrix=(
double **) RelinquishMagickMemory(matrix);
338 return((
double **) NULL);
340 for (j=0; j < (ssize_t) size; j++)
372 assert(matrix_info->signature == MagickCoreSignature);
373 LockSemaphoreInfo(matrix_info->semaphore);
374 switch (matrix_info->type)
378 if (matrix_info->mapped == MagickFalse)
379 matrix_info->elements=RelinquishMagickMemory(matrix_info->elements);
382 (void) UnmapBlob(matrix_info->elements,(
size_t) matrix_info->length);
383 matrix_info->elements=(
unsigned short *) NULL;
385 RelinquishMagickResource(MemoryResource,matrix_info->length);
390 (void) UnmapBlob(matrix_info->elements,(
size_t) matrix_info->length);
391 matrix_info->elements=NULL;
392 RelinquishMagickResource(MapResource,matrix_info->length);
397 if (matrix_info->file != -1)
398 (void) close(matrix_info->file);
399 (void) RelinquishUniqueFileResource(matrix_info->path);
400 RelinquishMagickResource(DiskResource,matrix_info->length);
406 UnlockSemaphoreInfo(matrix_info->semaphore);
407 RelinquishSemaphoreInfo(&matrix_info->semaphore);
408 return((
MatrixInfo *) RelinquishMagickMemory(matrix_info));
481MagickPrivate MagickBooleanType GaussJordanElimination(
double **matrix,
482 double **vectors,
const size_t rank,
const size_t number_vectors)
484#define GaussJordanSwap(x,y) \
510 columns=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*columns));
511 rows=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*rows));
512 pivots=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*pivots));
513 if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) ||
514 (pivots == (ssize_t *) NULL))
516 if (pivots != (ssize_t *) NULL)
517 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
518 if (columns != (ssize_t *) NULL)
519 columns=(ssize_t *) RelinquishMagickMemory(columns);
520 if (rows != (ssize_t *) NULL)
521 rows=(ssize_t *) RelinquishMagickMemory(rows);
524 (void) memset(columns,0,rank*
sizeof(*columns));
525 (void) memset(rows,0,rank*
sizeof(*rows));
526 (void) memset(pivots,0,rank*
sizeof(*pivots));
529 for (i=0; i < (ssize_t) rank; i++)
532 for (j=0; j < (ssize_t) rank; j++)
535 for (k=0; k < (ssize_t) rank; k++)
542 if (fabs(matrix[j][k]) >= max)
544 max=fabs(matrix[j][k]);
552 for (k=0; k < (ssize_t) rank; k++)
553 GaussJordanSwap(matrix[row][k],matrix[column][k]);
554 for (k=0; k < (ssize_t) number_vectors; k++)
555 GaussJordanSwap(vectors[k][row],vectors[k][column]);
559 if (matrix[column][column] == 0.0)
561 scale=PerceptibleReciprocal(matrix[column][column]);
562 matrix[column][column]=1.0;
563 for (j=0; j < (ssize_t) rank; j++)
564 matrix[column][j]*=scale;
565 for (j=0; j < (ssize_t) number_vectors; j++)
566 vectors[j][column]*=scale;
567 for (j=0; j < (ssize_t) rank; j++)
570 scale=matrix[j][column];
571 matrix[j][column]=0.0;
572 for (k=0; k < (ssize_t) rank; k++)
573 matrix[j][k]-=scale*matrix[column][k];
574 for (k=0; k < (ssize_t) number_vectors; k++)
575 vectors[k][j]-=scale*vectors[k][column];
578 for (j=(ssize_t) rank-1; j >= 0; j--)
579 if (columns[j] != rows[j])
580 for (i=0; i < (ssize_t) rank; i++)
581 GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]);
582 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
583 rows=(ssize_t *) RelinquishMagickMemory(rows);
584 columns=(ssize_t *) RelinquishMagickMemory(columns);
610MagickExport
size_t GetMatrixColumns(
const MatrixInfo *matrix_info)
613 assert(matrix_info->signature == MagickCoreSignature);
614 return(matrix_info->columns);
647static inline ssize_t EdgeX(
const ssize_t x,
const size_t columns)
651 if (x >= (ssize_t) columns)
652 return((ssize_t) (columns-1));
656static inline ssize_t EdgeY(
const ssize_t y,
const size_t rows)
660 if (y >= (ssize_t) rows)
661 return((ssize_t) (rows-1));
665static inline MagickOffsetType ReadMatrixElements(
666 const MatrixInfo *magick_restrict matrix_info,
const MagickOffsetType offset,
667 const MagickSizeType length,
unsigned char *magick_restrict buffer)
675#if !defined(MAGICKCORE_HAVE_PREAD)
676 LockSemaphoreInfo(matrix_info->semaphore);
677 if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
679 UnlockSemaphoreInfo(matrix_info->semaphore);
680 return((MagickOffsetType) -1);
684 for (i=0; i < (MagickOffsetType) length; i+=count)
686#if !defined(MAGICKCORE_HAVE_PREAD)
687 count=read(matrix_info->file,buffer+i,(
size_t) MagickMin(length-i,
688 (MagickSizeType) MAGICK_SSIZE_MAX));
690 count=pread(matrix_info->file,buffer+i,(
size_t) MagickMin(length-
691 (MagickSizeType) i,(MagickSizeType) MAGICK_SSIZE_MAX),offset+i);
700#if !defined(MAGICKCORE_HAVE_PREAD)
701 UnlockSemaphoreInfo(matrix_info->semaphore);
706MagickExport MagickBooleanType GetMatrixElement(
const MatrixInfo *matrix_info,
707 const ssize_t x,
const ssize_t y,
void *value)
713 assert(matrix_info != (
const MatrixInfo *) NULL);
714 assert(matrix_info->signature == MagickCoreSignature);
715 i=EdgeY(y,matrix_info->rows)*(MagickOffsetType) matrix_info->columns+
716 EdgeX(x,matrix_info->columns);
717 if (matrix_info->type != DiskCache)
719 (void) memcpy(value,(
unsigned char *) matrix_info->elements+i*
720 (MagickOffsetType) matrix_info->stride,matrix_info->stride);
723 count=ReadMatrixElements(matrix_info,i*(MagickOffsetType) matrix_info->stride,
724 matrix_info->stride,(
unsigned char *) value);
725 if (count != (MagickOffsetType) matrix_info->stride)
752MagickExport
size_t GetMatrixRows(
const MatrixInfo *matrix_info)
754 assert(matrix_info != (
const MatrixInfo *) NULL);
755 assert(matrix_info->signature == MagickCoreSignature);
756 return(matrix_info->rows);
830MagickPrivate
void LeastSquaresAddTerms(
double **matrix,
double **vectors,
831 const double *terms,
const double *results,
const size_t rank,
832 const size_t number_vectors)
838 for (j=0; j < (ssize_t) rank; j++)
840 for (i=0; i < (ssize_t) rank; i++)
841 matrix[i][j]+=terms[i]*terms[j];
842 for (i=0; i < (ssize_t) number_vectors; i++)
843 vectors[i][j]+=results[i]*terms[j];
893 assert(matrix_info != (
const MatrixInfo *) NULL);
894 assert(matrix_info->signature == MagickCoreSignature);
896 assert(exception->signature == MagickCoreSignature);
897 if (matrix_info->stride <
sizeof(
double))
898 return((
Image *) NULL);
902 (void) GetMatrixElement(matrix_info,0,0,&min_value);
904 for (y=0; y < (ssize_t) matrix_info->rows; y++)
909 for (x=0; x < (ssize_t) matrix_info->columns; x++)
914 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
916 if (value < min_value)
919 if (value > max_value)
923 if ((min_value == 0.0) && (max_value == 0.0))
926 if (min_value == max_value)
928 scale_factor=(double) QuantumRange/min_value;
932 scale_factor=(double) QuantumRange/(max_value-min_value);
936 image=AcquireImage((
ImageInfo *) NULL,exception);
937 image->columns=matrix_info->columns;
938 image->rows=matrix_info->rows;
939 image->colorspace=GRAYColorspace;
941 image_view=AcquireAuthenticCacheView(image,exception);
942#if defined(MAGICKCORE_OPENMP_SUPPORT)
943 #pragma omp parallel for schedule(static) shared(status) \
944 magick_number_threads(image,image,image->rows,2)
946 for (y=0; y < (ssize_t) image->rows; y++)
957 if (status == MagickFalse)
959 q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
960 if (q == (Quantum *) NULL)
965 for (x=0; x < (ssize_t) image->columns; x++)
967 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
969 value=scale_factor*(value-min_value);
970 *q=ClampToQuantum(value);
971 q+=GetPixelChannels(image);
973 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
976 image_view=DestroyCacheView(image_view);
977 if (status == MagickFalse)
978 image=DestroyImage(image);
1004MagickExport MagickBooleanType NullMatrix(
MatrixInfo *matrix_info)
1016 assert(matrix_info != (
const MatrixInfo *) NULL);
1017 assert(matrix_info->signature == MagickCoreSignature);
1018 if (matrix_info->type != DiskCache)
1020 (void) memset(matrix_info->elements,0,(
size_t)
1021 matrix_info->length);
1025 (void) lseek(matrix_info->file,0,SEEK_SET);
1026 for (y=0; y < (ssize_t) matrix_info->rows; y++)
1028 for (x=0; x < (ssize_t) matrix_info->length; x++)
1030 count=write(matrix_info->file,&value,
sizeof(value));
1031 if (count != (ssize_t)
sizeof(value))
1034 if (x < (ssize_t) matrix_info->length)
1037 return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue);
1067MagickExport
double **RelinquishMagickMatrix(
double **matrix,
1068 const size_t number_rows)
1073 if (matrix == (
double **) NULL )
1075 for (i=0; i < (ssize_t) number_rows; i++)
1076 matrix[i]=(
double *) RelinquishMagickMemory(matrix[i]);
1077 matrix=(
double **) RelinquishMagickMemory(matrix);
1111MagickExport MagickBooleanType SetMatrixElement(
const MatrixInfo *matrix_info,
1112 const ssize_t x,
const ssize_t y,
const void *value)
1118 assert(matrix_info != (
const MatrixInfo *) NULL);
1119 assert(matrix_info->signature == MagickCoreSignature);
1120 i=y*(MagickOffsetType) matrix_info->columns+x;
1122 (((MagickSizeType) i*matrix_info->stride) >= matrix_info->length))
1123 return(MagickFalse);
1124 if (matrix_info->type != DiskCache)
1126 (void) memcpy((
unsigned char *) matrix_info->elements+i*
1127 (MagickOffsetType) matrix_info->stride,value,matrix_info->stride);
1130 count=WriteMatrixElements(matrix_info,i*(MagickOffsetType)
1131 matrix_info->stride,matrix_info->stride,(
unsigned char *) value);
1132 if (count != (MagickOffsetType) matrix_info->stride)
1133 return(MagickFalse);