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