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