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