MagickCore 6.9.13-50
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 (number_arguments < 1) {
1289 coeff = (double *) RelinquishMagickMemory(coeff);
1290 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1291 "InvalidArgument", "%s : 'Needs at least 1 argument'",
1292 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1293 return((double *) NULL);
1294 }
1295 if ( arguments[0] < MagickEpsilon || arguments[0] > 160.0 ) {
1296 coeff=(double *) RelinquishMagickMemory(coeff);
1297 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1298 "InvalidArgument", "%s : Invalid FOV Angle",
1299 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1300 return((double *) NULL);
1301 }
1302 coeff[0] = DegreesToRadians(arguments[0]);
1303 if ( *method == Cylinder2PlaneDistortion )
1304 /* image is curved around cylinder, so FOV angle (in radians)
1305 * scales directly to image X coordinate, according to its radius.
1306 */
1307 coeff[1] = (double) image->columns/coeff[0];
1308 else
1309 /* radius is distance away from an image with this angular FOV */
1310 coeff[1] = (double) image->columns / ( 2 * tan(coeff[0]/2) );
1311
1312 coeff[2] = (double)(image->columns)/2.0+image->page.x;
1313 coeff[3] = (double)(image->rows)/2.0+image->page.y;
1314 coeff[4] = coeff[2];
1315 coeff[5] = coeff[3]; /* assuming image size is the same */
1316 return(coeff);
1317 }
1318 case BarrelDistortion:
1319 case BarrelInverseDistortion:
1320 {
1321 /* Barrel Distortion
1322 Rs=(A*Rd^3 + B*Rd^2 + C*Rd + D)*Rd
1323 BarrelInv Distortion
1324 Rs=Rd/(A*Rd^3 + B*Rd^2 + C*Rd + D)
1325
1326 Where Rd is the normalized radius from corner to middle of image
1327 Input Arguments are one of the following forms (number of arguments)...
1328 3: A,B,C
1329 4: A,B,C,D
1330 5: A,B,C X,Y
1331 6: A,B,C,D X,Y
1332 8: Ax,Bx,Cx,Dx Ay,By,Cy,Dy
1333 10: Ax,Bx,Cx,Dx Ay,By,Cy,Dy X,Y
1334
1335 Returns 10 coefficient values, which are de-normalized (pixel scale)
1336 Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Xc, Yc
1337 */
1338 /* Radius de-normalization scaling factor */
1339 double
1340 rscale = 2.0/MagickMin((double) image->columns,(double) image->rows);
1341
1342 /* sanity check number of args must = 3,4,5,6,8,10 or error */
1343 if ( (number_arguments < 3) || (number_arguments == 7) ||
1344 (number_arguments == 9) || (number_arguments > 10) )
1345 {
1346 coeff=(double *) RelinquishMagickMemory(coeff);
1347 (void) ThrowMagickException(exception,GetMagickModule(),
1348 OptionError,"InvalidArgument", "%s : number of arguments",
1349 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1350 return((double *) NULL);
1351 }
1352 /* A,B,C,D coefficients */
1353 coeff[0] = arguments[0];
1354 coeff[1] = arguments[1];
1355 coeff[2] = arguments[2];
1356 if ((number_arguments == 3) || (number_arguments == 5) )
1357 coeff[3] = 1.0 - coeff[0] - coeff[1] - coeff[2];
1358 else
1359 coeff[3] = arguments[3];
1360 /* de-normalize the coefficients */
1361 coeff[0] *= pow(rscale,3.0);
1362 coeff[1] *= rscale*rscale;
1363 coeff[2] *= rscale;
1364 /* Y coefficients: as given OR same as X coefficients */
1365 if ( number_arguments >= 8 ) {
1366 coeff[4] = arguments[4] * pow(rscale,3.0);
1367 coeff[5] = arguments[5] * rscale*rscale;
1368 coeff[6] = arguments[6] * rscale;
1369 coeff[7] = arguments[7];
1370 }
1371 else {
1372 coeff[4] = coeff[0];
1373 coeff[5] = coeff[1];
1374 coeff[6] = coeff[2];
1375 coeff[7] = coeff[3];
1376 }
1377 /* X,Y Center of Distortion (image coordinates) */
1378 if ( number_arguments == 5 ) {
1379 coeff[8] = arguments[3];
1380 coeff[9] = arguments[4];
1381 }
1382 else if ( number_arguments == 6 ) {
1383 coeff[8] = arguments[4];
1384 coeff[9] = arguments[5];
1385 }
1386 else if ( number_arguments == 10 ) {
1387 coeff[8] = arguments[8];
1388 coeff[9] = arguments[9];
1389 }
1390 else {
1391 /* center of the image provided (image coordinates) */
1392 coeff[8] = (double)image->columns/2.0 + image->page.x;
1393 coeff[9] = (double)image->rows/2.0 + image->page.y;
1394 }
1395 return(coeff);
1396 }
1397 case ShepardsDistortion:
1398 {
1399 /* Shepards Distortion input arguments are the coefficients!
1400 Just check the number of arguments is valid!
1401 Args: u1,v1, x1,y1, ...
1402 OR : u1,v1, r1,g1,c1, ...
1403 */
1404 if ( number_arguments%cp_size != 0 ||
1405 number_arguments < cp_size ) {
1406 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1407 "InvalidArgument", "%s : 'requires CP's (4 numbers each)'",
1408 CommandOptionToMnemonic(MagickDistortOptions, *method));
1409 coeff=(double *) RelinquishMagickMemory(coeff);
1410 return((double *) NULL);
1411 }
1412 /* User defined weighting power for Shepard's Method */
1413 { const char *artifact=GetImageArtifact(image,"shepards:power");
1414 if ( artifact != (const char *) NULL ) {
1415 coeff[0]=StringToDouble(artifact,(char **) NULL) / 2.0;
1416 if ( coeff[0] < MagickEpsilon ) {
1417 (void) ThrowMagickException(exception,GetMagickModule(),
1418 OptionError,"InvalidArgument","%s", "-define shepards:power" );
1419 coeff=(double *) RelinquishMagickMemory(coeff);
1420 return((double *) NULL);
1421 }
1422 }
1423 else
1424 coeff[0]=1.0; /* Default power of 2 (Inverse Squared) */
1425 }
1426 return(coeff);
1427 }
1428 default:
1429 break;
1430 }
1431 /* you should never reach this point */
1432 perror("no method handler"); /* just fail assertion */
1433 return((double *) NULL);
1434}
1435␌
1436/*
1437%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1438% %
1439% %
1440% %
1441+ D i s t o r t R e s i z e I m a g e %
1442% %
1443% %
1444% %
1445%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1446%
1447% DistortResizeImage() resize image using the equivalent but slower image
1448% distortion operator. The filter is applied using a EWA cylindrical
1449% resampling. But like resize the final image size is limited to whole pixels
1450% with no effects by virtual-pixels on the result.
1451%
1452% Note that images containing a transparency channel will be twice as slow to
1453% resize as images one without transparency.
1454%
1455% The format of the DistortResizeImage method is:
1456%
1457% Image *DistortResizeImage(const Image *image,const size_t columns,
1458% const size_t rows,ExceptionInfo *exception)
1459%
1460% A description of each parameter follows:
1461%
1462% o image: the image.
1463%
1464% o columns: the number of columns in the resized image.
1465%
1466% o rows: the number of rows in the resized image.
1467%
1468% o exception: return any errors or warnings in this structure.
1469%
1470*/
1471MagickExport Image *DistortResizeImage(const Image *image,
1472 const size_t columns,const size_t rows,ExceptionInfo *exception)
1473{
1474#define DistortResizeImageTag "Distort/Image"
1475
1476 Image
1477 *resize_image,
1478 *tmp_image;
1479
1480 RectangleInfo
1481 crop_area;
1482
1483 double
1484 distort_args[12];
1485
1486 VirtualPixelMethod
1487 vp_save;
1488
1489 /*
1490 Distort resize image.
1491 */
1492 assert(image != (const Image *) NULL);
1493 assert(image->signature == MagickCoreSignature);
1494 assert(exception != (ExceptionInfo *) NULL);
1495 assert(exception->signature == MagickCoreSignature);
1496 if (IsEventLogging() != MagickFalse)
1497 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1498 if ((columns == 0) || (rows == 0))
1499 return((Image *) NULL);
1500 /* Do not short-circuit this resize if final image size is unchanged */
1501
1502 (void) memset(distort_args,0,12*sizeof(double));
1503 distort_args[4]=(double) image->columns;
1504 distort_args[6]=(double) columns;
1505 distort_args[9]=(double) image->rows;
1506 distort_args[11]=(double) rows;
1507
1508 vp_save=GetImageVirtualPixelMethod(image);
1509
1510 tmp_image=CloneImage(image,0,0,MagickTrue,exception);
1511 if ( tmp_image == (Image *) NULL )
1512 return((Image *) NULL);
1513 (void) SetImageVirtualPixelMethod(tmp_image,TransparentVirtualPixelMethod);
1514
1515 if (image->matte == MagickFalse)
1516 {
1517 /*
1518 Image has not transparency channel, so we free to use it
1519 */
1520 (void) SetImageAlphaChannel(tmp_image,SetAlphaChannel);
1521 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1522 MagickTrue,exception),
1523
1524 tmp_image=DestroyImage(tmp_image);
1525 if ( resize_image == (Image *) NULL )
1526 return((Image *) NULL);
1527
1528 (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel);
1529 InheritException(exception,&image->exception);
1530 }
1531 else
1532 {
1533 /*
1534 Image has transparency so handle colors and alpha separately.
1535 Basically we need to separate Virtual-Pixel alpha in the resized
1536 image, so only the actual original images alpha channel is used.
1537 */
1538 Image
1539 *resize_alpha;
1540
1541 /* distort alpha channel separately */
1542 (void) SeparateImageChannel(tmp_image,TrueAlphaChannel);
1543 (void) SetImageAlphaChannel(tmp_image,OpaqueAlphaChannel);
1544 resize_alpha=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1545 MagickTrue,exception),
1546 tmp_image=DestroyImage(tmp_image);
1547 if ( resize_alpha == (Image *) NULL )
1548 return((Image *) NULL);
1549
1550 /* distort the actual image containing alpha + VP alpha */
1551 tmp_image=CloneImage(image,0,0,MagickTrue,exception);
1552 if ( tmp_image == (Image *) NULL )
1553 return((Image *) NULL);
1554 (void) SetImageVirtualPixelMethod(tmp_image,
1555 TransparentVirtualPixelMethod);
1556 (void) SetImageVirtualPixelMethod(tmp_image,
1557 TransparentVirtualPixelMethod);
1558 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1559 MagickTrue,exception),
1560 tmp_image=DestroyImage(tmp_image);
1561 if ( resize_image == (Image *) NULL)
1562 {
1563 resize_alpha=DestroyImage(resize_alpha);
1564 return((Image *) NULL);
1565 }
1566
1567 /* replace resize images alpha with the separally distorted alpha */
1568 (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel);
1569 (void) SetImageAlphaChannel(resize_alpha,DeactivateAlphaChannel);
1570 (void) CompositeImage(resize_image,CopyOpacityCompositeOp,resize_alpha,
1571 0,0);
1572 InheritException(exception,&resize_image->exception);
1573 resize_image->matte=image->matte;
1574 resize_image->compose=image->compose;
1575 resize_alpha=DestroyImage(resize_alpha);
1576 }
1577 (void) SetImageVirtualPixelMethod(resize_image,vp_save);
1578
1579 /*
1580 Clean up the results of the Distortion
1581 */
1582 crop_area.width=columns;
1583 crop_area.height=rows;
1584 crop_area.x=0;
1585 crop_area.y=0;
1586
1587 tmp_image=resize_image;
1588 resize_image=CropImage(tmp_image,&crop_area,exception);
1589 tmp_image=DestroyImage(tmp_image);
1590 if (resize_image != (Image *) NULL)
1591 {
1592 resize_image->page.width=0;
1593 resize_image->page.height=0;
1594 }
1595 return(resize_image);
1596}
1597␌
1598/*
1599%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1600% %
1601% %
1602% %
1603% D i s t o r t I m a g e %
1604% %
1605% %
1606% %
1607%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1608%
1609% DistortImage() distorts an image using various distortion methods, by
1610% mapping color lookups of the source image to a new destination image
1611% usually of the same size as the source image, unless 'bestfit' is set to
1612% true.
1613%
1614% If 'bestfit' is enabled, and distortion allows it, the destination image is
1615% adjusted to ensure the whole source 'image' will just fit within the final
1616% destination image, which will be sized and offset accordingly. Also in
1617% many cases the virtual offset of the source image will be taken into
1618% account in the mapping.
1619%
1620% If the '-verbose' control option has been set print to standard error the
1621% equivalent '-fx' formula with coefficients for the function, if practical.
1622%
1623% The format of the DistortImage() method is:
1624%
1625% Image *DistortImage(const Image *image,const DistortImageMethod method,
1626% const size_t number_arguments,const double *arguments,
1627% MagickBooleanType bestfit, ExceptionInfo *exception)
1628%
1629% A description of each parameter follows:
1630%
1631% o image: the image to be distorted.
1632%
1633% o method: the method of image distortion.
1634%
1635% ArcDistortion always ignores source image offset, and always
1636% 'bestfit' the destination image with the top left corner offset
1637% relative to the polar mapping center.
1638%
1639% Affine, Perspective, and Bilinear, do least squares fitting of the
1640% distortion when more than the minimum number of control point pairs
1641% are provided.
1642%
1643% Perspective, and Bilinear, fall back to a Affine distortion when less
1644% than 4 control point pairs are provided. While Affine distortions
1645% let you use any number of control point pairs, that is Zero pairs is
1646% a No-Op (viewport only) distortion, one pair is a translation and
1647% two pairs of control points do a scale-rotate-translate, without any
1648% shearing.
1649%
1650% o number_arguments: the number of arguments given.
1651%
1652% o arguments: an array of floating point arguments for this method.
1653%
1654% o bestfit: Attempt to 'bestfit' the size of the resulting image.
1655% This also forces the resulting image to be a 'layered' virtual
1656% canvas image. Can be overridden using 'distort:viewport' setting.
1657%
1658% o exception: return any errors or warnings in this structure
1659%
1660% Extra Controls from Image meta-data (artifacts)...
1661%
1662% o "verbose"
1663% Output to stderr alternatives, internal coefficients, and FX
1664% equivalents for the distortion operation (if feasible).
1665% This forms an extra check of the distortion method, and allows users
1666% access to the internal constants IM calculates for the distortion.
1667%
1668% o "distort:viewport"
1669% Directly set the output image canvas area and offset to use for the
1670% resulting image, rather than use the original images canvas, or a
1671% calculated 'bestfit' canvas.
1672%
1673% o "distort:scale"
1674% Scale the size of the output canvas by this amount to provide a
1675% method of Zooming, and for super-sampling the results.
1676%
1677% Other settings that can effect results include
1678%
1679% o 'interpolate' For source image lookups (scale enlargements)
1680%
1681% o 'filter' Set filter to use for area-resampling (scale shrinking).
1682% Set to 'point' to turn off and use 'interpolate' lookup
1683% instead
1684%
1685*/
1686
1687MagickExport Image *DistortImage(const Image *image,DistortImageMethod method,
1688 const size_t number_arguments,const double *arguments,
1689 MagickBooleanType bestfit,ExceptionInfo *exception)
1690{
1691#define DistortImageTag "Distort/Image"
1692
1693 double
1694 *coeff,
1695 output_scaling;
1696
1697 Image
1698 *distort_image;
1699
1700 RectangleInfo
1701 geometry; /* geometry of the distorted space viewport */
1702
1703 MagickBooleanType
1704 viewport_given;
1705
1706 assert(image != (Image *) NULL);
1707 assert(image->signature == MagickCoreSignature);
1708 assert(exception != (ExceptionInfo *) NULL);
1709 assert(exception->signature == MagickCoreSignature);
1710 if (IsEventLogging() != MagickFalse)
1711 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1712 /*
1713 Handle Special Compound Distortions
1714 */
1715 if (method == ResizeDistortion)
1716 {
1717 if (number_arguments != 2)
1718 {
1719 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1720 "InvalidArgument","%s : '%s'","Resize",
1721 "Invalid number of args: 2 only");
1722 return((Image *) NULL);
1723 }
1724 distort_image=DistortResizeImage(image,CastDoubleToSizeT(arguments[0]),
1725 CastDoubleToSizeT(arguments[1]),exception);
1726 return(distort_image);
1727 }
1728
1729 /*
1730 Convert input arguments (usually as control points for reverse mapping)
1731 into mapping coefficients to apply the distortion.
1732
1733 Note that some distortions are mapped to other distortions,
1734 and as such do not require specific code after this point.
1735 */
1736 coeff=GenerateCoefficients(image,&method,number_arguments,arguments,0,
1737 exception);
1738 if (coeff == (double *) NULL)
1739 return((Image *) NULL);
1740
1741 /*
1742 Determine the size and offset for a 'bestfit' destination.
1743 Usually the four corners of the source image is enough.
1744 */
1745
1746 /* default output image bounds, when no 'bestfit' is requested */
1747 geometry.width=image->columns;
1748 geometry.height=image->rows;
1749 geometry.x=0;
1750 geometry.y=0;
1751
1752 if ( method == ArcDistortion ) {
1753 bestfit = MagickTrue; /* always calculate a 'best fit' viewport */
1754 }
1755
1756 /* Work out the 'best fit', (required for ArcDistortion) */
1757 if ( bestfit ) {
1758 PointInfo
1759 s,d,min,max; /* source, dest coords --mapping--> min, max coords */
1760
1761 MagickBooleanType
1762 fix_bounds = MagickTrue; /* enlarge bounds for VP handling */
1763
1764 s.x=s.y=min.x=max.x=min.y=max.y=0.0; /* keep compiler happy */
1765
1766/* defines to figure out the bounds of the distorted image */
1767#define InitalBounds(p) \
1768{ \
1769 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1770 min.x = max.x = p.x; \
1771 min.y = max.y = p.y; \
1772}
1773#define ExpandBounds(p) \
1774{ \
1775 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1776 min.x = MagickMin(min.x,p.x); \
1777 max.x = MagickMax(max.x,p.x); \
1778 min.y = MagickMin(min.y,p.y); \
1779 max.y = MagickMax(max.y,p.y); \
1780}
1781
1782 switch (method)
1783 {
1784 case AffineDistortion:
1785 { double inverse[6];
1786 InvertAffineCoefficients(coeff, inverse);
1787 s.x = (double) image->page.x;
1788 s.y = (double) image->page.y;
1789 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1790 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1791 InitalBounds(d);
1792 s.x = (double) image->page.x+image->columns;
1793 s.y = (double) image->page.y;
1794 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1795 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1796 ExpandBounds(d);
1797 s.x = (double) image->page.x;
1798 s.y = (double) image->page.y+image->rows;
1799 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1800 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1801 ExpandBounds(d);
1802 s.x = (double) image->page.x+image->columns;
1803 s.y = (double) image->page.y+image->rows;
1804 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1805 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1806 ExpandBounds(d);
1807 break;
1808 }
1809 case PerspectiveDistortion:
1810 { double inverse[8], scale;
1811 InvertPerspectiveCoefficients(coeff, inverse);
1812 s.x = (double) image->page.x;
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 InitalBounds(d);
1819 s.x = (double) image->page.x+image->columns;
1820 s.y = (double) image->page.y;
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;
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 s.x = (double) image->page.x+image->columns;
1834 s.y = (double) image->page.y+image->rows;
1835 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1836 scale=MagickSafeReciprocal(scale);
1837 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1838 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1839 ExpandBounds(d);
1840 break;
1841 }
1842 case ArcDistortion:
1843 { double a, ca, sa;
1844 /* Forward Map Corners */
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 InitalBounds(d);
1849 d.x = (coeff[2]-coeff[3])*ca;
1850 d.y = (coeff[2]-coeff[3])*sa;
1851 ExpandBounds(d);
1852 a = coeff[0]+coeff[1]/2; ca = cos(a); sa = sin(a);
1853 d.x = coeff[2]*ca;
1854 d.y = coeff[2]*sa;
1855 ExpandBounds(d);
1856 d.x = (coeff[2]-coeff[3])*ca;
1857 d.y = (coeff[2]-coeff[3])*sa;
1858 ExpandBounds(d);
1859 /* Orthogonal points along top of arc */
1860 for( a=(double) (ceil((double) ((coeff[0]-coeff[1]/2.0)/MagickPI2))*MagickPI2);
1861 a<(coeff[0]+coeff[1]/2.0); a+=MagickPI2 ) {
1862 ca = cos(a); sa = sin(a);
1863 d.x = coeff[2]*ca;
1864 d.y = coeff[2]*sa;
1865 ExpandBounds(d);
1866 }
1867 /*
1868 Convert the angle_to_width and radius_to_height
1869 to appropriate scaling factors, to allow faster processing
1870 in the mapping function.
1871 */
1872 coeff[1] = (double) (Magick2PI*image->columns/coeff[1]);
1873 coeff[3] = (double)image->rows/coeff[3];
1874 break;
1875 }
1876 case PolarDistortion:
1877 {
1878 if (number_arguments < 2)
1879 coeff[2] = coeff[3] = 0.0;
1880 min.x = coeff[2]-coeff[0];
1881 max.x = coeff[2]+coeff[0];
1882 min.y = coeff[3]-coeff[0];
1883 max.y = coeff[3]+coeff[0];
1884 /* should be about 1.0 if Rmin = 0 */
1885 coeff[7]=(double) geometry.height/(coeff[0]-coeff[1]);
1886 break;
1887 }
1888 case DePolarDistortion:
1889 {
1890 /* direct calculation as it needs to tile correctly
1891 * for reversibility in a DePolar-Polar cycle */
1892 fix_bounds = MagickFalse;
1893 geometry.x = geometry.y = 0;
1894 geometry.height = CastDoubleToSizeT(ceil(coeff[0]-coeff[1]));
1895 geometry.width = CastDoubleToSizeT(ceil((coeff[0]-coeff[1])*
1896 (coeff[5]-coeff[4])*0.5));
1897 /* correct scaling factors relative to new size */
1898 coeff[6]=(coeff[5]-coeff[4])*MagickSafeReciprocal(geometry.width); /* changed width */
1899 coeff[7]=(coeff[0]-coeff[1])*MagickSafeReciprocal(geometry.height); /* should be about 1.0 */
1900 break;
1901 }
1902 case Cylinder2PlaneDistortion:
1903 {
1904 /* direct calculation so center of distortion is either a pixel
1905 * center, or pixel edge. This allows for reversibility of the
1906 * distortion */
1907 geometry.x = geometry.y = 0;
1908 geometry.width = CastDoubleToSizeT(ceil( 2.0*coeff[1]*tan(coeff[0]/2.0) ));
1909 geometry.height = CastDoubleToSizeT(ceil( 2.0*coeff[3]/cos(coeff[0]/2.0) ));
1910 /* correct center of distortion relative to new size */
1911 coeff[4] = (double) geometry.width/2.0;
1912 coeff[5] = (double) geometry.height/2.0;
1913 fix_bounds = MagickFalse;
1914 break;
1915 }
1916 case Plane2CylinderDistortion:
1917 {
1918 /* direct calculation center is either pixel center, or pixel edge
1919 * so as to allow reversibility of the image distortion */
1920 geometry.x = geometry.y = 0;
1921 geometry.width = CastDoubleToSizeT(ceil(coeff[0]*coeff[1])); /* FOV * radius */
1922 geometry.height = CastDoubleToSizeT(2.0*coeff[3]); /* input image height */
1923 /* correct center of distortion relative to new size */
1924 coeff[4] = (double) geometry.width/2.0;
1925 coeff[5] = (double) geometry.height/2.0;
1926 fix_bounds = MagickFalse;
1927 break;
1928 }
1929
1930 case ShepardsDistortion:
1931 case BilinearForwardDistortion:
1932 case BilinearReverseDistortion:
1933#if 0
1934 case QuadrilateralDistortion:
1935#endif
1936 case PolynomialDistortion:
1937 case BarrelDistortion:
1938 case BarrelInverseDistortion:
1939 default:
1940 /* no calculated bestfit available for these distortions */
1941 bestfit = MagickFalse;
1942 fix_bounds = MagickFalse;
1943 break;
1944 }
1945
1946 /* Set the output image geometry to calculated 'bestfit'.
1947 Yes this tends to 'over do' the file image size, ON PURPOSE!
1948 Do not do this for DePolar which needs to be exact for virtual tiling.
1949 */
1950 if ( fix_bounds ) {
1951 geometry.x = (ssize_t) floor(min.x-0.5);
1952 geometry.y = (ssize_t) floor(min.y-0.5);
1953 geometry.width=CastDoubleToSizeT(ceil(max.x-geometry.x+0.5));
1954 geometry.height=CastDoubleToSizeT(ceil(max.y-geometry.y+0.5));
1955 }
1956
1957 } /* end bestfit destination image calculations */
1958
1959 /* The user provided a 'viewport' expert option which may
1960 overrides some parts of the current output image geometry.
1961 This also overrides its default 'bestfit' setting.
1962 */
1963 { const char *artifact=GetImageArtifact(image,"distort:viewport");
1964 viewport_given = MagickFalse;
1965 if ( artifact != (const char *) NULL ) {
1966 MagickStatusType flags=ParseAbsoluteGeometry(artifact,&geometry);
1967 if (flags==NoValue)
1968 (void) ThrowMagickException(exception,GetMagickModule(),
1969 OptionWarning,"InvalidGeometry","`%s' `%s'",
1970 "distort:viewport",artifact);
1971 else
1972 viewport_given = MagickTrue;
1973 }
1974 }
1975
1976 /* Verbose output */
1977 if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
1978 ssize_t
1979 i;
1980 char image_gen[MaxTextExtent];
1981 const char *lookup;
1982
1983 /* Set destination image size and virtual offset */
1984 if ( bestfit || viewport_given ) {
1985 (void) FormatLocaleString(image_gen, MaxTextExtent," -size %.20gx%.20g "
1986 "-page %+.20g%+.20g xc: +insert \\\n",(double) geometry.width,
1987 (double) geometry.height,(double) geometry.x,(double) geometry.y);
1988 lookup="v.p{ xx-v.page.x-.5, yy-v.page.y-.5 }";
1989 }
1990 else {
1991 image_gen[0] = '\0'; /* no destination to generate */
1992 lookup = "p{ xx-page.x-.5, yy-page.y-.5 }"; /* simplify lookup */
1993 }
1994
1995 switch (method)
1996 {
1997 case AffineDistortion:
1998 {
1999 double
2000 *inverse;
2001
2002 inverse=(double *) AcquireQuantumMemory(6,sizeof(*inverse));
2003 if (inverse == (double *) NULL)
2004 {
2005 coeff=(double *) RelinquishMagickMemory(coeff);
2006 (void) ThrowMagickException(exception,GetMagickModule(),
2007 ResourceLimitError,"MemoryAllocationFailed","%s","DistortImages");
2008 return((Image *) NULL);
2009 }
2010 InvertAffineCoefficients(coeff, inverse);
2011 CoefficientsToAffineArgs(inverse);
2012 (void) FormatLocaleFile(stderr, "Affine Projection:\n");
2013 (void) FormatLocaleFile(stderr,
2014 " -distort AffineProjection \\\n '");
2015 for (i=0; i < 5; i++)
2016 (void) FormatLocaleFile(stderr, "%lf,", inverse[i]);
2017 (void) FormatLocaleFile(stderr, "%lf'\n", inverse[5]);
2018 inverse=(double *) RelinquishMagickMemory(inverse);
2019 (void) FormatLocaleFile(stderr, "Affine Distort, FX Equivelent:\n");
2020 (void) FormatLocaleFile(stderr, "%s", image_gen);
2021 (void) FormatLocaleFile(stderr,
2022 " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2023 (void) FormatLocaleFile(stderr," xx=%+lf*ii %+lf*jj %+lf;\n",
2024 coeff[0],coeff[1],coeff[2]);
2025 (void) FormatLocaleFile(stderr," yy=%+lf*ii %+lf*jj %+lf;\n",
2026 coeff[3],coeff[4],coeff[5]);
2027 (void) FormatLocaleFile(stderr," %s' \\\n",lookup);
2028 break;
2029 }
2030 case PerspectiveDistortion:
2031 {
2032 double
2033 *inverse;
2034
2035 inverse=(double *) AcquireQuantumMemory(8,sizeof(*inverse));
2036 if (inverse == (double *) NULL)
2037 {
2038 coeff=(double *) RelinquishMagickMemory(coeff);
2039 (void) ThrowMagickException(exception,GetMagickModule(),
2040 ResourceLimitError,"MemoryAllocationFailed","%s",
2041 "DistortCoefficients");
2042 return((Image *) NULL);
2043 }
2044 InvertPerspectiveCoefficients(coeff, inverse);
2045 (void) FormatLocaleFile(stderr,"Perspective Projection:\n");
2046 (void) FormatLocaleFile(stderr,
2047 " -distort PerspectiveProjection \\\n '");
2048 for (i=0; i < 4; i++)
2049 (void) FormatLocaleFile(stderr, "%.*g, ",GetMagickPrecision(),
2050 inverse[i]);
2051 (void) FormatLocaleFile(stderr, "\n ");
2052 for ( ; i < 7; i++)
2053 (void) FormatLocaleFile(stderr, "%.*g, ",GetMagickPrecision(),
2054 inverse[i]);
2055 (void) FormatLocaleFile(stderr, "%.*g'\n",GetMagickPrecision(),
2056 inverse[7]);
2057 inverse=(double *) RelinquishMagickMemory(inverse);
2058 (void) FormatLocaleFile(stderr,"Perspective Distort, FX Equivalent:\n");
2059 (void) FormatLocaleFile(stderr,"%.1024s",image_gen);
2060 (void) FormatLocaleFile(stderr,
2061 " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2062 (void) FormatLocaleFile(stderr," rr=%+.*g*ii %+.*g*jj + 1;\n",
2063 GetMagickPrecision(),coeff[6],GetMagickPrecision(),coeff[7]);
2064 (void) FormatLocaleFile(stderr,
2065 " xx=(%+.*g*ii %+.*g*jj %+.*g)/rr;\n",
2066 GetMagickPrecision(),coeff[0],GetMagickPrecision(),coeff[1],
2067 GetMagickPrecision(),coeff[2]);
2068 (void) FormatLocaleFile(stderr,
2069 " yy=(%+.*g*ii %+.*g*jj %+.*g)/rr;\n",
2070 GetMagickPrecision(),coeff[3],GetMagickPrecision(),coeff[4],
2071 GetMagickPrecision(),coeff[5]);
2072 (void) FormatLocaleFile(stderr," rr%s0 ? %s : blue' \\\n",
2073 coeff[8] < 0.0 ? "<" : ">", lookup);
2074 break;
2075 }
2076 case BilinearForwardDistortion:
2077 {
2078 (void) FormatLocaleFile(stderr,"BilinearForward Mapping Equations:\n");
2079 (void) FormatLocaleFile(stderr,"%s", image_gen);
2080 (void) FormatLocaleFile(stderr," i = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
2081 coeff[0],coeff[1],coeff[2],coeff[3]);
2082 (void) FormatLocaleFile(stderr," j = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
2083 coeff[4],coeff[5],coeff[6],coeff[7]);
2084#if 0
2085 /* for debugging */
2086 (void) FormatLocaleFile(stderr, " c8 = %+lf c9 = 2*a = %+lf;\n",
2087 coeff[8], coeff[9]);
2088#endif
2089 (void) FormatLocaleFile(stderr,
2090 "BilinearForward Distort, FX Equivalent:\n");
2091 (void) FormatLocaleFile(stderr,"%s", image_gen);
2092 (void) FormatLocaleFile(stderr,
2093 " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",0.5-coeff[3],0.5-
2094 coeff[7]);
2095 (void) FormatLocaleFile(stderr," bb=%lf*ii %+lf*jj %+lf;\n",
2096 coeff[6], -coeff[2], coeff[8]);
2097 /* Handle Special degenerate (non-quadratic) or trapezoidal case */
2098 if (coeff[9] != 0)
2099 {
2100 (void) FormatLocaleFile(stderr,
2101 " rt=bb*bb %+lf*(%lf*ii%+lf*jj);\n",-2*coeff[9],coeff[4],
2102 -coeff[0]);
2103 (void) FormatLocaleFile(stderr,
2104 " yy=( -bb + sqrt(rt) ) / %lf;\n",coeff[9]);
2105 }
2106 else
2107 (void) FormatLocaleFile(stderr," yy=(%lf*ii%+lf*jj)/bb;\n",
2108 -coeff[4],coeff[0]);
2109 (void) FormatLocaleFile(stderr,
2110 " xx=(ii %+lf*yy)/(%lf %+lf*yy);\n",-coeff[1],coeff[0],
2111 coeff[2]);
2112 if ( coeff[9] != 0 )
2113 (void) FormatLocaleFile(stderr," (rt < 0 ) ? red : %s'\n",
2114 lookup);
2115 else
2116 (void) FormatLocaleFile(stderr," %s' \\\n", lookup);
2117 break;
2118 }
2119 case BilinearReverseDistortion:
2120 {
2121#if 0
2122 (void) FormatLocaleFile(stderr, "Polynomial Projection Distort:\n");
2123 (void) FormatLocaleFile(stderr, " -distort PolynomialProjection \\\n");
2124 (void) FormatLocaleFile(stderr, " '1.5, %lf, %lf, %lf, %lf,\n",
2125 coeff[3], coeff[0], coeff[1], coeff[2]);
2126 (void) FormatLocaleFile(stderr, " %lf, %lf, %lf, %lf'\n",
2127 coeff[7], coeff[4], coeff[5], coeff[6]);
2128#endif
2129 (void) FormatLocaleFile(stderr,
2130 "BilinearReverse Distort, FX Equivalent:\n");
2131 (void) FormatLocaleFile(stderr,"%s", image_gen);
2132 (void) FormatLocaleFile(stderr,
2133 " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2134 (void) FormatLocaleFile(stderr,
2135 " xx=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",coeff[0],coeff[1],
2136 coeff[2], coeff[3]);
2137 (void) FormatLocaleFile(stderr,
2138 " yy=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",coeff[4],coeff[5],
2139 coeff[6], coeff[7]);
2140 (void) FormatLocaleFile(stderr," %s' \\\n", lookup);
2141 break;
2142 }
2143 case PolynomialDistortion:
2144 {
2145 size_t nterms = CastDoubleToSizeT(coeff[1]);
2146 (void) FormatLocaleFile(stderr,
2147 "Polynomial (order %lg, terms %lu), FX Equivalent\n",coeff[0],
2148 (unsigned long) nterms);
2149 (void) FormatLocaleFile(stderr,"%s", image_gen);
2150 (void) FormatLocaleFile(stderr,
2151 " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2152 (void) FormatLocaleFile(stderr, " xx =");
2153 for (i=0; i < (ssize_t) nterms; i++)
2154 {
2155 if ((i != 0) && (i%4 == 0))
2156 (void) FormatLocaleFile(stderr, "\n ");
2157 (void) FormatLocaleFile(stderr," %+lf%s",coeff[2+i],
2158 poly_basis_str(i));
2159 }
2160 (void) FormatLocaleFile(stderr,";\n yy =");
2161 for (i=0; i < (ssize_t) nterms; i++)
2162 {
2163 if ((i != 0) && (i%4 == 0))
2164 (void) FormatLocaleFile(stderr,"\n ");
2165 (void) FormatLocaleFile(stderr," %+lf%s",coeff[2+i+nterms],
2166 poly_basis_str(i));
2167 }
2168 (void) FormatLocaleFile(stderr,";\n %s' \\\n", lookup);
2169 break;
2170 }
2171 case ArcDistortion:
2172 {
2173 (void) FormatLocaleFile(stderr,"Arc Distort, Internal Coefficients:\n");
2174 for (i=0; i < 5; i++)
2175 (void) FormatLocaleFile(stderr,
2176 " c%.20g = %+lf\n",(double) i,coeff[i]);
2177 (void) FormatLocaleFile(stderr,"Arc Distort, FX Equivalent:\n");
2178 (void) FormatLocaleFile(stderr,"%s", image_gen);
2179 (void) FormatLocaleFile(stderr," -fx 'ii=i+page.x; jj=j+page.y;\n");
2180 (void) FormatLocaleFile(stderr," xx=(atan2(jj,ii)%+lf)/(2*pi);\n",
2181 -coeff[0]);
2182 (void) FormatLocaleFile(stderr," xx=xx-round(xx);\n");
2183 (void) FormatLocaleFile(stderr," xx=xx*%lf %+lf;\n",coeff[1],
2184 coeff[4]);
2185 (void) FormatLocaleFile(stderr,
2186 " yy=(%lf - hypot(ii,jj)) * %lf;\n",coeff[2],coeff[3]);
2187 (void) FormatLocaleFile(stderr," v.p{xx-.5,yy-.5}' \\\n");
2188 break;
2189 }
2190 case PolarDistortion:
2191 {
2192 (void) FormatLocaleFile(stderr,"Polar Distort, Internal Coefficients\n");
2193 for (i=0; i < 8; i++)
2194 (void) FormatLocaleFile(stderr," c%.20g = %+lf\n",(double) i,
2195 coeff[i]);
2196 (void) FormatLocaleFile(stderr,"Polar Distort, FX Equivalent:\n");
2197 (void) FormatLocaleFile(stderr,"%s", image_gen);
2198 (void) FormatLocaleFile(stderr,
2199 " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",-coeff[2],-coeff[3]);
2200 (void) FormatLocaleFile(stderr," xx=(atan2(ii,jj)%+lf)/(2*pi);\n",
2201 -(coeff[4]+coeff[5])/2 );
2202 (void) FormatLocaleFile(stderr," xx=xx-round(xx);\n");
2203 (void) FormatLocaleFile(stderr," xx=xx*2*pi*%lf + v.w/2;\n",
2204 coeff[6] );
2205 (void) FormatLocaleFile(stderr," yy=(hypot(ii,jj)%+lf)*%lf;\n",
2206 -coeff[1],coeff[7] );
2207 (void) FormatLocaleFile(stderr," v.p{xx-.5,yy-.5}' \\\n");
2208 break;
2209 }
2210 case DePolarDistortion:
2211 {
2212 (void) FormatLocaleFile(stderr,
2213 "DePolar Distort, Internal Coefficients\n");
2214 for (i=0; i < 8; i++)
2215 (void) FormatLocaleFile(stderr," c%.20g = %+lf\n",(double) i,
2216 coeff[i]);
2217 (void) FormatLocaleFile(stderr,"DePolar Distort, FX Equivalent:\n");
2218 (void) FormatLocaleFile(stderr,"%s", image_gen);
2219 (void) FormatLocaleFile(stderr," -fx 'aa=(i+.5)*%lf %+lf;\n",
2220 coeff[6],+coeff[4]);
2221 (void) FormatLocaleFile(stderr," rr=(j+.5)*%lf %+lf;\n",
2222 coeff[7],+coeff[1]);
2223 (void) FormatLocaleFile(stderr," xx=rr*sin(aa) %+lf;\n",
2224 coeff[2]);
2225 (void) FormatLocaleFile(stderr," yy=rr*cos(aa) %+lf;\n",
2226 coeff[3]);
2227 (void) FormatLocaleFile(stderr," v.p{xx-.5,yy-.5}' \\\n");
2228 break;
2229 }
2230 case Cylinder2PlaneDistortion:
2231 {
2232 (void) FormatLocaleFile(stderr,
2233 "Cylinder to Plane Distort, Internal Coefficients\n");
2234 (void) FormatLocaleFile(stderr," cylinder_radius = %+lf\n",coeff[1]);
2235 (void) FormatLocaleFile(stderr,
2236 "Cylinder to Plane Distort, FX Equivalent:\n");
2237 (void) FormatLocaleFile(stderr, "%s", image_gen);
2238 (void) FormatLocaleFile(stderr,
2239 " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",-coeff[4],
2240 -coeff[5]);
2241 (void) FormatLocaleFile(stderr," aa=atan(ii/%+lf);\n",coeff[1]);
2242 (void) FormatLocaleFile(stderr," xx=%lf*aa%+lf;\n",
2243 coeff[1],coeff[2]);
2244 (void) FormatLocaleFile(stderr," yy=jj*cos(aa)%+lf;\n",coeff[3]);
2245 (void) FormatLocaleFile(stderr," %s' \\\n", lookup);
2246 break;
2247 }
2248 case Plane2CylinderDistortion:
2249 {
2250 (void) FormatLocaleFile(stderr,
2251 "Plane to Cylinder Distort, Internal Coefficients\n");
2252 (void) FormatLocaleFile(stderr," cylinder_radius = %+lf\n",coeff[1]);
2253 (void) FormatLocaleFile(stderr,
2254 "Plane to Cylinder Distort, FX Equivalent:\n");
2255 (void) FormatLocaleFile(stderr,"%s", image_gen);
2256 (void) FormatLocaleFile(stderr,
2257 " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",-coeff[4],
2258 -coeff[5]);
2259 (void) FormatLocaleFile(stderr," ii=ii/%+lf;\n",coeff[1]);
2260 (void) FormatLocaleFile(stderr," xx=%lf*tan(ii)%+lf;\n",coeff[1],
2261 coeff[2] );
2262 (void) FormatLocaleFile(stderr," yy=jj/cos(ii)%+lf;\n",coeff[3]);
2263 (void) FormatLocaleFile(stderr," %s' \\\n", lookup);
2264 break;
2265 }
2266 case BarrelDistortion:
2267 case BarrelInverseDistortion:
2268 {
2269 double
2270 xc,
2271 yc;
2272
2273 /*
2274 NOTE: This does the barrel roll in pixel coords not image coords
2275 The internal distortion must do it in image coordinates, so that is
2276 what the center coeff (8,9) is given in.
2277 */
2278 xc=((double)image->columns-1.0)/2.0+image->page.x;
2279 yc=((double)image->rows-1.0)/2.0+image->page.y;
2280 (void) FormatLocaleFile(stderr, "Barrel%s Distort, FX Equivalent:\n",
2281 method == BarrelDistortion ? "" : "Inv");
2282 (void) FormatLocaleFile(stderr, "%s", image_gen);
2283 if ( fabs(coeff[8]-xc-0.5) < 0.1 && fabs(coeff[9]-yc-0.5) < 0.1 )
2284 (void) FormatLocaleFile(stderr," -fx 'xc=(w-1)/2; yc=(h-1)/2;\n");
2285 else
2286 (void) FormatLocaleFile(stderr," -fx 'xc=%lf; yc=%lf;\n",coeff[8]-
2287 0.5,coeff[9]-0.5);
2288 (void) FormatLocaleFile(stderr,
2289 " ii=i-xc; jj=j-yc; rr=hypot(ii,jj);\n");
2290 (void) FormatLocaleFile(stderr,
2291 " ii=ii%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2292 method == BarrelDistortion ? "*" : "/",coeff[0],coeff[1],coeff[2],
2293 coeff[3]);
2294 (void) FormatLocaleFile(stderr,
2295 " jj=jj%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2296 method == BarrelDistortion ? "*" : "/",coeff[4],coeff[5],coeff[6],
2297 coeff[7]);
2298 (void) FormatLocaleFile(stderr," p{ii+xc,jj+yc}' \\\n");
2299 }
2300 default:
2301 break;
2302 }
2303 }
2304 /*
2305 The user provided a 'scale' expert option will scale the output image size,
2306 by the factor given allowing for super-sampling of the distorted image
2307 space. Any scaling factors must naturally be halved as a result.
2308 */
2309 { const char *artifact;
2310 artifact=GetImageArtifact(image,"distort:scale");
2311 output_scaling = 1.0;
2312 if (artifact != (const char *) NULL) {
2313 output_scaling = fabs(StringToDouble(artifact,(char **) NULL));
2314 geometry.width=CastDoubleToSizeT(output_scaling*geometry.width+0.5);
2315 geometry.height=CastDoubleToSizeT(output_scaling*geometry.height+0.5);
2316 geometry.x=CastDoubleToSsizeT(output_scaling*geometry.x+0.5);
2317 geometry.y=CastDoubleToSsizeT(output_scaling*geometry.y+0.5);
2318 if ( output_scaling < 0.1 ) {
2319 coeff = (double *) RelinquishMagickMemory(coeff);
2320 (void) ThrowMagickException(exception,GetMagickModule(),
2321 OptionError,"InvalidArgument","%s","-define distort:scale" );
2322 return((Image *) NULL);
2323 }
2324 output_scaling = 1/output_scaling;
2325 }
2326 }
2327#define ScaleFilter(F,A,B,C,D) \
2328 ScaleResampleFilter( (F), \
2329 output_scaling*(A), output_scaling*(B), \
2330 output_scaling*(C), output_scaling*(D) )
2331
2332 /*
2333 Initialize the distort image attributes.
2334 */
2335 distort_image=CloneImage(image,geometry.width,geometry.height,MagickTrue,
2336 exception);
2337 if (distort_image == (Image *) NULL)
2338 {
2339 coeff=(double *) RelinquishMagickMemory(coeff);
2340 return((Image *) NULL);
2341 }
2342 /* if image is ColorMapped - change it to DirectClass */
2343 if (SetImageStorageClass(distort_image,DirectClass) == MagickFalse)
2344 {
2345 coeff=(double *) RelinquishMagickMemory(coeff);
2346 InheritException(exception,&distort_image->exception);
2347 distort_image=DestroyImage(distort_image);
2348 return((Image *) NULL);
2349 }
2350 if ((IsPixelGray(&distort_image->background_color) == MagickFalse) &&
2351 (IsGrayColorspace(distort_image->colorspace) != MagickFalse))
2352 (void) SetImageColorspace(distort_image,sRGBColorspace);
2353 if (distort_image->background_color.opacity != OpaqueOpacity)
2354 distort_image->matte=MagickTrue;
2355 distort_image->page.x=geometry.x;
2356 distort_image->page.y=geometry.y;
2357
2358 { /* ----- MAIN CODE -----
2359 Sample the source image to each pixel in the distort image.
2360 */
2361 CacheView
2362 *distort_view;
2363
2364 MagickBooleanType
2365 status;
2366
2367 MagickOffsetType
2368 progress;
2369
2370 MagickPixelPacket
2371 zero;
2372
2373 ResampleFilter
2374 **magick_restrict resample_filter;
2375
2376 ssize_t
2377 j;
2378
2379 status=MagickTrue;
2380 progress=0;
2381 GetMagickPixelPacket(distort_image,&zero);
2382 resample_filter=AcquireResampleFilterTLS(image,UndefinedVirtualPixelMethod,
2383 MagickFalse,exception);
2384 distort_view=AcquireAuthenticCacheView(distort_image,exception);
2385#if defined(MAGICKCORE_OPENMP_SUPPORT)
2386 #pragma omp parallel for schedule(static) shared(progress,status) \
2387 magick_number_threads(image,distort_image,distort_image->rows,1)
2388#endif
2389 for (j=0; j < (ssize_t) distort_image->rows; j++)
2390 {
2391 const int
2392 id = GetOpenMPThreadId();
2393
2394 double
2395 validity; /* how mathematically valid is this the mapping */
2396
2397 MagickBooleanType
2398 sync;
2399
2400 MagickPixelPacket
2401 pixel, /* pixel color to assign to distorted image */
2402 invalid; /* the color to assign when distort result is invalid */
2403
2404 PointInfo
2405 d,
2406 s; /* transform destination image x,y to source image x,y */
2407
2408 IndexPacket
2409 *magick_restrict indexes;
2410
2411 ssize_t
2412 i;
2413
2414 PixelPacket
2415 *magick_restrict q;
2416
2417 q=QueueCacheViewAuthenticPixels(distort_view,0,j,distort_image->columns,1,
2418 exception);
2419 if (q == (PixelPacket *) NULL)
2420 {
2421 status=MagickFalse;
2422 continue;
2423 }
2424 indexes=GetCacheViewAuthenticIndexQueue(distort_view);
2425 pixel=zero;
2426
2427 /* Define constant scaling vectors for Affine Distortions
2428 Other methods are either variable, or use interpolated lookup
2429 */
2430 switch (method)
2431 {
2432 case AffineDistortion:
2433 ScaleFilter( resample_filter[id],
2434 coeff[0], coeff[1],
2435 coeff[3], coeff[4] );
2436 break;
2437 default:
2438 break;
2439 }
2440
2441 /* Initialize default pixel validity
2442 * negative: pixel is invalid output 'matte_color'
2443 * 0.0 to 1.0: antialiased, mix with resample output
2444 * 1.0 or greater: use resampled output.
2445 */
2446 validity = 1.0;
2447
2448 GetMagickPixelPacket(distort_image,&invalid);
2449 SetMagickPixelPacket(distort_image,&distort_image->matte_color,
2450 (IndexPacket *) NULL, &invalid);
2451 if (distort_image->colorspace == CMYKColorspace)
2452 ConvertRGBToCMYK(&invalid); /* what about other color spaces? */
2453
2454 for (i=0; i < (ssize_t) distort_image->columns; i++)
2455 {
2456 /* map pixel coordinate to distortion space coordinate */
2457 d.x = (double) (geometry.x+i+0.5)*output_scaling;
2458 d.y = (double) (geometry.y+j+0.5)*output_scaling;
2459 s = d; /* default is a no-op mapping */
2460 switch (method)
2461 {
2462 case AffineDistortion:
2463 {
2464 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2465 s.y=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2466 /* Affine partial derivatives are constant -- set above */
2467 break;
2468 }
2469 case PerspectiveDistortion:
2470 {
2471 double
2472 p,q,r,abs_r,abs_c6,abs_c7,scale;
2473 /* perspective is a ratio of affines */
2474 p=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2475 q=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2476 r=coeff[6]*d.x+coeff[7]*d.y+1.0;
2477 /* Pixel Validity -- is it a 'sky' or 'ground' pixel */
2478 validity = (r*coeff[8] < 0.0) ? 0.0 : 1.0;
2479 /* Determine horizon anti-alias blending */
2480 abs_r = fabs(r)*2;
2481 abs_c6 = fabs(coeff[6]);
2482 abs_c7 = fabs(coeff[7]);
2483 if ( abs_c6 > abs_c7 ) {
2484 if ( abs_r < abs_c6*output_scaling )
2485 validity = 0.5 - coeff[8]*r/(coeff[6]*output_scaling);
2486 }
2487 else if ( abs_r < abs_c7*output_scaling )
2488 validity = 0.5 - coeff[8]*r/(coeff[7]*output_scaling);
2489 /* Perspective Sampling Point (if valid) */
2490 if ( validity > 0.0 ) {
2491 /* divide by r affine, for perspective scaling */
2492 scale = 1.0/r;
2493 s.x = p*scale;
2494 s.y = q*scale;
2495 /* Perspective Partial Derivatives or Scaling Vectors */
2496 scale *= scale;
2497 ScaleFilter( resample_filter[id],
2498 (r*coeff[0] - p*coeff[6])*scale,
2499 (r*coeff[1] - p*coeff[7])*scale,
2500 (r*coeff[3] - q*coeff[6])*scale,
2501 (r*coeff[4] - q*coeff[7])*scale );
2502 }
2503 break;
2504 }
2505 case BilinearReverseDistortion:
2506 {
2507 /* Reversed Mapped is just a simple polynomial */
2508 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2]*d.x*d.y+coeff[3];
2509 s.y=coeff[4]*d.x+coeff[5]*d.y
2510 +coeff[6]*d.x*d.y+coeff[7];
2511 /* Bilinear partial derivatives of scaling vectors */
2512 ScaleFilter( resample_filter[id],
2513 coeff[0] + coeff[2]*d.y,
2514 coeff[1] + coeff[2]*d.x,
2515 coeff[4] + coeff[6]*d.y,
2516 coeff[5] + coeff[6]*d.x );
2517 break;
2518 }
2519 case BilinearForwardDistortion:
2520 {
2521 /* Forward mapped needs reversed polynomial equations
2522 * which unfortunately requires a square root! */
2523 double b,c;
2524 d.x -= coeff[3]; d.y -= coeff[7];
2525 b = coeff[6]*d.x - coeff[2]*d.y + coeff[8];
2526 c = coeff[4]*d.x - coeff[0]*d.y;
2527
2528 validity = 1.0;
2529 /* Handle Special degenerate (non-quadratic) case
2530 * Currently without horizon anti-aliasing */
2531 if ( fabs(coeff[9]) < MagickEpsilon )
2532 s.y = -c/b;
2533 else {
2534 c = b*b - 2*coeff[9]*c;
2535 if ( c < 0.0 )
2536 validity = 0.0;
2537 else
2538 s.y = ( -b + sqrt(c) )/coeff[9];
2539 }
2540 if ( validity > 0.0 )
2541 s.x = ( d.x - coeff[1]*s.y) / ( coeff[0] + coeff[2]*s.y );
2542
2543 /* NOTE: the sign of the square root should be -ve for parts
2544 where the source image becomes 'flipped' or 'mirrored'.
2545 FUTURE: Horizon handling
2546 FUTURE: Scaling factors or Derivatives (how?)
2547 */
2548 break;
2549 }
2550#if 0
2551 case BilinearDistortion:
2552 /* Bilinear mapping of any Quadrilateral to any Quadrilateral */
2553 /* UNDER DEVELOPMENT */
2554 break;
2555#endif
2556 case PolynomialDistortion:
2557 {
2558 /* multi-ordered polynomial */
2559 ssize_t
2560 k;
2561
2562 ssize_t
2563 nterms=(ssize_t)coeff[1];
2564
2565 PointInfo
2566 du,dv; /* the du,dv vectors from unit dx,dy -- derivatives */
2567
2568 s.x=s.y=du.x=du.y=dv.x=dv.y=0.0;
2569 for(k=0; k < nterms; k++) {
2570 s.x += poly_basis_fn(k,d.x,d.y)*coeff[2+k];
2571 du.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k];
2572 du.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k];
2573 s.y += poly_basis_fn(k,d.x,d.y)*coeff[2+k+nterms];
2574 dv.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k+nterms];
2575 dv.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k+nterms];
2576 }
2577 ScaleFilter( resample_filter[id], du.x,du.y,dv.x,dv.y );
2578 break;
2579 }
2580 case ArcDistortion:
2581 {
2582 /* what is the angle and radius in the destination image */
2583 s.x = (double) ((atan2(d.y,d.x) - coeff[0])/Magick2PI);
2584 s.x -= MagickRound(s.x); /* angle */
2585 s.y = hypot(d.x,d.y); /* radius */
2586
2587 /* Arc Distortion Partial Scaling Vectors
2588 Are derived by mapping the perpendicular unit vectors
2589 dR and dA*R*2PI rather than trying to map dx and dy
2590 The results is a very simple orthogonal aligned ellipse.
2591 */
2592 if ( s.y > MagickEpsilon )
2593 ScaleFilter( resample_filter[id],
2594 (double) (coeff[1]/(Magick2PI*s.y)), 0, 0, coeff[3] );
2595 else
2596 ScaleFilter( resample_filter[id],
2597 distort_image->columns*2, 0, 0, coeff[3] );
2598
2599 /* now scale the angle and radius for source image lookup point */
2600 s.x = s.x*coeff[1] + coeff[4] + image->page.x +0.5;
2601 s.y = (coeff[2] - s.y) * coeff[3] + image->page.y;
2602 break;
2603 }
2604 case PolarDistortion:
2605 { /* 2D Cartesian to Polar View */
2606 d.x -= coeff[2];
2607 d.y -= coeff[3];
2608 s.x = atan2(d.x,d.y) - (coeff[4]+coeff[5])/2;
2609 s.x /= Magick2PI;
2610 s.x -= MagickRound(s.x);
2611 s.x *= Magick2PI; /* angle - relative to centerline */
2612 s.y = hypot(d.x,d.y); /* radius */
2613
2614 /* Polar Scaling vectors are based on mapping dR and dA vectors
2615 This results in very simple orthogonal scaling vectors
2616 */
2617 if ( s.y > MagickEpsilon )
2618 ScaleFilter( resample_filter[id],
2619 (double) (coeff[6]/(Magick2PI*s.y)), 0, 0, coeff[7] );
2620 else
2621 ScaleFilter( resample_filter[id],
2622 distort_image->columns*2, 0, 0, coeff[7] );
2623
2624 /* now finish mapping radius/angle to source x,y coords */
2625 s.x = s.x*coeff[6] + (double)image->columns/2.0 + image->page.x;
2626 s.y = (s.y-coeff[1])*coeff[7] + image->page.y;
2627 break;
2628 }
2629 case DePolarDistortion:
2630 { /* @D Polar to Cartesian */
2631 /* ignore all destination virtual offsets */
2632 d.x = ((double)i+0.5)*output_scaling*coeff[6]+coeff[4];
2633 d.y = ((double)j+0.5)*output_scaling*coeff[7]+coeff[1];
2634 s.x = d.y*sin(d.x) + coeff[2];
2635 s.y = d.y*cos(d.x) + coeff[3];
2636 /* derivatives are useless - better to use SuperSampling */
2637 break;
2638 }
2639 case Cylinder2PlaneDistortion:
2640 { /* 3D Cylinder to Tangential Plane */
2641 double ax, cx;
2642 /* relative to center of distortion */
2643 d.x -= coeff[4]; d.y -= coeff[5];
2644 d.x /= coeff[1]; /* x' = x/r */
2645 ax=atan(d.x); /* aa = atan(x/r) = u/r */
2646 cx=cos(ax); /* cx = cos(atan(x/r)) = 1/sqrt(x^2+u^2) */
2647 s.x = coeff[1]*ax; /* u = r*atan(x/r) */
2648 s.y = d.y*cx; /* v = y*cos(u/r) */
2649 /* derivatives... (see personal notes) */
2650 ScaleFilter( resample_filter[id],
2651 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2652#if 0
2653if ( i == 0 && j == 0 ) {
2654 fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2655 fprintf(stderr, "phi = %lf\n", (double)(ax * 180.0/MagickPI) );
2656 fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n",
2657 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2658 fflush(stderr); }
2659#endif
2660 /* add center of distortion in source */
2661 s.x += coeff[2]; s.y += coeff[3];
2662 break;
2663 }
2664 case Plane2CylinderDistortion:
2665 { /* 3D Cylinder to Tangential Plane */
2666 /* relative to center of distortion */
2667 d.x -= coeff[4]; d.y -= coeff[5];
2668
2669 /* is pixel valid - horizon of a infinite Virtual-Pixel Plane
2670 * (see Anthony Thyssen's personal note) */
2671 validity = (double) ((coeff[1]*MagickPI2 - fabs(d.x))/output_scaling + 0.5);
2672
2673 if ( validity > 0.0 ) {
2674 double cx,tx;
2675 d.x /= coeff[1]; /* x'= x/r */
2676 cx = 1/cos(d.x); /* cx = 1/cos(x/r) */
2677 tx = tan(d.x); /* tx = tan(x/r) */
2678 s.x = coeff[1]*tx; /* u = r * tan(x/r) */
2679 s.y = d.y*cx; /* v = y / cos(x/r) */
2680 /* derivatives... (see Anthony Thyssen's personal notes) */
2681 ScaleFilter( resample_filter[id],
2682 cx*cx, 0.0, s.y*cx/coeff[1], cx );
2683#if 0
2684/*if ( i == 0 && j == 0 ) {*/
2685if ( d.x == 0.5 && d.y == 0.5 ) {
2686 fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2687 fprintf(stderr, "radius = %lf phi = %lf validity = %lf\n",
2688 coeff[1], (double)(d.x * 180.0/MagickPI), validity );
2689 fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n",
2690 cx*cx, 0.0, s.y*cx/coeff[1], cx);
2691 fflush(stderr); }
2692#endif
2693 }
2694 /* add center of distortion in source */
2695 s.x += coeff[2]; s.y += coeff[3];
2696 break;
2697 }
2698 case BarrelDistortion:
2699 case BarrelInverseDistortion:
2700 { /* Lens Barrel Distortion Correction */
2701 double r,fx,fy,gx,gy;
2702 /* Radial Polynomial Distortion (de-normalized) */
2703 d.x -= coeff[8];
2704 d.y -= coeff[9];
2705 r = sqrt(d.x*d.x+d.y*d.y);
2706 if ( r > MagickEpsilon ) {
2707 fx = ((coeff[0]*r + coeff[1])*r + coeff[2])*r + coeff[3];
2708 fy = ((coeff[4]*r + coeff[5])*r + coeff[6])*r + coeff[7];
2709 gx = ((3*coeff[0]*r + 2*coeff[1])*r + coeff[2])/r;
2710 gy = ((3*coeff[4]*r + 2*coeff[5])*r + coeff[6])/r;
2711 /* adjust functions and scaling for 'inverse' form */
2712 if ( method == BarrelInverseDistortion ) {
2713 fx = 1/fx; fy = 1/fy;
2714 gx *= -fx*fx; gy *= -fy*fy;
2715 }
2716 /* Set the source pixel to lookup and EWA derivative vectors */
2717 s.x = d.x*fx + coeff[8];
2718 s.y = d.y*fy + coeff[9];
2719 ScaleFilter( resample_filter[id],
2720 gx*d.x*d.x + fx, gx*d.x*d.y,
2721 gy*d.x*d.y, gy*d.y*d.y + fy );
2722 }
2723 else {
2724 /* Special handling to avoid divide by zero when r==0
2725 **
2726 ** The source and destination pixels match in this case
2727 ** which was set at the top of the loop using s = d;
2728 ** otherwise... s.x=coeff[8]; s.y=coeff[9];
2729 */
2730 if ( method == BarrelDistortion )
2731 ScaleFilter( resample_filter[id],
2732 coeff[3], 0, 0, coeff[7] );
2733 else /* method == BarrelInverseDistortion */
2734 /* FUTURE, trap for D==0 causing division by zero */
2735 ScaleFilter( resample_filter[id],
2736 1.0/coeff[3], 0, 0, 1.0/coeff[7] );
2737 }
2738 break;
2739 }
2740 case ShepardsDistortion:
2741 { /* Shepards Method, or Inverse Weighted Distance for
2742 displacement around the destination image control points
2743 The input arguments are the coefficients to the function.
2744 This is more of a 'displacement' function rather than an
2745 absolute distortion function.
2746
2747 Note: We can not determine derivatives using shepards method
2748 so only a point sample interpolation can be used.
2749 */
2750 size_t
2751 i;
2752 double
2753 denominator;
2754
2755 denominator = s.x = s.y = 0;
2756 for(i=0; i<number_arguments; i+=4) {
2757 double weight =
2758 ((double)d.x-arguments[i+2])*((double)d.x-arguments[i+2])
2759 + ((double)d.y-arguments[i+3])*((double)d.y-arguments[i+3]);
2760 weight = pow(weight,coeff[0]); /* shepards power factor */
2761 weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
2762
2763 s.x += (arguments[ i ]-arguments[i+2])*weight;
2764 s.y += (arguments[i+1]-arguments[i+3])*weight;
2765 denominator += weight;
2766 }
2767 s.x /= denominator;
2768 s.y /= denominator;
2769 s.x += d.x; /* make it as relative displacement */
2770 s.y += d.y;
2771 break;
2772 }
2773 default:
2774 break; /* use the default no-op given above */
2775 }
2776 /* map virtual canvas location back to real image coordinate */
2777 if ( bestfit && method != ArcDistortion ) {
2778 s.x -= image->page.x;
2779 s.y -= image->page.y;
2780 }
2781 s.x -= 0.5;
2782 s.y -= 0.5;
2783
2784 if ( validity <= 0.0 ) {
2785 /* result of distortion is an invalid pixel - don't resample */
2786 SetPixelPacket(distort_image,&invalid,q,indexes);
2787 }
2788 else {
2789 /* resample the source image to find its correct color */
2790 status=ResamplePixelColor(resample_filter[id],s.x,s.y,&pixel);
2791 if (status == MagickFalse)
2792 SetPixelPacket(distort_image,&invalid,q,indexes);
2793 else
2794 {
2795 /* if validity between 0.0 & 1.0 mix result with invalid pixel */
2796 if ( validity < 1.0 ) {
2797 /* Do a blend of sample color and invalid pixel */
2798 /* should this be a 'Blend', or an 'Over' compose */
2799 MagickPixelCompositeBlend(&pixel,validity,&invalid,(1.0-validity),
2800 &pixel);
2801 }
2802 }
2803 SetPixelPacket(distort_image,&pixel,q,indexes);
2804 }
2805 q++;
2806 indexes++;
2807 }
2808 sync=SyncCacheViewAuthenticPixels(distort_view,exception);
2809 if (sync == MagickFalse)
2810 status=MagickFalse;
2811 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2812 {
2813 MagickBooleanType
2814 proceed;
2815
2816#if defined(MAGICKCORE_OPENMP_SUPPORT)
2817 #pragma omp atomic
2818#endif
2819 progress++;
2820 proceed=SetImageProgress(image,DistortImageTag,progress,image->rows);
2821 if (proceed == MagickFalse)
2822 status=MagickFalse;
2823 }
2824 }
2825 distort_view=DestroyCacheView(distort_view);
2826 resample_filter=DestroyResampleFilterTLS(resample_filter);
2827
2828 if (status == MagickFalse)
2829 distort_image=DestroyImage(distort_image);
2830 }
2831
2832 /* Arc does not return an offset unless 'bestfit' is in effect
2833 And the user has not provided an overriding 'viewport'.
2834 */
2835 if ( method == ArcDistortion && !bestfit && !viewport_given ) {
2836 distort_image->page.x = 0;
2837 distort_image->page.y = 0;
2838 }
2839 coeff=(double *) RelinquishMagickMemory(coeff);
2840 return(distort_image);
2841}
2842␌
2843/*
2844%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2845% %
2846% %
2847% %
2848% R o t a t e I m a g e %
2849% %
2850% %
2851% %
2852%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2853%
2854% RotateImage() creates a new image that is a rotated copy of an existing
2855% one. Positive angles rotate counter-clockwise (right-hand rule), while
2856% negative angles rotate clockwise. Rotated images are usually larger than
2857% the originals and have 'empty' triangular corners. X axis. Empty
2858% triangles left over from shearing the image are filled with the background
2859% color defined by member 'background_color' of the image. RotateImage
2860% allocates the memory necessary for the new Image structure and returns a
2861% pointer to the new image.
2862%
2863% The format of the RotateImage method is:
2864%
2865% Image *RotateImage(const Image *image,const double degrees,
2866% ExceptionInfo *exception)
2867%
2868% A description of each parameter follows.
2869%
2870% o image: the image.
2871%
2872% o degrees: Specifies the number of degrees to rotate the image.
2873%
2874% o exception: return any errors or warnings in this structure.
2875%
2876*/
2877MagickExport Image *RotateImage(const Image *image,const double degrees,
2878 ExceptionInfo *exception)
2879{
2880 Image
2881 *distort_image,
2882 *rotate_image;
2883
2884 MagickRealType
2885 angle;
2886
2887 PointInfo
2888 shear;
2889
2890 size_t
2891 rotations;
2892
2893 /*
2894 Adjust rotation angle.
2895 */
2896 assert(image != (Image *) NULL);
2897 assert(image->signature == MagickCoreSignature);
2898 assert(exception != (ExceptionInfo *) NULL);
2899 assert(exception->signature == MagickCoreSignature);
2900 if (IsEventLogging() != MagickFalse)
2901 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2902 angle=fmod(degrees,360.0);
2903 while (angle < -45.0)
2904 angle+=360.0;
2905 for (rotations=0; angle > 45.0; rotations++)
2906 angle-=90.0;
2907 rotations%=4;
2908 shear.x=(-tan((double) DegreesToRadians(angle)/2.0));
2909 shear.y=sin((double) DegreesToRadians(angle));
2910 if ((fabs(shear.x) < MagickEpsilon) && (fabs(shear.y) < MagickEpsilon))
2911 return(IntegralRotateImage(image,rotations,exception));
2912 distort_image=CloneImage(image,0,0,MagickTrue,exception);
2913 if (distort_image == (Image *) NULL)
2914 return((Image *) NULL);
2915 (void) SetImageVirtualPixelMethod(distort_image,BackgroundVirtualPixelMethod);
2916 rotate_image=DistortImage(distort_image,ScaleRotateTranslateDistortion,1,
2917 &degrees,MagickTrue,exception);
2918 distort_image=DestroyImage(distort_image);
2919 return(rotate_image);
2920}
2921␌
2922/*
2923%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2924% %
2925% %
2926% %
2927% S p a r s e C o l o r I m a g e %
2928% %
2929% %
2930% %
2931%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2932%
2933% SparseColorImage(), given a set of coordinates, interpolates the colors
2934% found at those coordinates, across the whole image, using various methods.
2935%
2936% The format of the SparseColorImage() method is:
2937%
2938% Image *SparseColorImage(const Image *image,const ChannelType channel,
2939% const SparseColorMethod method,const size_t number_arguments,
2940% const double *arguments,ExceptionInfo *exception)
2941%
2942% A description of each parameter follows:
2943%
2944% o image: the image to be filled in.
2945%
2946% o channel: Specify which color values (in RGBKA sequence) are being set.
2947% This also determines the number of color_values in above.
2948%
2949% o method: the method to fill in the gradient between the control points.
2950%
2951% The methods used for SparseColor() are often simular to methods
2952% used for DistortImage(), and even share the same code for determination
2953% of the function coefficients, though with more dimensions (or resulting
2954% values).
2955%
2956% o number_arguments: the number of arguments given.
2957%
2958% o arguments: array of floating point arguments for this method--
2959% x,y,color_values-- with color_values given as normalized values.
2960%
2961% o exception: return any errors or warnings in this structure
2962%
2963*/
2964MagickExport Image *SparseColorImage(const Image *image,
2965 const ChannelType channel,const SparseColorMethod method,
2966 const size_t number_arguments,const double *arguments,
2967 ExceptionInfo *exception)
2968{
2969#define SparseColorTag "Distort/SparseColor"
2970
2971 SparseColorMethod
2972 sparse_method;
2973
2974 double
2975 *coeff;
2976
2977 Image
2978 *sparse_image;
2979
2980 size_t
2981 number_colors;
2982
2983 assert(image != (Image *) NULL);
2984 assert(image->signature == MagickCoreSignature);
2985 assert(exception != (ExceptionInfo *) NULL);
2986 assert(exception->signature == MagickCoreSignature);
2987 if (IsEventLogging() != MagickFalse)
2988 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2989 /*
2990 Determine number of color values needed per control point.
2991 */
2992 number_colors=0;
2993 if ( channel & RedChannel ) number_colors++;
2994 if ( channel & GreenChannel ) number_colors++;
2995 if ( channel & BlueChannel ) number_colors++;
2996 if ( channel & IndexChannel ) number_colors++;
2997 if ( channel & OpacityChannel ) number_colors++;
2998
2999 /*
3000 Convert input arguments into mapping coefficients, in this case
3001 we are mapping (distorting) colors, rather than coordinates.
3002 */
3003 { DistortImageMethod
3004 distort_method;
3005
3006 distort_method=(DistortImageMethod) method;
3007 if ( distort_method >= SentinelDistortion )
3008 distort_method = ShepardsDistortion; /* Pretend to be Shepards */
3009 coeff = GenerateCoefficients(image, &distort_method, number_arguments,
3010 arguments, number_colors, exception);
3011 if ( coeff == (double *) NULL )
3012 return((Image *) NULL);
3013 /*
3014 Note some Distort Methods may fall back to other simpler methods,
3015 Currently the only fallback of concern is Bilinear to Affine
3016 (Barycentric), which is also sparse_colr method. This also ensures
3017 correct two and one color Barycentric handling.
3018 */
3019 sparse_method = (SparseColorMethod) distort_method;
3020 if ( distort_method == ShepardsDistortion )
3021 sparse_method = method; /* return non-distort methods to normal */
3022 if ( sparse_method == InverseColorInterpolate )
3023 coeff[0]=0.5; /* sqrt() the squared distance for inverse */
3024 }
3025
3026 /* Verbose output */
3027 if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
3028
3029 switch (sparse_method) {
3030 case BarycentricColorInterpolate:
3031 {
3032 ssize_t x=0;
3033 (void) FormatLocaleFile(stderr, "Barycentric Sparse Color:\n");
3034 if ( channel & RedChannel )
3035 (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf' \\\n",
3036 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
3037 if ( channel & GreenChannel )
3038 (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf' \\\n",
3039 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
3040 if ( channel & BlueChannel )
3041 (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf' \\\n",
3042 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
3043 if ( channel & IndexChannel )
3044 (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf' \\\n",
3045 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
3046 if ( channel & OpacityChannel )
3047 (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf' \\\n",
3048 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
3049 break;
3050 }
3051 case BilinearColorInterpolate:
3052 {
3053 ssize_t x=0;
3054 (void) FormatLocaleFile(stderr, "Bilinear Sparse Color\n");
3055 if ( channel & RedChannel )
3056 (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3057 coeff[ x ], coeff[x+1],
3058 coeff[x+2], coeff[x+3]),x+=4;
3059 if ( channel & GreenChannel )
3060 (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3061 coeff[ x ], coeff[x+1],
3062 coeff[x+2], coeff[x+3]),x+=4;
3063 if ( channel & BlueChannel )
3064 (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3065 coeff[ x ], coeff[x+1],
3066 coeff[x+2], coeff[x+3]),x+=4;
3067 if ( channel & IndexChannel )
3068 (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3069 coeff[ x ], coeff[x+1],
3070 coeff[x+2], coeff[x+3]),x+=4;
3071 if ( channel & OpacityChannel )
3072 (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3073 coeff[ x ], coeff[x+1],
3074 coeff[x+2], coeff[x+3]),x+=4;
3075 break;
3076 }
3077 default:
3078 /* sparse color method is too complex for FX emulation */
3079 break;
3080 }
3081 }
3082
3083 /* Generate new image for generated interpolated gradient.
3084 * ASIDE: Actually we could have just replaced the colors of the original
3085 * image, but IM Core policy, is if storage class could change then clone
3086 * the image.
3087 */
3088
3089 sparse_image=CloneImage(image,0,0,MagickTrue,exception);
3090 if (sparse_image == (Image *) NULL)
3091 return((Image *) NULL);
3092 if (SetImageStorageClass(sparse_image,DirectClass) == MagickFalse)
3093 { /* if image is ColorMapped - change it to DirectClass */
3094 InheritException(exception,&image->exception);
3095 sparse_image=DestroyImage(sparse_image);
3096 return((Image *) NULL);
3097 }
3098 { /* ----- MAIN CODE ----- */
3099 CacheView
3100 *sparse_view;
3101
3102 MagickBooleanType
3103 status;
3104
3105 MagickOffsetType
3106 progress;
3107
3108 ssize_t
3109 j;
3110
3111 status=MagickTrue;
3112 progress=0;
3113 sparse_view=AcquireAuthenticCacheView(sparse_image,exception);
3114#if defined(MAGICKCORE_OPENMP_SUPPORT)
3115 #pragma omp parallel for schedule(static) shared(progress,status) \
3116 magick_number_threads(image,sparse_image,sparse_image->rows,1)
3117#endif
3118 for (j=0; j < (ssize_t) sparse_image->rows; j++)
3119 {
3120 MagickBooleanType
3121 sync;
3122
3123 MagickPixelPacket
3124 pixel; /* pixel to assign to distorted image */
3125
3126 IndexPacket
3127 *magick_restrict indexes;
3128
3129 ssize_t
3130 i;
3131
3132 PixelPacket
3133 *magick_restrict q;
3134
3135 q=GetCacheViewAuthenticPixels(sparse_view,0,j,sparse_image->columns,
3136 1,exception);
3137 if (q == (PixelPacket *) NULL)
3138 {
3139 status=MagickFalse;
3140 continue;
3141 }
3142 indexes=GetCacheViewAuthenticIndexQueue(sparse_view);
3143 GetMagickPixelPacket(sparse_image,&pixel);
3144 for (i=0; i < (ssize_t) image->columns; i++)
3145 {
3146 SetMagickPixelPacket(image,q,indexes,&pixel);
3147 switch (sparse_method)
3148 {
3149 case BarycentricColorInterpolate:
3150 {
3151 ssize_t x=0;
3152 if ( channel & RedChannel )
3153 pixel.red = coeff[x]*i +coeff[x+1]*j
3154 +coeff[x+2], x+=3;
3155 if ( channel & GreenChannel )
3156 pixel.green = coeff[x]*i +coeff[x+1]*j
3157 +coeff[x+2], x+=3;
3158 if ( channel & BlueChannel )
3159 pixel.blue = coeff[x]*i +coeff[x+1]*j
3160 +coeff[x+2], x+=3;
3161 if ( channel & IndexChannel )
3162 pixel.index = coeff[x]*i +coeff[x+1]*j
3163 +coeff[x+2], x+=3;
3164 if ( channel & OpacityChannel )
3165 pixel.opacity = coeff[x]*i +coeff[x+1]*j
3166 +coeff[x+2], x+=3;
3167 break;
3168 }
3169 case BilinearColorInterpolate:
3170 {
3171 ssize_t x=0;
3172 if ( channel & RedChannel )
3173 pixel.red = coeff[x]*i + coeff[x+1]*j +
3174 coeff[x+2]*i*j + coeff[x+3], x+=4;
3175 if ( channel & GreenChannel )
3176 pixel.green = coeff[x]*i + coeff[x+1]*j +
3177 coeff[x+2]*i*j + coeff[x+3], x+=4;
3178 if ( channel & BlueChannel )
3179 pixel.blue = coeff[x]*i + coeff[x+1]*j +
3180 coeff[x+2]*i*j + coeff[x+3], x+=4;
3181 if ( channel & IndexChannel )
3182 pixel.index = coeff[x]*i + coeff[x+1]*j +
3183 coeff[x+2]*i*j + coeff[x+3], x+=4;
3184 if ( channel & OpacityChannel )
3185 pixel.opacity = coeff[x]*i + coeff[x+1]*j +
3186 coeff[x+2]*i*j + coeff[x+3], x+=4;
3187 break;
3188 }
3189 case InverseColorInterpolate:
3190 case ShepardsColorInterpolate:
3191 { /* Inverse (Squared) Distance weights average (IDW) */
3192 size_t
3193 k;
3194 double
3195 denominator;
3196
3197 if ( channel & RedChannel ) pixel.red = 0.0;
3198 if ( channel & GreenChannel ) pixel.green = 0.0;
3199 if ( channel & BlueChannel ) pixel.blue = 0.0;
3200 if ( channel & IndexChannel ) pixel.index = 0.0;
3201 if ( channel & OpacityChannel ) pixel.opacity = 0.0;
3202 denominator = 0.0;
3203 for(k=0; k<number_arguments; k+=2+number_colors) {
3204 ssize_t x=(ssize_t) k+2;
3205 double weight =
3206 ((double)i-arguments[ k ])*((double)i-arguments[ k ])
3207 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
3208 weight = pow(weight,coeff[0]); /* inverse of power factor */
3209 weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
3210 if ( channel & RedChannel )
3211 pixel.red += arguments[x++]*weight;
3212 if ( channel & GreenChannel )
3213 pixel.green += arguments[x++]*weight;
3214 if ( channel & BlueChannel )
3215 pixel.blue += arguments[x++]*weight;
3216 if ( channel & IndexChannel )
3217 pixel.index += arguments[x++]*weight;
3218 if ( channel & OpacityChannel )
3219 pixel.opacity += arguments[x++]*weight;
3220 denominator += weight;
3221 }
3222 if ( channel & RedChannel ) pixel.red /= denominator;
3223 if ( channel & GreenChannel ) pixel.green /= denominator;
3224 if ( channel & BlueChannel ) pixel.blue /= denominator;
3225 if ( channel & IndexChannel ) pixel.index /= denominator;
3226 if ( channel & OpacityChannel ) pixel.opacity /= denominator;
3227 break;
3228 }
3229 case ManhattanColorInterpolate:
3230 {
3231 size_t
3232 k;
3233
3234 double
3235 minimum = MagickMaximumValue;
3236
3237 /*
3238 Just use the closest control point you can find!
3239 */
3240 for(k=0; k<number_arguments; k+=2+number_colors) {
3241 double distance =
3242 fabs((double)i-arguments[ k ])
3243 + fabs((double)j-arguments[k+1]);
3244 if ( distance < minimum ) {
3245 ssize_t x=(ssize_t) k+2;
3246 if ( channel & RedChannel ) pixel.red = arguments[x++];
3247 if ( channel & GreenChannel ) pixel.green = arguments[x++];
3248 if ( channel & BlueChannel ) pixel.blue = arguments[x++];
3249 if ( channel & IndexChannel ) pixel.index = arguments[x++];
3250 if ( channel & OpacityChannel ) pixel.opacity = arguments[x++];
3251 minimum = distance;
3252 }
3253 }
3254 break;
3255 }
3256 case VoronoiColorInterpolate:
3257 default:
3258 {
3259 size_t
3260 k;
3261
3262 double
3263 minimum = MagickMaximumValue;
3264
3265 /*
3266 Just use the closest control point you can find!
3267 */
3268 for(k=0; k<number_arguments; k+=2+number_colors) {
3269 double distance =
3270 ((double)i-arguments[ k ])*((double)i-arguments[ k ])
3271 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
3272 if ( distance < minimum ) {
3273 ssize_t x=(ssize_t) k+2;
3274 if ( channel & RedChannel ) pixel.red = arguments[x++];
3275 if ( channel & GreenChannel ) pixel.green = arguments[x++];
3276 if ( channel & BlueChannel ) pixel.blue = arguments[x++];
3277 if ( channel & IndexChannel ) pixel.index = arguments[x++];
3278 if ( channel & OpacityChannel ) pixel.opacity = arguments[x++];
3279 minimum = distance;
3280 }
3281 }
3282 break;
3283 }
3284 }
3285 /* set the color directly back into the source image */
3286 if ( channel & RedChannel )
3287 pixel.red=ClampPixel((MagickRealType) QuantumRange*pixel.red);
3288 if ( channel & GreenChannel )
3289 pixel.green=ClampPixel((MagickRealType) QuantumRange*pixel.green);
3290 if ( channel & BlueChannel )
3291 pixel.blue=ClampPixel((MagickRealType) QuantumRange*pixel.blue);
3292 if ( channel & IndexChannel )
3293 pixel.index=ClampPixel((MagickRealType) QuantumRange*pixel.index);
3294 if ( channel & OpacityChannel )
3295 pixel.opacity=ClampPixel((MagickRealType) QuantumRange*pixel.opacity);
3296 SetPixelPacket(sparse_image,&pixel,q,indexes);
3297 q++;
3298 indexes++;
3299 }
3300 sync=SyncCacheViewAuthenticPixels(sparse_view,exception);
3301 if (sync == MagickFalse)
3302 status=MagickFalse;
3303 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3304 {
3305 MagickBooleanType
3306 proceed;
3307
3308#if defined(MAGICKCORE_OPENMP_SUPPORT)
3309 #pragma omp atomic
3310#endif
3311 progress++;
3312 proceed=SetImageProgress(image,SparseColorTag,progress,image->rows);
3313 if (proceed == MagickFalse)
3314 status=MagickFalse;
3315 }
3316 }
3317 sparse_view=DestroyCacheView(sparse_view);
3318 if (status == MagickFalse)
3319 sparse_image=DestroyImage(sparse_image);
3320 }
3321 coeff = (double *) RelinquishMagickMemory(coeff);
3322 return(sparse_image);
3323}