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