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