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