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"
123static void MatrixSignalHandler(
int magick_unused(status))
125 magick_unreferenced(status);
126 ThrowFatalException(CacheFatalError,
"UnableToExtendMatrixCache");
130static inline MagickOffsetType WriteMatrixElements(
131 const MatrixInfo *magick_restrict matrix_info,
const MagickOffsetType offset,
132 const MagickSizeType length,
const unsigned char *magick_restrict buffer)
140#if !defined(MAGICKCORE_HAVE_PWRITE)
141 LockSemaphoreInfo(matrix_info->semaphore);
142 if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
144 UnlockSemaphoreInfo(matrix_info->semaphore);
145 return((MagickOffsetType) -1);
149 for (i=0; i < (MagickOffsetType) length; i+=count)
151#if !defined(MAGICKCORE_HAVE_PWRITE)
152 count=write(matrix_info->file,buffer+i,(
size_t) MagickMin(length-i,
153 (MagickSizeType) MagickMaxBufferExtent));
155 count=pwrite(matrix_info->file,buffer+i,(
size_t) MagickMin(length-i,
156 (MagickSizeType) MagickMaxBufferExtent),(off_t) (offset+i));
165#if !defined(MAGICKCORE_HAVE_PWRITE)
166 UnlockSemaphoreInfo(matrix_info->semaphore);
171static MagickBooleanType SetMatrixExtent(
172 MatrixInfo *magick_restrict matrix_info,MagickSizeType length)
179 if (length != (MagickSizeType) ((MagickOffsetType) length))
181 offset=(MagickOffsetType) lseek(matrix_info->file,0,SEEK_END);
184 if ((MagickSizeType) offset >= length)
186 extent=(MagickOffsetType) length-1;
187 count=WriteMatrixElements(matrix_info,extent,1,(
const unsigned char *)
"");
188#if defined(MAGICKCORE_HAVE_POSIX_FALLOCATE)
189 if (matrix_info->synchronize != MagickFalse)
190 (void) posix_fallocate(matrix_info->file,offset+1,extent-offset);
193 (void) signal(SIGBUS,MatrixSignalHandler);
195 return(count != (MagickOffsetType) 1 ? MagickFalse : MagickTrue);
198MagickExport MatrixInfo *AcquireMatrixInfo(
const size_t columns,
199 const size_t rows,
const size_t stride,ExceptionInfo *exception)
210 matrix_info=(MatrixInfo *) AcquireMagickMemory(
sizeof(*matrix_info));
211 if (matrix_info == (MatrixInfo *) NULL)
212 return((MatrixInfo *) NULL);
213 (void) memset(matrix_info,0,
sizeof(*matrix_info));
214 matrix_info->signature=MagickCoreSignature;
215 matrix_info->columns=columns;
216 matrix_info->rows=rows;
217 matrix_info->stride=stride;
218 matrix_info->semaphore=AllocateSemaphoreInfo();
219 synchronize=GetEnvironmentValue(
"MAGICK_SYNCHRONIZE");
220 if (synchronize != (
const char *) NULL)
222 matrix_info->synchronize=IsStringTrue(synchronize);
223 synchronize=DestroyString(synchronize);
225 matrix_info->length=(MagickSizeType) columns*rows*stride;
226 if (matrix_info->columns != (
size_t) (matrix_info->length/rows/stride))
228 (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
229 "CacheResourcesExhausted",
"`%s'",
"matrix cache");
230 return(DestroyMatrixInfo(matrix_info));
232 matrix_info->type=MemoryCache;
233 status=AcquireMagickResource(AreaResource,matrix_info->length);
234 if ((status != MagickFalse) &&
235 (matrix_info->length == (MagickSizeType) ((
size_t) matrix_info->length)))
237 status=AcquireMagickResource(MemoryResource,matrix_info->length);
238 if (status != MagickFalse)
240 matrix_info->mapped=MagickFalse;
241 matrix_info->elements=AcquireMagickMemory((
size_t)
242 matrix_info->length);
243 if (matrix_info->elements == NULL)
245 matrix_info->mapped=MagickTrue;
246 matrix_info->elements=MapBlob(-1,IOMode,0,(
size_t)
247 matrix_info->length);
249 if (matrix_info->elements == (
unsigned short *) NULL)
250 RelinquishMagickResource(MemoryResource,matrix_info->length);
253 matrix_info->file=(-1);
254 if (matrix_info->elements == (
unsigned short *) NULL)
256 status=AcquireMagickResource(DiskResource,matrix_info->length);
257 if (status == MagickFalse)
259 (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
260 "CacheResourcesExhausted",
"`%s'",
"matrix cache");
261 return(DestroyMatrixInfo(matrix_info));
263 matrix_info->type=DiskCache;
264 matrix_info->file=AcquireUniqueFileResource(matrix_info->path);
265 if (matrix_info->file == -1)
266 return(DestroyMatrixInfo(matrix_info));
267 status=AcquireMagickResource(MapResource,matrix_info->length);
268 if (status != MagickFalse)
270 status=SetMatrixExtent(matrix_info,matrix_info->length);
271 if (status != MagickFalse)
272 matrix_info->elements=(
void *) MapBlob(matrix_info->file,IOMode,0,
273 (
size_t) matrix_info->length);
274 if (matrix_info->elements != NULL)
275 matrix_info->type=MapCache;
277 RelinquishMagickResource(MapResource,matrix_info->length);
315MagickExport
double **AcquireMagickMatrix(
const size_t number_rows,
325 matrix=(
double **) AcquireQuantumMemory(number_rows,
sizeof(*matrix));
326 if (matrix == (
double **) NULL)
327 return((
double **) NULL);
328 for (i=0; i < (ssize_t) number_rows; i++)
330 matrix[i]=(
double *) AcquireQuantumMemory(size,
sizeof(*matrix[i]));
331 if (matrix[i] == (
double *) NULL)
333 for (j=0; j < i; j++)
334 matrix[j]=(
double *) RelinquishMagickMemory(matrix[j]);
335 matrix=(
double **) RelinquishMagickMemory(matrix);
336 return((
double **) NULL);
338 for (j=0; j < (ssize_t) size; j++)
367MagickExport MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info)
369 assert(matrix_info != (MatrixInfo *) NULL);
370 assert(matrix_info->signature == MagickCoreSignature);
371 LockSemaphoreInfo(matrix_info->semaphore);
372 switch (matrix_info->type)
376 if (matrix_info->mapped == MagickFalse)
377 matrix_info->elements=RelinquishMagickMemory(matrix_info->elements);
380 (void) UnmapBlob(matrix_info->elements,(
size_t) matrix_info->length);
381 matrix_info->elements=(
unsigned short *) NULL;
383 RelinquishMagickResource(MemoryResource,matrix_info->length);
388 (void) UnmapBlob(matrix_info->elements,(
size_t) matrix_info->length);
389 matrix_info->elements=NULL;
390 RelinquishMagickResource(MapResource,matrix_info->length);
395 if (matrix_info->file != -1)
396 (void) close(matrix_info->file);
397 (void) RelinquishUniqueFileResource(matrix_info->path);
398 RelinquishMagickResource(DiskResource,matrix_info->length);
404 UnlockSemaphoreInfo(matrix_info->semaphore);
405 DestroySemaphoreInfo(&matrix_info->semaphore);
406 return((MatrixInfo *) RelinquishMagickMemory(matrix_info));
479MagickExport MagickBooleanType GaussJordanElimination(
double **matrix,
480 double **vectors,
const size_t rank,
const size_t number_vectors)
482#define GaussJordanSwap(x,y) \
488#define GaussJordanSwapLD(x,y) \
490 long double temp = (x); \
494#define ThrowGaussJordanException() \
496 for (i=0; i < (ssize_t) rank; i++) \
497 hp_matrix[i]=(long double *) RelinquishMagickMemory(hp_matrix[i]); \
498 hp_matrix=(long double **) RelinquishMagickMemory(hp_matrix); \
499 if (pivots != (ssize_t *) NULL) \
500 pivots=(ssize_t *) RelinquishMagickMemory(pivots); \
501 if (rows != (ssize_t *) NULL) \
502 rows=(ssize_t *) RelinquishMagickMemory(rows); \
503 if (columns != (ssize_t *) NULL) \
504 columns=(ssize_t *) RelinquishMagickMemory(columns); \
505 return(MagickFalse); \
509 **hp_matrix = (
long double **) NULL,
514 *columns = (ssize_t *) NULL,
517 *pivots = (ssize_t *) NULL,
519 *rows = (ssize_t *) NULL;
524 hp_matrix=(
long double **) AcquireQuantumMemory(rank,
sizeof(*hp_matrix));
525 if (hp_matrix == (
long double **) NULL)
527 for (i=0; i < (ssize_t) rank; i++)
529 hp_matrix[i]=(
long double *) AcquireQuantumMemory(rank,
530 sizeof(*hp_matrix[i]));
531 if (hp_matrix[i] == (
long double *) NULL)
532 ThrowGaussJordanException();
533 for (j=0; j < (ssize_t) rank; j++)
534 hp_matrix[i][j]=(
long double)matrix[i][j];
536 columns=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*columns));
537 rows=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*rows));
538 pivots=(ssize_t *) AcquireQuantumMemory(rank,
sizeof(*pivots));
539 if ((columns == (ssize_t *) NULL) || (rows == (ssize_t *) NULL) ||
540 (pivots == (ssize_t *) NULL))
541 ThrowGaussJordanException();
542 (void) memset(columns,0,rank*
sizeof(*columns));
543 (void) memset(rows,0,rank*
sizeof(*rows));
544 (void) memset(pivots,0,rank*
sizeof(*pivots));
545 for (i=0; i < (ssize_t) rank; i++)
559 for (j=0; j < (ssize_t) rank; j++)
561 for (k=0; k < (ssize_t) rank; k++)
562 if ((pivots[k] == 0) && (fabsl(hp_matrix[j][k]) > max))
564 max=fabsl(hp_matrix[j][k]);
568 if ((column == -1) || (row == -1) || (fabsl(max) < LDBL_MIN))
569 ThrowGaussJordanException();
573 for (k=0; k < (ssize_t) rank; k++)
574 GaussJordanSwapLD(hp_matrix[row][k],hp_matrix[column][k]);
575 for (k=0; k < (ssize_t) number_vectors; k++)
576 GaussJordanSwap(vectors[k][row],vectors[k][column]);
580 if (fabsl(hp_matrix[column][column]) < LDBL_MIN)
581 ThrowGaussJordanException();
582 scale=1.0L/hp_matrix[column][column];
583 hp_matrix[column][column]=1.0;
584 for (j=0; j < (ssize_t) rank; j++)
585 hp_matrix[column][j]*=scale;
586 for (j=0; j < (ssize_t) number_vectors; j++)
587 vectors[j][column]*=(
double) scale;
588 for (j=0; j < (ssize_t) rank; j++)
591 scale=hp_matrix[j][column];
592 hp_matrix[j][column]=0.0;
593 for (k=0; k < (ssize_t) rank; k++)
594 hp_matrix[j][k]-=scale*hp_matrix[column][k];
595 for (k=0; k < (ssize_t) number_vectors; k++)
596 vectors[k][j]-=(
double)(scale*(
long double) vectors[k][column]);
599 for (j=(ssize_t) rank-1; j >= 0; j--)
600 if (columns[j] != rows[j])
601 for (i=0; i < (ssize_t) rank; i++)
602 GaussJordanSwapLD(hp_matrix[i][columns[j]],hp_matrix[i][rows[j]]);
606 for (i=0; i < (ssize_t) rank; i++)
607 for (j=0; j < (ssize_t) rank; j++)
608 matrix[i][j]=(
double)hp_matrix[i][j];
612 for (i=0; i < (ssize_t) rank; i++)
613 hp_matrix[i]=(
long double *) RelinquishMagickMemory(hp_matrix[i]);
614 hp_matrix=(
long double **) RelinquishMagickMemory(hp_matrix);
615 pivots=(ssize_t *) RelinquishMagickMemory(pivots);
616 rows=(ssize_t *) RelinquishMagickMemory(rows);
617 columns=(ssize_t *) RelinquishMagickMemory(columns);
643MagickExport
size_t GetMatrixColumns(
const MatrixInfo *matrix_info)
645 assert(matrix_info != (MatrixInfo *) NULL);
646 assert(matrix_info->signature == MagickCoreSignature);
647 return(matrix_info->columns);
680static inline ssize_t EdgeX(
const ssize_t x,
const size_t columns)
684 if (x >= (ssize_t) columns)
685 return((ssize_t) (columns-1));
689static inline ssize_t EdgeY(
const ssize_t y,
const size_t rows)
693 if (y >= (ssize_t) rows)
694 return((ssize_t) (rows-1));
698static inline MagickOffsetType ReadMatrixElements(
699 const MatrixInfo *magick_restrict matrix_info,
const MagickOffsetType offset,
700 const MagickSizeType length,
unsigned char *magick_restrict buffer)
708#if !defined(MAGICKCORE_HAVE_PREAD)
709 LockSemaphoreInfo(matrix_info->semaphore);
710 if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
712 UnlockSemaphoreInfo(matrix_info->semaphore);
713 return((MagickOffsetType) -1);
717 for (i=0; i < (MagickOffsetType) length; i+=count)
719#if !defined(MAGICKCORE_HAVE_PREAD)
720 count=read(matrix_info->file,buffer+i,(
size_t) MagickMin(length-i,
721 (MagickSizeType) MagickMaxBufferExtent));
723 count=pread(matrix_info->file,buffer+i,(
size_t) MagickMin(length-i,
724 (MagickSizeType) MagickMaxBufferExtent),(off_t) (offset+i));
733#if !defined(MAGICKCORE_HAVE_PREAD)
734 UnlockSemaphoreInfo(matrix_info->semaphore);
739MagickExport MagickBooleanType GetMatrixElement(
const MatrixInfo *matrix_info,
740 const ssize_t x,
const ssize_t y,
void *value)
746 assert(matrix_info != (
const MatrixInfo *) NULL);
747 assert(matrix_info->signature == MagickCoreSignature);
748 i=(MagickOffsetType) EdgeY(y,matrix_info->rows)*matrix_info->columns+
749 EdgeX(x,matrix_info->columns);
750 if (matrix_info->type != DiskCache)
752 (void) memcpy(value,(
unsigned char *) matrix_info->elements+i*
753 matrix_info->stride,matrix_info->stride);
756 count=ReadMatrixElements(matrix_info,i*matrix_info->stride,
757 matrix_info->stride,(
unsigned char *) value);
758 if (count != (MagickOffsetType) matrix_info->stride)
785MagickExport
size_t GetMatrixRows(
const MatrixInfo *matrix_info)
787 assert(matrix_info != (
const MatrixInfo *) NULL);
788 assert(matrix_info->signature == MagickCoreSignature);
789 return(matrix_info->rows);
863MagickExport
void LeastSquaresAddTerms(
double **matrix,
double **vectors,
864 const double *terms,
const double *results,
const size_t rank,
865 const size_t number_vectors)
871 for (j=0; j < (ssize_t) rank; j++)
873 for (i=0; i < (ssize_t) rank; i++)
874 matrix[i][j]+=terms[i]*terms[j];
875 for (i=0; i < (ssize_t) number_vectors; i++)
876 vectors[i][j]+=results[i]*terms[j];
906MagickExport Image *MatrixToImage(
const MatrixInfo *matrix_info,
907 ExceptionInfo *exception)
927 assert(matrix_info != (
const MatrixInfo *) NULL);
928 assert(matrix_info->signature == MagickCoreSignature);
929 assert(exception != (ExceptionInfo *) NULL);
930 assert(exception->signature == MagickCoreSignature);
931 if (matrix_info->stride <
sizeof(
double))
932 return((Image *) NULL);
936 (void) GetMatrixElement(matrix_info,0,0,&value);
939 for (y=0; y < (ssize_t) matrix_info->rows; y++)
944 for (x=0; x < (ssize_t) matrix_info->columns; x++)
946 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
948 if (value < min_value)
951 if (value > max_value)
955 if ((min_value == 0.0) && (max_value == 0.0))
958 if (min_value == max_value)
960 scale_factor=(double) QuantumRange/min_value;
964 scale_factor=(double) QuantumRange/(max_value-min_value);
968 image=AcquireImage((ImageInfo *) NULL);
969 image->columns=matrix_info->columns;
970 image->rows=matrix_info->rows;
971 image->colorspace=GRAYColorspace;
973 image_view=AcquireAuthenticCacheView(image,exception);
974#if defined(MAGICKCORE_OPENMP_SUPPORT)
975 #pragma omp parallel for schedule(static) shared(status) \
976 magick_number_threads(image,image,image->rows,2)
978 for (y=0; y < (ssize_t) image->rows; y++)
989 if (status == MagickFalse)
991 q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
992 if (q == (PixelPacket *) NULL)
997 for (x=0; x < (ssize_t) image->columns; x++)
999 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
1001 value=scale_factor*(value-min_value);
1002 q->red=ClampToQuantum(value);
1007 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1010 image_view=DestroyCacheView(image_view);
1011 if (status == MagickFalse)
1012 image=DestroyImage(image);
1038MagickExport MagickBooleanType NullMatrix(MatrixInfo *matrix_info)
1050 assert(matrix_info != (
const MatrixInfo *) NULL);
1051 assert(matrix_info->signature == MagickCoreSignature);
1052 if (matrix_info->type != DiskCache)
1054 (void) memset(matrix_info->elements,0,(
size_t)
1055 matrix_info->length);
1059 (void) lseek(matrix_info->file,0,SEEK_SET);
1060 for (y=0; y < (ssize_t) matrix_info->rows; y++)
1062 for (x=0; x < (ssize_t) matrix_info->length; x++)
1064 count=write(matrix_info->file,&value,
sizeof(value));
1065 if (count != (ssize_t)
sizeof(value))
1068 if (x < (ssize_t) matrix_info->length)
1071 return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue);
1101MagickExport
double **RelinquishMagickMatrix(
double **matrix,
1102 const size_t number_rows)
1107 if (matrix == (
double **) NULL )
1109 for (i=0; i < (ssize_t) number_rows; i++)
1110 matrix[i]=(
double *) RelinquishMagickMemory(matrix[i]);
1111 matrix=(
double **) RelinquishMagickMemory(matrix);
1145MagickExport MagickBooleanType SetMatrixElement(
const MatrixInfo *matrix_info,
1146 const ssize_t x,
const ssize_t y,
const void *value)
1152 assert(matrix_info != (
const MatrixInfo *) NULL);
1153 assert(matrix_info->signature == MagickCoreSignature);
1154 i=(MagickOffsetType) y*matrix_info->columns+x;
1156 ((MagickSizeType) (i*matrix_info->stride) >= matrix_info->length))
1157 return(MagickFalse);
1158 if (matrix_info->type != DiskCache)
1160 (void) memcpy((
unsigned char *) matrix_info->elements+i*
1161 matrix_info->stride,value,matrix_info->stride);
1164 count=WriteMatrixElements(matrix_info,i*matrix_info->stride,
1165 matrix_info->stride,(
unsigned char *) value);
1166 if (count != (MagickOffsetType) matrix_info->stride)
1167 return(MagickFalse);