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