1 /* 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % % 4 % % 5 % % 6 % DDDD IIIII SSSSS TTTTT OOO RRRR TTTTT % 7 % D D I SS T O O R R T % 8 % D D I SSS T O O RRRR T % 9 % D D I SS T O O R R T % 10 % DDDD IIIII SSSSS T OOO R R T % 11 % % 12 % % 13 % MagickCore Image Distortion Methods % 14 % % 15 % Software Design % 16 % Cristy % 17 % Anthony Thyssen % 18 % June 2007 % 19 % % 20 % % 21 % Copyright 1999-2016 ImageMagick Studio LLC, a non-profit organization % 22 % dedicated to making software imaging solutions freely available. % 23 % % 24 % You may not use this file except in compliance with the License. You may % 25 % obtain a copy of the License at % 26 % % 27 % http://www.imagemagick.org/script/license.php % 28 % % 29 % Unless required by applicable law or agreed to in writing, software % 30 % distributed under the License is distributed on an "AS IS" BASIS, % 31 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. % 32 % See the License for the specific language governing permissions and % 33 % limitations under the License. % 34 % % 35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 36 % 37 % 38 */ 39 40 /* 42 Include declarations. 43 */ 44 #include "MagickCore/studio.h" 45 #include "MagickCore/artifact.h" 46 #include "MagickCore/cache.h" 47 #include "MagickCore/cache-view.h" 48 #include "MagickCore/channel.h" 49 #include "MagickCore/colorspace-private.h" 50 #include "MagickCore/composite-private.h" 51 #include "MagickCore/distort.h" 52 #include "MagickCore/exception.h" 53 #include "MagickCore/exception-private.h" 54 #include "MagickCore/gem.h" 55 #include "MagickCore/image.h" 56 #include "MagickCore/linked-list.h" 57 #include "MagickCore/list.h" 58 #include "MagickCore/matrix.h" 59 #include "MagickCore/matrix-private.h" 60 #include "MagickCore/memory_.h" 61 #include "MagickCore/monitor-private.h" 62 #include "MagickCore/option.h" 63 #include "MagickCore/pixel.h" 64 #include "MagickCore/pixel-accessor.h" 65 #include "MagickCore/pixel-private.h" 66 #include "MagickCore/resample.h" 67 #include "MagickCore/resample-private.h" 68 #include "MagickCore/registry.h" 69 #include "MagickCore/resource_.h" 70 #include "MagickCore/semaphore.h" 71 #include "MagickCore/shear.h" 72 #include "MagickCore/string_.h" 73 #include "MagickCore/string-private.h" 74 #include "MagickCore/thread-private.h" 75 #include "MagickCore/token.h" 76 #include "MagickCore/transform.h" 77 78 /* 80 Numerous internal routines for image distortions. 81 */ 82 static inline void AffineArgsToCoefficients(double *affine) 83 { 84 /* map external sx,ry,rx,sy,tx,ty to internal c0,c2,c4,c1,c3,c5 */ 85 double tmp[4]; /* note indexes 0 and 5 remain unchanged */ 86 tmp[0]=affine[1]; tmp[1]=affine[2]; tmp[2]=affine[3]; tmp[3]=affine[4]; 87 affine[3]=tmp[0]; affine[1]=tmp[1]; affine[4]=tmp[2]; affine[2]=tmp[3]; 88 } 89 90 static inline void CoefficientsToAffineArgs(double *coeff) 91 { 92 /* map internal c0,c1,c2,c3,c4,c5 to external sx,ry,rx,sy,tx,ty */ 93 double tmp[4]; /* note indexes 0 and 5 remain unchanged */ 94 tmp[0]=coeff[3]; tmp[1]=coeff[1]; tmp[2]=coeff[4]; tmp[3]=coeff[2]; 95 coeff[1]=tmp[0]; coeff[2]=tmp[1]; coeff[3]=tmp[2]; coeff[4]=tmp[3]; 96 } 97 static void InvertAffineCoefficients(const double *coeff,double *inverse) 98 { 99 /* From "Digital Image Warping" by George Wolberg, page 50 */ 100 double determinant; 101 102 determinant=PerceptibleReciprocal(coeff[0]*coeff[4]-coeff[1]*coeff[3]); 103 inverse[0]=determinant*coeff[4]; 104 inverse[1]=determinant*(-coeff[1]); 105 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[2]*coeff[4]); 106 inverse[3]=determinant*(-coeff[3]); 107 inverse[4]=determinant*coeff[0]; 108 inverse[5]=determinant*(coeff[2]*coeff[3]-coeff[0]*coeff[5]); 109 } 110 111 static void InvertPerspectiveCoefficients(const double *coeff, 112 double *inverse) 113 { 114 /* From "Digital Image Warping" by George Wolberg, page 53 */ 115 double determinant; 116 117 determinant=PerceptibleReciprocal(coeff[0]*coeff[4]-coeff[3]*coeff[1]); 118 inverse[0]=determinant*(coeff[4]-coeff[7]*coeff[5]); 119 inverse[1]=determinant*(coeff[7]*coeff[2]-coeff[1]); 120 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[4]*coeff[2]); 121 inverse[3]=determinant*(coeff[6]*coeff[5]-coeff[3]); 122 inverse[4]=determinant*(coeff[0]-coeff[6]*coeff[2]); 123 inverse[5]=determinant*(coeff[3]*coeff[2]-coeff[0]*coeff[5]); 124 inverse[6]=determinant*(coeff[3]*coeff[7]-coeff[6]*coeff[4]); 125 inverse[7]=determinant*(coeff[6]*coeff[1]-coeff[0]*coeff[7]); 126 } 127 128 /* 129 * Polynomial Term Defining Functions 130 * 131 * Order must either be an integer, or 1.5 to produce 132 * the 2 number_valuesal polynomial function... 133 * affine 1 (3) u = c0 + c1*x + c2*y 134 * bilinear 1.5 (4) u = '' + c3*x*y 135 * quadratic 2 (6) u = '' + c4*x*x + c5*y*y 136 * cubic 3 (10) u = '' + c6*x^3 + c7*x*x*y + c8*x*y*y + c9*y^3 137 * quartic 4 (15) u = '' + c10*x^4 + ... + c14*y^4 138 * quintic 5 (21) u = '' + c15*x^5 + ... + c20*y^5 139 * number in parenthesis minimum number of points needed. 140 * Anything beyond quintic, has not been implemented until 141 * a more automated way of determining terms is found. 142 143 * Note the slight re-ordering of the terms for a quadratic polynomial 144 * which is to allow the use of a bi-linear (order=1.5) polynomial. 145 * All the later polynomials are ordered simply from x^N to y^N 146 */ 147 static size_t poly_number_terms(double order) 148 { 149 /* Return the number of terms for a 2d polynomial */ 150 if ( order < 1 || order > 5 || 151 ( order != floor(order) && (order-1.5) > MagickEpsilon) ) 152 return 0; /* invalid polynomial order */ 153 return((size_t) floor((order+1)*(order+2)/2)); 154 } 155 156 static double poly_basis_fn(ssize_t n, double x, double y) 157 { 158 /* Return the result for this polynomial term */ 159 switch(n) { 160 case 0: return( 1.0 ); /* constant */ 161 case 1: return( x ); 162 case 2: return( y ); /* affine order = 1 terms = 3 */ 163 case 3: return( x*y ); /* bilinear order = 1.5 terms = 4 */ 164 case 4: return( x*x ); 165 case 5: return( y*y ); /* quadratic order = 2 terms = 6 */ 166 case 6: return( x*x*x ); 167 case 7: return( x*x*y ); 168 case 8: return( x*y*y ); 169 case 9: return( y*y*y ); /* cubic order = 3 terms = 10 */ 170 case 10: return( x*x*x*x ); 171 case 11: return( x*x*x*y ); 172 case 12: return( x*x*y*y ); 173 case 13: return( x*y*y*y ); 174 case 14: return( y*y*y*y ); /* quartic order = 4 terms = 15 */ 175 case 15: return( x*x*x*x*x ); 176 case 16: return( x*x*x*x*y ); 177 case 17: return( x*x*x*y*y ); 178 case 18: return( x*x*y*y*y ); 179 case 19: return( x*y*y*y*y ); 180 case 20: return( y*y*y*y*y ); /* quintic order = 5 terms = 21 */ 181 } 182 return( 0 ); /* should never happen */ 183 } 184 static const char *poly_basis_str(ssize_t n) 185 { 186 /* return the result for this polynomial term */ 187 switch(n) { 188 case 0: return(""); /* constant */ 189 case 1: return("*ii"); 190 case 2: return("*jj"); /* affine order = 1 terms = 3 */ 191 case 3: return("*ii*jj"); /* bilinear order = 1.5 terms = 4 */ 192 case 4: return("*ii*ii"); 193 case 5: return("*jj*jj"); /* quadratic order = 2 terms = 6 */ 194 case 6: return("*ii*ii*ii"); 195 case 7: return("*ii*ii*jj"); 196 case 8: return("*ii*jj*jj"); 197 case 9: return("*jj*jj*jj"); /* cubic order = 3 terms = 10 */ 198 case 10: return("*ii*ii*ii*ii"); 199 case 11: return("*ii*ii*ii*jj"); 200 case 12: return("*ii*ii*jj*jj"); 201 case 13: return("*ii*jj*jj*jj"); 202 case 14: return("*jj*jj*jj*jj"); /* quartic order = 4 terms = 15 */ 203 case 15: return("*ii*ii*ii*ii*ii"); 204 case 16: return("*ii*ii*ii*ii*jj"); 205 case 17: return("*ii*ii*ii*jj*jj"); 206 case 18: return("*ii*ii*jj*jj*jj"); 207 case 19: return("*ii*jj*jj*jj*jj"); 208 case 20: return("*jj*jj*jj*jj*jj"); /* quintic order = 5 terms = 21 */ 209 } 210 return( "UNKNOWN" ); /* should never happen */ 211 } 212 static double poly_basis_dx(ssize_t n, double x, double y) 213 { 214 /* polynomial term for x derivative */ 215 switch(n) { 216 case 0: return( 0.0 ); /* constant */ 217 case 1: return( 1.0 ); 218 case 2: return( 0.0 ); /* affine order = 1 terms = 3 */ 219 case 3: return( y ); /* bilinear order = 1.5 terms = 4 */ 220 case 4: return( x ); 221 case 5: return( 0.0 ); /* quadratic order = 2 terms = 6 */ 222 case 6: return( x*x ); 223 case 7: return( x*y ); 224 case 8: return( y*y ); 225 case 9: return( 0.0 ); /* cubic order = 3 terms = 10 */ 226 case 10: return( x*x*x ); 227 case 11: return( x*x*y ); 228 case 12: return( x*y*y ); 229 case 13: return( y*y*y ); 230 case 14: return( 0.0 ); /* quartic order = 4 terms = 15 */ 231 case 15: return( x*x*x*x ); 232 case 16: return( x*x*x*y ); 233 case 17: return( x*x*y*y ); 234 case 18: return( x*y*y*y ); 235 case 19: return( y*y*y*y ); 236 case 20: return( 0.0 ); /* quintic order = 5 terms = 21 */ 237 } 238 return( 0.0 ); /* should never happen */ 239 } 240 static double poly_basis_dy(ssize_t n, double x, double y) 241 { 242 /* polynomial term for y derivative */ 243 switch(n) { 244 case 0: return( 0.0 ); /* constant */ 245 case 1: return( 0.0 ); 246 case 2: return( 1.0 ); /* affine order = 1 terms = 3 */ 247 case 3: return( x ); /* bilinear order = 1.5 terms = 4 */ 248 case 4: return( 0.0 ); 249 case 5: return( y ); /* quadratic order = 2 terms = 6 */ 250 default: return( poly_basis_dx(n-1,x,y) ); /* weird but true */ 251 } 252 /* NOTE: the only reason that last is not true for 'quadratic' 253 is due to the re-arrangement of terms to allow for 'bilinear' 254 */ 255 } 256 257 /* 259 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 260 % % 261 % % 262 % % 263 % A f f i n e T r a n s f o r m I m a g e % 264 % % 265 % % 266 % % 267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 268 % 269 % AffineTransformImage() transforms an image as dictated by the affine matrix. 270 % It allocates the memory necessary for the new Image structure and returns 271 % a pointer to the new image. 272 % 273 % The format of the AffineTransformImage method is: 274 % 275 % Image *AffineTransformImage(const Image *image, 276 % AffineMatrix *affine_matrix,ExceptionInfo *exception) 277 % 278 % A description of each parameter follows: 279 % 280 % o image: the image. 281 % 282 % o affine_matrix: the affine matrix. 283 % 284 % o exception: return any errors or warnings in this structure. 285 % 286 */ 287 MagickExport Image *AffineTransformImage(const Image *image, 288 const AffineMatrix *affine_matrix,ExceptionInfo *exception) 289 { 290 double 291 distort[6]; 292 293 Image 294 *deskew_image; 295 296 /* 297 Affine transform image. 298 */ 299 assert(image->signature == MagickCoreSignature); 300 if (image->debug != MagickFalse) 301 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 302 assert(affine_matrix != (AffineMatrix *) NULL); 303 assert(exception != (ExceptionInfo *) NULL); 304 assert(exception->signature == MagickCoreSignature); 305 distort[0]=affine_matrix->sx; 306 distort[1]=affine_matrix->rx; 307 distort[2]=affine_matrix->ry; 308 distort[3]=affine_matrix->sy; 309 distort[4]=affine_matrix->tx; 310 distort[5]=affine_matrix->ty; 311 deskew_image=DistortImage(image,AffineProjectionDistortion,6,distort, 312 MagickTrue,exception); 313 return(deskew_image); 314 } 315 316 /* 318 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 319 % % 320 % % 321 % % 322 + G e n e r a t e C o e f f i c i e n t s % 323 % % 324 % % 325 % % 326 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 327 % 328 % GenerateCoefficients() takes user provided input arguments and generates 329 % the coefficients, needed to apply the specific distortion for either 330 % distorting images (generally using control points) or generating a color 331 % gradient from sparsely separated color points. 332 % 333 % The format of the GenerateCoefficients() method is: 334 % 335 % Image *GenerateCoefficients(const Image *image,DistortMethod method, 336 % const size_t number_arguments,const double *arguments, 337 % size_t number_values, ExceptionInfo *exception) 338 % 339 % A description of each parameter follows: 340 % 341 % o image: the image to be distorted. 342 % 343 % o method: the method of image distortion/ sparse gradient 344 % 345 % o number_arguments: the number of arguments given. 346 % 347 % o arguments: the arguments for this distortion method. 348 % 349 % o number_values: the style and format of given control points, (caller type) 350 % 0: 2 dimensional mapping of control points (Distort) 351 % Format: u,v,x,y where u,v is the 'source' of the 352 % the color to be plotted, for DistortImage() 353 % N: Interpolation of control points with N values (usally r,g,b) 354 % Format: x,y,r,g,b mapping x,y to color values r,g,b 355 % IN future, variable number of values may be given (1 to N) 356 % 357 % o exception: return any errors or warnings in this structure 358 % 359 % Note that the returned array of double values must be freed by the 360 % calling method using RelinquishMagickMemory(). This however may change in 361 % the future to require a more 'method' specific method. 362 % 363 % Because of this this method should not be classed as stable or used 364 % outside other MagickCore library methods. 365 */ 366 367 static inline double MagickRound(double x) 368 { 369 /* 370 Round the fraction to nearest integer. 371 */ 372 if ((x-floor(x)) < (ceil(x)-x)) 373 return(floor(x)); 374 return(ceil(x)); 375 } 376 377 static double *GenerateCoefficients(const Image *image, 378 DistortMethod *method,const size_t number_arguments,const double *arguments, 379 size_t number_values,ExceptionInfo *exception) 380 { 381 double 382 *coeff; 383 384 register size_t 385 i; 386 387 size_t 388 number_coeff, /* number of coefficients to return (array size) */ 389 cp_size, /* number floating point numbers per control point */ 390 cp_x,cp_y, /* the x,y indexes for control point */ 391 cp_values; /* index of values for this control point */ 392 /* number_values Number of values given per control point */ 393 394 if ( number_values == 0 ) { 395 /* Image distortion using control points (or other distortion) 396 That is generate a mapping so that x,y->u,v given u,v,x,y 397 */ 398 number_values = 2; /* special case: two values of u,v */ 399 cp_values = 0; /* the values i,j are BEFORE the destination CP x,y */ 400 cp_x = 2; /* location of x,y in input control values */ 401 cp_y = 3; 402 /* NOTE: cp_values, also used for later 'reverse map distort' tests */ 403 } 404 else { 405 cp_x = 0; /* location of x,y in input control values */ 406 cp_y = 1; 407 cp_values = 2; /* and the other values are after x,y */ 408 /* Typically in this case the values are R,G,B color values */ 409 } 410 cp_size = number_values+2; /* each CP defintion involves this many numbers */ 411 412 /* If not enough control point pairs are found for specific distortions 413 fall back to Affine distortion (allowing 0 to 3 point pairs) 414 */ 415 if ( number_arguments < 4*cp_size && 416 ( *method == BilinearForwardDistortion 417 || *method == BilinearReverseDistortion 418 || *method == PerspectiveDistortion 419 ) ) 420 *method = AffineDistortion; 421 422 number_coeff=0; 423 switch (*method) { 424 case AffineDistortion: 425 /* also BarycentricColorInterpolate: */ 426 number_coeff=3*number_values; 427 break; 428 case PolynomialDistortion: 429 /* number of coefficents depend on the given polynomal 'order' */ 430 i = poly_number_terms(arguments[0]); 431 number_coeff = 2 + i*number_values; 432 if ( i == 0 ) { 433 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 434 "InvalidArgument","%s : '%s'","Polynomial", 435 "Invalid order, should be interger 1 to 5, or 1.5"); 436 return((double *) NULL); 437 } 438 if ( number_arguments < 1+i*cp_size ) { 439 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 440 "InvalidArgument", "%s : 'require at least %.20g CPs'", 441 "Polynomial", (double) i); 442 return((double *) NULL); 443 } 444 break; 445 case BilinearReverseDistortion: 446 number_coeff=4*number_values; 447 break; 448 /* 449 The rest are constants as they are only used for image distorts 450 */ 451 case BilinearForwardDistortion: 452 number_coeff=10; /* 2*4 coeff plus 2 constants */ 453 cp_x = 0; /* Reverse src/dest coords for forward mapping */ 454 cp_y = 1; 455 cp_values = 2; 456 break; 457 #if 0 458 case QuadraterialDistortion: 459 number_coeff=19; /* BilinearForward + BilinearReverse */ 460 #endif 461 break; 462 case ShepardsDistortion: 463 number_coeff=1; /* The power factor to use */ 464 break; 465 case ArcDistortion: 466 number_coeff=5; 467 break; 468 case ScaleRotateTranslateDistortion: 469 case AffineProjectionDistortion: 470 case Plane2CylinderDistortion: 471 case Cylinder2PlaneDistortion: 472 number_coeff=6; 473 break; 474 case PolarDistortion: 475 case DePolarDistortion: 476 number_coeff=8; 477 break; 478 case PerspectiveDistortion: 479 case PerspectiveProjectionDistortion: 480 number_coeff=9; 481 break; 482 case BarrelDistortion: 483 case BarrelInverseDistortion: 484 number_coeff=10; 485 break; 486 default: 487 perror("unknown method given"); /* just fail assertion */ 488 } 489 490 /* allocate the array of coefficients needed */ 491 coeff = (double *) AcquireQuantumMemory(number_coeff,sizeof(*coeff)); 492 if (coeff == (double *) NULL) { 493 (void) ThrowMagickException(exception,GetMagickModule(), 494 ResourceLimitError,"MemoryAllocationFailed", 495 "%s", "GenerateCoefficients"); 496 return((double *) NULL); 497 } 498 499 /* zero out coefficients array */ 500 for (i=0; i < number_coeff; i++) 501 coeff[i] = 0.0; 502 503 switch (*method) 504 { 505 case AffineDistortion: 506 { 507 /* Affine Distortion 508 v = c0*x + c1*y + c2 509 for each 'value' given 510 511 Input Arguments are sets of control points... 512 For Distort Images u,v, x,y ... 513 For Sparse Gradients x,y, r,g,b ... 514 */ 515 if ( number_arguments%cp_size != 0 || 516 number_arguments < cp_size ) { 517 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 518 "InvalidArgument", "%s : 'require at least %.20g CPs'", 519 "Affine", 1.0); 520 coeff=(double *) RelinquishMagickMemory(coeff); 521 return((double *) NULL); 522 } 523 /* handle special cases of not enough arguments */ 524 if ( number_arguments == cp_size ) { 525 /* Only 1 CP Set Given */ 526 if ( cp_values == 0 ) { 527 /* image distortion - translate the image */ 528 coeff[0] = 1.0; 529 coeff[2] = arguments[0] - arguments[2]; 530 coeff[4] = 1.0; 531 coeff[5] = arguments[1] - arguments[3]; 532 } 533 else { 534 /* sparse gradient - use the values directly */ 535 for (i=0; i<number_values; i++) 536 coeff[i*3+2] = arguments[cp_values+i]; 537 } 538 } 539 else { 540 /* 2 or more points (usally 3) given. 541 Solve a least squares simultaneous equation for coefficients. 542 */ 543 double 544 **matrix, 545 **vectors, 546 terms[3]; 547 548 MagickBooleanType 549 status; 550 551 /* create matrix, and a fake vectors matrix */ 552 matrix = AcquireMagickMatrix(3UL,3UL); 553 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors)); 554 if (matrix == (double **) NULL || vectors == (double **) NULL) 555 { 556 matrix = RelinquishMagickMatrix(matrix, 3UL); 557 vectors = (double **) RelinquishMagickMemory(vectors); 558 coeff = (double *) RelinquishMagickMemory(coeff); 559 (void) ThrowMagickException(exception,GetMagickModule(), 560 ResourceLimitError,"MemoryAllocationFailed", 561 "%s", "DistortCoefficients"); 562 return((double *) NULL); 563 } 564 /* fake a number_values x3 vectors matrix from coefficients array */ 565 for (i=0; i < number_values; i++) 566 vectors[i] = &(coeff[i*3]); 567 /* Add given control point pairs for least squares solving */ 568 for (i=0; i < number_arguments; i+=cp_size) { 569 terms[0] = arguments[i+cp_x]; /* x */ 570 terms[1] = arguments[i+cp_y]; /* y */ 571 terms[2] = 1; /* 1 */ 572 LeastSquaresAddTerms(matrix,vectors,terms, 573 &(arguments[i+cp_values]),3UL,number_values); 574 } 575 if ( number_arguments == 2*cp_size ) { 576 /* Only two pairs were given, but we need 3 to solve the affine. 577 Fake extra coordinates by rotating p1 around p0 by 90 degrees. 578 x2 = x0 - (y1-y0) y2 = y0 + (x1-x0) 579 */ 580 terms[0] = arguments[cp_x] 581 - ( arguments[cp_size+cp_y] - arguments[cp_y] ); /* x2 */ 582 terms[1] = arguments[cp_y] + 583 + ( arguments[cp_size+cp_x] - arguments[cp_x] ); /* y2 */ 584 terms[2] = 1; /* 1 */ 585 if ( cp_values == 0 ) { 586 /* Image Distortion - rotate the u,v coordients too */ 587 double 588 uv2[2]; 589 uv2[0] = arguments[0] - arguments[5] + arguments[1]; /* u2 */ 590 uv2[1] = arguments[1] + arguments[4] - arguments[0]; /* v2 */ 591 LeastSquaresAddTerms(matrix,vectors,terms,uv2,3UL,2UL); 592 } 593 else { 594 /* Sparse Gradient - use values of p0 for linear gradient */ 595 LeastSquaresAddTerms(matrix,vectors,terms, 596 &(arguments[cp_values]),3UL,number_values); 597 } 598 } 599 /* Solve for LeastSquares Coefficients */ 600 status=GaussJordanElimination(matrix,vectors,3UL,number_values); 601 matrix = RelinquishMagickMatrix(matrix, 3UL); 602 vectors = (double **) RelinquishMagickMemory(vectors); 603 if ( status == MagickFalse ) { 604 coeff = (double *) RelinquishMagickMemory(coeff); 605 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 606 "InvalidArgument","%s : 'Unsolvable Matrix'", 607 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 608 return((double *) NULL); 609 } 610 } 611 return(coeff); 612 } 613 case AffineProjectionDistortion: 614 { 615 /* 616 Arguments: Affine Matrix (forward mapping) 617 Arguments sx, rx, ry, sy, tx, ty 618 Where u = sx*x + ry*y + tx 619 v = rx*x + sy*y + ty 620 621 Returns coefficients (in there inverse form) ordered as... 622 sx ry tx rx sy ty 623 624 AffineProjection Distortion Notes... 625 + Will only work with a 2 number_values for Image Distortion 626 + Can not be used for generating a sparse gradient (interpolation) 627 */ 628 double inverse[8]; 629 if (number_arguments != 6) { 630 coeff = (double *) RelinquishMagickMemory(coeff); 631 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 632 "InvalidArgument","%s : 'Needs 6 coeff values'", 633 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 634 return((double *) NULL); 635 } 636 /* FUTURE: trap test for sx*sy-rx*ry == 0 (determinant = 0, no inverse) */ 637 for(i=0; i<6UL; i++ ) 638 inverse[i] = arguments[i]; 639 AffineArgsToCoefficients(inverse); /* map into coefficents */ 640 InvertAffineCoefficients(inverse, coeff); /* invert */ 641 *method = AffineDistortion; 642 643 return(coeff); 644 } 645 case ScaleRotateTranslateDistortion: 646 { 647 /* Scale, Rotate and Translate Distortion 648 An alternative Affine Distortion 649 Argument options, by number of arguments given: 650 7: x,y, sx,sy, a, nx,ny 651 6: x,y, s, a, nx,ny 652 5: x,y, sx,sy, a 653 4: x,y, s, a 654 3: x,y, a 655 2: s, a 656 1: a 657 Where actions are (in order of application) 658 x,y 'center' of transforms (default = image center) 659 sx,sy scale image by this amount (default = 1) 660 a angle of rotation (argument required) 661 nx,ny move 'center' here (default = x,y or no movement) 662 And convert to affine mapping coefficients 663 664 ScaleRotateTranslate Distortion Notes... 665 + Does not use a set of CPs in any normal way 666 + Will only work with a 2 number_valuesal Image Distortion 667 + Cannot be used for generating a sparse gradient (interpolation) 668 */ 669 double 670 cosine, sine, 671 x,y,sx,sy,a,nx,ny; 672 673 /* set default center, and default scale */ 674 x = nx = (double)(image->columns)/2.0 + (double)image->page.x; 675 y = ny = (double)(image->rows)/2.0 + (double)image->page.y; 676 sx = sy = 1.0; 677 switch ( number_arguments ) { 678 case 0: 679 coeff = (double *) RelinquishMagickMemory(coeff); 680 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 681 "InvalidArgument","%s : 'Needs at least 1 argument'", 682 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 683 return((double *) NULL); 684 case 1: 685 a = arguments[0]; 686 break; 687 case 2: 688 sx = sy = arguments[0]; 689 a = arguments[1]; 690 break; 691 default: 692 x = nx = arguments[0]; 693 y = ny = arguments[1]; 694 switch ( number_arguments ) { 695 case 3: 696 a = arguments[2]; 697 break; 698 case 4: 699 sx = sy = arguments[2]; 700 a = arguments[3]; 701 break; 702 case 5: 703 sx = arguments[2]; 704 sy = arguments[3]; 705 a = arguments[4]; 706 break; 707 case 6: 708 sx = sy = arguments[2]; 709 a = arguments[3]; 710 nx = arguments[4]; 711 ny = arguments[5]; 712 break; 713 case 7: 714 sx = arguments[2]; 715 sy = arguments[3]; 716 a = arguments[4]; 717 nx = arguments[5]; 718 ny = arguments[6]; 719 break; 720 default: 721 coeff = (double *) RelinquishMagickMemory(coeff); 722 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 723 "InvalidArgument","%s : 'Too Many Arguments (7 or less)'", 724 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 725 return((double *) NULL); 726 } 727 break; 728 } 729 /* Trap if sx or sy == 0 -- image is scaled out of existance! */ 730 if ( fabs(sx) < MagickEpsilon || fabs(sy) < MagickEpsilon ) { 731 coeff = (double *) RelinquishMagickMemory(coeff); 732 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 733 "InvalidArgument","%s : 'Zero Scale Given'", 734 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 735 return((double *) NULL); 736 } 737 /* Save the given arguments as an affine distortion */ 738 a=DegreesToRadians(a); cosine=cos(a); sine=sin(a); 739 740 *method = AffineDistortion; 741 coeff[0]=cosine/sx; 742 coeff[1]=sine/sx; 743 coeff[2]=x-nx*coeff[0]-ny*coeff[1]; 744 coeff[3]=(-sine)/sy; 745 coeff[4]=cosine/sy; 746 coeff[5]=y-nx*coeff[3]-ny*coeff[4]; 747 return(coeff); 748 } 749 case PerspectiveDistortion: 750 { /* 751 Perspective Distortion (a ratio of affine distortions) 752 753 p(x,y) c0*x + c1*y + c2 754 u = ------ = ------------------ 755 r(x,y) c6*x + c7*y + 1 756 757 q(x,y) c3*x + c4*y + c5 758 v = ------ = ------------------ 759 r(x,y) c6*x + c7*y + 1 760 761 c8 = Sign of 'r', or the denominator affine, for the actual image. 762 This determines what part of the distorted image is 'ground' 763 side of the horizon, the other part is 'sky' or invalid. 764 Valid values are +1.0 or -1.0 only. 765 766 Input Arguments are sets of control points... 767 For Distort Images u,v, x,y ... 768 For Sparse Gradients x,y, r,g,b ... 769 770 Perspective Distortion Notes... 771 + Can be thought of as ratio of 3 affine transformations 772 + Not separatable: r() or c6 and c7 are used by both equations 773 + All 8 coefficients must be determined simultaniously 774 + Will only work with a 2 number_valuesal Image Distortion 775 + Can not be used for generating a sparse gradient (interpolation) 776 + It is not linear, but is simple to generate an inverse 777 + All lines within an image remain lines. 778 + but distances between points may vary. 779 */ 780 double 781 **matrix, 782 *vectors[1], 783 terms[8]; 784 785 size_t 786 cp_u = cp_values, 787 cp_v = cp_values+1; 788 789 MagickBooleanType 790 status; 791 792 if ( number_arguments%cp_size != 0 || 793 number_arguments < cp_size*4 ) { 794 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 795 "InvalidArgument", "%s : 'require at least %.20g CPs'", 796 CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0); 797 coeff=(double *) RelinquishMagickMemory(coeff); 798 return((double *) NULL); 799 } 800 /* fake 1x8 vectors matrix directly using the coefficients array */ 801 vectors[0] = &(coeff[0]); 802 /* 8x8 least-squares matrix (zeroed) */ 803 matrix = AcquireMagickMatrix(8UL,8UL); 804 if (matrix == (double **) NULL) { 805 (void) ThrowMagickException(exception,GetMagickModule(), 806 ResourceLimitError,"MemoryAllocationFailed", 807 "%s", "DistortCoefficients"); 808 return((double *) NULL); 809 } 810 /* Add control points for least squares solving */ 811 for (i=0; i < number_arguments; i+=4) { 812 terms[0]=arguments[i+cp_x]; /* c0*x */ 813 terms[1]=arguments[i+cp_y]; /* c1*y */ 814 terms[2]=1.0; /* c2*1 */ 815 terms[3]=0.0; 816 terms[4]=0.0; 817 terms[5]=0.0; 818 terms[6]=-terms[0]*arguments[i+cp_u]; /* 1/(c6*x) */ 819 terms[7]=-terms[1]*arguments[i+cp_u]; /* 1/(c7*y) */ 820 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_u]), 821 8UL,1UL); 822 823 terms[0]=0.0; 824 terms[1]=0.0; 825 terms[2]=0.0; 826 terms[3]=arguments[i+cp_x]; /* c3*x */ 827 terms[4]=arguments[i+cp_y]; /* c4*y */ 828 terms[5]=1.0; /* c5*1 */ 829 terms[6]=-terms[3]*arguments[i+cp_v]; /* 1/(c6*x) */ 830 terms[7]=-terms[4]*arguments[i+cp_v]; /* 1/(c7*y) */ 831 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_v]), 832 8UL,1UL); 833 } 834 /* Solve for LeastSquares Coefficients */ 835 status=GaussJordanElimination(matrix,vectors,8UL,1UL); 836 matrix = RelinquishMagickMatrix(matrix, 8UL); 837 if ( status == MagickFalse ) { 838 coeff = (double *) RelinquishMagickMemory(coeff); 839 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 840 "InvalidArgument","%s : 'Unsolvable Matrix'", 841 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 842 return((double *) NULL); 843 } 844 /* 845 Calculate 9'th coefficient! The ground-sky determination. 846 What is sign of the 'ground' in r() denominator affine function? 847 Just use any valid image coordinate (first control point) in 848 destination for determination of what part of view is 'ground'. 849 */ 850 coeff[8] = coeff[6]*arguments[cp_x] 851 + coeff[7]*arguments[cp_y] + 1.0; 852 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0; 853 854 return(coeff); 855 } 856 case PerspectiveProjectionDistortion: 857 { 858 /* 859 Arguments: Perspective Coefficents (forward mapping) 860 */ 861 if (number_arguments != 8) { 862 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 863 "InvalidArgument", "%s : 'Needs 8 coefficient values'", 864 CommandOptionToMnemonic(MagickDistortOptions, *method)); 865 return((double *) NULL); 866 } 867 /* FUTURE: trap test c0*c4-c3*c1 == 0 (determinate = 0, no inverse) */ 868 InvertPerspectiveCoefficients(arguments, coeff); 869 /* 870 Calculate 9'th coefficient! The ground-sky determination. 871 What is sign of the 'ground' in r() denominator affine function? 872 Just use any valid image cocodinate in destination for determination. 873 For a forward mapped perspective the images 0,0 coord will map to 874 c2,c5 in the distorted image, so set the sign of denominator of that. 875 */ 876 coeff[8] = coeff[6]*arguments[2] 877 + coeff[7]*arguments[5] + 1.0; 878 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0; 879 *method = PerspectiveDistortion; 880 881 return(coeff); 882 } 883 case BilinearForwardDistortion: 884 case BilinearReverseDistortion: 885 { 886 /* Bilinear Distortion (Forward mapping) 887 v = c0*x + c1*y + c2*x*y + c3; 888 for each 'value' given 889 890 This is actually a simple polynomial Distortion! The difference 891 however is when we need to reverse the above equation to generate a 892 BilinearForwardDistortion (see below). 893 894 Input Arguments are sets of control points... 895 For Distort Images u,v, x,y ... 896 For Sparse Gradients x,y, r,g,b ... 897 898 */ 899 double 900 **matrix, 901 **vectors, 902 terms[4]; 903 904 MagickBooleanType 905 status; 906 907 /* check the number of arguments */ 908 if ( number_arguments%cp_size != 0 || 909 number_arguments < cp_size*4 ) { 910 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 911 "InvalidArgument", "%s : 'require at least %.20g CPs'", 912 CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0); 913 coeff=(double *) RelinquishMagickMemory(coeff); 914 return((double *) NULL); 915 } 916 /* create matrix, and a fake vectors matrix */ 917 matrix = AcquireMagickMatrix(4UL,4UL); 918 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors)); 919 if (matrix == (double **) NULL || vectors == (double **) NULL) 920 { 921 matrix = RelinquishMagickMatrix(matrix, 4UL); 922 vectors = (double **) RelinquishMagickMemory(vectors); 923 coeff = (double *) RelinquishMagickMemory(coeff); 924 (void) ThrowMagickException(exception,GetMagickModule(), 925 ResourceLimitError,"MemoryAllocationFailed", 926 "%s", "DistortCoefficients"); 927 return((double *) NULL); 928 } 929 /* fake a number_values x4 vectors matrix from coefficients array */ 930 for (i=0; i < number_values; i++) 931 vectors[i] = &(coeff[i*4]); 932 /* Add given control point pairs for least squares solving */ 933 for (i=0; i < number_arguments; i+=cp_size) { 934 terms[0] = arguments[i+cp_x]; /* x */ 935 terms[1] = arguments[i+cp_y]; /* y */ 936 terms[2] = terms[0]*terms[1]; /* x*y */ 937 terms[3] = 1; /* 1 */ 938 LeastSquaresAddTerms(matrix,vectors,terms, 939 &(arguments[i+cp_values]),4UL,number_values); 940 } 941 /* Solve for LeastSquares Coefficients */ 942 status=GaussJordanElimination(matrix,vectors,4UL,number_values); 943 matrix = RelinquishMagickMatrix(matrix, 4UL); 944 vectors = (double **) RelinquishMagickMemory(vectors); 945 if ( status == MagickFalse ) { 946 coeff = (double *) RelinquishMagickMemory(coeff); 947 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 948 "InvalidArgument","%s : 'Unsolvable Matrix'", 949 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 950 return((double *) NULL); 951 } 952 if ( *method == BilinearForwardDistortion ) { 953 /* Bilinear Forward Mapped Distortion 954 955 The above least-squares solved for coefficents but in the forward 956 direction, due to changes to indexing constants. 957 958 i = c0*x + c1*y + c2*x*y + c3; 959 j = c4*x + c5*y + c6*x*y + c7; 960 961 where i,j are in the destination image, NOT the source. 962 963 Reverse Pixel mapping however needs to use reverse of these 964 functions. It required a full page of algbra to work out the 965 reversed mapping formula, but resolves down to the following... 966 967 c8 = c0*c5-c1*c4; 968 c9 = 2*(c2*c5-c1*c6); // '2*a' in the quadratic formula 969 970 i = i - c3; j = j - c7; 971 b = c6*i - c2*j + c8; // So that a*y^2 + b*y + c == 0 972 c = c4*i - c0*j; // y = ( -b +- sqrt(bb - 4ac) ) / (2*a) 973 974 r = b*b - c9*(c+c); 975 if ( c9 != 0 ) 976 y = ( -b + sqrt(r) ) / c9; 977 else 978 y = -c/b; 979 980 x = ( i - c1*y) / ( c1 - c2*y ); 981 982 NB: if 'r' is negative there is no solution! 983 NB: the sign of the sqrt() should be negative if image becomes 984 flipped or flopped, or crosses over itself. 985 NB: techniqually coefficient c5 is not needed, anymore, 986 but kept for completness. 987 988 See Anthony Thyssen <A.Thyssen (at) griffith.edu.au> 989 or Fred Weinhaus <fmw (at) alink.net> for more details. 990 991 */ 992 coeff[8] = coeff[0]*coeff[5] - coeff[1]*coeff[4]; 993 coeff[9] = 2*(coeff[2]*coeff[5] - coeff[1]*coeff[6]); 994 } 995 return(coeff); 996 } 997 #if 0 998 case QuadrilateralDistortion: 999 { 1000 /* Map a Quadrilateral to a unit square using BilinearReverse 1001 Then map that unit square back to the final Quadrilateral 1002 using BilinearForward. 1003 1004 Input Arguments are sets of control points... 1005 For Distort Images u,v, x,y ... 1006 For Sparse Gradients x,y, r,g,b ... 1007 1008 */ 1009 /* UNDER CONSTRUCTION */ 1010 return(coeff); 1011 } 1012 #endif 1013 1014 case PolynomialDistortion: 1015 { 1016 /* Polynomial Distortion 1017 1018 First two coefficents are used to hole global polynomal information 1019 c0 = Order of the polynimial being created 1020 c1 = number_of_terms in one polynomial equation 1021 1022 Rest of the coefficients map to the equations.... 1023 v = c0 + c1*x + c2*y + c3*x*y + c4*x^2 + c5*y^2 + c6*x^3 + ... 1024 for each 'value' (number_values of them) given. 1025 As such total coefficients = 2 + number_terms * number_values 1026 1027 Input Arguments are sets of control points... 1028 For Distort Images order [u,v, x,y] ... 1029 For Sparse Gradients order [x,y, r,g,b] ... 1030 1031 Polynomial Distortion Notes... 1032 + UNDER DEVELOPMENT -- Do not expect this to remain as is. 1033 + Currently polynomial is a reversed mapped distortion. 1034 + Order 1.5 is fudged to map into a bilinear distortion. 1035 though it is not the same order as that distortion. 1036 */ 1037 double 1038 **matrix, 1039 **vectors, 1040 *terms; 1041 1042 size_t 1043 nterms; /* number of polynomial terms per number_values */ 1044 1045 register ssize_t 1046 j; 1047 1048 MagickBooleanType 1049 status; 1050 1051 /* first two coefficients hold polynomial order information */ 1052 coeff[0] = arguments[0]; 1053 coeff[1] = (double) poly_number_terms(arguments[0]); 1054 nterms = (size_t) coeff[1]; 1055 1056 /* create matrix, a fake vectors matrix, and least sqs terms */ 1057 matrix = AcquireMagickMatrix(nterms,nterms); 1058 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors)); 1059 terms = (double *) AcquireQuantumMemory(nterms, sizeof(*terms)); 1060 if (matrix == (double **) NULL || 1061 vectors == (double **) NULL || 1062 terms == (double *) NULL ) 1063 { 1064 matrix = RelinquishMagickMatrix(matrix, nterms); 1065 vectors = (double **) RelinquishMagickMemory(vectors); 1066 terms = (double *) RelinquishMagickMemory(terms); 1067 coeff = (double *) RelinquishMagickMemory(coeff); 1068 (void) ThrowMagickException(exception,GetMagickModule(), 1069 ResourceLimitError,"MemoryAllocationFailed", 1070 "%s", "DistortCoefficients"); 1071 return((double *) NULL); 1072 } 1073 /* fake a number_values x3 vectors matrix from coefficients array */ 1074 for (i=0; i < number_values; i++) 1075 vectors[i] = &(coeff[2+i*nterms]); 1076 /* Add given control point pairs for least squares solving */ 1077 for (i=1; i < number_arguments; i+=cp_size) { /* NB: start = 1 not 0 */ 1078 for (j=0; j < (ssize_t) nterms; j++) 1079 terms[j] = poly_basis_fn(j,arguments[i+cp_x],arguments[i+cp_y]); 1080 LeastSquaresAddTerms(matrix,vectors,terms, 1081 &(arguments[i+cp_values]),nterms,number_values); 1082 } 1083 terms = (double *) RelinquishMagickMemory(terms); 1084 /* Solve for LeastSquares Coefficients */ 1085 status=GaussJordanElimination(matrix,vectors,nterms,number_values); 1086 matrix = RelinquishMagickMatrix(matrix, nterms); 1087 vectors = (double **) RelinquishMagickMemory(vectors); 1088 if ( status == MagickFalse ) { 1089 coeff = (double *) RelinquishMagickMemory(coeff); 1090 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 1091 "InvalidArgument","%s : 'Unsolvable Matrix'", 1092 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 1093 return((double *) NULL); 1094 } 1095 return(coeff); 1096 } 1097 case ArcDistortion: 1098 { 1099 /* Arc Distortion 1100 Args: arc_width rotate top_edge_radius bottom_edge_radius 1101 All but first argument are optional 1102 arc_width The angle over which to arc the image side-to-side 1103 rotate Angle to rotate image from vertical center 1104 top_radius Set top edge of source image at this radius 1105 bottom_radius Set bootom edge to this radius (radial scaling) 1106 1107 By default, if the radii arguments are nor provided the image radius 1108 is calculated so the horizontal center-line is fits the given arc 1109 without scaling. 1110 1111 The output image size is ALWAYS adjusted to contain the whole image, 1112 and an offset is given to position image relative to the 0,0 point of 1113 the origin, allowing users to use relative positioning onto larger 1114 background (via -flatten). 1115 1116 The arguments are converted to these coefficients 1117 c0: angle for center of source image 1118 c1: angle scale for mapping to source image 1119 c2: radius for top of source image 1120 c3: radius scale for mapping source image 1121 c4: centerline of arc within source image 1122 1123 Note the coefficients use a center angle, so asymptotic join is 1124 furthest from both sides of the source image. This also means that 1125 for arc angles greater than 360 the sides of the image will be 1126 trimmed equally. 1127 1128 Arc Distortion Notes... 1129 + Does not use a set of CPs 1130 + Will only work with Image Distortion 1131 + Can not be used for generating a sparse gradient (interpolation) 1132 */ 1133 if ( number_arguments >= 1 && arguments[0] < MagickEpsilon ) { 1134 coeff = (double *) RelinquishMagickMemory(coeff); 1135 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 1136 "InvalidArgument","%s : 'Arc Angle Too Small'", 1137 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 1138 return((double *) NULL); 1139 } 1140 if ( number_arguments >= 3 && arguments[2] < MagickEpsilon ) { 1141 coeff = (double *) RelinquishMagickMemory(coeff); 1142 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 1143 "InvalidArgument","%s : 'Outer Radius Too Small'", 1144 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 1145 return((double *) NULL); 1146 } 1147 coeff[0] = -MagickPI2; /* -90, place at top! */ 1148 if ( number_arguments >= 1 ) 1149 coeff[1] = DegreesToRadians(arguments[0]); 1150 else 1151 coeff[1] = MagickPI2; /* zero arguments - center is at top */ 1152 if ( number_arguments >= 2 ) 1153 coeff[0] += DegreesToRadians(arguments[1]); 1154 coeff[0] /= Magick2PI; /* normalize radians */ 1155 coeff[0] -= MagickRound(coeff[0]); 1156 coeff[0] *= Magick2PI; /* de-normalize back to radians */ 1157 coeff[3] = (double)image->rows-1; 1158 coeff[2] = (double)image->columns/coeff[1] + coeff[3]/2.0; 1159 if ( number_arguments >= 3 ) { 1160 if ( number_arguments >= 4 ) 1161 coeff[3] = arguments[2] - arguments[3]; 1162 else 1163 coeff[3] *= arguments[2]/coeff[2]; 1164 coeff[2] = arguments[2]; 1165 } 1166 coeff[4] = ((double)image->columns-1.0)/2.0; 1167 1168 return(coeff); 1169 } 1170 case PolarDistortion: 1171 case DePolarDistortion: 1172 { 1173 /* (De)Polar Distortion (same set of arguments) 1174 Args: Rmax, Rmin, Xcenter,Ycenter, Afrom,Ato 1175 DePolar can also have the extra arguments of Width, Height 1176 1177 Coefficients 0 to 5 is the sanatized version first 6 input args 1178 Coefficient 6 is the angle to coord ratio and visa-versa 1179 Coefficient 7 is the radius to coord ratio and visa-versa 1180 1181 WARNING: It is possible for Radius max<min and/or Angle from>to 1182 */ 1183 if ( number_arguments == 3 1184 || ( number_arguments > 6 && *method == PolarDistortion ) 1185 || number_arguments > 8 ) { 1186 (void) ThrowMagickException(exception,GetMagickModule(), 1187 OptionError,"InvalidArgument", "%s : number of arguments", 1188 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 1189 coeff=(double *) RelinquishMagickMemory(coeff); 1190 return((double *) NULL); 1191 } 1192 /* Rmax - if 0 calculate appropriate value */ 1193 if ( number_arguments >= 1 ) 1194 coeff[0] = arguments[0]; 1195 else 1196 coeff[0] = 0.0; 1197 /* Rmin - usally 0 */ 1198 coeff[1] = number_arguments >= 2 ? arguments[1] : 0.0; 1199 /* Center X,Y */ 1200 if ( number_arguments >= 4 ) { 1201 coeff[2] = arguments[2]; 1202 coeff[3] = arguments[3]; 1203 } 1204 else { /* center of actual image */ 1205 coeff[2] = (double)(image->columns)/2.0+image->page.x; 1206 coeff[3] = (double)(image->rows)/2.0+image->page.y; 1207 } 1208 /* Angle from,to - about polar center 0 is downward */ 1209 coeff[4] = -MagickPI; 1210 if ( number_arguments >= 5 ) 1211 coeff[4] = DegreesToRadians(arguments[4]); 1212 coeff[5] = coeff[4]; 1213 if ( number_arguments >= 6 ) 1214 coeff[5] = DegreesToRadians(arguments[5]); 1215 if ( fabs(coeff[4]-coeff[5]) < MagickEpsilon ) 1216 coeff[5] += Magick2PI; /* same angle is a full circle */ 1217 /* if radius 0 or negative, its a special value... */ 1218 if ( coeff[0] < MagickEpsilon ) { 1219 /* Use closest edge if radius == 0 */ 1220 if ( fabs(coeff[0]) < MagickEpsilon ) { 1221 coeff[0]=MagickMin(fabs(coeff[2]-image->page.x), 1222 fabs(coeff[3]-image->page.y)); 1223 coeff[0]=MagickMin(coeff[0], 1224 fabs(coeff[2]-image->page.x-image->columns)); 1225 coeff[0]=MagickMin(coeff[0], 1226 fabs(coeff[3]-image->page.y-image->rows)); 1227 } 1228 /* furthest diagonal if radius == -1 */ 1229 if ( fabs(-1.0-coeff[0]) < MagickEpsilon ) { 1230 double rx,ry; 1231 rx = coeff[2]-image->page.x; 1232 ry = coeff[3]-image->page.y; 1233 coeff[0] = rx*rx+ry*ry; 1234 ry = coeff[3]-image->page.y-image->rows; 1235 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry); 1236 rx = coeff[2]-image->page.x-image->columns; 1237 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry); 1238 ry = coeff[3]-image->page.y; 1239 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry); 1240 coeff[0] = sqrt(coeff[0]); 1241 } 1242 } 1243 /* IF Rmax <= 0 or Rmin < 0 OR Rmax < Rmin, THEN error */ 1244 if ( coeff[0] < MagickEpsilon || coeff[1] < -MagickEpsilon 1245 || (coeff[0]-coeff[1]) < MagickEpsilon ) { 1246 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 1247 "InvalidArgument", "%s : Invalid Radius", 1248 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 1249 coeff=(double *) RelinquishMagickMemory(coeff); 1250 return((double *) NULL); 1251 } 1252 /* converstion ratios */ 1253 if ( *method == PolarDistortion ) { 1254 coeff[6]=(double) image->columns/(coeff[5]-coeff[4]); 1255 coeff[7]=(double) image->rows/(coeff[0]-coeff[1]); 1256 } 1257 else { /* *method == DePolarDistortion */ 1258 coeff[6]=(coeff[5]-coeff[4])/image->columns; 1259 coeff[7]=(coeff[0]-coeff[1])/image->rows; 1260 } 1261 return(coeff); 1262 } 1263 case Cylinder2PlaneDistortion: 1264 case Plane2CylinderDistortion: 1265 { 1266 /* 3D Cylinder to/from a Tangential Plane 1267 1268 Projection between a clinder and flat plain from a point on the 1269 center line of the cylinder. 1270 1271 The two surfaces coincide in 3D space at the given centers of 1272 distortion (perpendicular to projection point) on both images. 1273 1274 Args: FOV_arc_width 1275 Coefficents: FOV(radians), Radius, center_x,y, dest_center_x,y 1276 1277 FOV (Field Of View) the angular field of view of the distortion, 1278 across the width of the image, in degrees. The centers are the 1279 points of least distortion in the input and resulting images. 1280 1281 These centers are however determined later. 1282 1283 Coeff 0 is the FOV angle of view of image width in radians 1284 Coeff 1 is calculated radius of cylinder. 1285 Coeff 2,3 center of distortion of input image 1286 Coefficents 4,5 Center of Distortion of dest (determined later) 1287 */ 1288 if ( arguments[0] < MagickEpsilon || arguments[0] > 160.0 ) { 1289 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 1290 "InvalidArgument", "%s : Invalid FOV Angle", 1291 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 1292 coeff=(double *) RelinquishMagickMemory(coeff); 1293 return((double *) NULL); 1294 } 1295 coeff[0] = DegreesToRadians(arguments[0]); 1296 if ( *method == Cylinder2PlaneDistortion ) 1297 /* image is curved around cylinder, so FOV angle (in radians) 1298 * scales directly to image X coordinate, according to its radius. 1299 */ 1300 coeff[1] = (double) image->columns/coeff[0]; 1301 else 1302 /* radius is distance away from an image with this angular FOV */ 1303 coeff[1] = (double) image->columns / ( 2 * tan(coeff[0]/2) ); 1304 1305 coeff[2] = (double)(image->columns)/2.0+image->page.x; 1306 coeff[3] = (double)(image->rows)/2.0+image->page.y; 1307 coeff[4] = coeff[2]; 1308 coeff[5] = coeff[3]; /* assuming image size is the same */ 1309 return(coeff); 1310 } 1311 case BarrelDistortion: 1312 case BarrelInverseDistortion: 1313 { 1314 /* Barrel Distortion 1315 Rs=(A*Rd^3 + B*Rd^2 + C*Rd + D)*Rd 1316 BarrelInv Distortion 1317 Rs=Rd/(A*Rd^3 + B*Rd^2 + C*Rd + D) 1318 1319 Where Rd is the normalized radius from corner to middle of image 1320 Input Arguments are one of the following forms (number of arguments)... 1321 3: A,B,C 1322 4: A,B,C,D 1323 5: A,B,C X,Y 1324 6: A,B,C,D X,Y 1325 8: Ax,Bx,Cx,Dx Ay,By,Cy,Dy 1326 10: Ax,Bx,Cx,Dx Ay,By,Cy,Dy X,Y 1327 1328 Returns 10 coefficent values, which are de-normalized (pixel scale) 1329 Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Xc, Yc 1330 */ 1331 /* Radius de-normalization scaling factor */ 1332 double 1333 rscale = 2.0/MagickMin((double) image->columns,(double) image->rows); 1334 1335 /* sanity check number of args must = 3,4,5,6,8,10 or error */ 1336 if ( (number_arguments < 3) || (number_arguments == 7) || 1337 (number_arguments == 9) || (number_arguments > 10) ) 1338 { 1339 coeff=(double *) RelinquishMagickMemory(coeff); 1340 (void) ThrowMagickException(exception,GetMagickModule(), 1341 OptionError,"InvalidArgument", "%s : number of arguments", 1342 CommandOptionToMnemonic(MagickDistortOptions, *method) ); 1343 return((double *) NULL); 1344 } 1345 /* A,B,C,D coefficients */ 1346 coeff[0] = arguments[0]; 1347 coeff[1] = arguments[1]; 1348 coeff[2] = arguments[2]; 1349 if ((number_arguments == 3) || (number_arguments == 5) ) 1350 coeff[3] = 1.0 - coeff[0] - coeff[1] - coeff[2]; 1351 else 1352 coeff[3] = arguments[3]; 1353 /* de-normalize the coefficients */ 1354 coeff[0] *= pow(rscale,3.0); 1355 coeff[1] *= rscale*rscale; 1356 coeff[2] *= rscale; 1357 /* Y coefficients: as given OR same as X coefficients */ 1358 if ( number_arguments >= 8 ) { 1359 coeff[4] = arguments[4] * pow(rscale,3.0); 1360 coeff[5] = arguments[5] * rscale*rscale; 1361 coeff[6] = arguments[6] * rscale; 1362 coeff[7] = arguments[7]; 1363 } 1364 else { 1365 coeff[4] = coeff[0]; 1366 coeff[5] = coeff[1]; 1367 coeff[6] = coeff[2]; 1368 coeff[7] = coeff[3]; 1369 } 1370 /* X,Y Center of Distortion (image coodinates) */ 1371 if ( number_arguments == 5 ) { 1372 coeff[8] = arguments[3]; 1373 coeff[9] = arguments[4]; 1374 } 1375 else if ( number_arguments == 6 ) { 1376 coeff[8] = arguments[4]; 1377 coeff[9] = arguments[5]; 1378 } 1379 else if ( number_arguments == 10 ) { 1380 coeff[8] = arguments[8]; 1381 coeff[9] = arguments[9]; 1382 } 1383 else { 1384 /* center of the image provided (image coodinates) */ 1385 coeff[8] = (double)image->columns/2.0 + image->page.x; 1386 coeff[9] = (double)image->rows/2.0 + image->page.y; 1387 } 1388 return(coeff); 1389 } 1390 case ShepardsDistortion: 1391 { 1392 /* Shepards Distortion input arguments are the coefficents! 1393 Just check the number of arguments is valid! 1394 Args: u1,v1, x1,y1, ... 1395 OR : u1,v1, r1,g1,c1, ... 1396 */ 1397 if ( number_arguments%cp_size != 0 || 1398 number_arguments < cp_size ) { 1399 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 1400 "InvalidArgument", "%s : 'requires CP's (4 numbers each)'", 1401 CommandOptionToMnemonic(MagickDistortOptions, *method)); 1402 coeff=(double *) RelinquishMagickMemory(coeff); 1403 return((double *) NULL); 1404 } 1405 /* User defined weighting power for Shepard's Method */ 1406 { const char *artifact=GetImageArtifact(image,"shepards:power"); 1407 if ( artifact != (const char *) NULL ) { 1408 coeff[0]=StringToDouble(artifact,(char **) NULL) / 2.0; 1409 if ( coeff[0] < MagickEpsilon ) { 1410 (void) ThrowMagickException(exception,GetMagickModule(), 1411 OptionError,"InvalidArgument","%s", "-define shepards:power" ); 1412 coeff=(double *) RelinquishMagickMemory(coeff); 1413 return((double *) NULL); 1414 } 1415 } 1416 else 1417 coeff[0]=1.0; /* Default power of 2 (Inverse Squared) */ 1418 } 1419 return(coeff); 1420 } 1421 default: 1422 break; 1423 } 1424 /* you should never reach this point */ 1425 perror("no method handler"); /* just fail assertion */ 1426 return((double *) NULL); 1427 } 1428 1429 /* 1431 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1432 % % 1433 % % 1434 % % 1435 + D i s t o r t R e s i z e I m a g e % 1436 % % 1437 % % 1438 % % 1439 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1440 % 1441 % DistortResizeImage() resize image using the equivalent but slower image 1442 % distortion operator. The filter is applied using a EWA cylindrical 1443 % resampling. But like resize the final image size is limited to whole pixels 1444 % with no effects by virtual-pixels on the result. 1445 % 1446 % Note that images containing a transparency channel will be twice as slow to 1447 % resize as images one without transparency. 1448 % 1449 % The format of the DistortResizeImage method is: 1450 % 1451 % Image *DistortResizeImage(const Image *image,const size_t columns, 1452 % const size_t rows,ExceptionInfo *exception) 1453 % 1454 % A description of each parameter follows: 1455 % 1456 % o image: the image. 1457 % 1458 % o columns: the number of columns in the resized image. 1459 % 1460 % o rows: the number of rows in the resized image. 1461 % 1462 % o exception: return any errors or warnings in this structure. 1463 % 1464 */ 1465 MagickExport Image *DistortResizeImage(const Image *image, 1466 const size_t columns,const size_t rows,ExceptionInfo *exception) 1467 { 1468 #define DistortResizeImageTag "Distort/Image" 1469 1470 Image 1471 *resize_image, 1472 *tmp_image; 1473 1474 RectangleInfo 1475 crop_area; 1476 1477 double 1478 distort_args[12]; 1479 1480 VirtualPixelMethod 1481 vp_save; 1482 1483 /* 1484 Distort resize image. 1485 */ 1486 assert(image != (const Image *) NULL); 1487 assert(image->signature == MagickCoreSignature); 1488 if (image->debug != MagickFalse) 1489 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1490 assert(exception != (ExceptionInfo *) NULL); 1491 assert(exception->signature == MagickCoreSignature); 1492 if ((columns == 0) || (rows == 0)) 1493 return((Image *) NULL); 1494 /* Do not short-circuit this resize if final image size is unchanged */ 1495 1496 (void) ResetMagickMemory(distort_args,0,12*sizeof(double)); 1497 distort_args[4]=(double) image->columns; 1498 distort_args[6]=(double) columns; 1499 distort_args[9]=(double) image->rows; 1500 distort_args[11]=(double) rows; 1501 1502 vp_save=GetImageVirtualPixelMethod(image); 1503 1504 tmp_image=CloneImage(image,0,0,MagickTrue,exception); 1505 if ( tmp_image == (Image *) NULL ) 1506 return((Image *) NULL); 1507 (void) SetImageVirtualPixelMethod(tmp_image,TransparentVirtualPixelMethod, 1508 exception); 1509 1510 if (image->alpha_trait == UndefinedPixelTrait) 1511 { 1512 /* 1513 Image has not transparency channel, so we free to use it 1514 */ 1515 (void) SetImageAlphaChannel(tmp_image,SetAlphaChannel,exception); 1516 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args, 1517 MagickTrue,exception), 1518 1519 tmp_image=DestroyImage(tmp_image); 1520 if ( resize_image == (Image *) NULL ) 1521 return((Image *) NULL); 1522 1523 (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel, 1524 exception); 1525 } 1526 else 1527 { 1528 /* 1529 Image has transparency so handle colors and alpha separatly. 1530 Basically we need to separate Virtual-Pixel alpha in the resized 1531 image, so only the actual original images alpha channel is used. 1532 1533 distort alpha channel separately 1534 */ 1535 Image 1536 *resize_alpha; 1537 1538 (void) SetImageAlphaChannel(tmp_image,ExtractAlphaChannel,exception); 1539 (void) SetImageAlphaChannel(tmp_image,OpaqueAlphaChannel,exception); 1540 resize_alpha=DistortImage(tmp_image,AffineDistortion,12,distort_args, 1541 MagickTrue,exception), 1542 tmp_image=DestroyImage(tmp_image); 1543 if (resize_alpha == (Image *) NULL) 1544 return((Image *) NULL); 1545 1546 /* distort the actual image containing alpha + VP alpha */ 1547 tmp_image=CloneImage(image,0,0,MagickTrue,exception); 1548 if ( tmp_image == (Image *) NULL ) 1549 return((Image *) NULL); 1550 (void) SetImageVirtualPixelMethod(tmp_image,TransparentVirtualPixelMethod, exception); 1551 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args, 1552 MagickTrue,exception), 1553 tmp_image=DestroyImage(tmp_image); 1554 if ( resize_image == (Image *) NULL) 1555 { 1556 resize_alpha=DestroyImage(resize_alpha); 1557 return((Image *) NULL); 1558 } 1559 /* replace resize images alpha with the separally distorted alpha */ 1560 (void) SetImageAlphaChannel(resize_image,OffAlphaChannel,exception); 1561 (void) SetImageAlphaChannel(resize_alpha,OffAlphaChannel,exception); 1562 (void) CompositeImage(resize_image,resize_alpha,CopyAlphaCompositeOp, 1563 MagickTrue,0,0,exception); 1564 resize_alpha=DestroyImage(resize_alpha); 1565 } 1566 (void) SetImageVirtualPixelMethod(resize_image,vp_save,exception); 1567 1568 /* 1569 Clean up the results of the Distortion 1570 */ 1571 crop_area.width=columns; 1572 crop_area.height=rows; 1573 crop_area.x=0; 1574 crop_area.y=0; 1575 1576 tmp_image=resize_image; 1577 resize_image=CropImage(tmp_image,&crop_area,exception); 1578 tmp_image=DestroyImage(tmp_image); 1579 if (resize_image != (Image *) NULL) 1580 { 1581 resize_image->alpha_trait=image->alpha_trait; 1582 resize_image->compose=image->compose; 1583 resize_image->page.width=0; 1584 resize_image->page.height=0; 1585 } 1586 return(resize_image); 1587 } 1588 1589 /* 1591 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1592 % % 1593 % % 1594 % % 1595 % D i s t o r t I m a g e % 1596 % % 1597 % % 1598 % % 1599 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1600 % 1601 % DistortImage() distorts an image using various distortion methods, by 1602 % mapping color lookups of the source image to a new destination image 1603 % usally of the same size as the source image, unless 'bestfit' is set to 1604 % true. 1605 % 1606 % If 'bestfit' is enabled, and distortion allows it, the destination image is 1607 % adjusted to ensure the whole source 'image' will just fit within the final 1608 % destination image, which will be sized and offset accordingly. Also in 1609 % many cases the virtual offset of the source image will be taken into 1610 % account in the mapping. 1611 % 1612 % If the '-verbose' control option has been set print to standard error the 1613 % equicelent '-fx' formula with coefficients for the function, if practical. 1614 % 1615 % The format of the DistortImage() method is: 1616 % 1617 % Image *DistortImage(const Image *image,const DistortMethod method, 1618 % const size_t number_arguments,const double *arguments, 1619 % MagickBooleanType bestfit, ExceptionInfo *exception) 1620 % 1621 % A description of each parameter follows: 1622 % 1623 % o image: the image to be distorted. 1624 % 1625 % o method: the method of image distortion. 1626 % 1627 % ArcDistortion always ignores source image offset, and always 1628 % 'bestfit' the destination image with the top left corner offset 1629 % relative to the polar mapping center. 1630 % 1631 % Affine, Perspective, and Bilinear, do least squares fitting of the 1632 % distrotion when more than the minimum number of control point pairs 1633 % are provided. 1634 % 1635 % Perspective, and Bilinear, fall back to a Affine distortion when less 1636 % than 4 control point pairs are provided. While Affine distortions 1637 % let you use any number of control point pairs, that is Zero pairs is 1638 % a No-Op (viewport only) distortion, one pair is a translation and 1639 % two pairs of control points do a scale-rotate-translate, without any 1640 % shearing. 1641 % 1642 % o number_arguments: the number of arguments given. 1643 % 1644 % o arguments: an array of floating point arguments for this method. 1645 % 1646 % o bestfit: Attempt to 'bestfit' the size of the resulting image. 1647 % This also forces the resulting image to be a 'layered' virtual 1648 % canvas image. Can be overridden using 'distort:viewport' setting. 1649 % 1650 % o exception: return any errors or warnings in this structure 1651 % 1652 % Extra Controls from Image meta-data (artifacts)... 1653 % 1654 % o "verbose" 1655 % Output to stderr alternatives, internal coefficents, and FX 1656 % equivalents for the distortion operation (if feasible). 1657 % This forms an extra check of the distortion method, and allows users 1658 % access to the internal constants IM calculates for the distortion. 1659 % 1660 % o "distort:viewport" 1661 % Directly set the output image canvas area and offest to use for the 1662 % resulting image, rather than use the original images canvas, or a 1663 % calculated 'bestfit' canvas. 1664 % 1665 % o "distort:scale" 1666 % Scale the size of the output canvas by this amount to provide a 1667 % method of Zooming, and for super-sampling the results. 1668 % 1669 % Other settings that can effect results include 1670 % 1671 % o 'interpolate' For source image lookups (scale enlargements) 1672 % 1673 % o 'filter' Set filter to use for area-resampling (scale shrinking). 1674 % Set to 'point' to turn off and use 'interpolate' lookup 1675 % instead 1676 % 1677 */ 1678 MagickExport Image *DistortImage(const Image *image, DistortMethod method, 1679 const size_t number_arguments,const double *arguments, 1680 MagickBooleanType bestfit,ExceptionInfo *exception) 1681 { 1682 #define DistortImageTag "Distort/Image" 1683 1684 double 1685 *coeff, 1686 output_scaling; 1687 1688 Image 1689 *distort_image; 1690 1691 RectangleInfo 1692 geometry; /* geometry of the distorted space viewport */ 1693 1694 MagickBooleanType 1695 viewport_given; 1696 1697 assert(image != (Image *) NULL); 1698 assert(image->signature == MagickCoreSignature); 1699 if (image->debug != MagickFalse) 1700 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1701 assert(exception != (ExceptionInfo *) NULL); 1702 assert(exception->signature == MagickCoreSignature); 1703 1704 /* 1705 Handle Special Compound Distortions 1706 */ 1707 if ( method == ResizeDistortion ) 1708 { 1709 if ( number_arguments != 2 ) 1710 { 1711 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 1712 "InvalidArgument","%s : '%s'","Resize", 1713 "Invalid number of args: 2 only"); 1714 return((Image *) NULL); 1715 } 1716 distort_image=DistortResizeImage(image,(size_t)arguments[0], 1717 (size_t)arguments[1], exception); 1718 return(distort_image); 1719 } 1720 1721 /* 1722 Convert input arguments (usually as control points for reverse mapping) 1723 into mapping coefficients to apply the distortion. 1724 1725 Note that some distortions are mapped to other distortions, 1726 and as such do not require specific code after this point. 1727 */ 1728 coeff = GenerateCoefficients(image, &method, number_arguments, 1729 arguments, 0, exception); 1730 if ( coeff == (double *) NULL ) 1731 return((Image *) NULL); 1732 1733 /* 1734 Determine the size and offset for a 'bestfit' destination. 1735 Usally the four corners of the source image is enough. 1736 */ 1737 1738 /* default output image bounds, when no 'bestfit' is requested */ 1739 geometry.width=image->columns; 1740 geometry.height=image->rows; 1741 geometry.x=0; 1742 geometry.y=0; 1743 1744 if ( method == ArcDistortion ) { 1745 bestfit = MagickTrue; /* always calculate a 'best fit' viewport */ 1746 } 1747 1748 /* Work out the 'best fit', (required for ArcDistortion) */ 1749 if ( bestfit ) { 1750 PointInfo 1751 s,d,min,max; /* source, dest coords --mapping--> min, max coords */ 1752 1753 MagickBooleanType 1754 fix_bounds = MagickTrue; /* enlarge bounds for VP handling */ 1755 1756 s.x=s.y=min.x=max.x=min.y=max.y=0.0; /* keep compiler happy */ 1757 1758 /* defines to figure out the bounds of the distorted image */ 1759 #define InitalBounds(p) \ 1760 { \ 1761 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \ 1762 min.x = max.x = p.x; \ 1763 min.y = max.y = p.y; \ 1764 } 1765 #define ExpandBounds(p) \ 1766 { \ 1767 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \ 1768 min.x = MagickMin(min.x,p.x); \ 1769 max.x = MagickMax(max.x,p.x); \ 1770 min.y = MagickMin(min.y,p.y); \ 1771 max.y = MagickMax(max.y,p.y); \ 1772 } 1773 1774 switch (method) 1775 { 1776 case AffineDistortion: 1777 { double inverse[6]; 1778 InvertAffineCoefficients(coeff, inverse); 1779 s.x = (double) image->page.x; 1780 s.y = (double) image->page.y; 1781 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2]; 1782 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5]; 1783 InitalBounds(d); 1784 s.x = (double) image->page.x+image->columns; 1785 s.y = (double) image->page.y; 1786 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2]; 1787 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5]; 1788 ExpandBounds(d); 1789 s.x = (double) image->page.x; 1790 s.y = (double) image->page.y+image->rows; 1791 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2]; 1792 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5]; 1793 ExpandBounds(d); 1794 s.x = (double) image->page.x+image->columns; 1795 s.y = (double) image->page.y+image->rows; 1796 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2]; 1797 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5]; 1798 ExpandBounds(d); 1799 break; 1800 } 1801 case PerspectiveDistortion: 1802 { double inverse[8], scale; 1803 InvertPerspectiveCoefficients(coeff, inverse); 1804 s.x = (double) image->page.x; 1805 s.y = (double) image->page.y; 1806 scale=inverse[6]*s.x+inverse[7]*s.y+1.0; 1807 scale=PerceptibleReciprocal(scale); 1808 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]); 1809 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]); 1810 InitalBounds(d); 1811 s.x = (double) image->page.x+image->columns; 1812 s.y = (double) image->page.y; 1813 scale=inverse[6]*s.x+inverse[7]*s.y+1.0; 1814 scale=PerceptibleReciprocal(scale); 1815 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]); 1816 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]); 1817 ExpandBounds(d); 1818 s.x = (double) image->page.x; 1819 s.y = (double) image->page.y+image->rows; 1820 scale=inverse[6]*s.x+inverse[7]*s.y+1.0; 1821 scale=PerceptibleReciprocal(scale); 1822 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]); 1823 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]); 1824 ExpandBounds(d); 1825 s.x = (double) image->page.x+image->columns; 1826 s.y = (double) image->page.y+image->rows; 1827 scale=inverse[6]*s.x+inverse[7]*s.y+1.0; 1828 scale=PerceptibleReciprocal(scale); 1829 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]); 1830 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]); 1831 ExpandBounds(d); 1832 break; 1833 } 1834 case ArcDistortion: 1835 { double a, ca, sa; 1836 /* Forward Map Corners */ 1837 a = coeff[0]-coeff[1]/2; ca = cos(a); sa = sin(a); 1838 d.x = coeff[2]*ca; 1839 d.y = coeff[2]*sa; 1840 InitalBounds(d); 1841 d.x = (coeff[2]-coeff[3])*ca; 1842 d.y = (coeff[2]-coeff[3])*sa; 1843 ExpandBounds(d); 1844 a = coeff[0]+coeff[1]/2; ca = cos(a); sa = sin(a); 1845 d.x = coeff[2]*ca; 1846 d.y = coeff[2]*sa; 1847 ExpandBounds(d); 1848 d.x = (coeff[2]-coeff[3])*ca; 1849 d.y = (coeff[2]-coeff[3])*sa; 1850 ExpandBounds(d); 1851 /* Orthogonal points along top of arc */ 1852 for( a=(double) (ceil((double) ((coeff[0]-coeff[1]/2.0)/MagickPI2))*MagickPI2); 1853 a<(coeff[0]+coeff[1]/2.0); a+=MagickPI2 ) { 1854 ca = cos(a); sa = sin(a); 1855 d.x = coeff[2]*ca; 1856 d.y = coeff[2]*sa; 1857 ExpandBounds(d); 1858 } 1859 /* 1860 Convert the angle_to_width and radius_to_height 1861 to appropriate scaling factors, to allow faster processing 1862 in the mapping function. 1863 */ 1864 coeff[1] = (double) (Magick2PI*image->columns/coeff[1]); 1865 coeff[3] = (double)image->rows/coeff[3]; 1866 break; 1867 } 1868 case PolarDistortion: 1869 { 1870 if (number_arguments < 2) 1871 coeff[2] = coeff[3] = 0.0; 1872 min.x = coeff[2]-coeff[0]; 1873 max.x = coeff[2]+coeff[0]; 1874 min.y = coeff[3]-coeff[0]; 1875 max.y = coeff[3]+coeff[0]; 1876 /* should be about 1.0 if Rmin = 0 */ 1877 coeff[7]=(double) geometry.height/(coeff[0]-coeff[1]); 1878 break; 1879 } 1880 case DePolarDistortion: 1881 { 1882 /* direct calculation as it needs to tile correctly 1883 * for reversibility in a DePolar-Polar cycle */ 1884 fix_bounds = MagickFalse; 1885 geometry.x = geometry.y = 0; 1886 geometry.height = (size_t) ceil(coeff[0]-coeff[1]); 1887 geometry.width = (size_t) 1888 ceil((coeff[0]-coeff[1])*(coeff[5]-coeff[4])*0.5); 1889 /* correct scaling factors relative to new size */ 1890 coeff[6]=(coeff[5]-coeff[4])/geometry.width; /* changed width */ 1891 coeff[7]=(coeff[0]-coeff[1])/geometry.height; /* should be about 1.0 */ 1892 break; 1893 } 1894 case Cylinder2PlaneDistortion: 1895 { 1896 /* direct calculation so center of distortion is either a pixel 1897 * center, or pixel edge. This allows for reversibility of the 1898 * distortion */ 1899 geometry.x = geometry.y = 0; 1900 geometry.width = (size_t) ceil( 2.0*coeff[1]*tan(coeff[0]/2.0) ); 1901 geometry.height = (size_t) ceil( 2.0*coeff[3]/cos(coeff[0]/2.0) ); 1902 /* correct center of distortion relative to new size */ 1903 coeff[4] = (double) geometry.width/2.0; 1904 coeff[5] = (double) geometry.height/2.0; 1905 fix_bounds = MagickFalse; 1906 break; 1907 } 1908 case Plane2CylinderDistortion: 1909 { 1910 /* direct calculation center is either pixel center, or pixel edge 1911 * so as to allow reversibility of the image distortion */ 1912 geometry.x = geometry.y = 0; 1913 geometry.width = (size_t) ceil(coeff[0]*coeff[1]); /* FOV * radius */ 1914 geometry.height = (size_t) (2*coeff[3]); /* input image height */ 1915 /* correct center of distortion relative to new size */ 1916 coeff[4] = (double) geometry.width/2.0; 1917 coeff[5] = (double) geometry.height/2.0; 1918 fix_bounds = MagickFalse; 1919 break; 1920 } 1921 case ShepardsDistortion: 1922 case BilinearForwardDistortion: 1923 case BilinearReverseDistortion: 1924 #if 0 1925 case QuadrilateralDistortion: 1926 #endif 1927 case PolynomialDistortion: 1928 case BarrelDistortion: 1929 case BarrelInverseDistortion: 1930 default: 1931 /* no calculated bestfit available for these distortions */ 1932 bestfit = MagickFalse; 1933 fix_bounds = MagickFalse; 1934 break; 1935 } 1936 1937 /* Set the output image geometry to calculated 'bestfit'. 1938 Yes this tends to 'over do' the file image size, ON PURPOSE! 1939 Do not do this for DePolar which needs to be exact for virtual tiling. 1940 */ 1941 if ( fix_bounds ) { 1942 geometry.x = (ssize_t) floor(min.x-0.5); 1943 geometry.y = (ssize_t) floor(min.y-0.5); 1944 geometry.width=(size_t) ceil(max.x-geometry.x+0.5); 1945 geometry.height=(size_t) ceil(max.y-geometry.y+0.5); 1946 } 1947 1948 } /* end bestfit destination image calculations */ 1949 1950 /* The user provided a 'viewport' expert option which may 1951 overrides some parts of the current output image geometry. 1952 This also overrides its default 'bestfit' setting. 1953 */ 1954 { const char *artifact=GetImageArtifact(image,"distort:viewport"); 1955 viewport_given = MagickFalse; 1956 if ( artifact != (const char *) NULL ) { 1957 MagickStatusType flags=ParseAbsoluteGeometry(artifact,&geometry); 1958 if (flags==NoValue) 1959 (void) ThrowMagickException(exception,GetMagickModule(), 1960 OptionWarning,"InvalidSetting","'%s' '%s'", 1961 "distort:viewport",artifact); 1962 else 1963 viewport_given = MagickTrue; 1964 } 1965 } 1966 1967 /* Verbose output */ 1968 if (IsStringTrue(GetImageArtifact(image,"verbose")) != MagickFalse) { 1969 register ssize_t 1970 i; 1971 char image_gen[MagickPathExtent]; 1972 const char *lookup; 1973 1974 /* Set destination image size and virtual offset */ 1975 if ( bestfit || viewport_given ) { 1976 (void) FormatLocaleString(image_gen, MagickPathExtent," -size %.20gx%.20g " 1977 "-page %+.20g%+.20g xc: +insert \\\n",(double) geometry.width, 1978 (double) geometry.height,(double) geometry.x,(double) geometry.y); 1979 lookup="v.p{ xx-v.page.x-.5, yy-v.page.y-.5 }"; 1980 } 1981 else { 1982 image_gen[0] = '\0'; /* no destination to generate */ 1983 lookup = "p{ xx-page.x-.5, yy-page.y-.5 }"; /* simplify lookup */ 1984 } 1985 1986 switch (method) { 1987 case AffineDistortion: 1988 { 1989 double *inverse; 1990 1991 inverse = (double *) AcquireQuantumMemory(6,sizeof(*inverse)); 1992 if (inverse == (double *) NULL) { 1993 coeff = (double *) RelinquishMagickMemory(coeff); 1994 (void) ThrowMagickException(exception,GetMagickModule(), 1995 ResourceLimitError,"MemoryAllocationFailed", 1996 "%s", "DistortImages"); 1997 return((Image *) NULL); 1998 } 1999 InvertAffineCoefficients(coeff, inverse); 2000 CoefficientsToAffineArgs(inverse); 2001 (void) FormatLocaleFile(stderr, "Affine Projection:\n"); 2002 (void) FormatLocaleFile(stderr, " -distort AffineProjection \\\n '"); 2003 for (i=0; i < 5; i++) 2004 (void) FormatLocaleFile(stderr, "%lf,", inverse[i]); 2005 (void) FormatLocaleFile(stderr, "%lf'\n", inverse[5]); 2006 inverse = (double *) RelinquishMagickMemory(inverse); 2007 2008 (void) FormatLocaleFile(stderr, "Affine Distort, FX Equivelent:\n"); 2009 (void) FormatLocaleFile(stderr, "%s", image_gen); 2010 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n"); 2011 (void) FormatLocaleFile(stderr, " xx=%+lf*ii %+lf*jj %+lf;\n", 2012 coeff[0], coeff[1], coeff[2]); 2013 (void) FormatLocaleFile(stderr, " yy=%+lf*ii %+lf*jj %+lf;\n", 2014 coeff[3], coeff[4], coeff[5]); 2015 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup); 2016 2017 break; 2018 } 2019 2020 case PerspectiveDistortion: 2021 { 2022 double *inverse; 2023 2024 inverse = (double *) AcquireQuantumMemory(8,sizeof(*inverse)); 2025 if (inverse == (double *) NULL) { 2026 coeff = (double *) RelinquishMagickMemory(coeff); 2027 (void) ThrowMagickException(exception,GetMagickModule(), 2028 ResourceLimitError,"MemoryAllocationFailed", 2029 "%s", "DistortCoefficients"); 2030 return((Image *) NULL); 2031 } 2032 InvertPerspectiveCoefficients(coeff, inverse); 2033 (void) FormatLocaleFile(stderr, "Perspective Projection:\n"); 2034 (void) FormatLocaleFile(stderr, " -distort PerspectiveProjection \\\n '"); 2035 for (i=0; i<4; i++) 2036 (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]); 2037 (void) FormatLocaleFile(stderr, "\n "); 2038 for (; i<7; i++) 2039 (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]); 2040 (void) FormatLocaleFile(stderr, "%lf'\n", inverse[7]); 2041 inverse = (double *) RelinquishMagickMemory(inverse); 2042 2043 (void) FormatLocaleFile(stderr, "Perspective Distort, FX Equivelent:\n"); 2044 (void) FormatLocaleFile(stderr, "%s", image_gen); 2045 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n"); 2046 (void) FormatLocaleFile(stderr, " rr=%+lf*ii %+lf*jj + 1;\n", 2047 coeff[6], coeff[7]); 2048 (void) FormatLocaleFile(stderr, " xx=(%+lf*ii %+lf*jj %+lf)/rr;\n", 2049 coeff[0], coeff[1], coeff[2]); 2050 (void) FormatLocaleFile(stderr, " yy=(%+lf*ii %+lf*jj %+lf)/rr;\n", 2051 coeff[3], coeff[4], coeff[5]); 2052 (void) FormatLocaleFile(stderr, " rr%s0 ? %s : blue' \\\n", 2053 coeff[8] < 0 ? "<" : ">", lookup); 2054 break; 2055 } 2056 2057 case BilinearForwardDistortion: 2058 (void) FormatLocaleFile(stderr, "BilinearForward Mapping Equations:\n"); 2059 (void) FormatLocaleFile(stderr, "%s", image_gen); 2060 (void) FormatLocaleFile(stderr, " i = %+lf*x %+lf*y %+lf*x*y %+lf;\n", 2061 coeff[0], coeff[1], coeff[2], coeff[3]); 2062 (void) FormatLocaleFile(stderr, " j = %+lf*x %+lf*y %+lf*x*y %+lf;\n", 2063 coeff[4], coeff[5], coeff[6], coeff[7]); 2064 #if 0 2065 /* for debugging */ 2066 (void) FormatLocaleFile(stderr, " c8 = %+lf c9 = 2*a = %+lf;\n", 2067 coeff[8], coeff[9]); 2068 #endif 2069 (void) FormatLocaleFile(stderr, "BilinearForward Distort, FX Equivelent:\n"); 2070 (void) FormatLocaleFile(stderr, "%s", image_gen); 2071 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n", 2072 0.5-coeff[3], 0.5-coeff[7]); 2073 (void) FormatLocaleFile(stderr, " bb=%lf*ii %+lf*jj %+lf;\n", 2074 coeff[6], -coeff[2], coeff[8]); 2075 /* Handle Special degenerate (non-quadratic) or trapezoidal case */ 2076 if ( coeff[9] != 0 ) { 2077 (void) FormatLocaleFile(stderr, " rt=bb*bb %+lf*(%lf*ii%+lf*jj);\n", 2078 -2*coeff[9], coeff[4], -coeff[0]); 2079 (void) FormatLocaleFile(stderr, " yy=( -bb + sqrt(rt) ) / %lf;\n", 2080 coeff[9]); 2081 } else 2082 (void) FormatLocaleFile(stderr, " yy=(%lf*ii%+lf*jj)/bb;\n", 2083 -coeff[4], coeff[0]); 2084 (void) FormatLocaleFile(stderr, " xx=(ii %+lf*yy)/(%lf %+lf*yy);\n", 2085 -coeff[1], coeff[0], coeff[2]); 2086 if ( coeff[9] != 0 ) 2087 (void) FormatLocaleFile(stderr, " (rt < 0 ) ? red : %s'\n", lookup); 2088 else 2089 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup); 2090 break; 2091 2092 case BilinearReverseDistortion: 2093 #if 0 2094 (void) FormatLocaleFile(stderr, "Polynomial Projection Distort:\n"); 2095 (void) FormatLocaleFile(stderr, " -distort PolynomialProjection \\\n"); 2096 (void) FormatLocaleFile(stderr, " '1.5, %lf, %lf, %lf, %lf,\n", 2097 coeff[3], coeff[0], coeff[1], coeff[2]); 2098 (void) FormatLocaleFile(stderr, " %lf, %lf, %lf, %lf'\n", 2099 coeff[7], coeff[4], coeff[5], coeff[6]); 2100 #endif 2101 (void) FormatLocaleFile(stderr, "BilinearReverse Distort, FX Equivelent:\n"); 2102 (void) FormatLocaleFile(stderr, "%s", image_gen); 2103 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n"); 2104 (void) FormatLocaleFile(stderr, " xx=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n", 2105 coeff[0], coeff[1], coeff[2], coeff[3]); 2106 (void) FormatLocaleFile(stderr, " yy=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n", 2107 coeff[4], coeff[5], coeff[6], coeff[7]); 2108 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup); 2109 break; 2110 2111 case PolynomialDistortion: 2112 { 2113 size_t nterms = (size_t) coeff[1]; 2114 (void) FormatLocaleFile(stderr, "Polynomial (order %lg, terms %lu), FX Equivelent\n", 2115 coeff[0],(unsigned long) nterms); 2116 (void) FormatLocaleFile(stderr, "%s", image_gen); 2117 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n"); 2118 (void) FormatLocaleFile(stderr, " xx ="); 2119 for (i=0; i<(ssize_t) nterms; i++) { 2120 if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n "); 2121 (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i], 2122 poly_basis_str(i)); 2123 } 2124 (void) FormatLocaleFile(stderr, ";\n yy ="); 2125 for (i=0; i<(ssize_t) nterms; i++) { 2126 if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n "); 2127 (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i+nterms], 2128 poly_basis_str(i)); 2129 } 2130 (void) FormatLocaleFile(stderr, ";\n %s' \\\n", lookup); 2131 break; 2132 } 2133 case ArcDistortion: 2134 { 2135 (void) FormatLocaleFile(stderr, "Arc Distort, Internal Coefficients:\n"); 2136 for ( i=0; i<5; i++ ) 2137 (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]); 2138 (void) FormatLocaleFile(stderr, "Arc Distort, FX Equivelent:\n"); 2139 (void) FormatLocaleFile(stderr, "%s", image_gen); 2140 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x; jj=j+page.y;\n"); 2141 (void) FormatLocaleFile(stderr, " xx=(atan2(jj,ii)%+lf)/(2*pi);\n", 2142 -coeff[0]); 2143 (void) FormatLocaleFile(stderr, " xx=xx-round(xx);\n"); 2144 (void) FormatLocaleFile(stderr, " xx=xx*%lf %+lf;\n", 2145 coeff[1], coeff[4]); 2146 (void) FormatLocaleFile(stderr, " yy=(%lf - hypot(ii,jj)) * %lf;\n", 2147 coeff[2], coeff[3]); 2148 (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n"); 2149 break; 2150 } 2151 case PolarDistortion: 2152 { 2153 (void) FormatLocaleFile(stderr, "Polar Distort, Internal Coefficents\n"); 2154 for ( i=0; i<8; i++ ) 2155 (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]); 2156 (void) FormatLocaleFile(stderr, "Polar Distort, FX Equivelent:\n"); 2157 (void) FormatLocaleFile(stderr, "%s", image_gen); 2158 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n", 2159 -coeff[2], -coeff[3]); 2160 (void) FormatLocaleFile(stderr, " xx=(atan2(ii,jj)%+lf)/(2*pi);\n", 2161 -(coeff[4]+coeff[5])/2 ); 2162 (void) FormatLocaleFile(stderr, " xx=xx-round(xx);\n"); 2163 (void) FormatLocaleFile(stderr, " xx=xx*2*pi*%lf + v.w/2;\n", 2164 coeff[6] ); 2165 (void) FormatLocaleFile(stderr, " yy=(hypot(ii,jj)%+lf)*%lf;\n", 2166 -coeff[1], coeff[7] ); 2167 (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n"); 2168 break; 2169 } 2170 case DePolarDistortion: 2171 { 2172 (void) FormatLocaleFile(stderr, "DePolar Distort, Internal Coefficents\n"); 2173 for ( i=0; i<8; i++ ) 2174 (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]); 2175 (void) FormatLocaleFile(stderr, "DePolar Distort, FX Equivelent:\n"); 2176 (void) FormatLocaleFile(stderr, "%s", image_gen); 2177 (void) FormatLocaleFile(stderr, " -fx 'aa=(i+.5)*%lf %+lf;\n", coeff[6], +coeff[4] ); 2178 (void) FormatLocaleFile(stderr, " rr=(j+.5)*%lf %+lf;\n", coeff[7], +coeff[1] ); 2179 (void) FormatLocaleFile(stderr, " xx=rr*sin(aa) %+lf;\n", coeff[2] ); 2180 (void) FormatLocaleFile(stderr, " yy=rr*cos(aa) %+lf;\n", coeff[3] ); 2181 (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n"); 2182 break; 2183 } 2184 case Cylinder2PlaneDistortion: 2185 { 2186 (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, Internal Coefficents\n"); 2187 (void) FormatLocaleFile(stderr, " cylinder_radius = %+lf\n", coeff[1]); 2188 (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, FX Equivelent:\n"); 2189 (void) FormatLocaleFile(stderr, "%s", image_gen); 2190 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n", 2191 -coeff[4], -coeff[5]); 2192 (void) FormatLocaleFile(stderr, " aa=atan(ii/%+lf);\n", coeff[1] ); 2193 (void) FormatLocaleFile(stderr, " xx=%lf*aa%+lf;\n", 2194 coeff[1], coeff[2] ); 2195 (void) FormatLocaleFile(stderr, " yy=jj*cos(aa)%+lf;\n", coeff[3] ); 2196 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup); 2197 break; 2198 } 2199 case Plane2CylinderDistortion: 2200 { 2201 (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, Internal Coefficents\n"); 2202 (void) FormatLocaleFile(stderr, " cylinder_radius = %+lf\n", coeff[1]); 2203 (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, FX Equivelent:\n"); 2204 (void) FormatLocaleFile(stderr, "%s", image_gen); 2205 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n", 2206 -coeff[4], -coeff[5]); 2207 (void) FormatLocaleFile(stderr, " ii=ii/%+lf;\n", coeff[1] ); 2208 (void) FormatLocaleFile(stderr, " xx=%lf*tan(ii)%+lf;\n", 2209 coeff[1], coeff[2] ); 2210 (void) FormatLocaleFile(stderr, " yy=jj/cos(ii)%+lf;\n", 2211 coeff[3] ); 2212 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup); 2213 break; 2214 } 2215 case BarrelDistortion: 2216 case BarrelInverseDistortion: 2217 { double xc,yc; 2218 /* NOTE: This does the barrel roll in pixel coords not image coords 2219 ** The internal distortion must do it in image coordinates, 2220 ** so that is what the center coeff (8,9) is given in. 2221 */ 2222 xc = ((double)image->columns-1.0)/2.0 + image->page.x; 2223 yc = ((double)image->rows-1.0)/2.0 + image->page.y; 2224 (void) FormatLocaleFile(stderr, "Barrel%s Distort, FX Equivelent:\n", 2225 method == BarrelDistortion ? "" : "Inv"); 2226 (void) FormatLocaleFile(stderr, "%s", image_gen); 2227 if ( fabs(coeff[8]-xc-0.5) < 0.1 && fabs(coeff[9]-yc-0.5) < 0.1 ) 2228 (void) FormatLocaleFile(stderr, " -fx 'xc=(w-1)/2; yc=(h-1)/2;\n"); 2229 else 2230 (void) FormatLocaleFile(stderr, " -fx 'xc=%lf; yc=%lf;\n", 2231 coeff[8]-0.5, coeff[9]-0.5); 2232 (void) FormatLocaleFile(stderr, 2233 " ii=i-xc; jj=j-yc; rr=hypot(ii,jj);\n"); 2234 (void) FormatLocaleFile(stderr, " ii=ii%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n", 2235 method == BarrelDistortion ? "*" : "/", 2236 coeff[0],coeff[1],coeff[2],coeff[3]); 2237 (void) FormatLocaleFile(stderr, " jj=jj%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n", 2238 method == BarrelDistortion ? "*" : "/", 2239 coeff[4],coeff[5],coeff[6],coeff[7]); 2240 (void) FormatLocaleFile(stderr, " v.p{fx*ii+xc,fy*jj+yc}' \\\n"); 2241 } 2242 default: 2243 break; 2244 } 2245 } 2246 2247 /* The user provided a 'scale' expert option will scale the 2248 output image size, by the factor given allowing for super-sampling 2249 of the distorted image space. Any scaling factors must naturally 2250 be halved as a result. 2251 */ 2252 { const char *artifact; 2253 artifact=GetImageArtifact(image,"distort:scale"); 2254 output_scaling = 1.0; 2255 if (artifact != (const char *) NULL) { 2256 output_scaling = fabs(StringToDouble(artifact,(char **) NULL)); 2257 geometry.width=(size_t) (output_scaling*geometry.width+0.5); 2258 geometry.height=(size_t) (output_scaling*geometry.height+0.5); 2259 geometry.x=(ssize_t) (output_scaling*geometry.x+0.5); 2260 geometry.y=(ssize_t) (output_scaling*geometry.y+0.5); 2261 if ( output_scaling < 0.1 ) { 2262 coeff = (double *) RelinquishMagickMemory(coeff); 2263 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 2264 "InvalidArgument","%s", "-set option:distort:scale" ); 2265 return((Image *) NULL); 2266 } 2267 output_scaling = 1/output_scaling; 2268 } 2269 } 2270 #define ScaleFilter(F,A,B,C,D) \ 2271 ScaleResampleFilter( (F), \ 2272 output_scaling*(A), output_scaling*(B), \ 2273 output_scaling*(C), output_scaling*(D) ) 2274 2275 /* 2276 Initialize the distort image attributes. 2277 */ 2278 distort_image=CloneImage(image,geometry.width,geometry.height,MagickTrue, 2279 exception); 2280 if (distort_image == (Image *) NULL) 2281 return((Image *) NULL); 2282 /* if image is ColorMapped - change it to DirectClass */ 2283 if (SetImageStorageClass(distort_image,DirectClass,exception) == MagickFalse) 2284 { 2285 distort_image=DestroyImage(distort_image); 2286 return((Image *) NULL); 2287 } 2288 if ((IsPixelInfoGray(&distort_image->background_color) == MagickFalse) && 2289 (IsGrayColorspace(distort_image->colorspace) != MagickFalse)) 2290 (void) SetImageColorspace(distort_image,sRGBColorspace,exception); 2291 if (distort_image->background_color.alpha_trait != UndefinedPixelTrait) 2292 distort_image->alpha_trait=BlendPixelTrait; 2293 distort_image->page.x=geometry.x; 2294 distort_image->page.y=geometry.y; 2295 2296 { /* ----- MAIN CODE ----- 2297 Sample the source image to each pixel in the distort image. 2298 */ 2299 CacheView 2300 *distort_view; 2301 2302 MagickBooleanType 2303 status; 2304 2305 MagickOffsetType 2306 progress; 2307 2308 PixelInfo 2309 zero; 2310 2311 ResampleFilter 2312 **magick_restrict resample_filter; 2313 2314 ssize_t 2315 j; 2316 2317 status=MagickTrue; 2318 progress=0; 2319 GetPixelInfo(distort_image,&zero); 2320 resample_filter=AcquireResampleFilterThreadSet(image, 2321 UndefinedVirtualPixelMethod,MagickFalse,exception); 2322 distort_view=AcquireAuthenticCacheView(distort_image,exception); 2323 #if defined(MAGICKCORE_OPENMP_SUPPORT) 2324 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 2325 magick_threads(image,distort_image,distort_image->rows,1) 2326 #endif 2327 for (j=0; j < (ssize_t) distort_image->rows; j++) 2328 { 2329 const int 2330 id = GetOpenMPThreadId(); 2331 2332 double 2333 validity; /* how mathematically valid is this the mapping */ 2334 2335 MagickBooleanType 2336 sync; 2337 2338 PixelInfo 2339 pixel, /* pixel color to assign to distorted image */ 2340 invalid; /* the color to assign when distort result is invalid */ 2341 2342 PointInfo 2343 d, 2344 s; /* transform destination image x,y to source image x,y */ 2345 2346 register ssize_t 2347 i; 2348 2349 register Quantum 2350 *magick_restrict q; 2351 2352 q=QueueCacheViewAuthenticPixels(distort_view,0,j,distort_image->columns,1, 2353 exception); 2354 if (q == (Quantum *) NULL) 2355 { 2356 status=MagickFalse; 2357 continue; 2358 } 2359 pixel=zero; 2360 2361 /* Define constant scaling vectors for Affine Distortions 2362 Other methods are either variable, or use interpolated lookup 2363 */ 2364 switch (method) 2365 { 2366 case AffineDistortion: 2367 ScaleFilter( resample_filter[id], 2368 coeff[0], coeff[1], 2369 coeff[3], coeff[4] ); 2370 break; 2371 default: 2372 break; 2373 } 2374 2375 /* Initialize default pixel validity 2376 * negative: pixel is invalid output 'alpha_color' 2377 * 0.0 to 1.0: antialiased, mix with resample output 2378 * 1.0 or greater: use resampled output. 2379 */ 2380 validity = 1.0; 2381 2382 ConformPixelInfo(distort_image,&distort_image->alpha_color,&invalid, 2383 exception); 2384 for (i=0; i < (ssize_t) distort_image->columns; i++) 2385 { 2386 /* map pixel coordinate to distortion space coordinate */ 2387 d.x = (double) (geometry.x+i+0.5)*output_scaling; 2388 d.y = (double) (geometry.y+j+0.5)*output_scaling; 2389 s = d; /* default is a no-op mapping */ 2390 switch (method) 2391 { 2392 case AffineDistortion: 2393 { 2394 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2]; 2395 s.y=coeff[3]*d.x+coeff[4]*d.y+coeff[5]; 2396 /* Affine partial derivitives are constant -- set above */ 2397 break; 2398 } 2399 case PerspectiveDistortion: 2400 { 2401 double 2402 p,q,r,abs_r,abs_c6,abs_c7,scale; 2403 /* perspective is a ratio of affines */ 2404 p=coeff[0]*d.x+coeff[1]*d.y+coeff[2]; 2405 q=coeff[3]*d.x+coeff[4]*d.y+coeff[5]; 2406 r=coeff[6]*d.x+coeff[7]*d.y+1.0; 2407 /* Pixel Validity -- is it a 'sky' or 'ground' pixel */ 2408 validity = (r*coeff[8] < 0.0) ? 0.0 : 1.0; 2409 /* Determine horizon anti-alias blending */ 2410 abs_r = fabs(r)*2; 2411 abs_c6 = fabs(coeff[6]); 2412 abs_c7 = fabs(coeff[7]); 2413 if ( abs_c6 > abs_c7 ) { 2414 if ( abs_r < abs_c6*output_scaling ) 2415 validity = 0.5 - coeff[8]*r/(coeff[6]*output_scaling); 2416 } 2417 else if ( abs_r < abs_c7*output_scaling ) 2418 validity = 0.5 - coeff[8]*r/(coeff[7]*output_scaling); 2419 /* Perspective Sampling Point (if valid) */ 2420 if ( validity > 0.0 ) { 2421 /* divide by r affine, for perspective scaling */ 2422 scale = 1.0/r; 2423 s.x = p*scale; 2424 s.y = q*scale; 2425 /* Perspective Partial Derivatives or Scaling Vectors */ 2426 scale *= scale; 2427 ScaleFilter( resample_filter[id], 2428 (r*coeff[0] - p*coeff[6])*scale, 2429 (r*coeff[1] - p*coeff[7])*scale, 2430 (r*coeff[3] - q*coeff[6])*scale, 2431 (r*coeff[4] - q*coeff[7])*scale ); 2432 } 2433 break; 2434 } 2435 case BilinearReverseDistortion: 2436 { 2437 /* Reversed Mapped is just a simple polynomial */ 2438 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2]*d.x*d.y+coeff[3]; 2439 s.y=coeff[4]*d.x+coeff[5]*d.y 2440 +coeff[6]*d.x*d.y+coeff[7]; 2441 /* Bilinear partial derivitives of scaling vectors */ 2442 ScaleFilter( resample_filter[id], 2443 coeff[0] + coeff[2]*d.y, 2444 coeff[1] + coeff[2]*d.x, 2445 coeff[4] + coeff[6]*d.y, 2446 coeff[5] + coeff[6]*d.x ); 2447 break; 2448 } 2449 case BilinearForwardDistortion: 2450 { 2451 /* Forward mapped needs reversed polynomial equations 2452 * which unfortunatally requires a square root! */ 2453 double b,c; 2454 d.x -= coeff[3]; d.y -= coeff[7]; 2455 b = coeff[6]*d.x - coeff[2]*d.y + coeff[8]; 2456 c = coeff[4]*d.x - coeff[0]*d.y; 2457 2458 validity = 1.0; 2459 /* Handle Special degenerate (non-quadratic) case 2460 * Currently without horizon anti-alising */ 2461 if ( fabs(coeff[9]) < MagickEpsilon ) 2462 s.y = -c/b; 2463 else { 2464 c = b*b - 2*coeff[9]*c; 2465 if ( c < 0.0 ) 2466 validity = 0.0; 2467 else 2468 s.y = ( -b + sqrt(c) )/coeff[9]; 2469 } 2470 if ( validity > 0.0 ) 2471 s.x = ( d.x - coeff[1]*s.y) / ( coeff[0] + coeff[2]*s.y ); 2472 2473 /* NOTE: the sign of the square root should be -ve for parts 2474 where the source image becomes 'flipped' or 'mirrored'. 2475 FUTURE: Horizon handling 2476 FUTURE: Scaling factors or Deritives (how?) 2477 */ 2478 break; 2479 } 2480 #if 0 2481 case BilinearDistortion: 2482 /* Bilinear mapping of any Quadrilateral to any Quadrilateral */ 2483 /* UNDER DEVELOPMENT */ 2484 break; 2485 #endif 2486 case PolynomialDistortion: 2487 { 2488 /* multi-ordered polynomial */ 2489 register ssize_t 2490 k; 2491 2492 ssize_t 2493 nterms=(ssize_t)coeff[1]; 2494 2495 PointInfo 2496 du,dv; /* the du,dv vectors from unit dx,dy -- derivatives */ 2497 2498 s.x=s.y=du.x=du.y=dv.x=dv.y=0.0; 2499 for(k=0; k < nterms; k++) { 2500 s.x += poly_basis_fn(k,d.x,d.y)*coeff[2+k]; 2501 du.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k]; 2502 du.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k]; 2503 s.y += poly_basis_fn(k,d.x,d.y)*coeff[2+k+nterms]; 2504 dv.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k+nterms]; 2505 dv.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k+nterms]; 2506 } 2507 ScaleFilter( resample_filter[id], du.x,du.y,dv.x,dv.y ); 2508 break; 2509 } 2510 case ArcDistortion: 2511 { 2512 /* what is the angle and radius in the destination image */ 2513 s.x = (double) ((atan2(d.y,d.x) - coeff[0])/Magick2PI); 2514 s.x -= MagickRound(s.x); /* angle */ 2515 s.y = hypot(d.x,d.y); /* radius */ 2516 2517 /* Arc Distortion Partial Scaling Vectors 2518 Are derived by mapping the perpendicular unit vectors 2519 dR and dA*R*2PI rather than trying to map dx and dy 2520 The results is a very simple orthogonal aligned ellipse. 2521 */ 2522 if ( s.y > MagickEpsilon ) 2523 ScaleFilter( resample_filter[id], 2524 (double) (coeff[1]/(Magick2PI*s.y)), 0, 0, coeff[3] ); 2525 else 2526 ScaleFilter( resample_filter[id], 2527 distort_image->columns*2, 0, 0, coeff[3] ); 2528 2529 /* now scale the angle and radius for source image lookup point */ 2530 s.x = s.x*coeff[1] + coeff[4] + image->page.x +0.5; 2531 s.y = (coeff[2] - s.y) * coeff[3] + image->page.y; 2532 break; 2533 } 2534 case PolarDistortion: 2535 { /* 2D Cartesain to Polar View */ 2536 d.x -= coeff[2]; 2537 d.y -= coeff[3]; 2538 s.x = atan2(d.x,d.y) - (coeff[4]+coeff[5])/2; 2539 s.x /= Magick2PI; 2540 s.x -= MagickRound(s.x); 2541 s.x *= Magick2PI; /* angle - relative to centerline */ 2542 s.y = hypot(d.x,d.y); /* radius */ 2543 2544 /* Polar Scaling vectors are based on mapping dR and dA vectors 2545 This results in very simple orthogonal scaling vectors 2546 */ 2547 if ( s.y > MagickEpsilon ) 2548 ScaleFilter( resample_filter[id], 2549 (double) (coeff[6]/(Magick2PI*s.y)), 0, 0, coeff[7] ); 2550 else 2551 ScaleFilter( resample_filter[id], 2552 distort_image->columns*2, 0, 0, coeff[7] ); 2553 2554 /* now finish mapping radius/angle to source x,y coords */ 2555 s.x = s.x*coeff[6] + (double)image->columns/2.0 + image->page.x; 2556 s.y = (s.y-coeff[1])*coeff[7] + image->page.y; 2557 break; 2558 } 2559 case DePolarDistortion: 2560 { /* @D Polar to Carteasain */ 2561 /* ignore all destination virtual offsets */ 2562 d.x = ((double)i+0.5)*output_scaling*coeff[6]+coeff[4]; 2563 d.y = ((double)j+0.5)*output_scaling*coeff[7]+coeff[1]; 2564 s.x = d.y*sin(d.x) + coeff[2]; 2565 s.y = d.y*cos(d.x) + coeff[3]; 2566 /* derivatives are usless - better to use SuperSampling */ 2567 break; 2568 } 2569 case Cylinder2PlaneDistortion: 2570 { /* 3D Cylinder to Tangential Plane */ 2571 double ax, cx; 2572 /* relative to center of distortion */ 2573 d.x -= coeff[4]; d.y -= coeff[5]; 2574 d.x /= coeff[1]; /* x' = x/r */ 2575 ax=atan(d.x); /* aa = atan(x/r) = u/r */ 2576 cx=cos(ax); /* cx = cos(atan(x/r)) = 1/sqrt(x^2+u^2) */ 2577 s.x = coeff[1]*ax; /* u = r*atan(x/r) */ 2578 s.y = d.y*cx; /* v = y*cos(u/r) */ 2579 /* derivatives... (see personnal notes) */ 2580 ScaleFilter( resample_filter[id], 2581 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y ); 2582 #if 0 2583 if ( i == 0 && j == 0 ) { 2584 fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y); 2585 fprintf(stderr, "phi = %lf\n", (double)(ax * 180.0/MagickPI) ); 2586 fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n", 2587 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y ); 2588 fflush(stderr); } 2589 #endif 2590 /* add center of distortion in source */ 2591 s.x += coeff[2]; s.y += coeff[3]; 2592 break; 2593 } 2594 case Plane2CylinderDistortion: 2595 { /* 3D Cylinder to Tangential Plane */ 2596 /* relative to center of distortion */ 2597 d.x -= coeff[4]; d.y -= coeff[5]; 2598 2599 /* is pixel valid - horizon of a infinite Virtual-Pixel Plane 2600 * (see Anthony Thyssen's personal note) */ 2601 validity = (double) (coeff[1]*MagickPI2 - fabs(d.x))/output_scaling + 0.5; 2602 2603 if ( validity > 0.0 ) { 2604 double cx,tx; 2605 d.x /= coeff[1]; /* x'= x/r */ 2606 cx = 1/cos(d.x); /* cx = 1/cos(x/r) */ 2607 tx = tan(d.x); /* tx = tan(x/r) */ 2608 s.x = coeff[1]*tx; /* u = r * tan(x/r) */ 2609 s.y = d.y*cx; /* v = y / cos(x/r) */ 2610 /* derivatives... (see Anthony Thyssen's personal notes) */ 2611 ScaleFilter( resample_filter[id], 2612 cx*cx, 0.0, s.y*cx/coeff[1], cx ); 2613 #if 0 2614 /*if ( i == 0 && j == 0 )*/ 2615 if ( d.x == 0.5 && d.y == 0.5 ) { 2616 fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y); 2617 fprintf(stderr, "radius = %lf phi = %lf validity = %lf\n", 2618 coeff[1], (double)(d.x * 180.0/MagickPI), validity ); 2619 fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n", 2620 cx*cx, 0.0, s.y*cx/coeff[1], cx); 2621 fflush(stderr); } 2622 #endif 2623 } 2624 /* add center of distortion in source */ 2625 s.x += coeff[2]; s.y += coeff[3]; 2626 break; 2627 } 2628 case BarrelDistortion: 2629 case BarrelInverseDistortion: 2630 { /* Lens Barrel Distionion Correction */ 2631 double r,fx,fy,gx,gy; 2632 /* Radial Polynomial Distortion (de-normalized) */ 2633 d.x -= coeff[8]; 2634 d.y -= coeff[9]; 2635 r = sqrt(d.x*d.x+d.y*d.y); 2636 if ( r > MagickEpsilon ) { 2637 fx = ((coeff[0]*r + coeff[1])*r + coeff[2])*r + coeff[3]; 2638 fy = ((coeff[4]*r + coeff[5])*r + coeff[6])*r + coeff[7]; 2639 gx = ((3*coeff[0]*r + 2*coeff[1])*r + coeff[2])/r; 2640 gy = ((3*coeff[4]*r + 2*coeff[5])*r + coeff[6])/r; 2641 /* adjust functions and scaling for 'inverse' form */ 2642 if ( method == BarrelInverseDistortion ) { 2643 fx = 1/fx; fy = 1/fy; 2644 gx *= -fx*fx; gy *= -fy*fy; 2645 } 2646 /* Set the source pixel to lookup and EWA derivative vectors */ 2647 s.x = d.x*fx + coeff[8]; 2648 s.y = d.y*fy + coeff[9]; 2649 ScaleFilter( resample_filter[id], 2650 gx*d.x*d.x + fx, gx*d.x*d.y, 2651 gy*d.x*d.y, gy*d.y*d.y + fy ); 2652 } 2653 else { 2654 /* Special handling to avoid divide by zero when r==0 2655 ** 2656 ** The source and destination pixels match in this case 2657 ** which was set at the top of the loop using s = d; 2658 ** otherwise... s.x=coeff[8]; s.y=coeff[9]; 2659 */ 2660 if ( method == BarrelDistortion ) 2661 ScaleFilter( resample_filter[id], 2662 coeff[3], 0, 0, coeff[7] ); 2663 else /* method == BarrelInverseDistortion */ 2664 /* FUTURE, trap for D==0 causing division by zero */ 2665 ScaleFilter( resample_filter[id], 2666 1.0/coeff[3], 0, 0, 1.0/coeff[7] ); 2667 } 2668 break; 2669 } 2670 case ShepardsDistortion: 2671 { /* Shepards Method, or Inverse Weighted Distance for 2672 displacement around the destination image control points 2673 The input arguments are the coefficents to the function. 2674 This is more of a 'displacement' function rather than an 2675 absolute distortion function. 2676 2677 Note: We can not determine derivatives using shepards method 2678 so only a point sample interpolatation can be used. 2679 */ 2680 size_t 2681 i; 2682 double 2683 denominator; 2684 2685 denominator = s.x = s.y = 0; 2686 for(i=0; i<number_arguments; i+=4) { 2687 double weight = 2688 ((double)d.x-arguments[i+2])*((double)d.x-arguments[i+2]) 2689 + ((double)d.y-arguments[i+3])*((double)d.y-arguments[i+3]); 2690 weight = pow(weight,coeff[0]); /* shepards power factor */ 2691 weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight; 2692 2693 s.x += (arguments[ i ]-arguments[i+2])*weight; 2694 s.y += (arguments[i+1]-arguments[i+3])*weight; 2695 denominator += weight; 2696 } 2697 s.x /= denominator; 2698 s.y /= denominator; 2699 s.x += d.x; /* make it as relative displacement */ 2700 s.y += d.y; 2701 break; 2702 } 2703 default: 2704 break; /* use the default no-op given above */ 2705 } 2706 /* map virtual canvas location back to real image coordinate */ 2707 if ( bestfit && method != ArcDistortion ) { 2708 s.x -= image->page.x; 2709 s.y -= image->page.y; 2710 } 2711 s.x -= 0.5; 2712 s.y -= 0.5; 2713 2714 if ( validity <= 0.0 ) { 2715 /* result of distortion is an invalid pixel - don't resample */ 2716 SetPixelViaPixelInfo(distort_image,&invalid,q); 2717 } 2718 else { 2719 /* resample the source image to find its correct color */ 2720 (void) ResamplePixelColor(resample_filter[id],s.x,s.y,&pixel, 2721 exception); 2722 /* if validity between 0.0 and 1.0 mix result with invalid pixel */ 2723 if ( validity < 1.0 ) { 2724 /* Do a blend of sample color and invalid pixel */ 2725 /* should this be a 'Blend', or an 'Over' compose */ 2726 CompositePixelInfoBlend(&pixel,validity,&invalid,(1.0-validity), 2727 &pixel); 2728 } 2729 SetPixelViaPixelInfo(distort_image,&pixel,q); 2730 } 2731 q+=GetPixelChannels(distort_image); 2732 } 2733 sync=SyncCacheViewAuthenticPixels(distort_view,exception); 2734 if (sync == MagickFalse) 2735 status=MagickFalse; 2736 if (image->progress_monitor != (MagickProgressMonitor) NULL) 2737 { 2738 MagickBooleanType 2739 proceed; 2740 2741 #if defined(MAGICKCORE_OPENMP_SUPPORT) 2742 #pragma omp critical (MagickCore_DistortImage) 2743 #endif 2744 proceed=SetImageProgress(image,DistortImageTag,progress++, 2745 image->rows); 2746 if (proceed == MagickFalse) 2747 status=MagickFalse; 2748 } 2749 } 2750 distort_view=DestroyCacheView(distort_view); 2751 resample_filter=DestroyResampleFilterThreadSet(resample_filter); 2752 2753 if (status == MagickFalse) 2754 distort_image=DestroyImage(distort_image); 2755 } 2756 2757 /* Arc does not return an offset unless 'bestfit' is in effect 2758 And the user has not provided an overriding 'viewport'. 2759 */ 2760 if ( method == ArcDistortion && !bestfit && !viewport_given ) { 2761 distort_image->page.x = 0; 2762 distort_image->page.y = 0; 2763 } 2764 coeff = (double *) RelinquishMagickMemory(coeff); 2765 return(distort_image); 2766 } 2767 2768 /* 2770 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2771 % % 2772 % % 2773 % % 2774 % R o t a t e I m a g e % 2775 % % 2776 % % 2777 % % 2778 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2779 % 2780 % RotateImage() creates a new image that is a rotated copy of an existing 2781 % one. Positive angles rotate counter-clockwise (right-hand rule), while 2782 % negative angles rotate clockwise. Rotated images are usually larger than 2783 % the originals and have 'empty' triangular corners. X axis. Empty 2784 % triangles left over from shearing the image are filled with the background 2785 % color defined by member 'background_color' of the image. RotateImage 2786 % allocates the memory necessary for the new Image structure and returns a 2787 % pointer to the new image. 2788 % 2789 % The format of the RotateImage method is: 2790 % 2791 % Image *RotateImage(const Image *image,const double degrees, 2792 % ExceptionInfo *exception) 2793 % 2794 % A description of each parameter follows. 2795 % 2796 % o image: the image. 2797 % 2798 % o degrees: Specifies the number of degrees to rotate the image. 2799 % 2800 % o exception: return any errors or warnings in this structure. 2801 % 2802 */ 2803 MagickExport Image *RotateImage(const Image *image,const double degrees, 2804 ExceptionInfo *exception) 2805 { 2806 Image 2807 *distort_image, 2808 *rotate_image; 2809 2810 double 2811 angle; 2812 2813 PointInfo 2814 shear; 2815 2816 size_t 2817 rotations; 2818 2819 /* 2820 Adjust rotation angle. 2821 */ 2822 assert(image != (Image *) NULL); 2823 assert(image->signature == MagickCoreSignature); 2824 if (image->debug != MagickFalse) 2825 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 2826 assert(exception != (ExceptionInfo *) NULL); 2827 assert(exception->signature == MagickCoreSignature); 2828 angle=degrees; 2829 while (angle < -45.0) 2830 angle+=360.0; 2831 for (rotations=0; angle > 45.0; rotations++) 2832 angle-=90.0; 2833 rotations%=4; 2834 shear.x=(-tan((double) DegreesToRadians(angle)/2.0)); 2835 shear.y=sin((double) DegreesToRadians(angle)); 2836 if ((fabs(shear.x) < MagickEpsilon) && (fabs(shear.y) < MagickEpsilon)) 2837 return(IntegralRotateImage(image,rotations,exception)); 2838 distort_image=CloneImage(image,0,0,MagickTrue,exception); 2839 if (distort_image == (Image *) NULL) 2840 return((Image *) NULL); 2841 (void) SetImageVirtualPixelMethod(distort_image,BackgroundVirtualPixelMethod, 2842 exception); 2843 rotate_image=DistortImage(distort_image,ScaleRotateTranslateDistortion,1, 2844 °rees,MagickTrue,exception); 2845 distort_image=DestroyImage(distort_image); 2846 return(rotate_image); 2847 } 2848 2849 /* 2851 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2852 % % 2853 % % 2854 % % 2855 % S p a r s e C o l o r I m a g e % 2856 % % 2857 % % 2858 % % 2859 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2860 % 2861 % SparseColorImage(), given a set of coordinates, interpolates the colors 2862 % found at those coordinates, across the whole image, using various methods. 2863 % 2864 % The format of the SparseColorImage() method is: 2865 % 2866 % Image *SparseColorImage(const Image *image, 2867 % const SparseColorMethod method,const size_t number_arguments, 2868 % const double *arguments,ExceptionInfo *exception) 2869 % 2870 % A description of each parameter follows: 2871 % 2872 % o image: the image to be filled in. 2873 % 2874 % o method: the method to fill in the gradient between the control points. 2875 % 2876 % The methods used for SparseColor() are often simular to methods 2877 % used for DistortImage(), and even share the same code for determination 2878 % of the function coefficents, though with more dimensions (or resulting 2879 % values). 2880 % 2881 % o number_arguments: the number of arguments given. 2882 % 2883 % o arguments: array of floating point arguments for this method-- 2884 % x,y,color_values-- with color_values given as normalized values. 2885 % 2886 % o exception: return any errors or warnings in this structure 2887 % 2888 */ 2889 MagickExport Image *SparseColorImage(const Image *image, 2890 const SparseColorMethod method,const size_t number_arguments, 2891 const double *arguments,ExceptionInfo *exception) 2892 { 2893 #define SparseColorTag "Distort/SparseColor" 2894 2895 SparseColorMethod 2896 sparse_method; 2897 2898 double 2899 *coeff; 2900 2901 Image 2902 *sparse_image; 2903 2904 size_t 2905 number_colors; 2906 2907 assert(image != (Image *) NULL); 2908 assert(image->signature == MagickCoreSignature); 2909 if (image->debug != MagickFalse) 2910 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 2911 assert(exception != (ExceptionInfo *) NULL); 2912 assert(exception->signature == MagickCoreSignature); 2913 2914 /* Determine number of color values needed per control point */ 2915 number_colors=0; 2916 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0) 2917 number_colors++; 2918 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0) 2919 number_colors++; 2920 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0) 2921 number_colors++; 2922 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) && 2923 (image->colorspace == CMYKColorspace)) 2924 number_colors++; 2925 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) && 2926 (image->alpha_trait != UndefinedPixelTrait)) 2927 number_colors++; 2928 2929 /* 2930 Convert input arguments into mapping coefficients, this this case 2931 we are mapping (distorting) colors, rather than coordinates. 2932 */ 2933 { DistortMethod 2934 distort_method; 2935 2936 distort_method=(DistortMethod) method; 2937 if ( distort_method >= SentinelDistortion ) 2938 distort_method = ShepardsDistortion; /* Pretend to be Shepards */ 2939 coeff = GenerateCoefficients(image, &distort_method, number_arguments, 2940 arguments, number_colors, exception); 2941 if ( coeff == (double *) NULL ) 2942 return((Image *) NULL); 2943 /* 2944 Note some Distort Methods may fall back to other simpler methods, 2945 Currently the only fallback of concern is Bilinear to Affine 2946 (Barycentric), which is alaso sparse_colr method. This also ensures 2947 correct two and one color Barycentric handling. 2948 */ 2949 sparse_method = (SparseColorMethod) distort_method; 2950 if ( distort_method == ShepardsDistortion ) 2951 sparse_method = method; /* return non-distort methods to normal */ 2952 if ( sparse_method == InverseColorInterpolate ) 2953 coeff[0]=0.5; /* sqrt() the squared distance for inverse */ 2954 } 2955 2956 /* Verbose output */ 2957 if (IsStringTrue(GetImageArtifact(image,"verbose")) != MagickFalse) { 2958 2959 switch (sparse_method) { 2960 case BarycentricColorInterpolate: 2961 { 2962 register ssize_t x=0; 2963 (void) FormatLocaleFile(stderr, "Barycentric Sparse Color:\n"); 2964 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0) 2965 (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf' \\\n", 2966 coeff[x], coeff[x+1], coeff[x+2]),x+=3; 2967 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0) 2968 (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf' \\\n", 2969 coeff[x], coeff[x+1], coeff[x+2]),x+=3; 2970 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0) 2971 (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf' \\\n", 2972 coeff[x], coeff[x+1], coeff[x+2]),x+=3; 2973 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) && 2974 (image->colorspace == CMYKColorspace)) 2975 (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf' \\\n", 2976 coeff[x], coeff[x+1], coeff[x+2]),x+=3; 2977 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) && 2978 (image->alpha_trait != UndefinedPixelTrait)) 2979 (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf' \\\n", 2980 coeff[x], coeff[x+1], coeff[x+2]),x+=3; 2981 break; 2982 } 2983 case BilinearColorInterpolate: 2984 { 2985 register ssize_t x=0; 2986 (void) FormatLocaleFile(stderr, "Bilinear Sparse Color\n"); 2987 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0) 2988 (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n", 2989 coeff[ x ], coeff[x+1], 2990 coeff[x+2], coeff[x+3]),x+=4; 2991 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0) 2992 (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n", 2993 coeff[ x ], coeff[x+1], 2994 coeff[x+2], coeff[x+3]),x+=4; 2995 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0) 2996 (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n", 2997 coeff[ x ], coeff[x+1], 2998 coeff[x+2], coeff[x+3]),x+=4; 2999 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) && 3000 (image->colorspace == CMYKColorspace)) 3001 (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n", 3002 coeff[ x ], coeff[x+1], 3003 coeff[x+2], coeff[x+3]),x+=4; 3004 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) && 3005 (image->alpha_trait != UndefinedPixelTrait)) 3006 (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n", 3007 coeff[ x ], coeff[x+1], 3008 coeff[x+2], coeff[x+3]),x+=4; 3009 break; 3010 } 3011 default: 3012 /* sparse color method is too complex for FX emulation */ 3013 break; 3014 } 3015 } 3016 3017 /* Generate new image for generated interpolated gradient. 3018 * ASIDE: Actually we could have just replaced the colors of the original 3019 * image, but IM Core policy, is if storage class could change then clone 3020 * the image. 3021 */ 3022 3023 sparse_image=CloneImage(image,0,0,MagickTrue,exception); 3024 if (sparse_image == (Image *) NULL) 3025 return((Image *) NULL); 3026 if (SetImageStorageClass(sparse_image,DirectClass,exception) == MagickFalse) 3027 { /* if image is ColorMapped - change it to DirectClass */ 3028 sparse_image=DestroyImage(sparse_image); 3029 return((Image *) NULL); 3030 } 3031 { /* ----- MAIN CODE ----- */ 3032 CacheView 3033 *sparse_view; 3034 3035 MagickBooleanType 3036 status; 3037 3038 MagickOffsetType 3039 progress; 3040 3041 ssize_t 3042 j; 3043 3044 status=MagickTrue; 3045 progress=0; 3046 sparse_view=AcquireAuthenticCacheView(sparse_image,exception); 3047 #if defined(MAGICKCORE_OPENMP_SUPPORT) 3048 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 3049 magick_threads(image,sparse_image,sparse_image->rows,1) 3050 #endif 3051 for (j=0; j < (ssize_t) sparse_image->rows; j++) 3052 { 3053 MagickBooleanType 3054 sync; 3055 3056 PixelInfo 3057 pixel; /* pixel to assign to distorted image */ 3058 3059 register ssize_t 3060 i; 3061 3062 register Quantum 3063 *magick_restrict q; 3064 3065 q=GetCacheViewAuthenticPixels(sparse_view,0,j,sparse_image->columns, 3066 1,exception); 3067 if (q == (Quantum *) NULL) 3068 { 3069 status=MagickFalse; 3070 continue; 3071 } 3072 GetPixelInfo(sparse_image,&pixel); 3073 for (i=0; i < (ssize_t) image->columns; i++) 3074 { 3075 GetPixelInfoPixel(image,q,&pixel); 3076 switch (sparse_method) 3077 { 3078 case BarycentricColorInterpolate: 3079 { 3080 register ssize_t x=0; 3081 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0) 3082 pixel.red = coeff[x]*i +coeff[x+1]*j 3083 +coeff[x+2], x+=3; 3084 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0) 3085 pixel.green = coeff[x]*i +coeff[x+1]*j 3086 +coeff[x+2], x+=3; 3087 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0) 3088 pixel.blue = coeff[x]*i +coeff[x+1]*j 3089 +coeff[x+2], x+=3; 3090 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) && 3091 (image->colorspace == CMYKColorspace)) 3092 pixel.black = coeff[x]*i +coeff[x+1]*j 3093 +coeff[x+2], x+=3; 3094 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) && 3095 (image->alpha_trait != UndefinedPixelTrait)) 3096 pixel.alpha = coeff[x]*i +coeff[x+1]*j 3097 +coeff[x+2], x+=3; 3098 break; 3099 } 3100 case BilinearColorInterpolate: 3101 { 3102 register ssize_t x=0; 3103 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0) 3104 pixel.red = coeff[x]*i + coeff[x+1]*j + 3105 coeff[x+2]*i*j + coeff[x+3], x+=4; 3106 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0) 3107 pixel.green = coeff[x]*i + coeff[x+1]*j + 3108 coeff[x+2]*i*j + coeff[x+3], x+=4; 3109 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0) 3110 pixel.blue = coeff[x]*i + coeff[x+1]*j + 3111 coeff[x+2]*i*j + coeff[x+3], x+=4; 3112 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) && 3113 (image->colorspace == CMYKColorspace)) 3114 pixel.black = coeff[x]*i + coeff[x+1]*j + 3115 coeff[x+2]*i*j + coeff[x+3], x+=4; 3116 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) && 3117 (image->alpha_trait != UndefinedPixelTrait)) 3118 pixel.alpha = coeff[x]*i + coeff[x+1]*j + 3119 coeff[x+2]*i*j + coeff[x+3], x+=4; 3120 break; 3121 } 3122 case InverseColorInterpolate: 3123 case ShepardsColorInterpolate: 3124 { /* Inverse (Squared) Distance weights average (IDW) */ 3125 size_t 3126 k; 3127 double 3128 denominator; 3129 3130 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0) 3131 pixel.red=0.0; 3132 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0) 3133 pixel.green=0.0; 3134 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0) 3135 pixel.blue=0.0; 3136 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) && 3137 (image->colorspace == CMYKColorspace)) 3138 pixel.black=0.0; 3139 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) && 3140 (image->alpha_trait != UndefinedPixelTrait)) 3141 pixel.alpha=0.0; 3142 denominator = 0.0; 3143 for(k=0; k<number_arguments; k+=2+number_colors) { 3144 register ssize_t x=(ssize_t) k+2; 3145 double weight = 3146 ((double)i-arguments[ k ])*((double)i-arguments[ k ]) 3147 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]); 3148 weight = pow(weight,coeff[0]); /* inverse of power factor */ 3149 weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight; 3150 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0) 3151 pixel.red += arguments[x++]*weight; 3152 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0) 3153 pixel.green += arguments[x++]*weight; 3154 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0) 3155 pixel.blue += arguments[x++]*weight; 3156 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) && 3157 (image->colorspace == CMYKColorspace)) 3158 pixel.black += arguments[x++]*weight; 3159 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) && 3160 (image->alpha_trait != UndefinedPixelTrait)) 3161 pixel.alpha += arguments[x++]*weight; 3162 denominator += weight; 3163 } 3164 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0) 3165 pixel.red/=denominator; 3166 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0) 3167 pixel.green/=denominator; 3168 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0) 3169 pixel.blue/=denominator; 3170 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) && 3171 (image->colorspace == CMYKColorspace)) 3172 pixel.black/=denominator; 3173 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) && 3174 (image->alpha_trait != UndefinedPixelTrait)) 3175 pixel.alpha/=denominator; 3176 break; 3177 } 3178 case ManhattanColorInterpolate: 3179 { 3180 size_t 3181 k; 3182 3183 double 3184 minimum = MagickMaximumValue; 3185 3186 /* 3187 Just use the closest control point you can find! 3188 */ 3189 for(k=0; k<number_arguments; k+=2+number_colors) { 3190 double distance = 3191 fabs((double)i-arguments[ k ]) 3192 + fabs((double)j-arguments[k+1]); 3193 if ( distance < minimum ) { 3194 register ssize_t x=(ssize_t) k+2; 3195 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0) 3196 pixel.red=arguments[x++]; 3197 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0) 3198 pixel.green=arguments[x++]; 3199 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0) 3200 pixel.blue=arguments[x++]; 3201 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) && 3202 (image->colorspace == CMYKColorspace)) 3203 pixel.black=arguments[x++]; 3204 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) && 3205 (image->alpha_trait != UndefinedPixelTrait)) 3206 pixel.alpha=arguments[x++]; 3207 minimum = distance; 3208 } 3209 } 3210 break; 3211 } 3212 case VoronoiColorInterpolate: 3213 default: 3214 { 3215 size_t 3216 k; 3217 3218 double 3219 minimum = MagickMaximumValue; 3220 3221 /* 3222 Just use the closest control point you can find! 3223 */ 3224 for (k=0; k<number_arguments; k+=2+number_colors) { 3225 double distance = 3226 ((double)i-arguments[ k ])*((double)i-arguments[ k ]) 3227 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]); 3228 if ( distance < minimum ) { 3229 register ssize_t x=(ssize_t) k+2; 3230 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0) 3231 pixel.red=arguments[x++]; 3232 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0) 3233 pixel.green=arguments[x++]; 3234 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0) 3235 pixel.blue=arguments[x++]; 3236 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) && 3237 (image->colorspace == CMYKColorspace)) 3238 pixel.black=arguments[x++]; 3239 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) && 3240 (image->alpha_trait != UndefinedPixelTrait)) 3241 pixel.alpha=arguments[x++]; 3242 minimum = distance; 3243 } 3244 } 3245 break; 3246 } 3247 } 3248 /* set the color directly back into the source image */ 3249 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0) 3250 pixel.red=ClampPixel(QuantumRange*pixel.red); 3251 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0) 3252 pixel.green=ClampPixel(QuantumRange*pixel.green); 3253 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0) 3254 pixel.blue=ClampPixel(QuantumRange*pixel.blue); 3255 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) && 3256 (image->colorspace == CMYKColorspace)) 3257 pixel.black=ClampPixel(QuantumRange*pixel.black); 3258 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) && 3259 (image->alpha_trait != UndefinedPixelTrait)) 3260 pixel.alpha=ClampPixel(QuantumRange*pixel.alpha); 3261 SetPixelViaPixelInfo(sparse_image,&pixel,q); 3262 q+=GetPixelChannels(sparse_image); 3263 } 3264 sync=SyncCacheViewAuthenticPixels(sparse_view,exception); 3265 if (sync == MagickFalse) 3266 status=MagickFalse; 3267 if (image->progress_monitor != (MagickProgressMonitor) NULL) 3268 { 3269 MagickBooleanType 3270 proceed; 3271 3272 #if defined(MAGICKCORE_OPENMP_SUPPORT) 3273 #pragma omp critical (MagickCore_SparseColorImage) 3274 #endif 3275 proceed=SetImageProgress(image,SparseColorTag,progress++,image->rows); 3276 if (proceed == MagickFalse) 3277 status=MagickFalse; 3278 } 3279 } 3280 sparse_view=DestroyCacheView(sparse_view); 3281 if (status == MagickFalse) 3282 sparse_image=DestroyImage(sparse_image); 3283 } 3284 coeff = (double *) RelinquishMagickMemory(coeff); 3285 return(sparse_image); 3286 } 3287