granicus.if.org Git - imagemagick/blob - MagickCore/resize.c

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