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