]> granicus.if.org Git - imagemagick/blob - MagickCore/resize.c
(no commit message)
[imagemagick] / MagickCore / resize.c
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %                                                                             %
4 %                                                                             %
5 %                                                                             %
6 %                 RRRR   EEEEE  SSSSS  IIIII  ZZZZZ  EEEEE                    %
7 %                 R   R  E      SS       I       ZZ  E                        %
8 %                 RRRR   EEE     SSS     I     ZZZ   EEE                      %
9 %                 R R    E         SS    I    ZZ     E                        %
10 %                 R  R   EEEEE  SSSSS  IIIII  ZZZZZ  EEEEE                    %
11 %                                                                             %
12 %                                                                             %
13 %                      MagickCore Image Resize Methods                        %
14 %                                                                             %
15 %                              Software Design                                %
16 %                                   Cristy                                    %
17 %                                 July 1992                                   %
18 %                                                                             %
19 %                                                                             %
20 %  Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization      %
21 %  dedicated to making software imaging solutions freely available.           %
22 %                                                                             %
23 %  You may not use this file except in compliance with the License.  You may  %
24 %  obtain a copy of the License at                                            %
25 %                                                                             %
26 %    http://www.imagemagick.org/script/license.php                            %
27 %                                                                             %
28 %  Unless required by applicable law or agreed to in writing, software        %
29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
31 %  See the License for the specific language governing permissions and        %
32 %  limitations under the License.                                             %
33 %                                                                             %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 */
38 \f
39 /*
40   Include declarations.
41 */
42 #include "MagickCore/studio.h"
43 #include "MagickCore/artifact.h"
44 #include "MagickCore/blob.h"
45 #include "MagickCore/cache.h"
46 #include "MagickCore/cache-view.h"
47 #include "MagickCore/channel.h"
48 #include "MagickCore/color.h"
49 #include "MagickCore/color-private.h"
50 #include "MagickCore/draw.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/image-private.h"
56 #include "MagickCore/list.h"
57 #include "MagickCore/memory_.h"
58 #include "MagickCore/memory-private.h"
59 #include "MagickCore/magick.h"
60 #include "MagickCore/pixel-accessor.h"
61 #include "MagickCore/property.h"
62 #include "MagickCore/monitor.h"
63 #include "MagickCore/monitor-private.h"
64 #include "MagickCore/nt-base-private.h"
65 #include "MagickCore/option.h"
66 #include "MagickCore/pixel.h"
67 #include "MagickCore/pixel-private.h"
68 #include "MagickCore/quantum-private.h"
69 #include "MagickCore/resample.h"
70 #include "MagickCore/resample-private.h"
71 #include "MagickCore/resize.h"
72 #include "MagickCore/resize-private.h"
73 #include "MagickCore/resource_.h"
74 #include "MagickCore/string_.h"
75 #include "MagickCore/string-private.h"
76 #include "MagickCore/thread-private.h"
77 #include "MagickCore/token.h"
78 #include "MagickCore/utility.h"
79 #include "MagickCore/utility-private.h"
80 #include "MagickCore/version.h"
81 #if defined(MAGICKCORE_LQR_DELEGATE)
82 #include <lqr.h>
83 #endif
84 \f
85 /*
86   Typedef declarations.
87 */
88 struct _ResizeFilter
89 {
90   double
91     (*filter)(const double,const ResizeFilter *),
92     (*window)(const double,const ResizeFilter *),
93     support,        /* filter region of support - the filter support limit */
94     window_support, /* window support, usally equal to support (expert only) */
95     scale,          /* dimension scaling to fit window support (usally 1.0) */
96     blur,           /* x-scale (blur-sharpen) */
97     coefficient[7]; /* cubic coefficents for BC-cubic filters */
98
99   ResizeWeightingFunctionType
100     filterWeightingType,
101     windowWeightingType;
102
103   size_t
104     signature;
105 };
106 \f
107 /*
108   Forward declaractions.
109 */
110 static double
111   I0(double x),
112   BesselOrderOne(double),
113   Sinc(const double, const ResizeFilter *),
114   SincFast(const double, const ResizeFilter *);
115 \f
116 /*
117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
118 %                                                                             %
119 %                                                                             %
120 %                                                                             %
121 +   F i l t e r F u n c t i o n s                                             %
122 %                                                                             %
123 %                                                                             %
124 %                                                                             %
125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
126 %
127 %  These are the various filter and windowing functions that are provided.
128 %
129 %  They are internal to this module only.  See AcquireResizeFilterInfo() for
130 %  details of the access to these functions, via the GetResizeFilterSupport()
131 %  and GetResizeFilterWeight() API interface.
132 %
133 %  The individual filter functions have this format...
134 %
135 %     static MagickRealtype *FilterName(const double x,const double support)
136 %
137 %  A description of each parameter follows:
138 %
139 %    o x: the distance from the sampling point generally in the range of 0 to
140 %      support.  The GetResizeFilterWeight() ensures this a positive value.
141 %
142 %    o resize_filter: current filter information.  This allows function to
143 %      access support, and possibly other pre-calculated information defining
144 %      the functions.
145 %
146 */
147
148 static double Blackman(const double x,
149   const ResizeFilter *magick_unused(resize_filter))
150 {
151   /*
152     Blackman: 2nd order cosine windowing function:
153       0.42 + 0.5 cos(pi x) + 0.08 cos(2pi x)
154
155     Refactored by Chantal Racette and Nicolas Robidoux to one trig call and
156     five flops.
157   */
158   const double cosine=cos((double) (MagickPI*x));
159   return(0.34+cosine*(0.5+cosine*0.16));
160 }
161
162 static double Bohman(const double x,
163   const ResizeFilter *magick_unused(resize_filter))
164 {
165   /*
166     Bohman: 2rd Order cosine windowing function:
167       (1-x) cos(pi x) + sin(pi x) / pi.
168
169     Refactored by Nicolas Robidoux to one trig call, one sqrt call, and 7 flops,
170     taking advantage of the fact that the support of Bohman is 1.0 (so that we
171     know that sin(pi x) >= 0).
172   */
173   const double cosine=cos((double) (MagickPI*x));
174   const double sine=sqrt(1.0-cosine*cosine);
175   return((1.0-x)*cosine+(1.0/MagickPI)*sine);
176 }
177
178 static double Box(const double magick_unused(x),
179   const ResizeFilter *magick_unused(resize_filter))
180 {
181   /*
182     A Box filter is a equal weighting function (all weights equal).
183     DO NOT LIMIT results by support or resize point sampling will work
184     as it requests points beyond its normal 0.0 support size.
185   */
186   return(1.0);
187 }
188
189 static double Cosine(const double x,
190   const ResizeFilter *magick_unused(resize_filter))
191 {
192   /*
193     Cosine window function:
194       cos((pi/2)*x).
195   */
196   return((double)cos((double) (MagickPI2*x)));
197 }
198
199 static double CubicBC(const double x,const ResizeFilter *resize_filter)
200 {
201   /*
202     Cubic Filters using B,C determined values:
203        Mitchell-Netravali  B = 1/3 C = 1/3  "Balanced" cubic spline filter
204        Catmull-Rom         B = 0   C = 1/2  Interpolatory and exact on linears
205        Spline              B = 1   C = 0    B-Spline Gaussian approximation
206        Hermite             B = 0   C = 0    B-Spline interpolator
207
208     See paper by Mitchell and Netravali, Reconstruction Filters in Computer
209     Graphics Computer Graphics, Volume 22, Number 4, August 1988
210     http://www.cs.utexas.edu/users/fussell/courses/cs384g/lectures/mitchell/
211     Mitchell.pdf.
212
213     Coefficents are determined from B,C values:
214        P0 = (  6 - 2*B       )/6 = coeff[0]
215        P1 =         0
216        P2 = (-18 +12*B + 6*C )/6 = coeff[1]
217        P3 = ( 12 - 9*B - 6*C )/6 = coeff[2]
218        Q0 = (      8*B +24*C )/6 = coeff[3]
219        Q1 = (    -12*B -48*C )/6 = coeff[4]
220        Q2 = (      6*B +30*C )/6 = coeff[5]
221        Q3 = (    - 1*B - 6*C )/6 = coeff[6]
222
223     which are used to define the filter:
224
225        P0 + P1*x + P2*x^2 + P3*x^3      0 <= x < 1
226        Q0 + Q1*x + Q2*x^2 + Q3*x^3      1 <= x < 2
227
228     which ensures function is continuous in value and derivative (slope).
229   */
230   if (x < 1.0)
231     return(resize_filter->coefficient[0]+x*(x*
232       (resize_filter->coefficient[1]+x*resize_filter->coefficient[2])));
233   if (x < 2.0)
234     return(resize_filter->coefficient[3]+x*(resize_filter->coefficient[4]+x*
235       (resize_filter->coefficient[5]+x*resize_filter->coefficient[6])));
236   return(0.0);
237 }
238
239 static double Gaussian(const double x,const ResizeFilter *resize_filter)
240 {
241   /*
242     Gaussian with a sigma = 1/2 (or as user specified)
243
244     Gaussian Formula (1D) ...
245         exp( -(x^2)/((2.0*sigma^2) ) / (sqrt(2*PI)*sigma^2))
246
247     Gaussian Formula (2D) ...
248         exp( -(x^2+y^2)/(2.0*sigma^2) ) / (PI*sigma^2) )
249     or for radius
250         exp( -(r^2)/(2.0*sigma^2) ) / (PI*sigma^2) )
251
252     Note that it is only a change from 1-d to radial form is in the
253     normalization multiplier which is not needed or used when Gaussian is used
254     as a filter.
255
256     The constants are pre-calculated...
257
258         coeff[0]=sigma;
259         coeff[1]=1.0/(2.0*sigma^2);
260         coeff[2]=1.0/(sqrt(2*PI)*sigma^2);
261
262         exp( -coeff[1]*(x^2)) ) * coeff[2];
263
264     However the multiplier coeff[1] is need, the others are informative only.
265
266     This separates the gaussian 'sigma' value from the 'blur/support'
267     settings allowing for its use in special 'small sigma' gaussians,
268     without the filter 'missing' pixels because the support becomes too
269     small.
270   */
271   return(exp((double)(-resize_filter->coefficient[1]*x*x)));
272 }
273
274 static double Hann(const double x,
275   const ResizeFilter *magick_unused(resize_filter))
276 {
277   /*
278     Cosine window function:
279       0.5+0.5*cos(pi*x).
280   */
281   const double cosine=cos((double) (MagickPI*x));
282   return(0.5+0.5*cosine);
283 }
284
285 static double Hamming(const double x,
286   const ResizeFilter *magick_unused(resize_filter))
287 {
288   /*
289     Offset cosine window function:
290      .54 + .46 cos(pi x).
291   */
292   const double cosine=cos((double) (MagickPI*x));
293   return(0.54+0.46*cosine);
294 }
295
296 static double Jinc(const double x,
297   const ResizeFilter *magick_unused(resize_filter))
298 {
299   /*
300     See Pratt "Digital Image Processing" p.97 for Jinc/Bessel functions.
301     http://mathworld.wolfram.com/JincFunction.html and page 11 of
302     http://www.ph.ed.ac.uk/%7ewjh/teaching/mo/slides/lens/lens.pdf
303
304     The original "zoom" program by Paul Heckbert called this "Bessel".  But
305     really it is more accurately named "Jinc".
306   */
307   if (x == 0.0)
308     return(0.5*MagickPI);
309   return(BesselOrderOne(MagickPI*x)/x);
310 }
311
312 static double Kaiser(const double x,const ResizeFilter *resize_filter)
313 {
314   /*
315     Kaiser Windowing Function (bessel windowing)
316
317        I0( beta * sqrt( 1-x^2) ) / IO(0)
318
319     Beta (coeff[0]) is a free value from 5 to 8 (defaults to 6.5).
320     However it is typically defined in terms of Alpha*PI
321
322     The normalization factor (coeff[1]) is not actually needed,
323     but without it the filters has a large value at x=0 making it
324     difficult to compare the function with other windowing functions.
325   */
326   return(resize_filter->coefficient[1]*I0(resize_filter->coefficient[0]*
327     sqrt((double) (1.0-x*x))));
328 }
329
330 static double Lagrange(const double x,const ResizeFilter *resize_filter)
331 {
332   double
333     value;
334
335   register ssize_t
336     i;
337
338   ssize_t
339     n,
340     order;
341
342   /*
343     Lagrange piecewise polynomial fit of sinc: N is the 'order' of the lagrange
344     function and depends on the overall support window size of the filter. That
345     is: for a support of 2, it gives a lagrange-4 (piecewise cubic function).
346
347     "n" identifies the piece of the piecewise polynomial.
348
349     See Survey: Interpolation Methods, IEEE Transactions on Medical Imaging,
350     Vol 18, No 11, November 1999, p1049-1075, -- Equation 27 on p1064.
351   */
352   if (x > resize_filter->support)
353     return(0.0);
354   order=(ssize_t) (2.0*resize_filter->window_support);  /* number of pieces */
355   n=(ssize_t) (resize_filter->window_support+x);
356   value=1.0f;
357   for (i=0; i < order; i++)
358     if (i != n)
359       value*=(n-i-x)/(n-i);
360   return(value);
361 }
362
363 static double Quadratic(const double x,
364   const ResizeFilter *magick_unused(resize_filter))
365 {
366   /*
367     2rd order (quadratic) B-Spline approximation of Gaussian.
368   */
369   if (x < 0.5)
370     return(0.75-x*x);
371   if (x < 1.5)
372     return(0.5*(x-1.5)*(x-1.5));
373   return(0.0);
374 }
375
376 static double Sinc(const double x,
377   const ResizeFilter *magick_unused(resize_filter))
378 {
379   /*
380     Scaled sinc(x) function using a trig call:
381       sinc(x) == sin(pi x)/(pi x).
382   */
383   if (x != 0.0)
384     {
385       const double alpha=(double) (MagickPI*x);
386       return(sin((double) alpha)/alpha);
387     }
388   return((double) 1.0);
389 }
390
391 static double SincFast(const double x,
392   const ResizeFilter *magick_unused(resize_filter))
393 {
394   /*
395     Approximations of the sinc function sin(pi x)/(pi x) over the interval
396     [-4,4] constructed by Nicolas Robidoux and Chantal Racette with funding
397     from the Natural Sciences and Engineering Research Council of Canada.
398
399     Although the approximations are polynomials (for low order of
400     approximation) and quotients of polynomials (for higher order of
401     approximation) and consequently are similar in form to Taylor polynomials /
402     Pade approximants, the approximations are computed with a completely
403     different technique.
404
405     Summary: These approximations are "the best" in terms of bang (accuracy)
406     for the buck (flops). More specifically: Among the polynomial quotients
407     that can be computed using a fixed number of flops (with a given "+ - * /
408     budget"), the chosen polynomial quotient is the one closest to the
409     approximated function with respect to maximum absolute relative error over
410     the given interval.
411
412     The Remez algorithm, as implemented in the boost library's minimax package,
413     is the key to the construction: http://www.boost.org/doc/libs/1_36_0/libs/
414     math/doc/sf_and_dist/html/math_toolkit/backgrounders/remez.html
415
416     If outside of the interval of approximation, use the standard trig formula.
417   */
418   if (x > 4.0)
419     {
420       const double alpha=(double) (MagickPI*x);
421       return(sin((double) alpha)/alpha);
422     }
423   {
424     /*
425       The approximations only depend on x^2 (sinc is an even function).
426     */
427     const double xx = x*x;
428 #if MAGICKCORE_QUANTUM_DEPTH <= 8
429     /*
430       Maximum absolute relative error 6.3e-6 < 1/2^17.
431     */
432     const double c0 = 0.173610016489197553621906385078711564924e-2L;
433     const double c1 = -0.384186115075660162081071290162149315834e-3L;
434     const double c2 = 0.393684603287860108352720146121813443561e-4L;
435     const double c3 = -0.248947210682259168029030370205389323899e-5L;
436     const double c4 = 0.107791837839662283066379987646635416692e-6L;
437     const double c5 = -0.324874073895735800961260474028013982211e-8L;
438     const double c6 = 0.628155216606695311524920882748052490116e-10L;
439     const double c7 = -0.586110644039348333520104379959307242711e-12L;
440     const double p =
441       c0+xx*(c1+xx*(c2+xx*(c3+xx*(c4+xx*(c5+xx*(c6+xx*c7))))));
442     return((xx-1.0)*(xx-4.0)*(xx-9.0)*(xx-16.0)*p);
443 #elif MAGICKCORE_QUANTUM_DEPTH <= 16
444     /*
445       Max. abs. rel. error 2.2e-8 < 1/2^25.
446     */
447     const double c0 = 0.173611107357320220183368594093166520811e-2L;
448     const double c1 = -0.384240921114946632192116762889211361285e-3L;
449     const double c2 = 0.394201182359318128221229891724947048771e-4L;
450     const double c3 = -0.250963301609117217660068889165550534856e-5L;
451     const double c4 = 0.111902032818095784414237782071368805120e-6L;
452     const double c5 = -0.372895101408779549368465614321137048875e-8L;
453     const double c6 = 0.957694196677572570319816780188718518330e-10L;
454     const double c7 = -0.187208577776590710853865174371617338991e-11L;
455     const double c8 = 0.253524321426864752676094495396308636823e-13L;
456     const double c9 = -0.177084805010701112639035485248501049364e-15L;
457     const double p =
458       c0+xx*(c1+xx*(c2+xx*(c3+xx*(c4+xx*(c5+xx*(c6+xx*(c7+xx*(c8+xx*c9))))))));
459     return((xx-1.0)*(xx-4.0)*(xx-9.0)*(xx-16.0)*p);
460 #else
461     /*
462       Max. abs. rel. error 1.2e-12 < 1/2^39.
463     */
464     const double c0 = 0.173611111110910715186413700076827593074e-2L;
465     const double c1 = -0.289105544717893415815859968653611245425e-3L;
466     const double c2 = 0.206952161241815727624413291940849294025e-4L;
467     const double c3 = -0.834446180169727178193268528095341741698e-6L;
468     const double c4 = 0.207010104171026718629622453275917944941e-7L;
469     const double c5 = -0.319724784938507108101517564300855542655e-9L;
470     const double c6 = 0.288101675249103266147006509214934493930e-11L;
471     const double c7 = -0.118218971804934245819960233886876537953e-13L;
472     const double p =
473       c0+xx*(c1+xx*(c2+xx*(c3+xx*(c4+xx*(c5+xx*(c6+xx*c7))))));
474     const double d0 = 1.0L;
475     const double d1 = 0.547981619622284827495856984100563583948e-1L;
476     const double d2 = 0.134226268835357312626304688047086921806e-2L;
477     const double d3 = 0.178994697503371051002463656833597608689e-4L;
478     const double d4 = 0.114633394140438168641246022557689759090e-6L;
479     const double q = d0+xx*(d1+xx*(d2+xx*(d3+xx*d4)));
480     return((xx-1.0)*(xx-4.0)*(xx-9.0)*(xx-16.0)/q*p);
481 #endif
482   }
483 }
484
485 static double Triangle(const double x,
486   const ResizeFilter *magick_unused(resize_filter))
487 {
488   /*
489     1st order (linear) B-Spline, bilinear interpolation, Tent 1D filter, or
490     a Bartlett 2D Cone filter.  Also used as a Bartlett Windowing function
491     for Sinc().
492   */
493   if (x < 1.0)
494     return(1.0-x);
495   return(0.0);
496 }
497
498 static double Welch(const double x,
499   const ResizeFilter *magick_unused(resize_filter))
500 {
501   /*
502     Welch parabolic windowing filter.
503   */
504   if (x < 1.0)
505     return(1.0-x*x);
506   return(0.0);
507 }
508 \f
509 /*
510 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
511 %                                                                             %
512 %                                                                             %
513 %                                                                             %
514 +   A c q u i r e R e s i z e F i l t e r                                     %
515 %                                                                             %
516 %                                                                             %
517 %                                                                             %
518 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
519 %
520 %  AcquireResizeFilter() allocates the ResizeFilter structure.  Choose from
521 %  these filters:
522 %
523 %  FIR (Finite impulse Response) Filters
524 %      Box         Triangle   Quadratic
525 %      Spline      Hermite    Catrom
526 %      Mitchell
527 %
528 %  IIR (Infinite impulse Response) Filters
529 %      Gaussian     Sinc        Jinc (Bessel)
530 %
531 %  Windowed Sinc/Jinc Filters
532 %      Blackman    Bohman     Lanczos
533 %      Hann        Hamming    Cosine
534 %      Kaiser      Welch      Parzen
535 %      Bartlett
536 %
537 %  Special Purpose Filters
538 %      Cubic  SincFast  LanczosSharp  Lanczos2  Lanczos2Sharp
539 %      Robidoux RobidouxSharp
540 %
541 %  The users "-filter" selection is used to lookup the default 'expert'
542 %  settings for that filter from a internal table.  However any provided
543 %  'expert' settings (see below) may override this selection.
544 %
545 %  FIR filters are used as is, and are limited to that filters support window
546 %  (unless over-ridden).  'Gaussian' while classed as an IIR filter, is also
547 %  simply clipped by its support size (currently 1.5 or approximately 3*sigma
548 %  as recommended by many references)
549 %
550 %  The special a 'cylindrical' filter flag will promote the default 4-lobed
551 %  Windowed Sinc filter to a 3-lobed Windowed Jinc equivalent, which is better
552 %  suited to this style of image resampling. This typically happens when using
553 %  such a filter for images distortions.
554 %
555 %  SPECIFIC FILTERS:
556 %
557 %  Directly requesting 'Sinc', 'Jinc' function as a filter will force the use
558 %  of function without any windowing, or promotion for cylindrical usage.  This
559 %  is not recommended, except by image processing experts, especially as part
560 %  of expert option filter function selection.
561 %
562 %  Two forms of the 'Sinc' function are available: Sinc and SincFast.  Sinc is
563 %  computed using the traditional sin(pi*x)/(pi*x); it is selected if the user
564 %  specifically specifies the use of a Sinc filter. SincFast uses highly
565 %  accurate (and fast) polynomial (low Q) and rational (high Q) approximations,
566 %  and will be used by default in most cases.
567 %
568 %  The Lanczos filter is a special 3-lobed Sinc-windowed Sinc filter (promoted
569 %  to Jinc-windowed Jinc for cylindrical (Elliptical Weighted Average) use).
570 %  The Sinc version is the most popular windowed filter.
571 %
572 %  LanczosSharp is a slightly sharpened (blur=0.9812505644269356 < 1) form of
573 %  the Lanczos filter, specifically designed for EWA distortion (as a
574 %  Jinc-Jinc); it can also be used as a slightly sharper orthogonal Lanczos
575 %  (Sinc-Sinc) filter. The chosen blur value comes as close as possible to
576 %  satisfying the following condition without changing the character of the
577 %  corresponding EWA filter:
578 %
579 %    'No-Op' Vertical and Horizontal Line Preservation Condition: Images with
580 %    only vertical or horizontal features are preserved when performing 'no-op"
581 %    with EWA distortion.
582 %
583 %  The Lanczos2 and Lanczos2Sharp filters are 2-lobe versions of the Lanczos
584 %  filters.  The 'sharp' version uses a blur factor of 0.9549963639785485,
585 %  again chosen because the resulting EWA filter comes as close as possible to
586 %  satisfying the above condition.
587 %
588 %  Robidoux is another filter tuned for EWA. It is the Keys cubic filter
589 %  defined by B=(228 - 108 sqrt(2))/199. Robidoux satisfies the "'No-Op'
590 %  Vertical and Horizontal Line Preservation Condition" exactly, and it
591 %  moderately blurs high frequency 'pixel-hash' patterns under no-op.  It turns
592 %  out to be close to both Mitchell and Lanczos2Sharp.  For example, its first
593 %  crossing is at (36 sqrt(2) + 123)/(72 sqrt(2) + 47), almost the same as the
594 %  first crossing of Mitchell and Lanczos2Sharp.
595 %
596 %  RodidouxSharp is a slightly sharper version of Rodidoux, some believe it
597 %  is too sharp.  It is designed to minimize the maximum possible change in
598 %  a pixel value which is at one of the extremes (e.g., 0 or 255) under no-op
599 %  conditions.  Amazingly Mitchell falls roughly between Rodidoux and
600 %  RodidouxSharp, though this seems to have been pure coincidence.
601 %
602 %  'EXPERT' OPTIONS:
603 %
604 %  These artifact "defines" are not recommended for production use without
605 %  expert knowledge of resampling, filtering, and the effects they have on the
606 %  resulting resampled (resized or distorted) image.
607 %
608 %  They can be used to override any and all filter default, and it is
609 %  recommended you make good use of "filter:verbose" to make sure that the
610 %  overall effect of your selection (before and after) is as expected.
611 %
612 %    "filter:verbose"  controls whether to output the exact results of the
613 %       filter selections made, as well as plotting data for graphing the
614 %       resulting filter over the filters support range.
615 %
616 %    "filter:filter"  select the main function associated with this filter
617 %       name, as the weighting function of the filter.  This can be used to
618 %       set a windowing function as a weighting function, for special
619 %       purposes, such as graphing.
620 %
621 %       If a "filter:window" operation has not been provided, a 'Box'
622 %       windowing function will be set to denote that no windowing function is
623 %       being used.
624 %
625 %    "filter:window"  Select this windowing function for the filter. While any
626 %       filter could be used as a windowing function, using the 'first lobe' of
627 %       that filter over the whole support window, using a non-windowing
628 %       function is not advisible. If no weighting filter function is specified
629 %       a 'SincFast' filter is used.
630 %
631 %    "filter:lobes"  Number of lobes to use for the Sinc/Jinc filter.  This a
632 %       simpler method of setting filter support size that will correctly
633 %       handle the Sinc/Jinc switch for an operators filtering requirements.
634 %       Only integers should be given.
635 %
636 %    "filter:support" Set the support size for filtering to the size given.
637 %       This not recommended for Sinc/Jinc windowed filters (lobes should be
638 %       used instead).  This will override any 'filter:lobes' option.
639 %
640 %    "filter:win-support" Scale windowing function to this size instead.  This
641 %       causes the windowing (or self-windowing Lagrange filter) to act is if
642 %       the support window it much much larger than what is actually supplied
643 %       to the calling operator.  The filter however is still clipped to the
644 %       real support size given, by the support range supplied to the caller.
645 %       If unset this will equal the normal filter support size.
646 %
647 %    "filter:blur" Scale the filter and support window by this amount.  A value
648 %       of > 1 will generally result in a more blurred image with more ringing
649 %       effects, while a value <1 will sharpen the resulting image with more
650 %       aliasing effects.
651 %
652 %    "filter:sigma" The sigma value to use for the Gaussian filter only.
653 %       Defaults to '1/2'.  Using a different sigma effectively provides a
654 %       method of using the filter as a 'blur' convolution.  Particularly when
655 %       using it for Distort.
656 %
657 %    "filter:b"
658 %    "filter:c" Override the preset B,C values for a Cubic filter.
659 %       If only one of these are given it is assumes to be a 'Keys' type of
660 %       filter such that B+2C=1, where Keys 'alpha' value = C.
661 %
662 %  Examples:
663 %
664 %  Set a true un-windowed Sinc filter with 10 lobes (very slow):
665 %     -define filter:filter=Sinc
666 %     -define filter:lobes=8
667 %
668 %  Set an 8 lobe Lanczos (Sinc or Jinc) filter:
669 %     -filter Lanczos
670 %     -define filter:lobes=8
671 %
672 %  The format of the AcquireResizeFilter method is:
673 %
674 %      ResizeFilter *AcquireResizeFilter(const Image *image,
675 %        const FilterTypes filter_type,const MagickBooleanType cylindrical,
676 %        ExceptionInfo *exception)
677 %
678 %  A description of each parameter follows:
679 %
680 %    o image: the image.
681 %
682 %    o filter: the filter type, defining a preset filter, window and support.
683 %      The artifact settings listed above will override those selections.
684 %
685 %    o blur: blur the filter by this amount, use 1.0 if unknown.  Image
686 %      artifact "filter:blur" will override this API call usage, including any
687 %      internal change (such as for cylindrical usage).
688 %
689 %    o radial: use a 1D orthogonal filter (Sinc) or 2D cylindrical (radial)
690 %      filter (Jinc).
691 %
692 %    o exception: return any errors or warnings in this structure.
693 %
694 */
695 MagickPrivate ResizeFilter *AcquireResizeFilter(const Image *image,
696   const FilterTypes filter,const MagickBooleanType cylindrical,
697   ExceptionInfo *exception)
698 {
699   const char
700     *artifact;
701
702   FilterTypes
703     filter_type,
704     window_type;
705
706   double
707     B,
708     C,
709     value;
710
711   register ResizeFilter
712     *resize_filter;
713
714   /*
715     Table Mapping given Filter, into Weighting and Windowing functions.
716     A 'Box' windowing function means its a simble non-windowed filter.
717     An 'SincFast' filter function could be upgraded to a 'Jinc' filter if a
718     "cylindrical" is requested, unless a 'Sinc' or 'SincFast' filter was
719     specifically requested by the user.
720
721     WARNING: The order of this table must match the order of the FilterTypes
722     enumeration specified in "resample.h", or the filter names will not match
723     the filter being setup.
724
725     You can check filter setups with the "filter:verbose" expert setting.
726   */
727   static struct
728   {
729     FilterTypes
730       filter,
731       window;
732   } const mapping[SentinelFilter] =
733   {
734     { UndefinedFilter,     BoxFilter      },  /* Undefined (default to Box)   */
735     { PointFilter,         BoxFilter      },  /* SPECIAL: Nearest neighbour   */
736     { BoxFilter,           BoxFilter      },  /* Box averaging filter         */
737     { TriangleFilter,      BoxFilter      },  /* Linear interpolation filter  */
738     { HermiteFilter,       BoxFilter      },  /* Hermite interpolation filter */
739     { SincFastFilter,      HannFilter     },  /* Hann -- cosine-sinc          */
740     { SincFastFilter,      HammingFilter  },  /* Hamming --   '' variation    */
741     { SincFastFilter,      BlackmanFilter },  /* Blackman -- 2*cosine-sinc    */
742     { GaussianFilter,      BoxFilter      },  /* Gaussian blur filter         */
743     { QuadraticFilter,     BoxFilter      },  /* Quadratic Gaussian approx    */
744     { CubicFilter,         BoxFilter      },  /* General Cubic Filter, Spline */
745     { CatromFilter,        BoxFilter      },  /* Cubic-Keys interpolator      */
746     { MitchellFilter,      BoxFilter      },  /* 'Ideal' Cubic-Keys filter    */
747     { JincFilter,          BoxFilter      },  /* Raw 3-lobed Jinc function    */
748     { SincFilter,          BoxFilter      },  /* Raw 4-lobed Sinc function    */
749     { SincFastFilter,      BoxFilter      },  /* Raw fast sinc ("Pade"-type)  */
750     { SincFastFilter,      KaiserFilter   },  /* Kaiser -- square root-sinc   */
751     { LanczosFilter,       WelchFilter    },  /* Welch -- parabolic (3 lobe)  */
752     { SincFastFilter,      CubicFilter    },  /* Parzen -- cubic-sinc         */
753     { SincFastFilter,      BohmanFilter   },  /* Bohman -- 2*cosine-sinc      */
754     { SincFastFilter,      TriangleFilter },  /* Bartlett -- triangle-sinc    */
755     { LagrangeFilter,      BoxFilter      },  /* Lagrange self-windowing      */
756     { LanczosFilter,       LanczosFilter  },  /* Lanczos Sinc-Sinc filters    */
757     { LanczosSharpFilter,  LanczosSharpFilter }, /* | these require */
758     { Lanczos2Filter,      Lanczos2Filter },     /* | special handling */
759     { Lanczos2SharpFilter, Lanczos2SharpFilter },
760     { RobidouxFilter,      BoxFilter      },  /* Cubic Keys tuned for EWA     */
761     { RobidouxSharpFilter, BoxFilter      },  /* Sharper Cubic Keys for EWA   */
762     { LanczosFilter,       CosineFilter   },  /* Cosine window (3 lobes)      */
763     { SplineFilter,        BoxFilter      },  /* Spline Cubic Filter          */
764     { LanczosRadiusFilter, LanczosFilter  },  /* Lanczos with integer radius  */
765   };
766   /*
767     Table mapping the filter/window from the above table to an actual function.
768     The default support size for that filter as a weighting function, the range
769     to scale with to use that function as a sinc windowing function, (typ 1.0).
770
771     Note that the filter_type -> function is 1 to 1 except for Sinc(),
772     SincFast(), and CubicBC() functions, which may have multiple filter to
773     function associations.
774
775     See "filter:verbose" handling below for the function -> filter mapping.
776   */
777   static struct
778   {
779     double
780       (*function)(const double,const ResizeFilter*),
781       support, /* Default lobes/support size of the weighting filter. */
782       scale,   /* Support when function used as a windowing function
783                  Typically equal to the location of the first zero crossing. */
784       B,C;     /* BC-spline coefficients, ignored if not a CubicBC filter. */
785   } const filters[SentinelFilter] =
786   {
787     /*            .---  support window (if used as a Weighting Function)
788                   |    .--- first crossing (if used as a Windowing Function)
789                   |    |    .--- B value for Cubic Function
790                   |    |    |    .---- C value for Cubic Function
791                   |    |    |    |                                    */
792     { Box,       0.5, 0.5, 0.0, 0.0 }, /* Undefined (default to Box)  */
793     { Box,       0.0, 0.5, 0.0, 0.0 }, /* Point (special handling)    */
794     { Box,       0.5, 0.5, 0.0, 0.0 }, /* Box                         */
795     { Triangle,  1.0, 1.0, 0.0, 0.0 }, /* Triangle                    */
796     { CubicBC,   1.0, 1.0, 0.0, 0.0 }, /* Hermite (cubic  B=C=0)      */
797     { Hann,      1.0, 1.0, 0.0, 0.0 }, /* Hann, cosine window         */
798     { Hamming,   1.0, 1.0, 0.0, 0.0 }, /* Hamming, '' variation       */
799     { Blackman,  1.0, 1.0, 0.0, 0.0 }, /* Blackman, 2*cosine window   */
800     { Gaussian,  2.0, 1.5, 0.0, 0.0 }, /* Gaussian                    */
801     { Quadratic, 1.5, 1.5, 0.0, 0.0 }, /* Quadratic gaussian          */
802     { CubicBC,   2.0, 2.0, 1.0, 0.0 }, /* General Cubic Filter        */
803     { CubicBC,   2.0, 1.0, 0.0, 0.5 }, /* Catmull-Rom    (B=0,C=1/2)  */
804     { CubicBC,   2.0, 8.0/7.0, 1./3., 1./3. }, /* Mitchell   (B=C=1/3)    */
805     { Jinc,      3.0, 1.2196698912665045, 0.0, 0.0 }, /* Raw 3-lobed Jinc */
806     { Sinc,      4.0, 1.0, 0.0, 0.0 }, /* Raw 4-lobed Sinc            */
807     { SincFast,  4.0, 1.0, 0.0, 0.0 }, /* Raw fast sinc ("Pade"-type) */
808     { Kaiser,    1.0, 1.0, 0.0, 0.0 }, /* Kaiser (square root window) */
809     { Welch,     1.0, 1.0, 0.0, 0.0 }, /* Welch (parabolic window)    */
810     { CubicBC,   2.0, 2.0, 1.0, 0.0 }, /* Parzen (B-Spline window)    */
811     { Bohman,    1.0, 1.0, 0.0, 0.0 }, /* Bohman, 2*Cosine window     */
812     { Triangle,  1.0, 1.0, 0.0, 0.0 }, /* Bartlett (triangle window)  */
813     { Lagrange,  2.0, 1.0, 0.0, 0.0 }, /* Lagrange sinc approximation */
814     { SincFast,  3.0, 1.0, 0.0, 0.0 }, /* Lanczos, 3-lobed Sinc-Sinc  */
815     { SincFast,  3.0, 1.0, 0.0, 0.0 }, /* Lanczos, Sharpened          */
816     { SincFast,  2.0, 1.0, 0.0, 0.0 }, /* Lanczos, 2-lobed            */
817     { SincFast,  2.0, 1.0, 0.0, 0.0 }, /* Lanczos2, sharpened         */
818     /* Robidoux: Keys cubic close to Lanczos2D sharpened */
819     { CubicBC,   2.0, 1.1685777620836932,
820                             0.37821575509399867, 0.31089212245300067 },
821     /* RobidouxSharp: Sharper version of Robidoux */
822     { CubicBC,   2.0, 1.105822933719019,
823                             0.2620145123990142,  0.3689927438004929  },
824     { Cosine,    1.0, 1.0, 0.0, 0.0 }, /* Low level cosine window     */
825     { CubicBC,   2.0, 2.0, 1.0, 0.0 }, /* Cubic B-Spline (B=1,C=0)    */
826     { SincFast,  3.0, 1.0, 0.0, 0.0 }, /* Lanczos, Interger Radius    */
827   };
828   /*
829     The known zero crossings of the Jinc() or more accurately the Jinc(x*PI)
830     function being used as a filter. It is used by the "filter:lobes" expert
831     setting and for 'lobes' for Jinc functions in the previous table. This way
832     users do not have to deal with the highly irrational lobe sizes of the Jinc
833     filter.
834
835     Values taken from
836     http://cose.math.bas.bg/webMathematica/webComputing/BesselZeros.jsp
837     using Jv-function with v=1, then dividing by PI.
838   */
839   static double
840     jinc_zeros[16] =
841     {
842       1.2196698912665045,
843       2.2331305943815286,
844       3.2383154841662362,
845       4.2410628637960699,
846       5.2427643768701817,
847       6.2439216898644877,
848       7.2447598687199570,
849       8.2453949139520427,
850       9.2458926849494673,
851       10.246293348754916,
852       11.246622794877883,
853       12.246898461138105,
854       13.247132522181061,
855       14.247333735806849,
856       15.247508563037300,
857       16.247661874700962
858    };
859
860   /*
861     Allocate resize filter.
862   */
863   assert(image != (const Image *) NULL);
864   assert(image->signature == MagickSignature);
865   if (image->debug != MagickFalse)
866     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
867   assert(UndefinedFilter < filter && filter < SentinelFilter);
868   assert(exception != (ExceptionInfo *) NULL);
869   assert(exception->signature == MagickSignature);
870   resize_filter=(ResizeFilter *) AcquireMagickMemory(sizeof(*resize_filter));
871   if (resize_filter == (ResizeFilter *) NULL)
872     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
873   (void) ResetMagickMemory(resize_filter,0,sizeof(*resize_filter));
874   /*
875     Defaults for the requested filter.
876   */
877   filter_type=mapping[filter].filter;
878   window_type=mapping[filter].window;
879   resize_filter->blur=1.0;
880   /* Promote 1D Windowed Sinc Filters to a 2D Windowed Jinc filters */
881   if ( IfMagickTrue(cylindrical) && (filter_type == SincFastFilter) &&
882       (filter != SincFastFilter))
883     filter_type=JincFilter;  /* 1D Windowed Sinc => 2D Windowed Jinc filters */
884
885   /* Expert filter setting override */
886   artifact=GetImageArtifact(image,"filter:filter");
887   if (artifact != (const char *) NULL)
888     {
889       ssize_t
890         option;
891
892       option=ParseCommandOption(MagickFilterOptions,MagickFalse,artifact);
893       if ((UndefinedFilter < option) && (option < SentinelFilter))
894         { /* Raw filter request - no window function. */
895           filter_type=(FilterTypes) option;
896           window_type=BoxFilter;
897         }
898       /* Filter override with a specific window function. */
899       artifact=GetImageArtifact(image,"filter:window");
900       if (artifact != (const char *) NULL)
901         {
902           option=ParseCommandOption(MagickFilterOptions,MagickFalse,artifact);
903           if ((UndefinedFilter < option) && (option < SentinelFilter))
904             window_type=(FilterTypes) option;
905         }
906     }
907   else
908     {
909       /* Window specified, but no filter function?  Assume Sinc/Jinc. */
910       artifact=GetImageArtifact(image,"filter:window");
911       if (artifact != (const char *) NULL)
912         {
913           ssize_t
914             option;
915
916           option=ParseCommandOption(MagickFilterOptions,MagickFalse,artifact);
917           if ((UndefinedFilter < option) && (option < SentinelFilter))
918             {
919               filter_type= IfMagickTrue(cylindrical) ? JincFilter
920                                                      : SincFastFilter;
921               window_type=(FilterTypes) option;
922             }
923         }
924     }
925
926   /* Assign the real functions to use for the filters selected. */
927   resize_filter->filter=filters[filter_type].function;
928   resize_filter->support=filters[filter_type].support;
929   resize_filter->window=filters[window_type].function;
930   resize_filter->scale=filters[window_type].scale;
931   resize_filter->signature=MagickSignature;
932
933   /* Filter Modifications for orthogonal/cylindrical usage */
934   if (cylindrical != MagickFalse)
935     switch (filter_type)
936     {
937       case BoxFilter:
938         /* Support for Cylindrical Box should be sqrt(2)/2 */
939         resize_filter->support=(double) MagickSQ1_2;
940         break;
941       case LanczosFilter:
942       case LanczosSharpFilter:
943       case Lanczos2Filter:
944       case Lanczos2SharpFilter:
945       case LanczosRadiusFilter:
946         resize_filter->filter=filters[JincFilter].function;
947         resize_filter->window=filters[JincFilter].function;
948         resize_filter->scale=filters[JincFilter].scale;
949         /* number of lobes (support window size) remain unchanged */
950         break;
951       default:
952         break;
953     }
954   /* Global Sharpening (regardless of orthoginal/cylindrical) */
955   switch (filter_type)
956   {
957     case LanczosSharpFilter:
958       resize_filter->blur *= 0.9812505644269356;
959       break;
960     case Lanczos2SharpFilter:
961       resize_filter->blur *= 0.9549963639785485;
962       break;
963     /* case LanczosRadius:  blur adjust is done after lobes */
964     default:
965       break;
966   }
967
968   /*
969     Expert Option Modifications.
970   */
971
972   /* User Gaussian Sigma Override - no support change */
973   if ((resize_filter->filter == Gaussian) ||
974       (resize_filter->window == Gaussian) ) {
975     value=0.5;    /* guassian sigma default, half pixel */
976     artifact=GetImageArtifact(image,"filter:sigma");
977     if (artifact != (const char *) NULL)
978       value=StringToDouble(artifact,(char **) NULL);
979     /* Define coefficents for Gaussian */
980     resize_filter->coefficient[0]=value;                 /* note sigma too */
981     resize_filter->coefficient[1]=PerceptibleReciprocal(2.0*value*value); /* sigma scaling */
982     resize_filter->coefficient[2]=PerceptibleReciprocal(Magick2PI*value*value);
983        /* normalization - not actually needed or used! */
984     if ( value > 0.5 )
985       resize_filter->support *= 2*value;  /* increase support linearly */
986   }
987
988   /* User Kaiser Alpha Override - no support change */
989   if ((resize_filter->filter == Kaiser) ||
990       (resize_filter->window == Kaiser) ) {
991     value=6.5; /* default beta value for Kaiser bessel windowing function */
992     artifact=GetImageArtifact(image,"filter:alpha");  /* FUTURE: depreciate */
993     if (artifact != (const char *) NULL)
994       value=StringToDouble(artifact,(char **) NULL);
995     artifact=GetImageArtifact(image,"filter:kaiser-beta");
996     if (artifact != (const char *) NULL)
997       value=StringToDouble(artifact,(char **) NULL);
998     artifact=GetImageArtifact(image,"filter:kaiser-alpha");
999     if (artifact != (const char *) NULL)
1000       value=StringToDouble(artifact,(char **) NULL)*MagickPI;
1001     /* Define coefficents for Kaiser Windowing Function */
1002     resize_filter->coefficient[0]=value;         /* alpha */
1003     resize_filter->coefficient[1]=PerceptibleReciprocal(I0(value));
1004       /* normalization */
1005   }
1006
1007   /* Support Overrides */
1008   artifact=GetImageArtifact(image,"filter:lobes");
1009   if (artifact != (const char *) NULL)
1010     {
1011       ssize_t
1012         lobes;
1013
1014       lobes=(ssize_t) StringToLong(artifact);
1015       if (lobes < 1)
1016         lobes=1;
1017       resize_filter->support=(double) lobes;
1018     }
1019   if (resize_filter->filter == Jinc)
1020     {
1021       /*
1022         Convert a Jinc function lobes value to a real support value.
1023       */
1024       if (resize_filter->support > 16)
1025         resize_filter->support=jinc_zeros[15];  /* largest entry in table */
1026       else
1027         resize_filter->support=jinc_zeros[((long) resize_filter->support)-1];
1028       /*
1029         Blur this filter so support is a integer value (lobes dependant).
1030       */
1031       if (filter_type == LanczosRadiusFilter)
1032         resize_filter->blur*=floor(resize_filter->support)/
1033           resize_filter->support;
1034     }
1035   /*
1036     Expert blur override.
1037   */
1038   artifact=GetImageArtifact(image,"filter:blur");
1039   if (artifact != (const char *) NULL)
1040     resize_filter->blur*=StringToDouble(artifact,(char **) NULL);
1041   if (resize_filter->blur < MagickEpsilon)
1042     resize_filter->blur=(double) MagickEpsilon;
1043   /*
1044     Expert override of the support setting.
1045   */
1046   artifact=GetImageArtifact(image,"filter:support");
1047   if (artifact != (const char *) NULL)
1048     resize_filter->support=fabs(StringToDouble(artifact,(char **) NULL));
1049   /*
1050     Scale windowing function separately to the support 'clipping' window
1051     that calling operator is planning to actually use. (Expert override)
1052   */
1053   resize_filter->window_support=resize_filter->support; /* default */
1054   artifact=GetImageArtifact(image,"filter:win-support");
1055   if (artifact != (const char *) NULL)
1056     resize_filter->window_support=fabs(StringToDouble(artifact,(char **) NULL));
1057   /*
1058     Adjust window function scaling to match windowing support for weighting
1059     function.  This avoids a division on every filter call.
1060   */
1061   resize_filter->scale/=resize_filter->window_support;
1062   /*
1063    * Set Cubic Spline B,C values, calculate Cubic coefficients.
1064   */
1065   B=0.0;
1066   C=0.0;
1067   if ((resize_filter->filter == CubicBC) ||
1068       (resize_filter->window == CubicBC) )
1069     {
1070       B=filters[filter_type].B;
1071       C=filters[filter_type].C;
1072       if (filters[window_type].function == CubicBC)
1073         {
1074           B=filters[window_type].B;
1075           C=filters[window_type].C;
1076         }
1077       artifact=GetImageArtifact(image,"filter:b");
1078       if (artifact != (const char *) NULL)
1079         {
1080           B=StringToDouble(artifact,(char **) NULL);
1081           C=(1.0-B)/2.0; /* Calculate C to get a Keys cubic filter. */
1082           artifact=GetImageArtifact(image,"filter:c"); /* user C override */
1083           if (artifact != (const char *) NULL)
1084             C=StringToDouble(artifact,(char **) NULL);
1085         }
1086       else
1087         {
1088           artifact=GetImageArtifact(image,"filter:c");
1089           if (artifact != (const char *) NULL)
1090             {
1091               C=StringToDouble(artifact,(char **) NULL);
1092               B=1.0-2.0*C; /* Calculate B to get a Keys cubic filter. */
1093             }
1094         }
1095       {
1096         const double
1097           twoB = B+B;
1098
1099         /*
1100           Convert B,C values into Cubic Coefficents. See CubicBC().
1101         */
1102         resize_filter->coefficient[0]=1.0-(1.0/3.0)*B;
1103         resize_filter->coefficient[1]=-3.0+twoB+C;
1104         resize_filter->coefficient[2]=2.0-1.5*B-C;
1105         resize_filter->coefficient[3]=(4.0/3.0)*B+4.0*C;
1106         resize_filter->coefficient[4]=-8.0*C-twoB;
1107         resize_filter->coefficient[5]=B+5.0*C;
1108         resize_filter->coefficient[6]=(-1.0/6.0)*B-C;
1109       }
1110     }
1111
1112   /*
1113     Expert Option Request for verbose details of the resulting filter.
1114   */
1115 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1116   #pragma omp master
1117   {
1118 #endif
1119     if (IfStringTrue(GetImageArtifact(image,"filter:verbose")))
1120       {
1121         double
1122           support,
1123           x;
1124
1125         /*
1126           Set the weighting function properly when the weighting function
1127           may not exactly match the filter of the same name.  EG: a Point
1128           filter is really uses a Box weighting function with a different
1129           support than is typically used.
1130         */
1131         if (resize_filter->filter == Box)       filter_type=BoxFilter;
1132         if (resize_filter->filter == Sinc)      filter_type=SincFilter;
1133         if (resize_filter->filter == SincFast)  filter_type=SincFastFilter;
1134         if (resize_filter->filter == Jinc)      filter_type=JincFilter;
1135         if (resize_filter->filter == CubicBC)   filter_type=CubicFilter;
1136         if (resize_filter->window == Box)       window_type=BoxFilter;
1137         if (resize_filter->window == Sinc)      window_type=SincFilter;
1138         if (resize_filter->window == SincFast)  window_type=SincFastFilter;
1139         if (resize_filter->window == Jinc)      window_type=JincFilter;
1140         if (resize_filter->window == CubicBC)   window_type=CubicFilter;
1141         /*
1142           Report Filter Details.
1143         */
1144         support=GetResizeFilterSupport(resize_filter);  /* practical_support */
1145         (void) FormatLocaleFile(stdout,
1146           "# Resampling Filter (for graphing)\n#\n");
1147         (void) FormatLocaleFile(stdout,"# filter = %s\n",
1148           CommandOptionToMnemonic(MagickFilterOptions,filter_type));
1149         (void) FormatLocaleFile(stdout,"# window = %s\n",
1150           CommandOptionToMnemonic(MagickFilterOptions,window_type));
1151         (void) FormatLocaleFile(stdout,"# support = %.*g\n",
1152           GetMagickPrecision(),(double) resize_filter->support);
1153         (void) FormatLocaleFile(stdout,"# window-support = %.*g\n",
1154           GetMagickPrecision(),(double) resize_filter->window_support);
1155         (void) FormatLocaleFile(stdout,"# scale-blur = %.*g\n",
1156           GetMagickPrecision(),(double)resize_filter->blur);
1157         if ((filter_type == GaussianFilter) || (window_type == GaussianFilter))
1158           (void) FormatLocaleFile(stdout,"# gaussian-sigma = %.*g\n",
1159             GetMagickPrecision(),(double)resize_filter->coefficient[0]);
1160         if ( filter_type == KaiserFilter || window_type == KaiserFilter )
1161           (void) FormatLocaleFile(stdout,"# kaiser-beta = %.*g\n",
1162             GetMagickPrecision(),(double)resize_filter->coefficient[0]);
1163         (void) FormatLocaleFile(stdout,"# practical-support = %.*g\n",
1164           GetMagickPrecision(), (double)support);
1165         if ( filter_type == CubicFilter || window_type == CubicFilter )
1166           (void) FormatLocaleFile(stdout,"# B,C = %.*g,%.*g\n",
1167             GetMagickPrecision(),(double)B, GetMagickPrecision(),(double)C);
1168         (void) FormatLocaleFile(stdout,"\n");
1169         /*
1170           Output values of resulting filter graph -- for graphing filter result.
1171         */
1172         for (x=0.0; x <= support; x+=0.01f)
1173           (void) FormatLocaleFile(stdout,"%5.2lf\t%.*g\n",x,
1174             GetMagickPrecision(),(double)
1175             GetResizeFilterWeight(resize_filter,x));
1176         /*
1177           A final value so gnuplot can graph the 'stop' properly.
1178         */
1179         (void) FormatLocaleFile(stdout,"%5.2lf\t%.*g\n",support,
1180           GetMagickPrecision(),0.0);
1181       }
1182       /* Output the above once only for each image - remove setting */
1183     (void) DeleteImageArtifact((Image *) image,"filter:verbose");
1184 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1185   }
1186 #endif
1187   return(resize_filter);
1188 }
1189 \f
1190 /*
1191 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1192 %                                                                             %
1193 %                                                                             %
1194 %                                                                             %
1195 %   A d a p t i v e R e s i z e I m a g e                                     %
1196 %                                                                             %
1197 %                                                                             %
1198 %                                                                             %
1199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1200 %
1201 %  AdaptiveResizeImage() adaptively resize image with pixel resampling.
1202 %
1203 %  This is shortcut function for a fast interpolative resize using mesh
1204 %  interpolation.  It works well for small resizes of less than +/- 50%
1205 %  of the original image size.  For larger resizing on images a full
1206 %  filtered and slower resize function should be used instead.
1207 %
1208 %  The format of the AdaptiveResizeImage method is:
1209 %
1210 %      Image *AdaptiveResizeImage(const Image *image,const size_t columns,
1211 %        const size_t rows,ExceptionInfo *exception)
1212 %
1213 %  A description of each parameter follows:
1214 %
1215 %    o image: the image.
1216 %
1217 %    o columns: the number of columns in the resized image.
1218 %
1219 %    o rows: the number of rows in the resized image.
1220 %
1221 %    o exception: return any errors or warnings in this structure.
1222 %
1223 */
1224 MagickExport Image *AdaptiveResizeImage(const Image *image,
1225   const size_t columns,const size_t rows,ExceptionInfo *exception)
1226 {
1227   Image
1228     *resize_image;
1229
1230   resize_image=InterpolativeResizeImage(image,columns,rows,MeshInterpolatePixel,
1231     exception);
1232   return(resize_image);
1233 }
1234 \f
1235 /*
1236 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1237 %                                                                             %
1238 %                                                                             %
1239 %                                                                             %
1240 +   B e s s e l O r d e r O n e                                               %
1241 %                                                                             %
1242 %                                                                             %
1243 %                                                                             %
1244 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1245 %
1246 %  BesselOrderOne() computes the Bessel function of x of the first kind of
1247 %  order 0.  This is used to create the Jinc() filter function below.
1248 %
1249 %    Reduce x to |x| since j1(x)= -j1(-x), and for x in (0,8]
1250 %
1251 %       j1(x) = x*j1(x);
1252 %
1253 %    For x in (8,inf)
1254 %
1255 %       j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
1256 %
1257 %    where x1 = x-3*pi/4. Compute sin(x1) and cos(x1) as follow:
1258 %
1259 %       cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
1260 %               =  1/sqrt(2) * (sin(x) - cos(x))
1261 %       sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
1262 %               = -1/sqrt(2) * (sin(x) + cos(x))
1263 %
1264 %  The format of the BesselOrderOne method is:
1265 %
1266 %      double BesselOrderOne(double x)
1267 %
1268 %  A description of each parameter follows:
1269 %
1270 %    o x: double value.
1271 %
1272 */
1273
1274 #undef I0
1275 static double I0(double x)
1276 {
1277   double
1278     sum,
1279     t,
1280     y;
1281
1282   register ssize_t
1283     i;
1284
1285   /*
1286     Zeroth order Bessel function of the first kind.
1287   */
1288   sum=1.0;
1289   y=x*x/4.0;
1290   t=y;
1291   for (i=2; t > MagickEpsilon; i++)
1292   {
1293     sum+=t;
1294     t*=y/((double) i*i);
1295   }
1296   return(sum);
1297 }
1298
1299 #undef J1
1300 static double J1(double x)
1301 {
1302   double
1303     p,
1304     q;
1305
1306   register ssize_t
1307     i;
1308
1309   static const double
1310     Pone[] =
1311     {
1312        0.581199354001606143928050809e+21,
1313       -0.6672106568924916298020941484e+20,
1314        0.2316433580634002297931815435e+19,
1315       -0.3588817569910106050743641413e+17,
1316        0.2908795263834775409737601689e+15,
1317       -0.1322983480332126453125473247e+13,
1318        0.3413234182301700539091292655e+10,
1319       -0.4695753530642995859767162166e+7,
1320        0.270112271089232341485679099e+4
1321     },
1322     Qone[] =
1323     {
1324       0.11623987080032122878585294e+22,
1325       0.1185770712190320999837113348e+20,
1326       0.6092061398917521746105196863e+17,
1327       0.2081661221307607351240184229e+15,
1328       0.5243710262167649715406728642e+12,
1329       0.1013863514358673989967045588e+10,
1330       0.1501793594998585505921097578e+7,
1331       0.1606931573481487801970916749e+4,
1332       0.1e+1
1333     };
1334
1335   p=Pone[8];
1336   q=Qone[8];
1337   for (i=7; i >= 0; i--)
1338   {
1339     p=p*x*x+Pone[i];
1340     q=q*x*x+Qone[i];
1341   }
1342   return(p/q);
1343 }
1344
1345 #undef P1
1346 static double P1(double x)
1347 {
1348   double
1349     p,
1350     q;
1351
1352   register ssize_t
1353     i;
1354
1355   static const double
1356     Pone[] =
1357     {
1358       0.352246649133679798341724373e+5,
1359       0.62758845247161281269005675e+5,
1360       0.313539631109159574238669888e+5,
1361       0.49854832060594338434500455e+4,
1362       0.2111529182853962382105718e+3,
1363       0.12571716929145341558495e+1
1364     },
1365     Qone[] =
1366     {
1367       0.352246649133679798068390431e+5,
1368       0.626943469593560511888833731e+5,
1369       0.312404063819041039923015703e+5,
1370       0.4930396490181088979386097e+4,
1371       0.2030775189134759322293574e+3,
1372       0.1e+1
1373     };
1374
1375   p=Pone[5];
1376   q=Qone[5];
1377   for (i=4; i >= 0; i--)
1378   {
1379     p=p*(8.0/x)*(8.0/x)+Pone[i];
1380     q=q*(8.0/x)*(8.0/x)+Qone[i];
1381   }
1382   return(p/q);
1383 }
1384
1385 #undef Q1
1386 static double Q1(double x)
1387 {
1388   double
1389     p,
1390     q;
1391
1392   register ssize_t
1393     i;
1394
1395   static const double
1396     Pone[] =
1397     {
1398       0.3511751914303552822533318e+3,
1399       0.7210391804904475039280863e+3,
1400       0.4259873011654442389886993e+3,
1401       0.831898957673850827325226e+2,
1402       0.45681716295512267064405e+1,
1403       0.3532840052740123642735e-1
1404     },
1405     Qone[] =
1406     {
1407       0.74917374171809127714519505e+4,
1408       0.154141773392650970499848051e+5,
1409       0.91522317015169922705904727e+4,
1410       0.18111867005523513506724158e+4,
1411       0.1038187585462133728776636e+3,
1412       0.1e+1
1413     };
1414
1415   p=Pone[5];
1416   q=Qone[5];
1417   for (i=4; i >= 0; i--)
1418   {
1419     p=p*(8.0/x)*(8.0/x)+Pone[i];
1420     q=q*(8.0/x)*(8.0/x)+Qone[i];
1421   }
1422   return(p/q);
1423 }
1424
1425 static double BesselOrderOne(double x)
1426 {
1427   double
1428     p,
1429     q;
1430
1431   if (x == 0.0)
1432     return(0.0);
1433   p=x;
1434   if (x < 0.0)
1435     x=(-x);
1436   if (x < 8.0)
1437     return(p*J1(x));
1438   q=sqrt((double) (2.0/(MagickPI*x)))*(P1(x)*(1.0/sqrt(2.0)*(sin((double) x)-
1439     cos((double) x)))-8.0/x*Q1(x)*(-1.0/sqrt(2.0)*(sin((double) x)+
1440     cos((double) x))));
1441   if (p < 0.0)
1442     q=(-q);
1443   return(q);
1444 }
1445 \f
1446 /*
1447 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1448 %                                                                             %
1449 %                                                                             %
1450 %                                                                             %
1451 +   D e s t r o y R e s i z e F i l t e r                                     %
1452 %                                                                             %
1453 %                                                                             %
1454 %                                                                             %
1455 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1456 %
1457 %  DestroyResizeFilter() destroy the resize filter.
1458 %
1459 %  The format of the DestroyResizeFilter method is:
1460 %
1461 %      ResizeFilter *DestroyResizeFilter(ResizeFilter *resize_filter)
1462 %
1463 %  A description of each parameter follows:
1464 %
1465 %    o resize_filter: the resize filter.
1466 %
1467 */
1468 MagickPrivate ResizeFilter *DestroyResizeFilter(ResizeFilter *resize_filter)
1469 {
1470   assert(resize_filter != (ResizeFilter *) NULL);
1471   assert(resize_filter->signature == MagickSignature);
1472   resize_filter->signature=(~MagickSignature);
1473   resize_filter=(ResizeFilter *) RelinquishMagickMemory(resize_filter);
1474   return(resize_filter);
1475 }
1476 \f
1477 /*
1478 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1479 %                                                                             %
1480 %                                                                             %
1481 %                                                                             %
1482 +   G e t R e s i z e F i l t e r S u p p o r t                               %
1483 %                                                                             %
1484 %                                                                             %
1485 %                                                                             %
1486 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1487 %
1488 %  GetResizeFilterSupport() return the current support window size for this
1489 %  filter.  Note that this may have been enlarged by filter:blur factor.
1490 %
1491 %  The format of the GetResizeFilterSupport method is:
1492 %
1493 %      double GetResizeFilterSupport(const ResizeFilter *resize_filter)
1494 %
1495 %  A description of each parameter follows:
1496 %
1497 %    o filter: Image filter to use.
1498 %
1499 */
1500
1501 MagickPrivate double *GetResizeFilterCoefficient(
1502   const ResizeFilter *resize_filter)
1503 {
1504   assert(resize_filter != (ResizeFilter *) NULL);
1505   assert(resize_filter->signature == MagickSignature);
1506   return((double *) resize_filter->coefficient);
1507 }
1508
1509 MagickPrivate double GetResizeFilterBlur(
1510   const ResizeFilter *resize_filter)
1511 {
1512   assert(resize_filter != (ResizeFilter *) NULL);
1513   assert(resize_filter->signature == MagickSignature);
1514   return(resize_filter->blur);
1515 }
1516
1517 MagickPrivate double GetResizeFilterScale(
1518   const ResizeFilter *resize_filter)
1519 {
1520   assert(resize_filter != (ResizeFilter *) NULL);
1521   assert(resize_filter->signature == MagickSignature);
1522   return(resize_filter->scale);
1523 }
1524
1525 MagickPrivate double GetResizeFilterWindowSupport(
1526   const ResizeFilter *resize_filter)
1527 {
1528   assert(resize_filter != (ResizeFilter *) NULL);
1529   assert(resize_filter->signature == MagickSignature);
1530   return(resize_filter->window_support);
1531 }
1532
1533 MagickPrivate ResizeWeightingFunctionType GetResizeFilterWeightingType(
1534   const ResizeFilter *resize_filter)
1535 {
1536   assert(resize_filter != (ResizeFilter *) NULL);
1537   assert(resize_filter->signature == MagickSignature);
1538   return(resize_filter->filterWeightingType);
1539 }
1540
1541 MagickPrivate ResizeWeightingFunctionType GetResizeFilterWindowWeightingType(
1542   const ResizeFilter *resize_filter)
1543 {
1544   assert(resize_filter != (ResizeFilter *) NULL);
1545   assert(resize_filter->signature == MagickSignature);
1546   return(resize_filter->windowWeightingType);
1547 }
1548
1549 MagickPrivate double GetResizeFilterSupport(
1550   const ResizeFilter *resize_filter)
1551 {
1552   assert(resize_filter != (ResizeFilter *) NULL);
1553   assert(resize_filter->signature == MagickSignature);
1554   return(resize_filter->support*resize_filter->blur);
1555 }
1556 \f
1557 /*
1558 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1559 %                                                                             %
1560 %                                                                             %
1561 %                                                                             %
1562 +   G e t R e s i z e F i l t e r W e i g h t                                 %
1563 %                                                                             %
1564 %                                                                             %
1565 %                                                                             %
1566 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1567 %
1568 %  GetResizeFilterWeight evaluates the specified resize filter at the point x
1569 %  which usally lies between zero and the filters current 'support' and
1570 %  returns the weight of the filter function at that point.
1571 %
1572 %  The format of the GetResizeFilterWeight method is:
1573 %
1574 %      double GetResizeFilterWeight(const ResizeFilter *resize_filter,
1575 %        const double x)
1576 %
1577 %  A description of each parameter follows:
1578 %
1579 %    o filter: the filter type.
1580 %
1581 %    o x: the point.
1582 %
1583 */
1584 MagickPrivate double GetResizeFilterWeight(const ResizeFilter *resize_filter,
1585   const double x)
1586 {
1587   double
1588     scale,
1589     weight,
1590     x_blur;
1591
1592   /*
1593     Windowing function - scale the weighting filter by this amount.
1594   */
1595   assert(resize_filter != (ResizeFilter *) NULL);
1596   assert(resize_filter->signature == MagickSignature);
1597   x_blur=fabs((double) x)/resize_filter->blur;  /* X offset with blur scaling */
1598   if ((resize_filter->window_support < MagickEpsilon) ||
1599       (resize_filter->window == Box))
1600     scale=1.0;  /* Point or Box Filter -- avoid division by zero */
1601   else
1602     {
1603       scale=resize_filter->scale;
1604       scale=resize_filter->window(x_blur*scale,resize_filter);
1605     }
1606   weight=scale*resize_filter->filter(x_blur,resize_filter);
1607   return(weight);
1608 }
1609 \f
1610 /*
1611 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1612 %                                                                             %
1613 %                                                                             %
1614 %                                                                             %
1615 %   I n t e r p o l a t i v e R e s i z e I m a g e                           %
1616 %                                                                             %
1617 %                                                                             %
1618 %                                                                             %
1619 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1620 %
1621 %  InterpolativeResizeImage() resizes an image using the specified
1622 %  interpolation method.
1623 %
1624 %  The format of the InterpolativeResizeImage method is:
1625 %
1626 %      Image *InterpolativeResizeImage(const Image *image,const size_t columns,
1627 %        const size_t rows,const PixelInterpolateMethod method,
1628 %        ExceptionInfo *exception)
1629 %
1630 %  A description of each parameter follows:
1631 %
1632 %    o image: the image.
1633 %
1634 %    o columns: the number of columns in the resized image.
1635 %
1636 %    o rows: the number of rows in the resized image.
1637 %
1638 %    o method: the pixel interpolation method.
1639 %
1640 %    o exception: return any errors or warnings in this structure.
1641 %
1642 */
1643 MagickExport Image *InterpolativeResizeImage(const Image *image,
1644   const size_t columns,const size_t rows,const PixelInterpolateMethod method,
1645   ExceptionInfo *exception)
1646 {
1647 #define InterpolativeResizeImageTag  "Resize/Image"
1648
1649   CacheView
1650     *image_view,
1651     *resize_view;
1652
1653   Image
1654     *resize_image;
1655
1656   MagickBooleanType
1657     status;
1658
1659   MagickOffsetType
1660     progress;
1661
1662   PointInfo
1663     scale;
1664
1665   ssize_t
1666     y;
1667
1668   /*
1669     Interpolatively resize image.
1670   */
1671   assert(image != (const Image *) NULL);
1672   assert(image->signature == MagickSignature);
1673   if (image->debug != MagickFalse)
1674     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1675   assert(exception != (ExceptionInfo *) NULL);
1676   assert(exception->signature == MagickSignature);
1677   if ((columns == 0) || (rows == 0))
1678     return((Image *) NULL);
1679   if ((columns == image->columns) && (rows == image->rows))
1680     return(CloneImage(image,0,0,MagickTrue,exception));
1681   resize_image=CloneImage(image,columns,rows,MagickTrue,exception);
1682   if (resize_image == (Image *) NULL)
1683     return((Image *) NULL);
1684   if (SetImageStorageClass(resize_image,DirectClass,exception) == MagickFalse)
1685     {
1686       resize_image=DestroyImage(resize_image);
1687       return((Image *) NULL);
1688     }
1689   status=MagickTrue;
1690   progress=0;
1691   image_view=AcquireVirtualCacheView(image,exception);
1692   resize_view=AcquireAuthenticCacheView(resize_image,exception);
1693   scale.x=(double) image->columns/resize_image->columns;
1694   scale.y=(double) image->rows/resize_image->rows;
1695 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1696   #pragma omp parallel for schedule(static,4) shared(progress,status) \
1697     magick_threads(image,resize_image,resize_image->rows,1)
1698 #endif
1699   for (y=0; y < (ssize_t) resize_image->rows; y++)
1700   {
1701     PointInfo
1702       offset;
1703
1704     register Quantum
1705       *restrict q;
1706
1707     register ssize_t
1708       x;
1709
1710     if (status == MagickFalse)
1711       continue;
1712     q=QueueCacheViewAuthenticPixels(resize_view,0,y,resize_image->columns,1,
1713       exception);
1714     if (q == (Quantum *) NULL)
1715       continue;
1716     offset.y=((double) y+0.5)*scale.y-0.5;
1717     for (x=0; x < (ssize_t) resize_image->columns; x++)
1718     {
1719       register ssize_t
1720         i;
1721
1722       if (GetPixelReadMask(resize_image,q) == 0)
1723         {
1724           q+=GetPixelChannels(resize_image);
1725           continue;
1726         }
1727       for (i=0; i < (ssize_t) GetPixelChannels(resize_image); i++)
1728       {
1729         PixelChannel
1730           channel;
1731
1732         PixelTrait
1733           resize_traits,
1734           traits;
1735
1736         channel=GetPixelChannelChannel(image,i);
1737         traits=GetPixelChannelTraits(image,channel);
1738         resize_traits=GetPixelChannelTraits(resize_image,channel);
1739         if ((traits == UndefinedPixelTrait) ||
1740             (resize_traits == UndefinedPixelTrait))
1741           continue;
1742         offset.x=((double) x+0.5)*scale.x-0.5;
1743         status=InterpolatePixelChannels(image,image_view,resize_image,method,
1744           offset.x,offset.y,q,exception);
1745       }
1746       q+=GetPixelChannels(resize_image);
1747     }
1748     if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
1749       status=MagickFalse;
1750     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1751       {
1752         MagickBooleanType
1753           proceed;
1754
1755 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1756         #pragma omp critical (MagickCore_InterpolativeResizeImage)
1757 #endif
1758         proceed=SetImageProgress(image,InterpolativeResizeImageTag,progress++,
1759           image->rows);
1760         if (proceed == MagickFalse)
1761           status=MagickFalse;
1762       }
1763   }
1764   resize_view=DestroyCacheView(resize_view);
1765   image_view=DestroyCacheView(image_view);
1766   if (status == MagickFalse)
1767     resize_image=DestroyImage(resize_image);
1768   return(resize_image);
1769 }
1770 #if defined(MAGICKCORE_LQR_DELEGATE)
1771 \f
1772 /*
1773 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1774 %                                                                             %
1775 %                                                                             %
1776 %                                                                             %
1777 %   L i q u i d R e s c a l e I m a g e                                       %
1778 %                                                                             %
1779 %                                                                             %
1780 %                                                                             %
1781 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1782 %
1783 %  LiquidRescaleImage() rescales image with seam carving.
1784 %
1785 %  The format of the LiquidRescaleImage method is:
1786 %
1787 %      Image *LiquidRescaleImage(const Image *image,const size_t columns,
1788 %        const size_t rows,const double delta_x,const double rigidity,
1789 %        ExceptionInfo *exception)
1790 %
1791 %  A description of each parameter follows:
1792 %
1793 %    o image: the image.
1794 %
1795 %    o columns: the number of columns in the rescaled image.
1796 %
1797 %    o rows: the number of rows in the rescaled image.
1798 %
1799 %    o delta_x: maximum seam transversal step (0 means straight seams).
1800 %
1801 %    o rigidity: introduce a bias for non-straight seams (typically 0).
1802 %
1803 %    o exception: return any errors or warnings in this structure.
1804 %
1805 */
1806 MagickExport Image *LiquidRescaleImage(const Image *image,const size_t columns,
1807   const size_t rows,const double delta_x,const double rigidity,
1808   ExceptionInfo *exception)
1809 {
1810 #define LiquidRescaleImageTag  "Rescale/Image"
1811
1812   CacheView
1813     *image_view,
1814     *rescale_view;
1815
1816   gfloat
1817     *packet,
1818     *pixels;
1819
1820   Image
1821     *rescale_image;
1822
1823   int
1824     x_offset,
1825     y_offset;
1826
1827   LqrCarver
1828     *carver;
1829
1830   LqrRetVal
1831     lqr_status;
1832
1833   MagickBooleanType
1834     status;
1835
1836   MemoryInfo
1837     *pixel_info;
1838
1839   register gfloat
1840     *q;
1841
1842   ssize_t
1843     y;
1844
1845   /*
1846     Liquid rescale image.
1847   */
1848   assert(image != (const Image *) NULL);
1849   assert(image->signature == MagickSignature);
1850   if (image->debug != MagickFalse)
1851     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1852   assert(exception != (ExceptionInfo *) NULL);
1853   assert(exception->signature == MagickSignature);
1854   if ((columns == 0) || (rows == 0))
1855     return((Image *) NULL);
1856   if ((columns == image->columns) && (rows == image->rows))
1857     return(CloneImage(image,0,0,MagickTrue,exception));
1858   if ((columns <= 2) || (rows <= 2))
1859     return(ResizeImage(image,columns,rows,image->filter,exception));
1860   if ((columns >= (2*image->columns)) || (rows >= (2*image->rows)))
1861     {
1862       Image
1863         *resize_image;
1864
1865       size_t
1866         height,
1867         width;
1868
1869       /*
1870         Honor liquid resize size limitations.
1871       */
1872       for (width=image->columns; columns >= (2*width-1); width*=2) ;
1873       for (height=image->rows; rows >= (2*height-1); height*=2) ;
1874       resize_image=ResizeImage(image,width,height,image->filter,exception);
1875       if (resize_image == (Image *) NULL)
1876         return((Image *) NULL);
1877       rescale_image=LiquidRescaleImage(resize_image,columns,rows,delta_x,
1878         rigidity,exception);
1879       resize_image=DestroyImage(resize_image);
1880       return(rescale_image);
1881     }
1882   pixel_info=AcquireVirtualMemory(image->columns,image->rows*
1883     GetPixelChannels(image)*sizeof(*pixels));
1884   if (pixel_info == (MemoryInfo *) NULL)
1885     return((Image *) NULL);
1886   pixels=(gfloat *) GetVirtualMemoryBlob(pixel_info);
1887   status=MagickTrue;
1888   q=pixels;
1889   image_view=AcquireVirtualCacheView(image,exception);
1890   for (y=0; y < (ssize_t) image->rows; y++)
1891   {
1892     register const Quantum
1893       *restrict p;
1894
1895     register ssize_t
1896       x;
1897
1898     if (status == MagickFalse)
1899       continue;
1900     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1901     if (p == (const Quantum *) NULL)
1902       {
1903         status=MagickFalse;
1904         continue;
1905       }
1906     for (x=0; x < (ssize_t) image->columns; x++)
1907     {
1908       register ssize_t
1909         i;
1910
1911       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1912         *q++=QuantumScale*p[i];
1913       p+=GetPixelChannels(image);
1914     }
1915   }
1916   image_view=DestroyCacheView(image_view);
1917   carver=lqr_carver_new_ext(pixels,(int) image->columns,(int) image->rows,
1918     (int) GetPixelChannels(image),LQR_COLDEPTH_32F);
1919   if (carver == (LqrCarver *) NULL)
1920     {
1921       pixel_info=RelinquishVirtualMemory(pixel_info);
1922       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1923     }
1924   lqr_carver_set_preserve_input_image(carver);
1925   lqr_status=lqr_carver_init(carver,(int) delta_x,rigidity);
1926   lqr_status=lqr_carver_resize(carver,(int) columns,(int) rows);
1927   (void) lqr_status;
1928   rescale_image=CloneImage(image,lqr_carver_get_width(carver),
1929     lqr_carver_get_height(carver),MagickTrue,exception);
1930   if (rescale_image == (Image *) NULL)
1931     {
1932       pixel_info=RelinquishVirtualMemory(pixel_info);
1933       return((Image *) NULL);
1934     }
1935   if (SetImageStorageClass(rescale_image,DirectClass,exception) == MagickFalse)
1936     {
1937       pixel_info=RelinquishVirtualMemory(pixel_info);
1938       rescale_image=DestroyImage(rescale_image);
1939       return((Image *) NULL);
1940     }
1941   rescale_view=AcquireAuthenticCacheView(rescale_image,exception);
1942   (void) lqr_carver_scan_reset(carver);
1943   while (lqr_carver_scan_ext(carver,&x_offset,&y_offset,(void **) &packet) != 0)
1944   {
1945     register Quantum
1946       *restrict q;
1947
1948     register ssize_t
1949       i;
1950
1951     q=QueueCacheViewAuthenticPixels(rescale_view,x_offset,y_offset,1,1,
1952       exception);
1953     if (q == (Quantum *) NULL)
1954       break;
1955     for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1956     {
1957       PixelChannel
1958         channel;
1959
1960       PixelTrait
1961         rescale_traits,
1962         traits;
1963
1964       channel=GetPixelChannelChannel(image,i);
1965       traits=GetPixelChannelTraits(image,channel);
1966       rescale_traits=GetPixelChannelTraits(rescale_image,channel);
1967       if ((traits == UndefinedPixelTrait) ||
1968           (rescale_traits == UndefinedPixelTrait))
1969         continue;
1970       SetPixelChannel(rescale_image,channel,ClampToQuantum(QuantumRange*
1971         packet[i]),q);
1972     }
1973     if (SyncCacheViewAuthenticPixels(rescale_view,exception) == MagickFalse)
1974       break;
1975   }
1976   rescale_view=DestroyCacheView(rescale_view);
1977   pixel_info=RelinquishVirtualMemory(pixel_info);
1978   lqr_carver_destroy(carver);
1979   return(rescale_image);
1980 }
1981 #else
1982 MagickExport Image *LiquidRescaleImage(const Image *image,
1983   const size_t magick_unused(columns),const size_t magick_unused(rows),
1984   const double magick_unused(delta_x),const double magick_unused(rigidity),
1985   ExceptionInfo *exception)
1986 {
1987   assert(image != (const Image *) NULL);
1988   assert(image->signature == MagickSignature);
1989   if (image->debug != MagickFalse)
1990     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1991   assert(exception != (ExceptionInfo *) NULL);
1992   assert(exception->signature == MagickSignature);
1993   (void) ThrowMagickException(exception,GetMagickModule(),MissingDelegateError,
1994     "DelegateLibrarySupportNotBuiltIn","'%s' (LQR)",image->filename);
1995   return((Image *) NULL);
1996 }
1997 #endif
1998 \f
1999 /*
2000 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2001 %                                                                             %
2002 %                                                                             %
2003 %                                                                             %
2004 %   M a g n i f y I m a g e                                                   %
2005 %                                                                             %
2006 %                                                                             %
2007 %                                                                             %
2008 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2009 %
2010 %  MagnifyImage() doubles the size of the image with a pixel art scaling
2011 %  algorithm.
2012 %
2013 %  The format of the MagnifyImage method is:
2014 %
2015 %      Image *MagnifyImage(const Image *image,ExceptionInfo *exception)
2016 %
2017 %  A description of each parameter follows:
2018 %
2019 %    o image: the image.
2020 %
2021 %    o exception: return any errors or warnings in this structure.
2022 %
2023 */
2024 MagickExport Image *MagnifyImage(const Image *image,ExceptionInfo *exception)
2025 {
2026 #define MagnifyImageTag  "Magnify/Image"
2027
2028   CacheView
2029     *image_view,
2030     *magnify_view;
2031
2032   Image
2033     *magnify_image;
2034
2035   MagickBooleanType
2036     status;
2037
2038   MagickOffsetType
2039     progress;
2040
2041   ssize_t
2042     y;
2043
2044   /*
2045     Initialize magnified image attributes.
2046   */
2047   assert(image != (const Image *) NULL);
2048   assert(image->signature == MagickSignature);
2049   if (image->debug != MagickFalse)
2050     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2051   assert(exception != (ExceptionInfo *) NULL);
2052   assert(exception->signature == MagickSignature);
2053   magnify_image=CloneImage(image,2*image->columns,2*image->rows,MagickTrue,
2054     exception);
2055   if (magnify_image == (Image *) NULL)
2056     return((Image *) NULL);
2057   /*
2058     Magnify image.
2059   */
2060   status=MagickTrue;
2061   progress=0;
2062   image_view=AcquireVirtualCacheView(image,exception);
2063   magnify_view=AcquireAuthenticCacheView(magnify_image,exception);
2064 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2065   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2066     magick_threads(image,magnify_image,image->rows,1)
2067 #endif
2068   for (y=0; y < (ssize_t) image->rows; y++)
2069   {
2070     register Quantum
2071       *restrict q;
2072
2073     register ssize_t
2074       x;
2075
2076     if (status == MagickFalse)
2077       continue;
2078     q=QueueCacheViewAuthenticPixels(magnify_view,0,2*y,magnify_image->columns,2,
2079       exception);
2080     if (q == (Quantum *) NULL)
2081       {
2082         status=MagickFalse;
2083         continue;
2084       }
2085     /*
2086       Magnify this row of pixels.
2087     */
2088     for (x=0; x < (ssize_t) image->columns; x++)
2089     {
2090       MagickRealType
2091         intensity[9];
2092
2093       register const Quantum
2094         *restrict p;
2095
2096       register Quantum
2097         *restrict r;
2098
2099       register ssize_t
2100         i;
2101
2102       size_t
2103         channels;
2104
2105       p=GetCacheViewVirtualPixels(image_view,x-1,y-1,3,3,exception);
2106       if (p == (const Quantum *) NULL)
2107         {
2108           status=MagickFalse;
2109           continue;
2110         }
2111       channels=GetPixelChannels(image);
2112       for (i=0; i < 9; i++)
2113         intensity[i]=GetPixelIntensity(image,p+i*channels);
2114       r=q;
2115       if ((fabs(intensity[1]-intensity[7]) < MagickEpsilon) ||
2116           (fabs(intensity[3]-intensity[5]) < MagickEpsilon))
2117         {
2118           /*
2119             Clone center pixel.
2120           */
2121           for (i=0; i < (ssize_t) channels; i++)
2122             r[i]=p[4*channels+i];
2123           r+=GetPixelChannels(magnify_image);
2124           for (i=0; i < (ssize_t) channels; i++)
2125             r[i]=p[4*channels+i];
2126           r+=(magnify_image->columns-1)*GetPixelChannels(magnify_image);
2127           for (i=0; i < (ssize_t) channels; i++)
2128             r[i]=p[4*channels+i];
2129           r+=GetPixelChannels(magnify_image);
2130           for (i=0; i < (ssize_t) channels; i++)
2131             r[i]=p[4*channels+i];
2132         }
2133       else
2134         {
2135           /*
2136             Selectively clone pixel.
2137           */
2138           if (fabs(intensity[1]-intensity[3]) < MagickEpsilon)
2139             for (i=0; i < (ssize_t) channels; i++)
2140               r[i]=p[3*channels+i];
2141           else
2142             for (i=0; i < (ssize_t) channels; i++)
2143               r[i]=p[4*channels+i];
2144           r+=GetPixelChannels(magnify_image);
2145           if (fabs(intensity[1]-intensity[5]) < MagickEpsilon)
2146             for (i=0; i < (ssize_t) channels; i++)
2147               r[i]=p[5*channels+i];
2148           else
2149             for (i=0; i < (ssize_t) channels; i++)
2150               r[i]=p[4*channels+i];
2151           r+=(magnify_image->columns-1)*GetPixelChannels(magnify_image);
2152           if (fabs(intensity[3]-intensity[7]) < MagickEpsilon)
2153             for (i=0; i < (ssize_t) channels; i++)
2154               r[i]=p[3*channels+i];
2155           else
2156             for (i=0; i < (ssize_t) channels; i++)
2157               r[i]=p[4*channels+i];
2158           r+=GetPixelChannels(magnify_image);
2159           if (fabs(intensity[5]-intensity[7]) < MagickEpsilon)
2160             for (i=0; i < (ssize_t) channels; i++)
2161               r[i]=p[5*channels+i];
2162           else
2163             for (i=0; i < (ssize_t) channels; i++)
2164               r[i]=p[4*channels+i];
2165         }
2166       q+=2*GetPixelChannels(magnify_image);
2167     }
2168     if (SyncCacheViewAuthenticPixels(magnify_view,exception) == MagickFalse)
2169       status=MagickFalse;
2170     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2171       {
2172         MagickBooleanType
2173           proceed;
2174
2175 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2176         #pragma omp critical (MagickCore_MagnifyImage)
2177 #endif
2178         proceed=SetImageProgress(image,MagnifyImageTag,progress++,image->rows);
2179         if (proceed == MagickFalse)
2180           status=MagickFalse;
2181       }
2182   }
2183   magnify_view=DestroyCacheView(magnify_view);
2184   image_view=DestroyCacheView(image_view);
2185   if (status == MagickFalse)
2186     magnify_image=DestroyImage(magnify_image);
2187   return(magnify_image);
2188 }
2189 \f
2190 /*
2191 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2192 %                                                                             %
2193 %                                                                             %
2194 %                                                                             %
2195 %   M i n i f y I m a g e                                                     %
2196 %                                                                             %
2197 %                                                                             %
2198 %                                                                             %
2199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2200 %
2201 %  MinifyImage() is a convenience method that scales an image proportionally to
2202 %  half its size.
2203 %
2204 %  The format of the MinifyImage method is:
2205 %
2206 %      Image *MinifyImage(const Image *image,ExceptionInfo *exception)
2207 %
2208 %  A description of each parameter follows:
2209 %
2210 %    o image: the image.
2211 %
2212 %    o exception: return any errors or warnings in this structure.
2213 %
2214 */
2215 MagickExport Image *MinifyImage(const Image *image,ExceptionInfo *exception)
2216 {
2217   Image
2218     *minify_image;
2219
2220   assert(image != (Image *) NULL);
2221   assert(image->signature == MagickSignature);
2222   if (image->debug != MagickFalse)
2223     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2224   assert(exception != (ExceptionInfo *) NULL);
2225   assert(exception->signature == MagickSignature);
2226   minify_image=ResizeImage(image,image->columns/2,image->rows/2,SplineFilter,
2227     exception);
2228   return(minify_image);
2229 }
2230 \f
2231 /*
2232 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2233 %                                                                             %
2234 %                                                                             %
2235 %                                                                             %
2236 %   R e s a m p l e I m a g e                                                 %
2237 %                                                                             %
2238 %                                                                             %
2239 %                                                                             %
2240 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2241 %
2242 %  ResampleImage() resize image in terms of its pixel size, so that when
2243 %  displayed at the given resolution it will be the same size in terms of
2244 %  real world units as the original image at the original resolution.
2245 %
2246 %  The format of the ResampleImage method is:
2247 %
2248 %      Image *ResampleImage(Image *image,const double x_resolution,
2249 %        const double y_resolution,const FilterTypes filter,
2250 %        ExceptionInfo *exception)
2251 %
2252 %  A description of each parameter follows:
2253 %
2254 %    o image: the image to be resized to fit the given resolution.
2255 %
2256 %    o x_resolution: the new image x resolution.
2257 %
2258 %    o y_resolution: the new image y resolution.
2259 %
2260 %    o filter: Image filter to use.
2261 %
2262 %    o exception: return any errors or warnings in this structure.
2263 %
2264 */
2265 MagickExport Image *ResampleImage(const Image *image,const double x_resolution,
2266   const double y_resolution,const FilterTypes filter,ExceptionInfo *exception)
2267 {
2268 #define ResampleImageTag  "Resample/Image"
2269
2270   Image
2271     *resample_image;
2272
2273   size_t
2274     height,
2275     width;
2276
2277   /*
2278     Initialize sampled image attributes.
2279   */
2280   assert(image != (const Image *) NULL);
2281   assert(image->signature == MagickSignature);
2282   if (image->debug != MagickFalse)
2283     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2284   assert(exception != (ExceptionInfo *) NULL);
2285   assert(exception->signature == MagickSignature);
2286   width=(size_t) (x_resolution*image->columns/(image->resolution.x == 0.0 ?
2287     72.0 : image->resolution.x)+0.5);
2288   height=(size_t) (y_resolution*image->rows/(image->resolution.y == 0.0 ?
2289     72.0 : image->resolution.y)+0.5);
2290   resample_image=ResizeImage(image,width,height,filter,exception);
2291   if (resample_image != (Image *) NULL)
2292     {
2293       resample_image->resolution.x=x_resolution;
2294       resample_image->resolution.y=y_resolution;
2295     }
2296   return(resample_image);
2297 }
2298 \f
2299 /*
2300 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2301 %                                                                             %
2302 %                                                                             %
2303 %                                                                             %
2304 %   R e s i z e I m a g e                                                     %
2305 %                                                                             %
2306 %                                                                             %
2307 %                                                                             %
2308 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2309 %
2310 %  ResizeImage() scales an image to the desired dimensions, using the given
2311 %  filter (see AcquireFilterInfo()).
2312 %
2313 %  If an undefined filter is given the filter defaults to Mitchell for a
2314 %  colormapped image, a image with a matte channel, or if the image is
2315 %  enlarged.  Otherwise the filter defaults to a Lanczos.
2316 %
2317 %  ResizeImage() was inspired by Paul Heckbert's "zoom" program.
2318 %
2319 %  The format of the ResizeImage method is:
2320 %
2321 %      Image *ResizeImage(Image *image,const size_t columns,const size_t rows,
2322 %        const FilterTypes filter,ExceptionInfo *exception)
2323 %
2324 %  A description of each parameter follows:
2325 %
2326 %    o image: the image.
2327 %
2328 %    o columns: the number of columns in the scaled image.
2329 %
2330 %    o rows: the number of rows in the scaled image.
2331 %
2332 %    o filter: Image filter to use.
2333 %
2334 %    o exception: return any errors or warnings in this structure.
2335 %
2336 */
2337
2338 typedef struct _ContributionInfo
2339 {
2340   double
2341     weight;
2342
2343   ssize_t
2344     pixel;
2345 } ContributionInfo;
2346
2347 static ContributionInfo **DestroyContributionThreadSet(
2348   ContributionInfo **contribution)
2349 {
2350   register ssize_t
2351     i;
2352
2353   assert(contribution != (ContributionInfo **) NULL);
2354   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2355     if (contribution[i] != (ContributionInfo *) NULL)
2356       contribution[i]=(ContributionInfo *) RelinquishAlignedMemory(
2357         contribution[i]);
2358   contribution=(ContributionInfo **) RelinquishMagickMemory(contribution);
2359   return(contribution);
2360 }
2361
2362 static ContributionInfo **AcquireContributionThreadSet(const size_t count)
2363 {
2364   register ssize_t
2365     i;
2366
2367   ContributionInfo
2368     **contribution;
2369
2370   size_t
2371     number_threads;
2372
2373   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2374   contribution=(ContributionInfo **) AcquireQuantumMemory(number_threads,
2375     sizeof(*contribution));
2376   if (contribution == (ContributionInfo **) NULL)
2377     return((ContributionInfo **) NULL);
2378   (void) ResetMagickMemory(contribution,0,number_threads*sizeof(*contribution));
2379   for (i=0; i < (ssize_t) number_threads; i++)
2380   {
2381     contribution[i]=(ContributionInfo *) MagickAssumeAligned(
2382       AcquireAlignedMemory(count,sizeof(**contribution)));
2383     if (contribution[i] == (ContributionInfo *) NULL)
2384       return(DestroyContributionThreadSet(contribution));
2385   }
2386   return(contribution);
2387 }
2388
2389 static inline double MagickMax(const double x,const double y)
2390 {
2391   if (x > y)
2392     return(x);
2393   return(y);
2394 }
2395
2396 static inline double MagickMin(const double x,const double y)
2397 {
2398   if (x < y)
2399     return(x);
2400   return(y);
2401 }
2402
2403 static MagickBooleanType HorizontalFilter(const ResizeFilter *resize_filter,
2404   const Image *image,Image *resize_image,const double x_factor,
2405   const MagickSizeType span,MagickOffsetType *offset,ExceptionInfo *exception)
2406 {
2407 #define ResizeImageTag  "Resize/Image"
2408
2409   CacheView
2410     *image_view,
2411     *resize_view;
2412
2413   ClassType
2414     storage_class;
2415
2416   ContributionInfo
2417     **restrict contributions;
2418
2419   MagickBooleanType
2420     status;
2421
2422   double
2423     scale,
2424     support;
2425
2426   ssize_t
2427     x;
2428
2429   /*
2430     Apply filter to resize horizontally from image to resize image.
2431   */
2432   scale=MagickMax(1.0/x_factor+MagickEpsilon,1.0);
2433   support=scale*GetResizeFilterSupport(resize_filter);
2434   storage_class=support > 0.5 ? DirectClass : image->storage_class;
2435   if (SetImageStorageClass(resize_image,storage_class,exception) == MagickFalse)
2436     return(MagickFalse);
2437   if (support < 0.5)
2438     {
2439       /*
2440         Support too small even for nearest neighbour: Reduce to point sampling.
2441       */
2442       support=(double) 0.5;
2443       scale=1.0;
2444     }
2445   contributions=AcquireContributionThreadSet((size_t) (2.0*support+3.0));
2446   if (contributions == (ContributionInfo **) NULL)
2447     {
2448       (void) ThrowMagickException(exception,GetMagickModule(),
2449         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2450       return(MagickFalse);
2451     }
2452   status=MagickTrue;
2453   scale=PerceptibleReciprocal(scale);
2454   image_view=AcquireVirtualCacheView(image,exception);
2455   resize_view=AcquireAuthenticCacheView(resize_image,exception);
2456 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2457   #pragma omp parallel for schedule(static,4) shared(status) \
2458     magick_threads(image,resize_image,resize_image->columns,1)
2459 #endif
2460   for (x=0; x < (ssize_t) resize_image->columns; x++)
2461   {
2462     const int
2463       id = GetOpenMPThreadId();
2464
2465     double
2466       bisect,
2467       density;
2468
2469     register const Quantum
2470       *restrict p;
2471
2472     register ContributionInfo
2473       *restrict contribution;
2474
2475     register Quantum
2476       *restrict q;
2477
2478     register ssize_t
2479       y;
2480
2481     ssize_t
2482       n,
2483       start,
2484       stop;
2485
2486     if (status == MagickFalse)
2487       continue;
2488     bisect=(double) (x+0.5)/x_factor+MagickEpsilon;
2489     start=(ssize_t) MagickMax(bisect-support+0.5,0.0);
2490     stop=(ssize_t) MagickMin(bisect+support+0.5,(double) image->columns);
2491     density=0.0;
2492     contribution=contributions[id];
2493     for (n=0; n < (stop-start); n++)
2494     {
2495       contribution[n].pixel=start+n;
2496       contribution[n].weight=GetResizeFilterWeight(resize_filter,scale*
2497         ((double) (start+n)-bisect+0.5));
2498       density+=contribution[n].weight;
2499     }
2500     if ((density != 0.0) && (density != 1.0))
2501       {
2502         register ssize_t
2503           i;
2504
2505         /*
2506           Normalize.
2507         */
2508         density=PerceptibleReciprocal(density);
2509         for (i=0; i < n; i++)
2510           contribution[i].weight*=density;
2511       }
2512     p=GetCacheViewVirtualPixels(image_view,contribution[0].pixel,0,(size_t)
2513       (contribution[n-1].pixel-contribution[0].pixel+1),image->rows,exception);
2514     q=QueueCacheViewAuthenticPixels(resize_view,x,0,1,resize_image->rows,
2515       exception);
2516     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2517       {
2518         status=MagickFalse;
2519         continue;
2520       }
2521     for (y=0; y < (ssize_t) resize_image->rows; y++)
2522     {
2523       register ssize_t
2524         i;
2525
2526       for (i=0; i < (ssize_t) GetPixelChannels(resize_image); i++)
2527       {
2528         double
2529           alpha,
2530           gamma,
2531           pixel;
2532
2533         PixelChannel
2534           channel;
2535
2536         PixelTrait
2537           resize_traits,
2538           traits;
2539
2540         register ssize_t
2541           j;
2542
2543         ssize_t
2544           k;
2545
2546         channel=GetPixelChannelChannel(image,i);
2547         traits=GetPixelChannelTraits(image,channel);
2548         resize_traits=GetPixelChannelTraits(resize_image,channel);
2549         if ((traits == UndefinedPixelTrait) ||
2550             (resize_traits == UndefinedPixelTrait))
2551           continue;
2552         if (((resize_traits & CopyPixelTrait) != 0) ||
2553             (GetPixelReadMask(resize_image,q) == 0))
2554           {
2555             j=(ssize_t) (MagickMin(MagickMax(bisect,(double) start),(double)
2556               stop-1.0)+0.5);
2557             k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2558               (contribution[j-start].pixel-contribution[0].pixel);
2559             SetPixelChannel(resize_image,channel,p[k*GetPixelChannels(image)+i],
2560               q);
2561             continue;
2562           }
2563         pixel=0.0;
2564         if ((resize_traits & BlendPixelTrait) == 0)
2565           {
2566             /*
2567               No alpha blending.
2568             */
2569             for (j=0; j < n; j++)
2570             {
2571               k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2572                 (contribution[j].pixel-contribution[0].pixel);
2573               alpha=contribution[j].weight;
2574               pixel+=alpha*p[k*GetPixelChannels(image)+i];
2575             }
2576             SetPixelChannel(resize_image,channel,ClampToQuantum(pixel),q);
2577             continue;
2578           }
2579         /*
2580           Alpha blending.
2581         */
2582         gamma=0.0;
2583         for (j=0; j < n; j++)
2584         {
2585           k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2586             (contribution[j].pixel-contribution[0].pixel);
2587           alpha=contribution[j].weight*QuantumScale*
2588             GetPixelAlpha(image,p+k*GetPixelChannels(image));
2589           pixel+=alpha*p[k*GetPixelChannels(image)+i];
2590           gamma+=alpha;
2591         }
2592         gamma=PerceptibleReciprocal(gamma);
2593         SetPixelChannel(resize_image,channel,ClampToQuantum(gamma*pixel),q);
2594       }
2595       q+=GetPixelChannels(resize_image);
2596     }
2597     if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
2598       status=MagickFalse;
2599     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2600       {
2601         MagickBooleanType
2602           proceed;
2603
2604 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2605         #pragma omp critical (MagickCore_HorizontalFilter)
2606 #endif
2607         proceed=SetImageProgress(image,ResizeImageTag,(*offset)++,span);
2608         if (proceed == MagickFalse)
2609           status=MagickFalse;
2610       }
2611   }
2612   resize_view=DestroyCacheView(resize_view);
2613   image_view=DestroyCacheView(image_view);
2614   contributions=DestroyContributionThreadSet(contributions);
2615   return(status);
2616 }
2617
2618 static MagickBooleanType VerticalFilter(const ResizeFilter *resize_filter,
2619   const Image *image,Image *resize_image,const double y_factor,
2620   const MagickSizeType span,MagickOffsetType *offset,ExceptionInfo *exception)
2621 {
2622   CacheView
2623     *image_view,
2624     *resize_view;
2625
2626   ClassType
2627     storage_class;
2628
2629   ContributionInfo
2630     **restrict contributions;
2631
2632   double
2633     scale,
2634     support;
2635
2636   MagickBooleanType
2637     status;
2638
2639   ssize_t
2640     y;
2641
2642   /*
2643     Apply filter to resize vertically from image to resize image.
2644   */
2645   scale=MagickMax(1.0/y_factor+MagickEpsilon,1.0);
2646   support=scale*GetResizeFilterSupport(resize_filter);
2647   storage_class=support > 0.5 ? DirectClass : image->storage_class;
2648   if (SetImageStorageClass(resize_image,storage_class,exception) == MagickFalse)
2649     return(MagickFalse);
2650   if (support < 0.5)
2651     {
2652       /*
2653         Support too small even for nearest neighbour: Reduce to point sampling.
2654       */
2655       support=(double) 0.5;
2656       scale=1.0;
2657     }
2658   contributions=AcquireContributionThreadSet((size_t) (2.0*support+3.0));
2659   if (contributions == (ContributionInfo **) NULL)
2660     {
2661       (void) ThrowMagickException(exception,GetMagickModule(),
2662         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2663       return(MagickFalse);
2664     }
2665   status=MagickTrue;
2666   scale=PerceptibleReciprocal(scale);
2667   image_view=AcquireVirtualCacheView(image,exception);
2668   resize_view=AcquireAuthenticCacheView(resize_image,exception);
2669 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2670   #pragma omp parallel for schedule(static,4) shared(status) \
2671     magick_threads(image,resize_image,resize_image->rows,1)
2672 #endif
2673   for (y=0; y < (ssize_t) resize_image->rows; y++)
2674   {
2675     const int
2676       id = GetOpenMPThreadId();
2677
2678     double
2679       bisect,
2680       density;
2681
2682     register const Quantum
2683       *restrict p;
2684
2685     register ContributionInfo
2686       *restrict contribution;
2687
2688     register Quantum
2689       *restrict q;
2690
2691     register ssize_t
2692       x;
2693
2694     ssize_t
2695       n,
2696       start,
2697       stop;
2698
2699     if (status == MagickFalse)
2700       continue;
2701     bisect=(double) (y+0.5)/y_factor+MagickEpsilon;
2702     start=(ssize_t) MagickMax(bisect-support+0.5,0.0);
2703     stop=(ssize_t) MagickMin(bisect+support+0.5,(double) image->rows);
2704     density=0.0;
2705     contribution=contributions[id];
2706     for (n=0; n < (stop-start); n++)
2707     {
2708       contribution[n].pixel=start+n;
2709       contribution[n].weight=GetResizeFilterWeight(resize_filter,scale*
2710         ((double) (start+n)-bisect+0.5));
2711       density+=contribution[n].weight;
2712     }
2713     if ((density != 0.0) && (density != 1.0))
2714       {
2715         register ssize_t
2716           i;
2717
2718         /*
2719           Normalize.
2720         */
2721         density=PerceptibleReciprocal(density);
2722         for (i=0; i < n; i++)
2723           contribution[i].weight*=density;
2724       }
2725     p=GetCacheViewVirtualPixels(image_view,0,contribution[0].pixel,
2726       image->columns,(size_t) (contribution[n-1].pixel-contribution[0].pixel+1),
2727       exception);
2728     q=QueueCacheViewAuthenticPixels(resize_view,0,y,resize_image->columns,1,
2729       exception);
2730     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2731       {
2732         status=MagickFalse;
2733         continue;
2734       }
2735     for (x=0; x < (ssize_t) resize_image->columns; x++)
2736     {
2737       register ssize_t
2738         i;
2739
2740       for (i=0; i < (ssize_t) GetPixelChannels(resize_image); i++)
2741       {
2742         double
2743           alpha,
2744           gamma,
2745           pixel;
2746
2747         PixelChannel
2748           channel;
2749
2750         PixelTrait
2751           resize_traits,
2752           traits;
2753
2754         register ssize_t
2755           j;
2756
2757         ssize_t
2758           k;
2759
2760         channel=GetPixelChannelChannel(image,i);
2761         traits=GetPixelChannelTraits(image,channel);
2762         resize_traits=GetPixelChannelTraits(resize_image,channel);
2763         if ((traits == UndefinedPixelTrait) ||
2764             (resize_traits == UndefinedPixelTrait))
2765           continue;
2766         if (((resize_traits & CopyPixelTrait) != 0) ||
2767             (GetPixelReadMask(resize_image,q) == 0))
2768           {
2769             j=(ssize_t) (MagickMin(MagickMax(bisect,(double) start),(double)
2770               stop-1.0)+0.5);
2771             k=(ssize_t) ((contribution[j-start].pixel-contribution[0].pixel)*
2772               image->columns+x);
2773             SetPixelChannel(resize_image,channel,p[k*GetPixelChannels(image)+i],
2774               q);
2775             continue;
2776           }
2777         pixel=0.0;
2778         if ((resize_traits & BlendPixelTrait) == 0)
2779           {
2780             /*
2781               No alpha blending.
2782             */
2783             for (j=0; j < n; j++)
2784             {
2785               k=(ssize_t) ((contribution[j].pixel-contribution[0].pixel)*
2786                 image->columns+x);
2787               alpha=contribution[j].weight;
2788               pixel+=alpha*p[k*GetPixelChannels(image)+i];
2789             }
2790             SetPixelChannel(resize_image,channel,ClampToQuantum(pixel),q);
2791             continue;
2792           }
2793         gamma=0.0;
2794         for (j=0; j < n; j++)
2795         {
2796           k=(ssize_t) ((contribution[j].pixel-contribution[0].pixel)*
2797             image->columns+x);
2798           alpha=contribution[j].weight*QuantumScale*GetPixelAlpha(image,p+k*
2799             GetPixelChannels(image));
2800           pixel+=alpha*p[k*GetPixelChannels(image)+i];
2801           gamma+=alpha;
2802         }
2803         gamma=PerceptibleReciprocal(gamma);
2804         SetPixelChannel(resize_image,channel,ClampToQuantum(gamma*pixel),q);
2805       }
2806       q+=GetPixelChannels(resize_image);
2807     }
2808     if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
2809       status=MagickFalse;
2810     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2811       {
2812         MagickBooleanType
2813           proceed;
2814
2815 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2816         #pragma omp critical (MagickCore_VerticalFilter)
2817 #endif
2818         proceed=SetImageProgress(image,ResizeImageTag,(*offset)++,span);
2819         if (proceed == MagickFalse)
2820           status=MagickFalse;
2821       }
2822   }
2823   resize_view=DestroyCacheView(resize_view);
2824   image_view=DestroyCacheView(image_view);
2825   contributions=DestroyContributionThreadSet(contributions);
2826   return(status);
2827 }
2828
2829 MagickExport Image *ResizeImage(const Image *image,const size_t columns,
2830   const size_t rows,const FilterTypes filter,ExceptionInfo *exception)
2831 {
2832   double
2833     x_factor,
2834     y_factor;
2835
2836   FilterTypes
2837     filter_type;
2838
2839   Image
2840     *filter_image,
2841     *resize_image;
2842
2843   MagickOffsetType
2844     offset;
2845
2846   MagickSizeType
2847     span;
2848
2849   MagickStatusType
2850     status;
2851
2852   ResizeFilter
2853     *resize_filter;
2854
2855   /*
2856     Acquire resize image.
2857   */
2858   assert(image != (Image *) NULL);
2859   assert(image->signature == MagickSignature);
2860   if (image->debug != MagickFalse)
2861     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2862   assert(exception != (ExceptionInfo *) NULL);
2863   assert(exception->signature == MagickSignature);
2864   if ((columns == 0) || (rows == 0))
2865     ThrowImageException(ImageError,"NegativeOrZeroImageSize");
2866   if ((columns == image->columns) && (rows == image->rows) &&
2867       (filter == UndefinedFilter))
2868     return(CloneImage(image,0,0,MagickTrue,exception));
2869   resize_image=CloneImage(image,columns,rows,MagickTrue,exception);
2870   if (resize_image == (Image *) NULL)
2871     return(resize_image);
2872   /*
2873     Acquire resize filter.
2874   */
2875   x_factor=(double) columns/(double) image->columns;
2876   y_factor=(double) rows/(double) image->rows;
2877   if (x_factor > y_factor)
2878     filter_image=CloneImage(image,columns,image->rows,MagickTrue,exception);
2879   else
2880     filter_image=CloneImage(image,image->columns,rows,MagickTrue,exception);
2881   if (filter_image == (Image *) NULL)
2882     return(DestroyImage(resize_image));
2883   filter_type=LanczosFilter;
2884   if (filter != UndefinedFilter)
2885     filter_type=filter;
2886   else
2887     if ((x_factor == 1.0) && (y_factor == 1.0))
2888       filter_type=PointFilter;
2889     else
2890       if ((image->storage_class == PseudoClass) ||
2891           (image->alpha_trait == BlendPixelTrait) ||
2892           ((x_factor*y_factor) > 1.0))
2893         filter_type=MitchellFilter;
2894   resize_filter=AcquireResizeFilter(image,filter_type,MagickFalse,exception);
2895   /*
2896     Resize image.
2897   */
2898   offset=0;
2899   if (x_factor > y_factor)
2900     {
2901       span=(MagickSizeType) (filter_image->columns+rows);
2902       status=HorizontalFilter(resize_filter,image,filter_image,x_factor,span,
2903         &offset,exception);
2904       status&=VerticalFilter(resize_filter,filter_image,resize_image,y_factor,
2905         span,&offset,exception);
2906     }
2907   else
2908     {
2909       span=(MagickSizeType) (filter_image->rows+columns);
2910       status=VerticalFilter(resize_filter,image,filter_image,y_factor,span,
2911         &offset,exception);
2912       status&=HorizontalFilter(resize_filter,filter_image,resize_image,x_factor,
2913         span,&offset,exception);
2914     }
2915   /*
2916     Free resources.
2917   */
2918   filter_image=DestroyImage(filter_image);
2919   resize_filter=DestroyResizeFilter(resize_filter);
2920   if (status == MagickFalse)
2921     {
2922       resize_image=DestroyImage(resize_image);
2923       return((Image *) NULL);
2924     }
2925   resize_image->type=image->type;
2926   return(resize_image);
2927 }
2928 \f
2929 /*
2930 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2931 %                                                                             %
2932 %                                                                             %
2933 %                                                                             %
2934 %   S a m p l e I m a g e                                                     %
2935 %                                                                             %
2936 %                                                                             %
2937 %                                                                             %
2938 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2939 %
2940 %  SampleImage() scales an image to the desired dimensions with pixel
2941 %  sampling.  Unlike other scaling methods, this method does not introduce
2942 %  any additional color into the scaled image.
2943 %
2944 %  The format of the SampleImage method is:
2945 %
2946 %      Image *SampleImage(const Image *image,const size_t columns,
2947 %        const size_t rows,ExceptionInfo *exception)
2948 %
2949 %  A description of each parameter follows:
2950 %
2951 %    o image: the image.
2952 %
2953 %    o columns: the number of columns in the sampled image.
2954 %
2955 %    o rows: the number of rows in the sampled image.
2956 %
2957 %    o exception: return any errors or warnings in this structure.
2958 %
2959 */
2960 MagickExport Image *SampleImage(const Image *image,const size_t columns,
2961   const size_t rows,ExceptionInfo *exception)
2962 {
2963 #define SampleImageTag  "Sample/Image"
2964
2965   CacheView
2966     *image_view,
2967     *sample_view;
2968
2969   Image
2970     *sample_image;
2971
2972   MagickBooleanType
2973     status;
2974
2975   MagickOffsetType
2976     progress;
2977
2978   register ssize_t
2979     x;
2980
2981   ssize_t
2982     *x_offset,
2983     y;
2984
2985   PointInfo
2986     sample_offset;
2987
2988   /*
2989     Initialize sampled image attributes.
2990   */
2991   assert(image != (const Image *) NULL);
2992   assert(image->signature == MagickSignature);
2993   if (image->debug != MagickFalse)
2994     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2995   assert(exception != (ExceptionInfo *) NULL);
2996   assert(exception->signature == MagickSignature);
2997   if ((columns == 0) || (rows == 0))
2998     ThrowImageException(ImageError,"NegativeOrZeroImageSize");
2999   if ((columns == image->columns) && (rows == image->rows))
3000     return(CloneImage(image,0,0,MagickTrue,exception));
3001   sample_image=CloneImage(image,columns,rows,MagickTrue,exception);
3002   if (sample_image == (Image *) NULL)
3003     return((Image *) NULL);
3004   /*
3005     Set the sampling offset, default is in the mid-point of sample regions.
3006   */
3007   sample_offset.x=sample_offset.y=0.5-MagickEpsilon;
3008   {
3009     const char
3010       *value;
3011
3012     value=GetImageArtifact(image,"sample:offset");
3013     if (value != (char *) NULL)
3014       {
3015         GeometryInfo
3016           geometry_info;
3017
3018         MagickStatusType
3019           flags;
3020
3021         (void) ParseGeometry(value,&geometry_info);
3022         flags=ParseGeometry(value,&geometry_info);
3023         sample_offset.x=sample_offset.y=geometry_info.rho/100.0-MagickEpsilon;
3024         if ((flags & SigmaValue) != 0)
3025           sample_offset.y=geometry_info.sigma/100.0-MagickEpsilon;
3026       }
3027   }
3028   /*
3029     Allocate scan line buffer and column offset buffers.
3030   */
3031   x_offset=(ssize_t *) AcquireQuantumMemory((size_t) sample_image->columns,
3032     sizeof(*x_offset));
3033   if (x_offset == (ssize_t *) NULL)
3034     {
3035       sample_image=DestroyImage(sample_image);
3036       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3037     }
3038   for (x=0; x < (ssize_t) sample_image->columns; x++)
3039     x_offset[x]=(ssize_t) ((((double) x+sample_offset.x)*image->columns)/
3040       sample_image->columns);
3041   /*
3042     Sample each row.
3043   */
3044   status=MagickTrue;
3045   progress=0;
3046   image_view=AcquireVirtualCacheView(image,exception);
3047   sample_view=AcquireAuthenticCacheView(sample_image,exception);
3048 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3049   #pragma omp parallel for schedule(static,4) shared(status) \
3050     magick_threads(image,sample_image,1,1)
3051 #endif
3052   for (y=0; y < (ssize_t) sample_image->rows; y++)
3053   {
3054     register const Quantum
3055       *restrict p;
3056
3057     register Quantum
3058       *restrict q;
3059
3060     register ssize_t
3061       x;
3062
3063     ssize_t
3064       y_offset;
3065
3066     if (status == MagickFalse)
3067       continue;
3068     y_offset=(ssize_t) ((((double) y+sample_offset.y)*image->rows)/
3069       sample_image->rows);
3070     p=GetCacheViewVirtualPixels(image_view,0,y_offset,image->columns,1,
3071       exception);
3072     q=QueueCacheViewAuthenticPixels(sample_view,0,y,sample_image->columns,1,
3073       exception);
3074     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3075       {
3076         status=MagickFalse;
3077         continue;
3078       }
3079     /*
3080       Sample each column.
3081     */
3082     for (x=0; x < (ssize_t) sample_image->columns; x++)
3083     {
3084       register ssize_t
3085         i;
3086
3087       if (GetPixelReadMask(sample_image,q) == 0)
3088         {
3089           q+=GetPixelChannels(sample_image);
3090           continue;
3091         }
3092       for (i=0; i < (ssize_t) GetPixelChannels(sample_image); i++)
3093       {
3094         PixelChannel
3095           channel;
3096
3097         PixelTrait
3098           sample_traits,
3099           traits;
3100
3101         channel=GetPixelChannelChannel(image,i);
3102         traits=GetPixelChannelTraits(image,channel);
3103         sample_traits=GetPixelChannelTraits(sample_image,channel);
3104         if ((traits == UndefinedPixelTrait) ||
3105             (sample_traits == UndefinedPixelTrait))
3106           continue;
3107         SetPixelChannel(sample_image,channel,p[x_offset[x]*GetPixelChannels(
3108           image)+i],q);
3109       }
3110       q+=GetPixelChannels(sample_image);
3111     }
3112     if (SyncCacheViewAuthenticPixels(sample_view,exception) == MagickFalse)
3113       status=MagickFalse;
3114     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3115       {
3116         MagickBooleanType
3117           proceed;
3118
3119 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3120         #pragma omp critical (MagickCore_SampleImage)
3121 #endif
3122         proceed=SetImageProgress(image,SampleImageTag,progress++,image->rows);
3123         if (proceed == MagickFalse)
3124           status=MagickFalse;
3125       }
3126   }
3127   image_view=DestroyCacheView(image_view);
3128   sample_view=DestroyCacheView(sample_view);
3129   x_offset=(ssize_t *) RelinquishMagickMemory(x_offset);
3130   sample_image->type=image->type;
3131   if (status == MagickFalse)
3132     sample_image=DestroyImage(sample_image);
3133   return(sample_image);
3134 }
3135 \f
3136 /*
3137 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3138 %                                                                             %
3139 %                                                                             %
3140 %                                                                             %
3141 %   S c a l e I m a g e                                                       %
3142 %                                                                             %
3143 %                                                                             %
3144 %                                                                             %
3145 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3146 %
3147 %  ScaleImage() changes the size of an image to the given dimensions.
3148 %
3149 %  The format of the ScaleImage method is:
3150 %
3151 %      Image *ScaleImage(const Image *image,const size_t columns,
3152 %        const size_t rows,ExceptionInfo *exception)
3153 %
3154 %  A description of each parameter follows:
3155 %
3156 %    o image: the image.
3157 %
3158 %    o columns: the number of columns in the scaled image.
3159 %
3160 %    o rows: the number of rows in the scaled image.
3161 %
3162 %    o exception: return any errors or warnings in this structure.
3163 %
3164 */
3165 MagickExport Image *ScaleImage(const Image *image,const size_t columns,
3166   const size_t rows,ExceptionInfo *exception)
3167 {
3168 #define ScaleImageTag  "Scale/Image"
3169
3170   CacheView
3171     *image_view,
3172     *scale_view;
3173
3174   double
3175     alpha,
3176     gamma,
3177     pixel[CompositePixelChannel],
3178     *scale_scanline,
3179     *scanline,
3180     *x_vector,
3181     *y_vector;
3182
3183   Image
3184     *scale_image;
3185
3186   MagickBooleanType
3187     next_column,
3188     next_row,
3189     proceed,
3190     status;
3191
3192   PixelChannel
3193     channel;
3194
3195   PixelTrait
3196     scale_traits,
3197     traits;
3198
3199   PointInfo
3200     scale,
3201     span;
3202
3203   register ssize_t
3204     i;
3205
3206   ssize_t
3207     n,
3208     number_rows,
3209     y;
3210
3211   /*
3212     Initialize scaled image attributes.
3213   */
3214   assert(image != (const Image *) NULL);
3215   assert(image->signature == MagickSignature);
3216   if (image->debug != MagickFalse)
3217     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3218   assert(exception != (ExceptionInfo *) NULL);
3219   assert(exception->signature == MagickSignature);
3220   if ((columns == 0) || (rows == 0))
3221     return((Image *) NULL);
3222   if ((columns == image->columns) && (rows == image->rows))
3223     return(CloneImage(image,0,0,MagickTrue,exception));
3224   scale_image=CloneImage(image,columns,rows,MagickTrue,exception);
3225   if (scale_image == (Image *) NULL)
3226     return((Image *) NULL);
3227   if (SetImageStorageClass(scale_image,DirectClass,exception) == MagickFalse)
3228     {
3229       scale_image=DestroyImage(scale_image);
3230       return((Image *) NULL);
3231     }
3232   /*
3233     Allocate memory.
3234   */
3235   x_vector=(double *) AcquireQuantumMemory((size_t) image->columns,
3236     GetPixelChannels(image)*sizeof(*x_vector));
3237   scanline=x_vector;
3238   if (image->rows != scale_image->rows)
3239     scanline=(double *) AcquireQuantumMemory((size_t) image->columns,
3240       GetPixelChannels(image)*sizeof(*scanline));
3241   scale_scanline=(double *) AcquireQuantumMemory((size_t)
3242     scale_image->columns,MaxPixelChannels*sizeof(*scale_scanline));
3243   y_vector=(double *) AcquireQuantumMemory((size_t) image->columns,
3244     GetPixelChannels(image)*sizeof(*y_vector));
3245   if ((scanline == (double *) NULL) ||
3246       (scale_scanline == (double *) NULL) ||
3247       (x_vector == (double *) NULL) ||
3248       (y_vector == (double *) NULL))
3249     {
3250       scale_image=DestroyImage(scale_image);
3251       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3252     }
3253   /*
3254     Scale image.
3255   */
3256   number_rows=0;
3257   next_row=MagickTrue;
3258   span.y=1.0;
3259   scale.y=(double) scale_image->rows/(double) image->rows;
3260   for (i=0; i < (ssize_t) (GetPixelChannels(image)*image->columns); i++)
3261     y_vector[i]=0.0;
3262   n=0;
3263   status=MagickTrue;
3264   image_view=AcquireVirtualCacheView(image,exception);
3265   scale_view=AcquireAuthenticCacheView(scale_image,exception);
3266   for (y=0; y < (ssize_t) scale_image->rows; y++)
3267   {
3268     register const Quantum
3269       *restrict p;
3270
3271     register Quantum
3272       *restrict q;
3273
3274     register ssize_t
3275       x;
3276
3277     if (status == MagickFalse)
3278       break;
3279     q=QueueCacheViewAuthenticPixels(scale_view,0,y,scale_image->columns,1,
3280       exception);
3281     if (q == (Quantum *) NULL)
3282       {
3283         status=MagickFalse;
3284         break;
3285       }
3286     alpha=1.0;
3287     if (scale_image->rows == image->rows)
3288       {
3289         /*
3290           Read a new scanline.
3291         */
3292         p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
3293           exception);
3294         if (p == (const Quantum *) NULL)
3295           {
3296             status=MagickFalse;
3297             break;
3298           }
3299         for (x=0; x < (ssize_t) image->columns; x++)
3300         {
3301           if (GetPixelReadMask(image,p) == 0)
3302             {
3303               p+=GetPixelChannels(image);
3304               continue;
3305             }
3306           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3307           {
3308             PixelChannel
3309               channel;
3310
3311             PixelTrait
3312               traits;
3313
3314             channel=GetPixelChannelChannel(image,i);
3315             traits=GetPixelChannelTraits(image,channel);
3316             if ((traits & BlendPixelTrait) == 0)
3317               {
3318                 x_vector[x*GetPixelChannels(image)+i]=(double) p[i];
3319                 continue;
3320               }
3321             alpha=QuantumScale*GetPixelAlpha(image,p);
3322             x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
3323           }
3324           p+=GetPixelChannels(image);
3325         }
3326       }
3327     else
3328       {
3329         /*
3330           Scale Y direction.
3331         */
3332         while (scale.y < span.y)
3333         {
3334           if ((next_row != MagickFalse) &&
3335               (number_rows < (ssize_t) image->rows))
3336             {
3337               /*
3338                 Read a new scanline.
3339               */
3340               p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
3341                 exception);
3342               if (p == (const Quantum *) NULL)
3343                 {
3344                   status=MagickFalse;
3345                   break;
3346                 }
3347               for (x=0; x < (ssize_t) image->columns; x++)
3348               {
3349                 if (GetPixelReadMask(image,p) == 0)
3350                   {
3351                     p+=GetPixelChannels(image);
3352                     continue;
3353                   }
3354                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3355                 {
3356                   PixelChannel
3357                     channel;
3358
3359                   PixelTrait
3360                     traits;
3361
3362                   channel=GetPixelChannelChannel(image,i);
3363                   traits=GetPixelChannelTraits(image,channel);
3364                   if ((traits & BlendPixelTrait) == 0)
3365                     {
3366                       x_vector[x*GetPixelChannels(image)+i]=(double)
3367                         p[i];
3368                       continue;
3369                     }
3370                   alpha=QuantumScale*GetPixelAlpha(image,p);
3371                   x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
3372                 }
3373                 p+=GetPixelChannels(image);
3374               }
3375               number_rows++;
3376             }
3377           for (x=0; x < (ssize_t) image->columns; x++)
3378             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3379               y_vector[x*GetPixelChannels(image)+i]+=scale.y*
3380                 x_vector[x*GetPixelChannels(image)+i];
3381           span.y-=scale.y;
3382           scale.y=(double) scale_image->rows/(double) image->rows;
3383           next_row=MagickTrue;
3384         }
3385         if ((next_row != MagickFalse) && (number_rows < (ssize_t) image->rows))
3386           {
3387             /*
3388               Read a new scanline.
3389             */
3390             p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
3391               exception);
3392             if (p == (const Quantum *) NULL)
3393               {
3394                 status=MagickFalse;
3395                 break;
3396               }
3397             for (x=0; x < (ssize_t) image->columns; x++)
3398             {
3399               if (GetPixelReadMask(image,p) == 0)
3400                 {
3401                   p+=GetPixelChannels(image);
3402                   continue;
3403                 }
3404               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3405               {
3406                 PixelChannel
3407                   channel;
3408
3409                 PixelTrait
3410                   traits;
3411
3412                 channel=GetPixelChannelChannel(image,i);
3413                 traits=GetPixelChannelTraits(image,channel);
3414                 if ((traits & BlendPixelTrait) == 0)
3415                   {
3416                     x_vector[x*GetPixelChannels(image)+i]=(double) p[i];
3417                     continue;
3418                   }
3419                 alpha=QuantumScale*GetPixelAlpha(image,p);
3420                 x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
3421               }
3422               p+=GetPixelChannels(image);
3423             }
3424             number_rows++;
3425             next_row=MagickFalse;
3426           }
3427         for (x=0; x < (ssize_t) image->columns; x++)
3428         {
3429           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3430           {
3431             pixel[i]=y_vector[x*GetPixelChannels(image)+i]+span.y*
3432               x_vector[x*GetPixelChannels(image)+i];
3433             scanline[x*GetPixelChannels(image)+i]=pixel[i];
3434             y_vector[x*GetPixelChannels(image)+i]=0.0;
3435           }
3436         }
3437         scale.y-=span.y;
3438         if (scale.y <= 0)
3439           {
3440             scale.y=(double) scale_image->rows/(double) image->rows;
3441             next_row=MagickTrue;
3442           }
3443         span.y=1.0;
3444       }
3445     if (scale_image->columns == image->columns)
3446       {
3447         /*
3448           Transfer scanline to scaled image.
3449         */
3450         for (x=0; x < (ssize_t) scale_image->columns; x++)
3451         {
3452           if (GetPixelReadMask(scale_image,q) == 0)
3453             {
3454               q+=GetPixelChannels(scale_image);
3455               continue;
3456             }
3457           for (i=0; i < (ssize_t) GetPixelChannels(scale_image); i++)
3458           {
3459             ssize_t
3460               offset;
3461
3462             channel=GetPixelChannelChannel(scale_image,i);
3463             traits=GetPixelChannelTraits(image,channel);
3464             scale_traits=GetPixelChannelTraits(scale_image,channel);
3465             if ((traits == UndefinedPixelTrait) ||
3466                 (scale_traits == UndefinedPixelTrait))
3467               continue;
3468             offset=GetPixelChannelOffset(image,channel);
3469             if ((traits & BlendPixelTrait) == 0)
3470               {
3471                 SetPixelChannel(scale_image,channel,ClampToQuantum(
3472                   scanline[x*GetPixelChannels(image)+offset]),q);
3473                 continue;
3474               }
3475             alpha=QuantumScale*scanline[x*GetPixelChannels(image)+
3476               GetPixelChannelChannel(image,AlphaPixelChannel)];
3477             gamma=PerceptibleReciprocal(alpha);
3478             SetPixelChannel(scale_image,channel,ClampToQuantum(gamma*scanline[
3479               x*GetPixelChannels(image)+offset]),q);
3480           }
3481           q+=GetPixelChannels(scale_image);
3482         }
3483       }
3484     else
3485       {
3486         ssize_t
3487           n;
3488
3489         /*
3490           Scale X direction.
3491         */
3492         next_column=MagickFalse;
3493         n=0;
3494         span.x=1.0;
3495         for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3496           pixel[i]=0.0;
3497         for (x=0; x < (ssize_t) image->columns; x++)
3498         {
3499           scale.x=(double) scale_image->columns/(double) image->columns;
3500           while (scale.x >= span.x)
3501           {
3502             if (next_column != MagickFalse)
3503               {
3504                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3505                   pixel[i]=0.0;
3506                 n++;
3507               }
3508             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3509             {
3510               PixelChannel
3511                 channel;
3512
3513               PixelTrait
3514                 traits;
3515
3516               channel=GetPixelChannelChannel(image,i);
3517               traits=GetPixelChannelTraits(image,channel);
3518               if (traits == UndefinedPixelTrait)
3519                 continue;
3520               pixel[i]+=span.x*scanline[x*GetPixelChannels(image)+i];
3521               scale_scanline[n*MaxPixelChannels+channel]=pixel[i];
3522             }
3523             scale.x-=span.x;
3524             span.x=1.0;
3525             next_column=MagickTrue;
3526           }
3527         if (scale.x > 0)
3528           {
3529             if (next_column != MagickFalse)
3530               {
3531                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3532                   pixel[i]=0.0;
3533                 n++;
3534                 next_column=MagickFalse;
3535               }
3536             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3537               pixel[i]+=scale.x*scanline[x*GetPixelChannels(image)+i];
3538             span.x-=scale.x;
3539           }
3540       }
3541       if (span.x > 0)
3542         {
3543           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3544             pixel[i]+=span.x*scanline[(x-1)*GetPixelChannels(image)+i];
3545         }
3546       if ((next_column == MagickFalse) &&
3547           ((ssize_t) n < (ssize_t) scale_image->columns))
3548         for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3549         {
3550           channel=GetPixelChannelChannel(image,i);
3551           scale_scanline[n*MaxPixelChannels+channel]=pixel[i];
3552         }
3553       /*
3554         Transfer scanline to scaled image.
3555       */
3556       for (x=0; x < (ssize_t) scale_image->columns; x++)
3557       {
3558         if (GetPixelReadMask(scale_image,q) == 0)
3559           {
3560             q+=GetPixelChannels(scale_image);
3561             continue;
3562           }
3563         for (i=0; i < (ssize_t) GetPixelChannels(scale_image); i++)
3564         {
3565           channel=GetPixelChannelChannel(scale_image,i);
3566           traits=GetPixelChannelTraits(image,channel);
3567           scale_traits=GetPixelChannelTraits(scale_image,channel);
3568           if ((traits == UndefinedPixelTrait) ||
3569               (scale_traits == UndefinedPixelTrait))
3570             continue;
3571           if ((traits & BlendPixelTrait) == 0)
3572             {
3573               SetPixelChannel(scale_image,channel,ClampToQuantum(
3574                 scale_scanline[x*MaxPixelChannels+channel]),q);
3575               continue;
3576             }
3577           SetPixelChannel(scale_image,channel,ClampToQuantum(scale_scanline[
3578             x*MaxPixelChannels+channel]),q);
3579         }
3580         q+=GetPixelChannels(scale_image);
3581       }
3582     }
3583     if (SyncCacheViewAuthenticPixels(scale_view,exception) == MagickFalse)
3584       {
3585         status=MagickFalse;
3586         break;
3587       }
3588     proceed=SetImageProgress(image,ScaleImageTag,(MagickOffsetType) y,
3589       image->rows);
3590     if (proceed == MagickFalse)
3591       {
3592         status=MagickFalse;
3593         break;
3594       }
3595   }
3596   scale_view=DestroyCacheView(scale_view);
3597   image_view=DestroyCacheView(image_view);
3598   /*
3599     Free allocated memory.
3600   */
3601   y_vector=(double *) RelinquishMagickMemory(y_vector);
3602   scale_scanline=(double *) RelinquishMagickMemory(scale_scanline);
3603   if (scale_image->rows != image->rows)
3604     scanline=(double *) RelinquishMagickMemory(scanline);
3605   x_vector=(double *) RelinquishMagickMemory(x_vector);
3606   scale_image->type=image->type;
3607   if (status == MagickFalse)
3608     scale_image=DestroyImage(scale_image);
3609   return(scale_image);
3610 }
3611 \f
3612 /*
3613 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3614 %                                                                             %
3615 %                                                                             %
3616 %                                                                             %
3617 %   T h u m b n a i l I m a g e                                               %
3618 %                                                                             %
3619 %                                                                             %
3620 %                                                                             %
3621 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3622 %
3623 %  ThumbnailImage() changes the size of an image to the given dimensions and
3624 %  removes any associated profiles.  The goal is to produce small low cost
3625 %  thumbnail images suited for display on the Web.
3626 %
3627 %  The format of the ThumbnailImage method is:
3628 %
3629 %      Image *ThumbnailImage(const Image *image,const size_t columns,
3630 %        const size_t rows,ExceptionInfo *exception)
3631 %
3632 %  A description of each parameter follows:
3633 %
3634 %    o image: the image.
3635 %
3636 %    o columns: the number of columns in the scaled image.
3637 %
3638 %    o rows: the number of rows in the scaled image.
3639 %
3640 %    o exception: return any errors or warnings in this structure.
3641 %
3642 */
3643 MagickExport Image *ThumbnailImage(const Image *image,const size_t columns,
3644   const size_t rows,ExceptionInfo *exception)
3645 {
3646 #define SampleFactor  5
3647
3648   char
3649     value[MaxTextExtent];
3650
3651   const char
3652     *name;
3653
3654   Image
3655     *thumbnail_image;
3656
3657   double
3658     x_factor,
3659     y_factor;
3660
3661   size_t
3662     version;
3663
3664   struct stat
3665     attributes;
3666
3667   assert(image != (Image *) NULL);
3668   assert(image->signature == MagickSignature);
3669   if (image->debug != MagickFalse)
3670     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3671   assert(exception != (ExceptionInfo *) NULL);
3672   assert(exception->signature == MagickSignature);
3673   x_factor=(double) columns/(double) image->columns;
3674   y_factor=(double) rows/(double) image->rows;
3675   if ((x_factor*y_factor) > 0.1)
3676     thumbnail_image=ResizeImage(image,columns,rows,image->filter,exception);
3677   else
3678     if (((SampleFactor*columns) < 128) || ((SampleFactor*rows) < 128))
3679       thumbnail_image=ResizeImage(image,columns,rows,image->filter,exception);
3680     else
3681       {
3682         Image
3683           *sample_image;
3684
3685         sample_image=SampleImage(image,SampleFactor*columns,SampleFactor*rows,
3686           exception);
3687         if (sample_image == (Image *) NULL)
3688           return((Image *) NULL);
3689         thumbnail_image=ResizeImage(sample_image,columns,rows,image->filter,
3690           exception);
3691         sample_image=DestroyImage(sample_image);
3692       }
3693   if (thumbnail_image == (Image *) NULL)
3694     return(thumbnail_image);
3695   (void) ParseAbsoluteGeometry("0x0+0+0",&thumbnail_image->page);
3696   if (thumbnail_image->alpha_trait != BlendPixelTrait)
3697     (void) SetImageAlphaChannel(thumbnail_image,OpaqueAlphaChannel,exception);
3698   thumbnail_image->depth=8;
3699   thumbnail_image->interlace=NoInterlace;
3700   /*
3701     Strip all profiles except color profiles.
3702   */
3703   ResetImageProfileIterator(thumbnail_image);
3704   for (name=GetNextImageProfile(thumbnail_image); name != (const char *) NULL; )
3705   {
3706     if ((LocaleCompare(name,"icc") != 0) && (LocaleCompare(name,"icm") != 0))
3707      {
3708        (void) DeleteImageProfile(thumbnail_image,name);
3709        ResetImageProfileIterator(thumbnail_image);
3710      }
3711     name=GetNextImageProfile(thumbnail_image);
3712   }
3713   (void) DeleteImageProperty(thumbnail_image,"comment");
3714   (void) CopyMagickString(value,image->magick_filename,MaxTextExtent);
3715   if (strstr(image->magick_filename,"//") == (char *) NULL)
3716     (void) FormatLocaleString(value,MaxTextExtent,"file://%s",
3717       image->magick_filename);
3718   (void) SetImageProperty(thumbnail_image,"Thumb::URI",value,exception);
3719   (void) CopyMagickString(value,image->magick_filename,MaxTextExtent);
3720   if ( IfMagickTrue(GetPathAttributes(image->filename,&attributes)) )
3721     {
3722       (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3723         attributes.st_mtime);
3724       (void) SetImageProperty(thumbnail_image,"Thumb::MTime",value,exception);
3725     }
3726   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3727     attributes.st_mtime);
3728   (void) FormatMagickSize(GetBlobSize(image),MagickFalse,value);
3729   (void) ConcatenateMagickString(value,"B",MaxTextExtent);
3730   (void) SetImageProperty(thumbnail_image,"Thumb::Size",value,exception);
3731   (void) FormatLocaleString(value,MaxTextExtent,"image/%s",image->magick);
3732   LocaleLower(value);
3733   (void) SetImageProperty(thumbnail_image,"Thumb::Mimetype",value,exception);
3734   (void) SetImageProperty(thumbnail_image,"software",GetMagickVersion(&version),
3735     exception);
3736   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3737     image->magick_columns);
3738   (void) SetImageProperty(thumbnail_image,"Thumb::Image::Width",value,
3739     exception);
3740   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3741     image->magick_rows);
3742   (void) SetImageProperty(thumbnail_image,"Thumb::Image::Height",value,
3743     exception);
3744   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3745     GetImageListLength(image));
3746   (void) SetImageProperty(thumbnail_image,"Thumb::Document::Pages",value,
3747     exception);
3748   return(thumbnail_image);
3749 }