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