1/* 2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3% % 4% % 5% % 6% M M AAA TTTTT RRRR IIIII X X % 7% MM MM A A T R R I X X % 8% M M M AAAAA T RRRR I X % 9% M M A A T R R I X X % 10% M M A A T R R IIIII X X % 11% % 12% % 13% MagickCore Matrix Methods % 14% % 15% Software Design % 16% Cristy % 17% August 2007 % 18% % 19% % 20% Copyright 1999-2016 ImageMagick Studio LLC, a non-profit organization % 21% dedicated to making software imaging solutions freely available. % 22% % 23% You may not use this file except in compliance with the License. You may % 24% obtain a copy of the License at % 25% % 26% http://www.imagemagick.org/script/license.php % 27% % 28% Unless required by applicable law or agreed to in writing, software % 29% distributed under the License is distributed on an "AS IS" BASIS, % 30% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. % 31% See the License for the specific language governing permissions and % 32% limitations under the License. % 33% % 34%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 35% 36% 37*/ 38 39/* 40 Include declarations. 41*/ 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/memory_.h" 51#include "MagickCore/pixel-accessor.h" 52#include "MagickCore/pixel-private.h" 53#include "MagickCore/resource_.h" 54#include "MagickCore/semaphore.h" 55#include "MagickCore/thread-private.h" 56#include "MagickCore/utility.h" 57 58/* 59 Typedef declaration. 60*/ 61struct _MatrixInfo 62{ 63 CacheType 64 type; 65 66 size_t 67 columns, 68 rows, 69 stride; 70 71 MagickSizeType 72 length; 73 74 MagickBooleanType 75 mapped, 76 synchronize; 77 78 char 79 path[MagickPathExtent]; 80 81 int 82 file; 83 84 void 85 *elements; 86 87 SemaphoreInfo 88 *semaphore; 89 90 size_t 91 signature; 92}; 93 94/* 95%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 96% % 97% % 98% % 99% A c q u i r e M a t r i x I n f o % 100% % 101% % 102% % 103%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 104% 105% AcquireMatrixInfo() allocates the ImageInfo structure. 106% 107% The format of the AcquireMatrixInfo method is: 108% 109% MatrixInfo *AcquireMatrixInfo(const size_t columns,const size_t rows, 110% const size_t stride,ExceptionInfo *exception) 111% 112% A description of each parameter follows: 113% 114% o columns: the matrix columns. 115% 116% o rows: the matrix rows. 117% 118% o stride: the matrix stride. 119% 120% o exception: return any errors or warnings in this structure. 121% 122*/ 123 124#if defined(SIGBUS) 125static void MatrixSignalHandler(int status) 126{ 127 ThrowFatalException(CacheFatalError,"UnableToExtendMatrixCache"); 128} 129#endif 130 131static inline MagickOffsetType WriteMatrixElements( 132 const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset, 133 const MagickSizeType length,const unsigned char *magick_restrict buffer) 134{ 135 register MagickOffsetType 136 i; 137 138 ssize_t 139 count; 140 141#if !defined(MAGICKCORE_HAVE_PWRITE) 142 LockSemaphoreInfo(matrix_info->semaphore); 143 if (lseek(matrix_info->file,offset,SEEK_SET) < 0) 144 { 145 UnlockSemaphoreInfo(matrix_info->semaphore); 146 return((MagickOffsetType) -1); 147 } 148#endif 149 count=0; 150 for (i=0; i < (MagickOffsetType) length; i+=count) 151 { 152#if !defined(MAGICKCORE_HAVE_PWRITE) 153 count=write(matrix_info->file,buffer+i,(size_t) MagickMin(length-i, 154 (MagickSizeType) SSIZE_MAX)); 155#else 156 count=pwrite(matrix_info->file,buffer+i,(size_t) MagickMin(length-i, 157 (MagickSizeType) SSIZE_MAX),(off_t) (offset+i)); 158#endif 159 if (count <= 0) 160 { 161 count=0; 162 if (errno != EINTR) 163 break; 164 } 165 } 166#if !defined(MAGICKCORE_HAVE_PWRITE) 167 UnlockSemaphoreInfo(matrix_info->semaphore); 168#endif 169 return(i); 170} 171 172static MagickBooleanType SetMatrixExtent( 173 MatrixInfo *magick_restrict matrix_info, 174 MagickSizeType length) 175{ 176 MagickOffsetType 177 count, 178 extent, 179 offset; 180 181 if (length != (MagickSizeType) ((MagickOffsetType) length)) 182 return(MagickFalse); 183 offset=(MagickOffsetType) lseek(matrix_info->file,0,SEEK_END); 184 if (offset < 0) 185 return(MagickFalse); 186 if ((MagickSizeType) offset >= length) 187 return(MagickTrue); 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); 193#endif 194#if defined(SIGBUS) 195 (void) signal(SIGBUS,MatrixSignalHandler); 196#endif 197 return(count != (MagickOffsetType) 1 ? MagickFalse : MagickTrue); 198} 199 200MagickExport MatrixInfo *AcquireMatrixInfo(const size_t columns, 201 const size_t rows,const size_t stride,ExceptionInfo *exception) 202{ 203 char 204 *synchronize; 205 206 MagickBooleanType 207 status; 208 209 MatrixInfo 210 *matrix_info; 211 212 matrix_info=(MatrixInfo *) AcquireMagickMemory(sizeof(*matrix_info)); 213 if (matrix_info == (MatrixInfo *) NULL) 214 return((MatrixInfo *) NULL); 215 (void) ResetMagickMemory(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) 223 { 224 matrix_info->synchronize=IsStringTrue(synchronize); 225 synchronize=DestroyString(synchronize); 226 } 227 matrix_info->length=(MagickSizeType) columns*rows*stride; 228 if (matrix_info->columns != (size_t) (matrix_info->length/rows/stride)) 229 { 230 (void) ThrowMagickException(exception,GetMagickModule(),CacheError, 231 "CacheResourcesExhausted","`%s'","matrix cache"); 232 return(DestroyMatrixInfo(matrix_info)); 233 } 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))) 238 { 239 status=AcquireMagickResource(MemoryResource,matrix_info->length); 240 if (status != MagickFalse) 241 { 242 matrix_info->mapped=MagickFalse; 243 matrix_info->elements=AcquireMagickMemory((size_t) 244 matrix_info->length); 245 if (matrix_info->elements == NULL) 246 { 247 matrix_info->mapped=MagickTrue; 248 matrix_info->elements=MapBlob(-1,IOMode,0,(size_t) 249 matrix_info->length); 250 } 251 if (matrix_info->elements == (unsigned short *) NULL) 252 RelinquishMagickResource(MemoryResource,matrix_info->length); 253 } 254 } 255 matrix_info->file=(-1); 256 if (matrix_info->elements == (unsigned short *) NULL) 257 { 258 status=AcquireMagickResource(DiskResource,matrix_info->length); 259 if (status == MagickFalse) 260 { 261 (void) ThrowMagickException(exception,GetMagickModule(),CacheError, 262 "CacheResourcesExhausted","`%s'","matrix cache"); 263 return(DestroyMatrixInfo(matrix_info)); 264 } 265 matrix_info->type=DiskCache; 266 (void) AcquireMagickResource(MemoryResource,matrix_info->length); 267 matrix_info->file=AcquireUniqueFileResource(matrix_info->path); 268 if (matrix_info->file == -1) 269 return(DestroyMatrixInfo(matrix_info)); 270 status=AcquireMagickResource(MapResource,matrix_info->length); 271 if (status != MagickFalse) 272 { 273 status=SetMatrixExtent(matrix_info,matrix_info->length); 274 if (status != MagickFalse) 275 { 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; 280 else 281 RelinquishMagickResource(MapResource,matrix_info->length); 282 } 283 } 284 } 285 return(matrix_info); 286} 287 288/* 289%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 290% % 291% % 292% % 293% A c q u i r e M a g i c k M a t r i x % 294% % 295% % 296% % 297%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 298% 299% AcquireMagickMatrix() allocates and returns a matrix in the form of an 300% array of pointers to an array of doubles, with all values pre-set to zero. 301% 302% This used to generate the two dimensional matrix, and vectors required 303% for the GaussJordanElimination() method below, solving some system of 304% simultanious equations. 305% 306% The format of the AcquireMagickMatrix method is: 307% 308% double **AcquireMagickMatrix(const size_t number_rows, 309% const size_t size) 310% 311% A description of each parameter follows: 312% 313% o number_rows: the number pointers for the array of pointers 314% (first dimension). 315% 316% o size: the size of the array of doubles each pointer points to 317% (second dimension). 318% 319*/ 320MagickExport double **AcquireMagickMatrix(const size_t number_rows, 321 const size_t size) 322{ 323 double 324 **matrix; 325 326 register ssize_t 327 i, 328 j; 329 330 matrix=(double **) AcquireQuantumMemory(number_rows,sizeof(*matrix)); 331 if (matrix == (double **) NULL) 332 return((double **) NULL); 333 for (i=0; i < (ssize_t) number_rows; i++) 334 { 335 matrix[i]=(double *) AcquireQuantumMemory(size,sizeof(*matrix[i])); 336 if (matrix[i] == (double *) NULL) 337 { 338 for (j=0; j < i; j++) 339 matrix[j]=(double *) RelinquishMagickMemory(matrix[j]); 340 matrix=(double **) RelinquishMagickMemory(matrix); 341 return((double **) NULL); 342 } 343 for (j=0; j < (ssize_t) size; j++) 344 matrix[i][j]=0.0; 345 } 346 return(matrix); 347} 348 349/* 350%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 351% % 352% % 353% % 354% D e s t r o y M a t r i x I n f o % 355% % 356% % 357% % 358%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 359% 360% DestroyMatrixInfo() dereferences a matrix, deallocating memory associated 361% with the matrix. 362% 363% The format of the DestroyImage method is: 364% 365% MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info) 366% 367% A description of each parameter follows: 368% 369% o matrix_info: the matrix. 370% 371*/ 372MagickExport MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info) 373{ 374 assert(matrix_info != (MatrixInfo *) NULL); 375 assert(matrix_info->signature == MagickCoreSignature); 376 LockSemaphoreInfo(matrix_info->semaphore); 377 switch (matrix_info->type) 378 { 379 case MemoryCache: 380 { 381 if (matrix_info->mapped == MagickFalse) 382 matrix_info->elements=RelinquishMagickMemory(matrix_info->elements); 383 else 384 { 385 (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length); 386 matrix_info->elements=(unsigned short *) NULL; 387 } 388 RelinquishMagickResource(MemoryResource,matrix_info->length); 389 break; 390 } 391 case MapCache: 392 { 393 (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length); 394 matrix_info->elements=NULL; 395 RelinquishMagickResource(MapResource,matrix_info->length); 396 } 397 case DiskCache: 398 { 399 if (matrix_info->file != -1) 400 (void) close(matrix_info->file); 401 (void) RelinquishUniqueFileResource(matrix_info->path); 402 RelinquishMagickResource(DiskResource,matrix_info->length); 403 break; 404 } 405 default: 406 break; 407 } 408 UnlockSemaphoreInfo(matrix_info->semaphore); 409 RelinquishSemaphoreInfo(&matrix_info->semaphore); 410 return((MatrixInfo *) RelinquishMagickMemory(matrix_info)); 411} 412 413/* 414%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 415% % 416% % 417% % 418+ G a u s s J o r d a n E l i m i n a t i o n % 419% % 420% % 421% % 422%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 423% 424% GaussJordanElimination() returns a matrix in reduced row echelon form, 425% while simultaneously reducing and thus solving the augumented results 426% matrix. 427% 428% See also http://en.wikipedia.org/wiki/Gauss-Jordan_elimination 429% 430% The format of the GaussJordanElimination method is: 431% 432% MagickBooleanType GaussJordanElimination(double **matrix, 433% double **vectors,const size_t rank,const size_t number_vectors) 434% 435% A description of each parameter follows: 436% 437% o matrix: the matrix to be reduced, as an 'array of row pointers'. 438% 439% o vectors: the additional matrix argumenting the matrix for row reduction. 440% Producing an 'array of column vectors'. 441% 442% o rank: The size of the matrix (both rows and columns). 443% Also represents the number terms that need to be solved. 444% 445% o number_vectors: Number of vectors columns, argumenting the above matrix. 446% Usally 1, but can be more for more complex equation solving. 447% 448% Note that the 'matrix' is given as a 'array of row pointers' of rank size. 449% That is values can be assigned as matrix[row][column] where 'row' is 450% typically the equation, and 'column' is the term of the equation. 451% That is the matrix is in the form of a 'row first array'. 452% 453% However 'vectors' is a 'array of column pointers' which can have any number 454% of columns, with each column array the same 'rank' size as 'matrix'. 455% 456% This allows for simpler handling of the results, especially is only one 457% column 'vector' is all that is required to produce the desired solution. 458% 459% For example, the 'vectors' can consist of a pointer to a simple array of 460% doubles. when only one set of simultanious equations is to be solved from 461% the given set of coefficient weighted terms. 462% 463% double **matrix = AcquireMagickMatrix(8UL,8UL); 464% double coefficents[8]; 465% ... 466% GaussJordanElimination(matrix, &coefficents, 8UL, 1UL); 467% 468% However by specifing more 'columns' (as an 'array of vector columns', 469% you can use this function to solve a set of 'separable' equations. 470% 471% For example a distortion function where u = U(x,y) v = V(x,y) 472% And the functions U() and V() have separate coefficents, but are being 473% generated from a common x,y->u,v data set. 474% 475% Another example is generation of a color gradient from a set of colors at 476% specific coordients, such as a list x,y -> r,g,b,a. 477% 478% You can also use the 'vectors' to generate an inverse of the given 'matrix' 479% though as a 'column first array' rather than a 'row first array'. For 480% details see http://en.wikipedia.org/wiki/Gauss-Jordan_elimination 481% 482*/ 483MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix, 484 double **vectors,const size_t rank,const size_t number_vectors) 485{ 486#define GaussJordanSwap(x,y) \ 487{ \ 488 if ((x) != (y)) \ 489 { \ 490 (x)+=(y); \ 491 (y)=(x)-(y); \ 492 (x)=(x)-(y); \ 493 } \ 494} 495 496 double 497 max, 498 scale; 499 500 register ssize_t 501 i, 502 j, 503 k; 504 505 ssize_t 506 column, 507 *columns, 508 *pivots, 509 row, 510 *rows; 511 512 columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns)); 513 rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows)); 514 pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots)); 515 if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) || 516 (pivots == (ssize_t *) NULL)) 517 { 518 if (pivots != (ssize_t *) NULL) 519 pivots=(ssize_t *) RelinquishMagickMemory(pivots); 520 if (columns != (ssize_t *) NULL) 521 columns=(ssize_t *) RelinquishMagickMemory(columns); 522 if (rows != (ssize_t *) NULL) 523 rows=(ssize_t *) RelinquishMagickMemory(rows); 524 return(MagickFalse); 525 } 526 (void) ResetMagickMemory(columns,0,rank*sizeof(*columns)); 527 (void) ResetMagickMemory(rows,0,rank*sizeof(*rows)); 528 (void) ResetMagickMemory(pivots,0,rank*sizeof(*pivots)); 529 column=0; 530 row=0; 531 for (i=0; i < (ssize_t) rank; i++) 532 { 533 max=0.0; 534 for (j=0; j < (ssize_t) rank; j++) 535 if (pivots[j] != 1) 536 { 537 for (k=0; k < (ssize_t) rank; k++) 538 if (pivots[k] != 0) 539 { 540 if (pivots[k] > 1) 541 return(MagickFalse); 542 } 543 else 544 if (fabs(matrix[j][k]) >= max) 545 { 546 max=fabs(matrix[j][k]); 547 row=j; 548 column=k; 549 } 550 } 551 pivots[column]++; 552 if (row != column) 553 { 554 for (k=0; k < (ssize_t) rank; k++) 555 GaussJordanSwap(matrix[row][k],matrix[column][k]); 556 for (k=0; k < (ssize_t) number_vectors; k++) 557 GaussJordanSwap(vectors[k][row],vectors[k][column]); 558 } 559 rows[i]=row; 560 columns[i]=column; 561 if (matrix[column][column] == 0.0) 562 return(MagickFalse); /* sigularity */ 563 scale=PerceptibleReciprocal(matrix[column][column]); 564 matrix[column][column]=1.0; 565 for (j=0; j < (ssize_t) rank; j++) 566 matrix[column][j]*=scale; 567 for (j=0; j < (ssize_t) number_vectors; j++) 568 vectors[j][column]*=scale; 569 for (j=0; j < (ssize_t) rank; j++) 570 if (j != column) 571 { 572 scale=matrix[j][column]; 573 matrix[j][column]=0.0; 574 for (k=0; k < (ssize_t) rank; k++) 575 matrix[j][k]-=scale*matrix[column][k]; 576 for (k=0; k < (ssize_t) number_vectors; k++) 577 vectors[k][j]-=scale*vectors[k][column]; 578 } 579 } 580 for (j=(ssize_t) rank-1; j >= 0; j--) 581 if (columns[j] != rows[j]) 582 for (i=0; i < (ssize_t) rank; i++) 583 GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]); 584 pivots=(ssize_t *) RelinquishMagickMemory(pivots); 585 rows=(ssize_t *) RelinquishMagickMemory(rows); 586 columns=(ssize_t *) RelinquishMagickMemory(columns); 587 return(MagickTrue); 588} 589 590/* 591%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 592% % 593% % 594% % 595% G e t M a t r i x C o l u m n s % 596% % 597% % 598% % 599%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 600% 601% GetMatrixColumns() returns the number of columns in the matrix. 602% 603% The format of the GetMatrixColumns method is: 604% 605% size_t GetMatrixColumns(const MatrixInfo *matrix_info) 606% 607% A description of each parameter follows: 608% 609% o matrix_info: the matrix. 610% 611*/ 612MagickExport size_t GetMatrixColumns(const MatrixInfo *matrix_info) 613{ 614 assert(matrix_info != (MatrixInfo *) NULL); 615 assert(matrix_info->signature == MagickCoreSignature); 616 return(matrix_info->columns); 617} 618 619/* 620%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 621% % 622% % 623% % 624% G e t M a t r i x E l e m e n t % 625% % 626% % 627% % 628%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 629% 630% GetMatrixElement() returns the specifed element in the matrix. 631% 632% The format of the GetMatrixElement method is: 633% 634% MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info, 635% const ssize_t x,const ssize_t y,void *value) 636% 637% A description of each parameter follows: 638% 639% o matrix_info: the matrix columns. 640% 641% o x: the matrix x-offset. 642% 643% o y: the matrix y-offset. 644% 645% o value: return the matrix element in this buffer. 646% 647*/ 648 649static inline ssize_t EdgeX(const ssize_t x,const size_t columns) 650{ 651 if (x < 0L) 652 return(0L); 653 if (x >= (ssize_t) columns) 654 return((ssize_t) (columns-1)); 655 return(x); 656} 657 658static inline ssize_t EdgeY(const ssize_t y,const size_t rows) 659{ 660 if (y < 0L) 661 return(0L); 662 if (y >= (ssize_t) rows) 663 return((ssize_t) (rows-1)); 664 return(y); 665} 666 667static inline MagickOffsetType ReadMatrixElements( 668 const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset, 669 const MagickSizeType length,unsigned char *magick_restrict buffer) 670{ 671 register MagickOffsetType 672 i; 673 674 ssize_t 675 count; 676 677#if !defined(MAGICKCORE_HAVE_PREAD) 678 LockSemaphoreInfo(matrix_info->semaphore); 679 if (lseek(matrix_info->file,offset,SEEK_SET) < 0) 680 { 681 UnlockSemaphoreInfo(matrix_info->semaphore); 682 return((MagickOffsetType) -1); 683 } 684#endif 685 count=0; 686 for (i=0; i < (MagickOffsetType) length; i+=count) 687 { 688#if !defined(MAGICKCORE_HAVE_PREAD) 689 count=read(matrix_info->file,buffer+i,(size_t) MagickMin(length-i, 690 (MagickSizeType) SSIZE_MAX)); 691#else 692 count=pread(matrix_info->file,buffer+i,(size_t) MagickMin(length-i, 693 (MagickSizeType) SSIZE_MAX),(off_t) (offset+i)); 694#endif 695 if (count <= 0) 696 { 697 count=0; 698 if (errno != EINTR) 699 break; 700 } 701 } 702#if !defined(MAGICKCORE_HAVE_PREAD) 703 UnlockSemaphoreInfo(matrix_info->semaphore); 704#endif 705 return(i); 706} 707 708MagickExport MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info, 709 const ssize_t x,const ssize_t y,void *value) 710{ 711 MagickOffsetType 712 count, 713 i; 714 715 assert(matrix_info != (const MatrixInfo *) NULL); 716 assert(matrix_info->signature == MagickCoreSignature); 717 i=(MagickOffsetType) EdgeY(y,matrix_info->rows)*matrix_info->columns+ 718 EdgeX(x,matrix_info->columns); 719 if (matrix_info->type != DiskCache) 720 { 721 (void) memcpy(value,(unsigned char *) matrix_info->elements+i* 722 matrix_info->stride,matrix_info->stride); 723 return(MagickTrue); 724 } 725 count=ReadMatrixElements(matrix_info,i*matrix_info->stride, 726 matrix_info->stride,(unsigned char *) value); 727 if (count != (MagickOffsetType) matrix_info->stride) 728 return(MagickFalse); 729 return(MagickTrue); 730} 731 732/* 733%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 734% % 735% % 736% % 737% G e t M a t r i x R o w s % 738% % 739% % 740% % 741%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 742% 743% GetMatrixRows() returns the number of rows in the matrix. 744% 745% The format of the GetMatrixRows method is: 746% 747% size_t GetMatrixRows(const MatrixInfo *matrix_info) 748% 749% A description of each parameter follows: 750% 751% o matrix_info: the matrix. 752% 753*/ 754MagickExport size_t GetMatrixRows(const MatrixInfo *matrix_info) 755{ 756 assert(matrix_info != (const MatrixInfo *) NULL); 757 assert(matrix_info->signature == MagickCoreSignature); 758 return(matrix_info->rows); 759} 760 761/* 762%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 763% % 764% % 765% % 766+ L e a s t S q u a r e s A d d T e r m s % 767% % 768% % 769% % 770%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 771% 772% LeastSquaresAddTerms() adds one set of terms and associate results to the 773% given matrix and vectors for solving using least-squares function fitting. 774% 775% The format of the AcquireMagickMatrix method is: 776% 777% void LeastSquaresAddTerms(double **matrix,double **vectors, 778% const double *terms,const double *results,const size_t rank, 779% const size_t number_vectors); 780% 781% A description of each parameter follows: 782% 783% o matrix: the square matrix to add given terms/results to. 784% 785% o vectors: the result vectors to add terms/results to. 786% 787% o terms: the pre-calculated terms (without the unknown coefficent 788% weights) that forms the equation being added. 789% 790% o results: the result(s) that should be generated from the given terms 791% weighted by the yet-to-be-solved coefficents. 792% 793% o rank: the rank or size of the dimensions of the square matrix. 794% Also the length of vectors, and number of terms being added. 795% 796% o number_vectors: Number of result vectors, and number or results being 797% added. Also represents the number of separable systems of equations 798% that is being solved. 799% 800% Example of use... 801% 802% 2 dimensional Affine Equations (which are separable) 803% c0*x + c2*y + c4*1 => u 804% c1*x + c3*y + c5*1 => v 805% 806% double **matrix = AcquireMagickMatrix(3UL,3UL); 807% double **vectors = AcquireMagickMatrix(2UL,3UL); 808% double terms[3], results[2]; 809% ... 810% for each given x,y -> u,v 811% terms[0] = x; 812% terms[1] = y; 813% terms[2] = 1; 814% results[0] = u; 815% results[1] = v; 816% LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL); 817% ... 818% if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) { 819% c0 = vectors[0][0]; 820% c2 = vectors[0][1]; 821% c4 = vectors[0][2]; 822% c1 = vectors[1][0]; 823% c3 = vectors[1][1]; 824% c5 = vectors[1][2]; 825% } 826% else 827% printf("Matrix unsolvable\n); 828% RelinquishMagickMatrix(matrix,3UL); 829% RelinquishMagickMatrix(vectors,2UL); 830% 831*/ 832MagickPrivate void LeastSquaresAddTerms(double **matrix,double **vectors, 833 const double *terms,const double *results,const size_t rank, 834 const size_t number_vectors) 835{ 836 register ssize_t 837 i, 838 j; 839 840 for (j=0; j < (ssize_t) rank; j++) 841 { 842 for (i=0; i < (ssize_t) rank; i++) 843 matrix[i][j]+=terms[i]*terms[j]; 844 for (i=0; i < (ssize_t) number_vectors; i++) 845 vectors[i][j]+=results[i]*terms[j]; 846 } 847} 848 849/* 850%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 851% % 852% % 853% % 854% M a t r i x T o I m a g e % 855% % 856% % 857% % 858%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 859% 860% MatrixToImage() returns a matrix as an image. The matrix elements must be 861% of type double otherwise nonsense is returned. 862% 863% The format of the MatrixToImage method is: 864% 865% Image *MatrixToImage(const MatrixInfo *matrix_info, 866% ExceptionInfo *exception) 867% 868% A description of each parameter follows: 869% 870% o matrix_info: the matrix. 871% 872% o exception: return any errors or warnings in this structure. 873% 874*/ 875MagickExport Image *MatrixToImage(const MatrixInfo *matrix_info, 876 ExceptionInfo *exception) 877{ 878 CacheView 879 *image_view; 880 881 double 882 max_value, 883 min_value, 884 scale_factor, 885 value; 886 887 Image 888 *image; 889 890 MagickBooleanType 891 status; 892 893 ssize_t 894 y; 895 896 assert(matrix_info != (const MatrixInfo *) NULL); 897 assert(matrix_info->signature == MagickCoreSignature); 898 assert(exception != (ExceptionInfo *) NULL); 899 assert(exception->signature == MagickCoreSignature); 900 if (matrix_info->stride < sizeof(double)) 901 return((Image *) NULL); 902 /* 903 Determine range of matrix. 904 */ 905 (void) GetMatrixElement(matrix_info,0,0,&value); 906 min_value=value; 907 max_value=value; 908 for (y=0; y < (ssize_t) matrix_info->rows; y++) 909 { 910 register ssize_t 911 x; 912 913 for (x=0; x < (ssize_t) matrix_info->columns; x++) 914 { 915 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse) 916 continue; 917 if (value < min_value) 918 min_value=value; 919 else 920 if (value > max_value) 921 max_value=value; 922 } 923 } 924 if ((min_value == 0.0) && (max_value == 0.0)) 925 scale_factor=0; 926 else 927 if (min_value == max_value) 928 { 929 scale_factor=(double) QuantumRange/min_value; 930 min_value=0; 931 } 932 else 933 scale_factor=(double) QuantumRange/(max_value-min_value); 934 /* 935 Convert matrix to image. 936 */ 937 image=AcquireImage((ImageInfo *) NULL,exception); 938 image->columns=matrix_info->columns; 939 image->rows=matrix_info->rows; 940 image->colorspace=GRAYColorspace; 941 status=MagickTrue; 942 image_view=AcquireAuthenticCacheView(image,exception); 943#if defined(MAGICKCORE_OPENMP_SUPPORT) 944 #pragma omp parallel for schedule(static,4) shared(status) \ 945 magick_threads(image,image,image->rows,1) 946#endif 947 for (y=0; y < (ssize_t) image->rows; y++) 948 { 949 double 950 value; 951 952 register Quantum 953 *q; 954 955 register ssize_t 956 x; 957 958 if (status == MagickFalse) 959 continue; 960 q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception); 961 if (q == (Quantum *) NULL) 962 { 963 status=MagickFalse; 964 continue; 965 } 966 for (x=0; x < (ssize_t) image->columns; x++) 967 { 968 if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse) 969 continue; 970 value=scale_factor*(value-min_value); 971 *q=ClampToQuantum(value); 972 q+=GetPixelChannels(image); 973 } 974 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse) 975 status=MagickFalse; 976 } 977 image_view=DestroyCacheView(image_view); 978 if (status == MagickFalse) 979 image=DestroyImage(image); 980 return(image); 981} 982 983/* 984%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 985% % 986% % 987% % 988% N u l l M a t r i x % 989% % 990% % 991% % 992%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 993% 994% NullMatrix() sets all elements of the matrix to zero. 995% 996% The format of the ResetMagickMemory method is: 997% 998% MagickBooleanType *NullMatrix(MatrixInfo *matrix_info) 999% 1000% A description of each parameter follows: 1001% 1002% o matrix_info: the matrix. 1003% 1004*/ 1005MagickExport MagickBooleanType NullMatrix(MatrixInfo *matrix_info) 1006{ 1007 register ssize_t 1008 x; 1009 1010 ssize_t 1011 count, 1012 y; 1013 1014 unsigned char 1015 value; 1016 1017 assert(matrix_info != (const MatrixInfo *) NULL); 1018 assert(matrix_info->signature == MagickCoreSignature); 1019 if (matrix_info->type != DiskCache) 1020 { 1021 (void) ResetMagickMemory(matrix_info->elements,0,(size_t) 1022 matrix_info->length); 1023 return(MagickTrue); 1024 } 1025 value=0; 1026 (void) lseek(matrix_info->file,0,SEEK_SET); 1027 for (y=0; y < (ssize_t) matrix_info->rows; y++) 1028 { 1029 for (x=0; x < (ssize_t) matrix_info->length; x++) 1030 { 1031 count=write(matrix_info->file,&value,sizeof(value)); 1032 if (count != (ssize_t) sizeof(value)) 1033 break; 1034 } 1035 if (x < (ssize_t) matrix_info->length) 1036 break; 1037 } 1038 return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue); 1039} 1040 1041/* 1042%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1043% % 1044% % 1045% % 1046% R e l i n q u i s h M a g i c k M a t r i x % 1047% % 1048% % 1049% % 1050%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1051% 1052% RelinquishMagickMatrix() frees the previously acquired matrix (array of 1053% pointers to arrays of doubles). 1054% 1055% The format of the RelinquishMagickMatrix method is: 1056% 1057% double **RelinquishMagickMatrix(double **matrix, 1058% const size_t number_rows) 1059% 1060% A description of each parameter follows: 1061% 1062% o matrix: the matrix to relinquish 1063% 1064% o number_rows: the first dimension of the acquired matrix (number of 1065% pointers) 1066% 1067*/ 1068MagickExport double **RelinquishMagickMatrix(double **matrix, 1069 const size_t number_rows) 1070{ 1071 register ssize_t 1072 i; 1073 1074 if (matrix == (double **) NULL ) 1075 return(matrix); 1076 for (i=0; i < (ssize_t) number_rows; i++) 1077 matrix[i]=(double *) RelinquishMagickMemory(matrix[i]); 1078 matrix=(double **) RelinquishMagickMemory(matrix); 1079 return(matrix); 1080} 1081 1082/* 1083%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1084% % 1085% % 1086% % 1087% S e t M a t r i x E l e m e n t % 1088% % 1089% % 1090% % 1091%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1092% 1093% SetMatrixElement() sets the specifed element in the matrix. 1094% 1095% The format of the SetMatrixElement method is: 1096% 1097% MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info, 1098% const ssize_t x,const ssize_t y,void *value) 1099% 1100% A description of each parameter follows: 1101% 1102% o matrix_info: the matrix columns. 1103% 1104% o x: the matrix x-offset. 1105% 1106% o y: the matrix y-offset. 1107% 1108% o value: set the matrix element to this value. 1109% 1110*/ 1111 1112MagickExport MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info, 1113 const ssize_t x,const ssize_t y,const void *value) 1114{ 1115 MagickOffsetType 1116 count, 1117 i; 1118 1119 assert(matrix_info != (const MatrixInfo *) NULL); 1120 assert(matrix_info->signature == MagickCoreSignature); 1121 i=(MagickOffsetType) y*matrix_info->columns+x; 1122 if ((i < 0) || 1123 ((MagickSizeType) (i*matrix_info->stride) >= matrix_info->length)) 1124 return(MagickFalse); 1125 if (matrix_info->type != DiskCache) 1126 { 1127 (void) memcpy((unsigned char *) matrix_info->elements+i* 1128 matrix_info->stride,value,matrix_info->stride); 1129 return(MagickTrue); 1130 } 1131 count=WriteMatrixElements(matrix_info,i*matrix_info->stride, 1132 matrix_info->stride,(unsigned char *) value); 1133 if (count != (MagickOffsetType) matrix_info->stride) 1134 return(MagickFalse); 1135 return(MagickTrue); 1136} 1137