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