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