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