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 channel=GetPixelChannelChannel(image,i);
1663         PixelTrait traits=GetPixelChannelTraits(image,channel);
1664         PixelTrait resize_traits=GetPixelChannelTraits(resize_image,channel);
1665         if ((traits == UndefinedPixelTrait) ||
1666             (resize_traits == UndefinedPixelTrait))
1667           continue;
1668         offset.x=((double) x+0.5)*scale.x-0.5;
1669         status=InterpolatePixelChannels(image,image_view,resize_image,method,
1670           offset.x,offset.y,q,exception);
1671       }
1672       q+=GetPixelChannels(resize_image);
1673     }
1674     if( IfMagickFalse(SyncCacheViewAuthenticPixels(resize_view,exception)) )
1675       continue;
1676     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1677       {
1678         MagickBooleanType
1679           proceed;
1680
1681 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1682         #pragma omp critical (MagickCore_InterpolativeResizeImage)
1683 #endif
1684         proceed=SetImageProgress(image,InterpolativeResizeImageTag,progress++,
1685           image->rows);
1686         if( IfMagickFalse(proceed) )
1687           status=MagickFalse;
1688       }
1689   }
1690   resize_view=DestroyCacheView(resize_view);
1691   image_view=DestroyCacheView(image_view);
1692   if (status == MagickFalse)
1693     resize_image=DestroyImage(resize_image);
1694   return(resize_image);
1695 }
1696 #if defined(MAGICKCORE_LQR_DELEGATE)
1697 \f
1698 /*
1699 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1700 %                                                                             %
1701 %                                                                             %
1702 %                                                                             %
1703 %   L i q u i d R e s c a l e I m a g e                                       %
1704 %                                                                             %
1705 %                                                                             %
1706 %                                                                             %
1707 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1708 %
1709 %  LiquidRescaleImage() rescales image with seam carving.
1710 %
1711 %  The format of the LiquidRescaleImage method is:
1712 %
1713 %      Image *LiquidRescaleImage(const Image *image,const size_t columns,
1714 %        const size_t rows,const double delta_x,const double rigidity,
1715 %        ExceptionInfo *exception)
1716 %
1717 %  A description of each parameter follows:
1718 %
1719 %    o image: the image.
1720 %
1721 %    o columns: the number of columns in the rescaled image.
1722 %
1723 %    o rows: the number of rows in the rescaled image.
1724 %
1725 %    o delta_x: maximum seam transversal step (0 means straight seams).
1726 %
1727 %    o rigidity: introduce a bias for non-straight seams (typically 0).
1728 %
1729 %    o exception: return any errors or warnings in this structure.
1730 %
1731 */
1732 MagickExport Image *LiquidRescaleImage(const Image *image,const size_t columns,
1733   const size_t rows,const double delta_x,const double rigidity,
1734   ExceptionInfo *exception)
1735 {
1736 #define LiquidRescaleImageTag  "Rescale/Image"
1737
1738   CacheView
1739     *image_view,
1740     *rescale_view;
1741
1742   gfloat
1743     *packet,
1744     *pixels;
1745
1746   Image
1747     *rescale_image;
1748
1749   int
1750     x_offset,
1751     y_offset;
1752
1753   LqrCarver
1754     *carver;
1755
1756   LqrRetVal
1757     lqr_status;
1758
1759   MagickBooleanType
1760     status;
1761
1762   register gfloat
1763     *q;
1764
1765   ssize_t
1766     y;
1767
1768   /*
1769     Liquid rescale image.
1770   */
1771   assert(image != (const Image *) NULL);
1772   assert(image->signature == MagickSignature);
1773   if( IfMagickTrue(image->debug) )
1774     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1775   assert(exception != (ExceptionInfo *) NULL);
1776   assert(exception->signature == MagickSignature);
1777   if ((columns == 0) || (rows == 0))
1778     return((Image *) NULL);
1779   if ((columns == image->columns) && (rows == image->rows))
1780     return(CloneImage(image,0,0,MagickTrue,exception));
1781   if ((columns <= 2) || (rows <= 2))
1782     return(ResizeImage(image,columns,rows,image->filter,exception));
1783   if ((columns >= (2*image->columns)) || (rows >= (2*image->rows)))
1784     {
1785       Image
1786         *resize_image;
1787
1788       size_t
1789         height,
1790         width;
1791
1792       /*
1793         Honor liquid resize size limitations.
1794       */
1795       for (width=image->columns; columns >= (2*width-1); width*=2);
1796       for (height=image->rows; rows >= (2*height-1); height*=2);
1797       resize_image=ResizeImage(image,width,height,image->filter,exception);
1798       if (resize_image == (Image *) NULL)
1799         return((Image *) NULL);
1800       rescale_image=LiquidRescaleImage(resize_image,columns,rows,delta_x,
1801         rigidity,exception);
1802       resize_image=DestroyImage(resize_image);
1803       return(rescale_image);
1804     }
1805   pixels=(gfloat *) AcquireQuantumMemory(image->columns,image->rows*
1806     GetPixelChannels(image)*sizeof(*pixels));
1807   if (pixels == (gfloat *) NULL)
1808     return((Image *) NULL);
1809   status=MagickTrue;
1810   q=pixels;
1811   image_view=AcquireVirtualCacheView(image,exception);
1812   for (y=0; y < (ssize_t) image->rows; y++)
1813   {
1814     register const Quantum
1815       *restrict p;
1816
1817     register ssize_t
1818       x;
1819
1820     if (status == MagickFalse)
1821       continue;
1822     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1823     if (p == (const Quantum *) NULL)
1824       {
1825         status=MagickFalse;
1826         continue;
1827       }
1828     for (x=0; x < (ssize_t) image->columns; x++)
1829     {
1830       register ssize_t
1831         i;
1832
1833       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1834         *q++=QuantumScale*p[i];
1835       p+=GetPixelChannels(image);
1836     }
1837   }
1838   image_view=DestroyCacheView(image_view);
1839   carver=lqr_carver_new_ext(pixels,image->columns,image->rows,
1840     GetPixelChannels(image),LQR_COLDEPTH_32F);
1841   if (carver == (LqrCarver *) NULL)
1842     {
1843       pixels=(gfloat *) RelinquishMagickMemory(pixels);
1844       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1845     }
1846   lqr_status=lqr_carver_init(carver,(int) delta_x,rigidity);
1847   lqr_status=lqr_carver_resize(carver,columns,rows);
1848   (void) lqr_status;
1849   rescale_image=CloneImage(image,lqr_carver_get_width(carver),
1850     lqr_carver_get_height(carver),MagickTrue,exception);
1851   if (rescale_image == (Image *) NULL)
1852     {
1853       pixels=(gfloat *) RelinquishMagickMemory(pixels);
1854       return((Image *) NULL);
1855     }
1856   if( IfMagickFalse(SetImageStorageClass(rescale_image,DirectClass,exception)) )
1857     {
1858       pixels=(gfloat *) RelinquishMagickMemory(pixels);
1859       rescale_image=DestroyImage(rescale_image);
1860       return((Image *) NULL);
1861     }
1862   rescale_view=AcquireAuthenticCacheView(rescale_image,exception);
1863   (void) lqr_carver_scan_reset(carver);
1864   while (lqr_carver_scan_ext(carver,&x_offset,&y_offset,(void **) &packet) != 0)
1865   {
1866     register Quantum
1867       *restrict q;
1868
1869     register ssize_t
1870       i;
1871
1872     q=QueueCacheViewAuthenticPixels(rescale_view,x_offset,y_offset,1,1,
1873       exception);
1874     if (q == (Quantum *) NULL)
1875       break;
1876     for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1877     {
1878       PixelChannel channel=GetPixelChannelChannel(image,i);
1879       PixelTrait traits=GetPixelChannelTraits(image,channel);
1880       PixelTrait rescale_traits=GetPixelChannelTraits(rescale_image,channel);
1881       if ((traits == UndefinedPixelTrait) ||
1882           (rescale_traits == UndefinedPixelTrait))
1883         continue;
1884       SetPixelChannel(rescale_image,channel,ClampToQuantum(QuantumRange*
1885         packet[i]),q);
1886     }
1887     if( IfMagickFalse(SyncCacheViewAuthenticPixels(rescale_view,exception)) )
1888       break;
1889   }
1890   rescale_view=DestroyCacheView(rescale_view);
1891   lqr_carver_destroy(carver);
1892   return(rescale_image);
1893 }
1894 #else
1895 MagickExport Image *LiquidRescaleImage(const Image *image,
1896   const size_t magick_unused(columns),const size_t magick_unused(rows),
1897   const double magick_unused(delta_x),const double magick_unused(rigidity),
1898   ExceptionInfo *exception)
1899 {
1900   assert(image != (const Image *) NULL);
1901   assert(image->signature == MagickSignature);
1902   if( IfMagickTrue(image->debug) )
1903     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1904   assert(exception != (ExceptionInfo *) NULL);
1905   assert(exception->signature == MagickSignature);
1906   (void) ThrowMagickException(exception,GetMagickModule(),MissingDelegateError,
1907     "DelegateLibrarySupportNotBuiltIn","'%s' (LQR)",image->filename);
1908   return((Image *) NULL);
1909 }
1910 #endif
1911 \f
1912 /*
1913 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1914 %                                                                             %
1915 %                                                                             %
1916 %                                                                             %
1917 %   M a g n i f y I m a g e                                                   %
1918 %                                                                             %
1919 %                                                                             %
1920 %                                                                             %
1921 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1922 %
1923 %  MagnifyImage() is a convenience method that scales an image proportionally
1924 %  to twice its size.
1925 %
1926 %  The format of the MagnifyImage method is:
1927 %
1928 %      Image *MagnifyImage(const Image *image,ExceptionInfo *exception)
1929 %
1930 %  A description of each parameter follows:
1931 %
1932 %    o image: the image.
1933 %
1934 %    o exception: return any errors or warnings in this structure.
1935 %
1936 */
1937 MagickExport Image *MagnifyImage(const Image *image,ExceptionInfo *exception)
1938 {
1939   Image
1940     *magnify_image;
1941
1942   assert(image != (Image *) NULL);
1943   assert(image->signature == MagickSignature);
1944   if( IfMagickTrue(image->debug) )
1945     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1946   assert(exception != (ExceptionInfo *) NULL);
1947   assert(exception->signature == MagickSignature);
1948   magnify_image=ResizeImage(image,2*image->columns,2*image->rows,SplineFilter,
1949     exception);
1950   return(magnify_image);
1951 }
1952 \f
1953 /*
1954 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1955 %                                                                             %
1956 %                                                                             %
1957 %                                                                             %
1958 %   M i n i f y I m a g e                                                     %
1959 %                                                                             %
1960 %                                                                             %
1961 %                                                                             %
1962 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1963 %
1964 %  MinifyImage() is a convenience method that scales an image proportionally to
1965 %  half its size.
1966 %
1967 %  The format of the MinifyImage method is:
1968 %
1969 %      Image *MinifyImage(const Image *image,ExceptionInfo *exception)
1970 %
1971 %  A description of each parameter follows:
1972 %
1973 %    o image: the image.
1974 %
1975 %    o exception: return any errors or warnings in this structure.
1976 %
1977 */
1978 MagickExport Image *MinifyImage(const Image *image,ExceptionInfo *exception)
1979 {
1980   Image
1981     *minify_image;
1982
1983   assert(image != (Image *) NULL);
1984   assert(image->signature == MagickSignature);
1985   if( IfMagickTrue(image->debug) )
1986     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1987   assert(exception != (ExceptionInfo *) NULL);
1988   assert(exception->signature == MagickSignature);
1989   minify_image=ResizeImage(image,image->columns/2,image->rows/2,SplineFilter,
1990     exception);
1991   return(minify_image);
1992 }
1993 \f
1994 /*
1995 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1996 %                                                                             %
1997 %                                                                             %
1998 %                                                                             %
1999 %   R e s a m p l e I m a g e                                                 %
2000 %                                                                             %
2001 %                                                                             %
2002 %                                                                             %
2003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2004 %
2005 %  ResampleImage() resize image in terms of its pixel size, so that when
2006 %  displayed at the given resolution it will be the same size in terms of
2007 %  real world units as the original image at the original resolution.
2008 %
2009 %  The format of the ResampleImage method is:
2010 %
2011 %      Image *ResampleImage(Image *image,const double x_resolution,
2012 %        const double y_resolution,const FilterTypes filter,
2013 %        ExceptionInfo *exception)
2014 %
2015 %  A description of each parameter follows:
2016 %
2017 %    o image: the image to be resized to fit the given resolution.
2018 %
2019 %    o x_resolution: the new image x resolution.
2020 %
2021 %    o y_resolution: the new image y resolution.
2022 %
2023 %    o filter: Image filter to use.
2024 %
2025 %    o exception: return any errors or warnings in this structure.
2026 %
2027 */
2028 MagickExport Image *ResampleImage(const Image *image,const double x_resolution,
2029   const double y_resolution,const FilterTypes filter,ExceptionInfo *exception)
2030 {
2031 #define ResampleImageTag  "Resample/Image"
2032
2033   Image
2034     *resample_image;
2035
2036   size_t
2037     height,
2038     width;
2039
2040   /*
2041     Initialize sampled image attributes.
2042   */
2043   assert(image != (const Image *) NULL);
2044   assert(image->signature == MagickSignature);
2045   if( IfMagickTrue(image->debug) )
2046     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2047   assert(exception != (ExceptionInfo *) NULL);
2048   assert(exception->signature == MagickSignature);
2049   width=(size_t) (x_resolution*image->columns/(image->resolution.x == 0.0 ?
2050     72.0 : image->resolution.x)+0.5);
2051   height=(size_t) (y_resolution*image->rows/(image->resolution.y == 0.0 ?
2052     72.0 : image->resolution.y)+0.5);
2053   resample_image=ResizeImage(image,width,height,filter,exception);
2054   if (resample_image != (Image *) NULL)
2055     {
2056       resample_image->resolution.x=x_resolution;
2057       resample_image->resolution.y=y_resolution;
2058     }
2059   return(resample_image);
2060 }
2061 \f
2062 /*
2063 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2064 %                                                                             %
2065 %                                                                             %
2066 %                                                                             %
2067 %   R e s i z e I m a g e                                                     %
2068 %                                                                             %
2069 %                                                                             %
2070 %                                                                             %
2071 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2072 %
2073 %  ResizeImage() scales an image to the desired dimensions, using the given
2074 %  filter (see AcquireFilterInfo()).
2075 %
2076 %  If an undefined filter is given the filter defaults to Mitchell for a
2077 %  colormapped image, a image with a matte channel, or if the image is
2078 %  enlarged.  Otherwise the filter defaults to a Lanczos.
2079 %
2080 %  ResizeImage() was inspired by Paul Heckbert's "zoom" program.
2081 %
2082 %  The format of the ResizeImage method is:
2083 %
2084 %      Image *ResizeImage(Image *image,const size_t columns,
2085 %        const size_t rows,const FilterTypes filter,const double blur,
2086 %        ExceptionInfo *exception)
2087 %
2088 %  A description of each parameter follows:
2089 %
2090 %    o image: the image.
2091 %
2092 %    o columns: the number of columns in the scaled image.
2093 %
2094 %    o rows: the number of rows in the scaled image.
2095 %
2096 %    o filter: Image filter to use.
2097 %
2098 %    o blur: the blur factor where > 1 is blurry, < 1 is sharp.  Typically set
2099 %      this to 1.0.
2100 %
2101 %    o exception: return any errors or warnings in this structure.
2102 %
2103 */
2104
2105 typedef struct _ContributionInfo
2106 {
2107   double
2108     weight;
2109
2110   ssize_t
2111     pixel;
2112 } ContributionInfo;
2113
2114 static ContributionInfo **DestroyContributionThreadSet(
2115   ContributionInfo **contribution)
2116 {
2117   register ssize_t
2118     i;
2119
2120   assert(contribution != (ContributionInfo **) NULL);
2121   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2122     if (contribution[i] != (ContributionInfo *) NULL)
2123       contribution[i]=(ContributionInfo *) RelinquishAlignedMemory(
2124         contribution[i]);
2125   contribution=(ContributionInfo **) RelinquishMagickMemory(contribution);
2126   return(contribution);
2127 }
2128
2129 static ContributionInfo **AcquireContributionThreadSet(const size_t count)
2130 {
2131   register ssize_t
2132     i;
2133
2134   ContributionInfo
2135     **contribution;
2136
2137   size_t
2138     number_threads;
2139
2140   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2141   contribution=(ContributionInfo **) AcquireQuantumMemory(number_threads,
2142     sizeof(*contribution));
2143   if (contribution == (ContributionInfo **) NULL)
2144     return((ContributionInfo **) NULL);
2145   (void) ResetMagickMemory(contribution,0,number_threads*sizeof(*contribution));
2146   for (i=0; i < (ssize_t) number_threads; i++)
2147   {
2148     contribution[i]=(ContributionInfo *) MagickAssumeAligned(
2149       AcquireAlignedMemory(count,sizeof(**contribution)));
2150     if (contribution[i] == (ContributionInfo *) NULL)
2151       return(DestroyContributionThreadSet(contribution));
2152   }
2153   return(contribution);
2154 }
2155
2156 static inline double MagickMax(const double x,const double y)
2157 {
2158   if (x > y)
2159     return(x);
2160   return(y);
2161 }
2162
2163 static inline double MagickMin(const double x,const double y)
2164 {
2165   if (x < y)
2166     return(x);
2167   return(y);
2168 }
2169
2170 static MagickBooleanType HorizontalFilter(const ResizeFilter *resize_filter,
2171   const Image *image,Image *resize_image,const double x_factor,
2172   const MagickSizeType span,MagickOffsetType *offset,ExceptionInfo *exception)
2173 {
2174 #define ResizeImageTag  "Resize/Image"
2175
2176   CacheView
2177     *image_view,
2178     *resize_view;
2179
2180   ClassType
2181     storage_class;
2182
2183   ContributionInfo
2184     **restrict contributions;
2185
2186   double
2187     scale,
2188     support;
2189
2190   MagickBooleanType
2191     status;
2192
2193   ssize_t
2194     x;
2195
2196   /*
2197     Apply filter to resize horizontally from image to resize image.
2198   */
2199   scale=MagickMax(1.0/x_factor+MagickEpsilon,1.0);
2200   support=scale*GetResizeFilterSupport(resize_filter);
2201   storage_class=support > 0.5 ? DirectClass : image->storage_class;
2202   if (IfMagickFalse(SetImageStorageClass(resize_image,storage_class,exception)))
2203     return(MagickFalse);
2204   if (support < 0.5)
2205     {
2206       /*
2207         Support too small even for nearest neighbour: Reduce to point sampling.
2208       */
2209       support=(double) 0.5;
2210       scale=1.0;
2211     }
2212   contributions=AcquireContributionThreadSet((size_t) (2.0*support+3.0));
2213   if (contributions == (ContributionInfo **) NULL)
2214     {
2215       (void) ThrowMagickException(exception,GetMagickModule(),
2216         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2217       return(MagickFalse);
2218     }
2219   status=MagickTrue;
2220   scale=PerceptibleReciprocal(scale);
2221   image_view=AcquireVirtualCacheView(image,exception);
2222   resize_view=AcquireAuthenticCacheView(resize_image,exception);
2223 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2224   #pragma omp parallel for schedule(static,4) shared(status) \
2225     magick_threads(image,resize_image,resize_image->columns,1)
2226 #endif
2227   for (x=0; x < (ssize_t) resize_image->columns; x++)
2228   {
2229     double
2230       bisect,
2231       density;
2232
2233     register const Quantum
2234       *restrict p;
2235
2236     register ContributionInfo
2237       *restrict contribution;
2238
2239     register Quantum
2240       *restrict q;
2241
2242     register ssize_t
2243       y;
2244
2245     ssize_t
2246       n,
2247       start,
2248       stop;
2249
2250     if (status == MagickFalse)
2251       continue;
2252     bisect=(double) (x+0.5)/x_factor+MagickEpsilon;
2253     start=(ssize_t) MagickMax(bisect-support+0.5,0.0);
2254     stop=(ssize_t) MagickMin(bisect+support+0.5,(double) image->columns);
2255     density=0.0;
2256     contribution=contributions[GetOpenMPThreadId()];
2257     for (n=0; n < (stop-start); n++)
2258     {
2259       contribution[n].pixel=start+n;
2260       contribution[n].weight=GetResizeFilterWeight(resize_filter,scale*
2261         ((double) (start+n)-bisect+0.5));
2262       density+=contribution[n].weight;
2263     }
2264     if ((density != 0.0) && (density != 1.0))
2265       {
2266         register ssize_t
2267           i;
2268
2269         /*
2270           Normalize.
2271         */
2272         density=PerceptibleReciprocal(density);
2273         for (i=0; i < n; i++)
2274           contribution[i].weight*=density;
2275       }
2276     p=GetCacheViewVirtualPixels(image_view,contribution[0].pixel,0,(size_t)
2277       (contribution[n-1].pixel-contribution[0].pixel+1),image->rows,exception);
2278     q=QueueCacheViewAuthenticPixels(resize_view,x,0,1,resize_image->rows,
2279       exception);
2280     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2281       {
2282         status=MagickFalse;
2283         continue;
2284       }
2285     for (y=0; y < (ssize_t) resize_image->rows; y++)
2286     {
2287       register ssize_t
2288         i;
2289
2290       for (i=0; i < (ssize_t) GetPixelChannels(resize_image); i++)
2291       {
2292         double
2293           alpha,
2294           gamma,
2295           pixel;
2296
2297         register ssize_t
2298           j;
2299
2300         ssize_t
2301           k;
2302
2303         PixelChannel channel=GetPixelChannelChannel(image,i);
2304         PixelTrait traits=GetPixelChannelTraits(image,channel);
2305         PixelTrait resize_traits=GetPixelChannelTraits(resize_image,channel);
2306         if ((traits == UndefinedPixelTrait) ||
2307             (resize_traits == UndefinedPixelTrait))
2308           continue;
2309         if (((resize_traits & CopyPixelTrait) != 0) ||
2310             (GetPixelMask(resize_image,q) != 0))
2311           {
2312             j=(ssize_t) (MagickMin(MagickMax(bisect,(double) start),(double)
2313               stop-1.0)+0.5);
2314             k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2315               (contribution[j-start].pixel-contribution[0].pixel);
2316             SetPixelChannel(resize_image,channel,p[k*GetPixelChannels(image)+i],
2317               q);
2318             continue;
2319           }
2320         pixel=0.0;
2321         if ((GetPixelAlphaTraits(resize_image) & BlendPixelTrait) == 0)
2322           {
2323             /*
2324               No alpha blending.
2325             */
2326             for (j=0; j < n; j++)
2327             {
2328               k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2329                 (contribution[j].pixel-contribution[0].pixel);
2330               alpha=contribution[j].weight;
2331               pixel+=alpha*p[k*GetPixelChannels(image)+i];
2332             }
2333             SetPixelChannel(resize_image,channel,ClampToQuantum(pixel),q);
2334             continue;
2335           }
2336         /*
2337           Alpha blending.
2338         */
2339         gamma=0.0;
2340         for (j=0; j < n; j++)
2341         {
2342           k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2343             (contribution[j].pixel-contribution[0].pixel);
2344           alpha=contribution[j].weight*QuantumScale*
2345             GetPixelAlpha(image,p+k*GetPixelChannels(image));
2346           pixel+=alpha*p[k*GetPixelChannels(image)+i];
2347           gamma+=alpha;
2348         }
2349         gamma=PerceptibleReciprocal(gamma);
2350         SetPixelChannel(resize_image,channel,ClampToQuantum(gamma*pixel),q);
2351       }
2352       q+=GetPixelChannels(resize_image);
2353     }
2354     if( IfMagickFalse(SyncCacheViewAuthenticPixels(resize_view,exception)) )
2355       status=MagickFalse;
2356     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2357       {
2358         MagickBooleanType
2359           proceed;
2360
2361 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2362         #pragma omp critical (MagickCore_HorizontalFilter)
2363 #endif
2364         proceed=SetImageProgress(image,ResizeImageTag,(*offset)++,span);
2365         if( IfMagickFalse(proceed) )
2366           status=MagickFalse;
2367       }
2368   }
2369   resize_view=DestroyCacheView(resize_view);
2370   image_view=DestroyCacheView(image_view);
2371   contributions=DestroyContributionThreadSet(contributions);
2372   return(status);
2373 }
2374
2375 static MagickBooleanType VerticalFilter(const ResizeFilter *resize_filter,
2376   const Image *image,Image *resize_image,const double y_factor,
2377   const MagickSizeType span,MagickOffsetType *offset,ExceptionInfo *exception)
2378 {
2379   CacheView
2380     *image_view,
2381     *resize_view;
2382
2383   ClassType
2384     storage_class;
2385
2386   ContributionInfo
2387     **restrict contributions;
2388
2389   double
2390     scale,
2391     support;
2392
2393   MagickBooleanType
2394     status;
2395
2396   PixelInfo
2397     zero;
2398
2399   ssize_t
2400     y;
2401
2402   /*
2403     Apply filter to resize vertically from image to resize image.
2404   */
2405   scale=MagickMax(1.0/y_factor+MagickEpsilon,1.0);
2406   support=scale*GetResizeFilterSupport(resize_filter);
2407   storage_class=support > 0.5 ? DirectClass : image->storage_class;
2408   if (IfMagickFalse(SetImageStorageClass(resize_image,storage_class,exception)))
2409     return(MagickFalse);
2410   if (support < 0.5)
2411     {
2412       /*
2413         Support too small even for nearest neighbour: Reduce to point sampling.
2414       */
2415       support=(double) 0.5;
2416       scale=1.0;
2417     }
2418   contributions=AcquireContributionThreadSet((size_t) (2.0*support+3.0));
2419   if (contributions == (ContributionInfo **) NULL)
2420     {
2421       (void) ThrowMagickException(exception,GetMagickModule(),
2422         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2423       return(MagickFalse);
2424     }
2425   status=MagickTrue;
2426   scale=PerceptibleReciprocal(scale);
2427   (void) ResetMagickMemory(&zero,0,sizeof(zero));
2428   image_view=AcquireVirtualCacheView(image,exception);
2429   resize_view=AcquireAuthenticCacheView(resize_image,exception);
2430 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2431   #pragma omp parallel for schedule(static,4) shared(status) \
2432     magick_threads(image,resize_image,resize_image->rows,1)
2433 #endif
2434   for (y=0; y < (ssize_t) resize_image->rows; y++)
2435   {
2436     double
2437       bisect,
2438       density;
2439
2440     register const Quantum
2441       *restrict p;
2442
2443     register ContributionInfo
2444       *restrict contribution;
2445
2446     register Quantum
2447       *restrict q;
2448
2449     register ssize_t
2450       x;
2451
2452     ssize_t
2453       n,
2454       start,
2455       stop;
2456
2457     if (status == MagickFalse)
2458       continue;
2459     bisect=(double) (y+0.5)/y_factor+MagickEpsilon;
2460     start=(ssize_t) MagickMax(bisect-support+0.5,0.0);
2461     stop=(ssize_t) MagickMin(bisect+support+0.5,(double) image->rows);
2462     density=0.0;
2463     contribution=contributions[GetOpenMPThreadId()];
2464     for (n=0; n < (stop-start); n++)
2465     {
2466       contribution[n].pixel=start+n;
2467       contribution[n].weight=GetResizeFilterWeight(resize_filter,scale*
2468         ((double) (start+n)-bisect+0.5));
2469       density+=contribution[n].weight;
2470     }
2471     if ((density != 0.0) && (density != 1.0))
2472       {
2473         register ssize_t
2474           i;
2475
2476         /*
2477           Normalize.
2478         */
2479         density=PerceptibleReciprocal(density);
2480         for (i=0; i < n; i++)
2481           contribution[i].weight*=density;
2482       }
2483     p=GetCacheViewVirtualPixels(image_view,0,contribution[0].pixel,
2484       image->columns,(size_t) (contribution[n-1].pixel-contribution[0].pixel+1),
2485       exception);
2486     q=QueueCacheViewAuthenticPixels(resize_view,0,y,resize_image->columns,1,
2487       exception);
2488     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2489       {
2490         status=MagickFalse;
2491         continue;
2492       }
2493     for (x=0; x < (ssize_t) resize_image->columns; x++)
2494     {
2495       register ssize_t
2496         i;
2497
2498       for (i=0; i < (ssize_t) GetPixelChannels(resize_image); i++)
2499       {
2500         double
2501           alpha,
2502           gamma,
2503           pixel;
2504
2505         register ssize_t
2506           j;
2507
2508         ssize_t
2509           k;
2510
2511         PixelChannel channel=GetPixelChannelChannel(image,i);
2512         PixelTrait traits=GetPixelChannelTraits(image,channel);
2513         PixelTrait resize_traits=GetPixelChannelTraits(resize_image,channel);
2514         if ((traits == UndefinedPixelTrait) ||
2515             (resize_traits == UndefinedPixelTrait))
2516           continue;
2517         if (((resize_traits & CopyPixelTrait) != 0) ||
2518             (GetPixelMask(resize_image,q) != 0))
2519           {
2520             j=(ssize_t) (MagickMin(MagickMax(bisect,(double) start),(double)
2521               stop-1.0)+0.5);
2522             k=(ssize_t) ((contribution[j-start].pixel-contribution[0].pixel)*
2523               image->columns+x);
2524             SetPixelChannel(resize_image,channel,p[k*GetPixelChannels(image)+i],
2525               q);
2526             continue;
2527           }
2528         pixel=0.0;
2529         if ((GetPixelAlphaTraits(resize_image) & BlendPixelTrait) == 0)
2530           {
2531             /*
2532               No alpha blending.
2533             */
2534             for (j=0; j < n; j++)
2535             {
2536               k=(ssize_t) ((contribution[j].pixel-contribution[0].pixel)*
2537                 image->columns+x);
2538               alpha=contribution[j].weight;
2539               pixel+=alpha*p[k*GetPixelChannels(image)+i];
2540             }
2541             SetPixelChannel(resize_image,channel,ClampToQuantum(pixel),q);
2542             continue;
2543           }
2544         gamma=0.0;
2545         for (j=0; j < n; j++)
2546         {
2547           k=(ssize_t) ((contribution[j].pixel-contribution[0].pixel)*
2548             image->columns+x);
2549           alpha=contribution[j].weight*QuantumScale*GetPixelAlpha(image,p+k*
2550             GetPixelChannels(image));
2551           pixel+=alpha*p[k*GetPixelChannels(image)+i];
2552           gamma+=alpha;
2553         }
2554         gamma=PerceptibleReciprocal(gamma);
2555         SetPixelChannel(resize_image,channel,ClampToQuantum(gamma*pixel),q);
2556       }
2557       q+=GetPixelChannels(resize_image);
2558     }
2559     if( IfMagickFalse(SyncCacheViewAuthenticPixels(resize_view,exception)) )
2560       status=MagickFalse;
2561     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2562       {
2563         MagickBooleanType
2564           proceed;
2565
2566 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2567         #pragma omp critical (MagickCore_VerticalFilter)
2568 #endif
2569         proceed=SetImageProgress(image,ResizeImageTag,(*offset)++,span);
2570         if( IfMagickFalse(proceed) )
2571           status=MagickFalse;
2572       }
2573   }
2574   resize_view=DestroyCacheView(resize_view);
2575   image_view=DestroyCacheView(image_view);
2576   contributions=DestroyContributionThreadSet(contributions);
2577   return(status);
2578 }
2579
2580 MagickExport Image *ResizeImage(const Image *image,const size_t columns,
2581   const size_t rows,const FilterTypes filter,ExceptionInfo *exception)
2582 {
2583   FilterTypes
2584     filter_type;
2585
2586   Image
2587     *filter_image,
2588     *resize_image;
2589
2590   MagickOffsetType
2591     offset;
2592
2593   double
2594     x_factor,
2595     y_factor;
2596
2597   MagickSizeType
2598     span;
2599
2600   MagickStatusType
2601     status;
2602
2603   ResizeFilter
2604     *resize_filter;
2605
2606   /*
2607     Acquire resize image.
2608   */
2609   assert(image != (Image *) NULL);
2610   assert(image->signature == MagickSignature);
2611   if( IfMagickTrue(image->debug) )
2612     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2613   assert(exception != (ExceptionInfo *) NULL);
2614   assert(exception->signature == MagickSignature);
2615   if ((columns == 0) || (rows == 0))
2616     ThrowImageException(ImageError,"NegativeOrZeroImageSize");
2617   if ((columns == image->columns) && (rows == image->rows) &&
2618       (filter == UndefinedFilter))
2619     return(CloneImage(image,0,0,MagickTrue,exception));
2620   resize_image=CloneImage(image,columns,rows,MagickTrue,exception);
2621   if (resize_image == (Image *) NULL)
2622     return(resize_image);
2623   /*
2624     Acquire resize filter.
2625   */
2626   x_factor=(double) columns/(double) image->columns;
2627   y_factor=(double) rows/(double) image->rows;
2628   if (x_factor > y_factor)
2629     filter_image=CloneImage(image,columns,image->rows,MagickTrue,exception);
2630   else
2631     filter_image=CloneImage(image,image->columns,rows,MagickTrue,exception);
2632   if (filter_image == (Image *) NULL)
2633     return(DestroyImage(resize_image));
2634   filter_type=LanczosFilter;
2635   if (filter != UndefinedFilter)
2636     filter_type=filter;
2637   else
2638     if ((x_factor == 1.0) && (y_factor == 1.0))
2639       filter_type=PointFilter;
2640     else
2641       if ((image->storage_class == PseudoClass) ||
2642           (image->alpha_trait == BlendPixelTrait) ||
2643           ((x_factor*y_factor) > 1.0))
2644         filter_type=MitchellFilter;
2645   resize_filter=AcquireResizeFilter(image,filter_type,MagickFalse,exception);
2646   /*
2647     Resize image.
2648   */
2649   offset=0;
2650   if (x_factor > y_factor)
2651     {
2652       span=(MagickSizeType) (filter_image->columns+rows);
2653       status=HorizontalFilter(resize_filter,image,filter_image,x_factor,span,
2654         &offset,exception);
2655       status&=VerticalFilter(resize_filter,filter_image,resize_image,y_factor,
2656         span,&offset,exception);
2657     }
2658   else
2659     {
2660       span=(MagickSizeType) (filter_image->rows+columns);
2661       status=VerticalFilter(resize_filter,image,filter_image,y_factor,span,
2662         &offset,exception);
2663       status&=HorizontalFilter(resize_filter,filter_image,resize_image,x_factor,
2664         span,&offset,exception);
2665     }
2666   /*
2667     Free resources.
2668   */
2669   filter_image=DestroyImage(filter_image);
2670   resize_filter=DestroyResizeFilter(resize_filter);
2671   if (status == MagickFalse)
2672     {
2673       resize_image=DestroyImage(resize_image);
2674       return((Image *) NULL);
2675     }
2676   resize_image->type=image->type;
2677   return(resize_image);
2678 }
2679 \f
2680 /*
2681 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2682 %                                                                             %
2683 %                                                                             %
2684 %                                                                             %
2685 %   S a m p l e I m a g e                                                     %
2686 %                                                                             %
2687 %                                                                             %
2688 %                                                                             %
2689 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2690 %
2691 %  SampleImage() scales an image to the desired dimensions with pixel
2692 %  sampling.  Unlike other scaling methods, this method does not introduce
2693 %  any additional color into the scaled image.
2694 %
2695 %  The format of the SampleImage method is:
2696 %
2697 %      Image *SampleImage(const Image *image,const size_t columns,
2698 %        const size_t rows,ExceptionInfo *exception)
2699 %
2700 %  A description of each parameter follows:
2701 %
2702 %    o image: the image.
2703 %
2704 %    o columns: the number of columns in the sampled image.
2705 %
2706 %    o rows: the number of rows in the sampled image.
2707 %
2708 %    o exception: return any errors or warnings in this structure.
2709 %
2710 */
2711 MagickExport Image *SampleImage(const Image *image,const size_t columns,
2712   const size_t rows,ExceptionInfo *exception)
2713 {
2714 #define SampleImageTag  "Sample/Image"
2715
2716   CacheView
2717     *image_view,
2718     *sample_view;
2719
2720   Image
2721     *sample_image;
2722
2723   MagickBooleanType
2724     status;
2725
2726   MagickOffsetType
2727     progress;
2728
2729   register ssize_t
2730     x;
2731
2732   ssize_t
2733     *x_offset,
2734     y;
2735
2736   /*
2737     Initialize sampled image attributes.
2738   */
2739   assert(image != (const Image *) NULL);
2740   assert(image->signature == MagickSignature);
2741   if( IfMagickTrue(image->debug) )
2742     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2743   assert(exception != (ExceptionInfo *) NULL);
2744   assert(exception->signature == MagickSignature);
2745   if ((columns == 0) || (rows == 0))
2746     ThrowImageException(ImageError,"NegativeOrZeroImageSize");
2747   if ((columns == image->columns) && (rows == image->rows))
2748     return(CloneImage(image,0,0,MagickTrue,exception));
2749   sample_image=CloneImage(image,columns,rows,MagickTrue,exception);
2750   if (sample_image == (Image *) NULL)
2751     return((Image *) NULL);
2752   /*
2753     Allocate scan line buffer and column offset buffers.
2754   */
2755   x_offset=(ssize_t *) AcquireQuantumMemory((size_t) sample_image->columns,
2756     sizeof(*x_offset));
2757   if (x_offset == (ssize_t *) NULL)
2758     {
2759       sample_image=DestroyImage(sample_image);
2760       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2761     }
2762   for (x=0; x < (ssize_t) sample_image->columns; x++)
2763     x_offset[x]=(ssize_t) (((double) x*image->columns)/sample_image->columns+
2764       0.5);
2765   /*
2766     Sample each row.
2767   */
2768   status=MagickTrue;
2769   progress=0;
2770   image_view=AcquireVirtualCacheView(image,exception);
2771   sample_view=AcquireAuthenticCacheView(sample_image,exception);
2772 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2773   #pragma omp parallel for schedule(static,4) shared(status) \
2774     magick_threads(image,sample_image,1,1)
2775 #endif
2776   for (y=0; y < (ssize_t) sample_image->rows; y++)
2777   {
2778     register const Quantum
2779       *restrict p;
2780
2781     register Quantum
2782       *restrict q;
2783
2784     register ssize_t
2785       x;
2786
2787     ssize_t
2788       y_offset;
2789
2790     if (status == MagickFalse)
2791       continue;
2792     y_offset=(ssize_t) (((double) y*image->rows)/sample_image->rows+0.5);
2793     p=GetCacheViewVirtualPixels(image_view,0,y_offset,image->columns,1,
2794       exception);
2795     q=QueueCacheViewAuthenticPixels(sample_view,0,y,sample_image->columns,1,
2796       exception);
2797     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2798       {
2799         status=MagickFalse;
2800         continue;
2801       }
2802     /*
2803       Sample each column.
2804     */
2805     for (x=0; x < (ssize_t) sample_image->columns; x++)
2806     {
2807       register ssize_t
2808         i;
2809
2810       if (GetPixelMask(sample_image,q) != 0)
2811         {
2812           q+=GetPixelChannels(sample_image);
2813           continue;
2814         }
2815       for (i=0; i < (ssize_t) GetPixelChannels(sample_image); i++)
2816       {
2817         PixelChannel channel=GetPixelChannelChannel(image,i);
2818         PixelTrait traits=GetPixelChannelTraits(image,channel);
2819         PixelTrait sample_traits=GetPixelChannelTraits(sample_image,channel);
2820         if ((traits == UndefinedPixelTrait) ||
2821             (sample_traits == UndefinedPixelTrait))
2822           continue;
2823         SetPixelChannel(sample_image,channel,p[x_offset[x]*GetPixelChannels(
2824           image)+i],q);
2825       }
2826       q+=GetPixelChannels(sample_image);
2827     }
2828     if( IfMagickFalse(SyncCacheViewAuthenticPixels(sample_view,exception)) )
2829       status=MagickFalse;
2830     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2831       {
2832         MagickBooleanType
2833           proceed;
2834
2835 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2836         #pragma omp critical (MagickCore_SampleImage)
2837 #endif
2838         proceed=SetImageProgress(image,SampleImageTag,progress++,image->rows);
2839         if( IfMagickFalse(proceed) )
2840           status=MagickFalse;
2841       }
2842   }
2843   image_view=DestroyCacheView(image_view);
2844   sample_view=DestroyCacheView(sample_view);
2845   x_offset=(ssize_t *) RelinquishMagickMemory(x_offset);
2846   sample_image->type=image->type;
2847   if (status == MagickFalse)
2848     sample_image=DestroyImage(sample_image);
2849   return(sample_image);
2850 }
2851 \f
2852 /*
2853 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2854 %                                                                             %
2855 %                                                                             %
2856 %                                                                             %
2857 %   S c a l e I m a g e                                                       %
2858 %                                                                             %
2859 %                                                                             %
2860 %                                                                             %
2861 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2862 %
2863 %  ScaleImage() changes the size of an image to the given dimensions.
2864 %
2865 %  The format of the ScaleImage method is:
2866 %
2867 %      Image *ScaleImage(const Image *image,const size_t columns,
2868 %        const size_t rows,ExceptionInfo *exception)
2869 %
2870 %  A description of each parameter follows:
2871 %
2872 %    o image: the image.
2873 %
2874 %    o columns: the number of columns in the scaled image.
2875 %
2876 %    o rows: the number of rows in the scaled image.
2877 %
2878 %    o exception: return any errors or warnings in this structure.
2879 %
2880 */
2881 MagickExport Image *ScaleImage(const Image *image,const size_t columns,
2882   const size_t rows,ExceptionInfo *exception)
2883 {
2884 #define ScaleImageTag  "Scale/Image"
2885
2886   CacheView
2887     *image_view,
2888     *scale_view;
2889
2890   double
2891     alpha,
2892     gamma,
2893     pixel[CompositePixelChannel],
2894     *scale_scanline,
2895     *scanline,
2896     *x_vector,
2897     *y_vector;
2898
2899   Image
2900     *scale_image;
2901
2902   MagickBooleanType
2903     next_column,
2904     next_row,
2905     proceed,
2906     status;
2907
2908   PointInfo
2909     scale,
2910     span;
2911
2912   register ssize_t
2913     i;
2914
2915   ssize_t
2916     n,
2917     number_rows,
2918     y;
2919
2920   /*
2921     Initialize scaled image attributes.
2922   */
2923   assert(image != (const Image *) NULL);
2924   assert(image->signature == MagickSignature);
2925   if( IfMagickTrue(image->debug) )
2926     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2927   assert(exception != (ExceptionInfo *) NULL);
2928   assert(exception->signature == MagickSignature);
2929   if ((columns == 0) || (rows == 0))
2930     return((Image *) NULL);
2931   if ((columns == image->columns) && (rows == image->rows))
2932     return(CloneImage(image,0,0,MagickTrue,exception));
2933   scale_image=CloneImage(image,columns,rows,MagickTrue,exception);
2934   if (scale_image == (Image *) NULL)
2935     return((Image *) NULL);
2936   if( IfMagickFalse(SetImageStorageClass(scale_image,DirectClass,exception)) )
2937     {
2938       scale_image=DestroyImage(scale_image);
2939       return((Image *) NULL);
2940     }
2941   /*
2942     Allocate memory.
2943   */
2944   x_vector=(double *) AcquireQuantumMemory((size_t) image->columns,
2945     GetPixelChannels(image)*sizeof(*x_vector));
2946   scanline=x_vector;
2947   if (image->rows != scale_image->rows)
2948     scanline=(double *) AcquireQuantumMemory((size_t) image->columns,
2949       GetPixelChannels(image)*sizeof(*scanline));
2950   scale_scanline=(double *) AcquireQuantumMemory((size_t)
2951     scale_image->columns,MaxPixelChannels*sizeof(*scale_scanline));
2952   y_vector=(double *) AcquireQuantumMemory((size_t) image->columns,
2953     GetPixelChannels(image)*sizeof(*y_vector));
2954   if ((scanline == (double *) NULL) ||
2955       (scale_scanline == (double *) NULL) ||
2956       (x_vector == (double *) NULL) ||
2957       (y_vector == (double *) NULL))
2958     {
2959       scale_image=DestroyImage(scale_image);
2960       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2961     }
2962   /*
2963     Scale image.
2964   */
2965   number_rows=0;
2966   next_row=MagickTrue;
2967   span.y=1.0;
2968   scale.y=(double) scale_image->rows/(double) image->rows;
2969   for (i=0; i < (ssize_t) (GetPixelChannels(image)*image->columns); i++)
2970     y_vector[i]=0.0;
2971   n=0;
2972   status=MagickTrue;
2973   image_view=AcquireVirtualCacheView(image,exception);
2974   scale_view=AcquireAuthenticCacheView(scale_image,exception);
2975   for (y=0; y < (ssize_t) scale_image->rows; y++)
2976   {
2977     register const Quantum
2978       *restrict p;
2979
2980     register Quantum
2981       *restrict q;
2982
2983     register ssize_t
2984       x;
2985
2986     if (status == MagickFalse)
2987       break;
2988     q=QueueCacheViewAuthenticPixels(scale_view,0,y,scale_image->columns,1,
2989       exception);
2990     if (q == (Quantum *) NULL)
2991       {
2992         status=MagickFalse;
2993         break;
2994       }
2995     alpha=1.0;
2996     if (scale_image->rows == image->rows)
2997       {
2998         /*
2999           Read a new scanline.
3000         */
3001         p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
3002           exception);
3003         if (p == (const Quantum *) NULL)
3004           {
3005             status=MagickFalse;
3006             break;
3007           }
3008         for (x=0; x < (ssize_t) image->columns; x++)
3009         {
3010           if (GetPixelMask(image,p) != 0)
3011             {
3012               p+=GetPixelChannels(image);
3013               continue;
3014             }
3015           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3016           {
3017             PixelChannel channel=GetPixelChannelChannel(image,i);
3018             PixelTrait traits=GetPixelChannelTraits(image,channel);
3019             if ((traits & BlendPixelTrait) == 0)
3020               {
3021                 x_vector[x*GetPixelChannels(image)+i]=(double) p[i];
3022                 continue;
3023               }
3024             alpha=QuantumScale*GetPixelAlpha(image,p);
3025             x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
3026           }
3027           p+=GetPixelChannels(image);
3028         }
3029       }
3030     else
3031       {
3032         /*
3033           Scale Y direction.
3034         */
3035         while (scale.y < span.y)
3036         {
3037           if( IfMagickTrue(next_row) && (number_rows < (ssize_t) image->rows))
3038             {
3039               /*
3040                 Read a new scanline.
3041               */
3042               p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
3043                 exception);
3044               if (p == (const Quantum *) NULL)
3045                 {
3046                   status=MagickFalse;
3047                   break;
3048                 }
3049               for (x=0; x < (ssize_t) image->columns; x++)
3050               {
3051                 if (GetPixelMask(image,p) != 0)
3052                   {
3053                     p+=GetPixelChannels(image);
3054                     continue;
3055                   }
3056                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3057                 {
3058                   PixelChannel channel=GetPixelChannelChannel(image,i);
3059                   PixelTrait traits=GetPixelChannelTraits(image,channel);
3060                   if ((traits & BlendPixelTrait) == 0)
3061                     {
3062                       x_vector[x*GetPixelChannels(image)+i]=(double) p[i];
3063                       continue;
3064                     }
3065                   alpha=QuantumScale*GetPixelAlpha(image,p);
3066                   x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
3067                 }
3068                 p+=GetPixelChannels(image);
3069               }
3070               number_rows++;
3071             }
3072           for (x=0; x < (ssize_t) image->columns; x++)
3073             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3074               y_vector[x*GetPixelChannels(image)+i]+=scale.y*
3075                 x_vector[x*GetPixelChannels(image)+i];
3076           span.y-=scale.y;
3077           scale.y=(double) scale_image->rows/(double) image->rows;
3078           next_row=MagickTrue;
3079         }
3080         if( IfMagickTrue(next_row) && (number_rows < (ssize_t) image->rows))
3081           {
3082             /*
3083               Read a new scanline.
3084             */
3085             p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
3086               exception);
3087             if (p == (const Quantum *) NULL)
3088               {
3089                 status=MagickFalse;
3090                 break;
3091               }
3092             for (x=0; x < (ssize_t) image->columns; x++)
3093             {
3094               if (GetPixelMask(image,p) != 0)
3095                 {
3096                   p+=GetPixelChannels(image);
3097                   continue;
3098                 }
3099               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3100               {
3101                 PixelChannel channel=GetPixelChannelChannel(image,i);
3102                 PixelTrait traits=GetPixelChannelTraits(image,channel);
3103                 if ((traits & BlendPixelTrait) == 0)
3104                   {
3105                     x_vector[x*GetPixelChannels(image)+i]=(double) p[i];
3106                     continue;
3107                   }
3108                 alpha=QuantumScale*GetPixelAlpha(image,p);
3109                 x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
3110               }
3111               p+=GetPixelChannels(image);
3112             }
3113             number_rows++;
3114             next_row=MagickFalse;
3115           }
3116         for (x=0; x < (ssize_t) image->columns; x++)
3117         {
3118           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3119           {
3120             pixel[i]=y_vector[x*GetPixelChannels(image)+i]+span.y*
3121               x_vector[x*GetPixelChannels(image)+i];
3122             scanline[x*GetPixelChannels(image)+i]=pixel[i];
3123             y_vector[x*GetPixelChannels(image)+i]=0.0;
3124           }
3125         }
3126         scale.y-=span.y;
3127         if (scale.y <= 0)
3128           {
3129             scale.y=(double) scale_image->rows/(double) image->rows;
3130             next_row=MagickTrue;
3131           }
3132         span.y=1.0;
3133       }
3134     if (scale_image->columns == image->columns)
3135       {
3136         /*
3137           Transfer scanline to scaled image.
3138         */
3139         for (x=0; x < (ssize_t) scale_image->columns; x++)
3140         {
3141           if (GetPixelMask(scale_image,q) != 0)
3142             {
3143               q+=GetPixelChannels(scale_image);
3144               continue;
3145             }
3146           for (i=0; i < (ssize_t) GetPixelChannels(scale_image); i++)
3147           {
3148             ssize_t
3149               offset;
3150
3151             PixelChannel channel=GetPixelChannelChannel(scale_image,i);
3152             PixelTrait traits=GetPixelChannelTraits(image,channel);
3153             PixelTrait scale_traits=GetPixelChannelTraits(scale_image,channel);
3154             if ((traits == UndefinedPixelTrait) ||
3155                 (scale_traits == UndefinedPixelTrait))
3156               continue;
3157             offset=GetPixelChannelOffset(image,channel);
3158             if ((traits & BlendPixelTrait) == 0)
3159               {
3160                 SetPixelChannel(scale_image,channel,ClampToQuantum(
3161                   scanline[x*GetPixelChannels(image)+offset]),q);
3162                 continue;
3163               }
3164             alpha=QuantumScale*scanline[x*GetPixelChannels(image)+
3165               GetPixelChannelChannel(image,AlphaPixelChannel)];
3166             gamma=PerceptibleReciprocal(alpha);
3167             SetPixelChannel(scale_image,channel,ClampToQuantum(gamma*scanline[
3168               x*GetPixelChannels(image)+offset]),q);
3169           }
3170           q+=GetPixelChannels(scale_image);
3171         }
3172       }
3173     else
3174       {
3175         ssize_t
3176           n;
3177
3178         /*
3179           Scale X direction.
3180         */
3181         next_column=MagickFalse;
3182         n=0;
3183         span.x=1.0;
3184         for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3185           pixel[i]=0.0;
3186         for (x=0; x < (ssize_t) image->columns; x++)
3187         {
3188           scale.x=(double) scale_image->columns/(double) image->columns;
3189           while (scale.x >= span.x)
3190           {
3191             if( IfMagickTrue(next_column) )
3192               {
3193                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3194                   pixel[i]=0.0;
3195                 n++;
3196               }
3197             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3198             {
3199               PixelChannel channel=GetPixelChannelChannel(image,i);
3200               PixelTrait traits=GetPixelChannelTraits(image,channel);
3201               if (traits == UndefinedPixelTrait)
3202                 continue;
3203               pixel[i]+=span.x*scanline[x*GetPixelChannels(image)+i];
3204               scale_scanline[n*MaxPixelChannels+channel]=pixel[i];
3205             }
3206             scale.x-=span.x;
3207             span.x=1.0;
3208             next_column=MagickTrue;
3209           }
3210         if (scale.x > 0)
3211           {
3212             if( IfMagickTrue(next_column) )
3213               {
3214                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3215                   pixel[i]=0.0;
3216                 n++;
3217                 next_column=MagickFalse;
3218               }
3219             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3220               pixel[i]+=scale.x*scanline[x*GetPixelChannels(image)+i];
3221             span.x-=scale.x;
3222           }
3223       }
3224       if (span.x > 0)
3225         {
3226           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3227             pixel[i]+=span.x*scanline[(x-1)*GetPixelChannels(image)+i];
3228         }
3229       if( IfMagickFalse(next_column) &&
3230           ((ssize_t) n < (ssize_t) scale_image->columns))
3231         for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3232         {
3233           PixelChannel channel=GetPixelChannelChannel(image,i);
3234           scale_scanline[n*MaxPixelChannels+channel]=pixel[i];
3235         }
3236       /*
3237         Transfer scanline to scaled image.
3238       */
3239       for (x=0; x < (ssize_t) scale_image->columns; x++)
3240       {
3241         if (GetPixelMask(scale_image,q) != 0)
3242           {
3243             q+=GetPixelChannels(scale_image);
3244             continue;
3245           }
3246         for (i=0; i < (ssize_t) GetPixelChannels(scale_image); i++)
3247         {
3248           PixelChannel channel=GetPixelChannelChannel(scale_image,i);
3249           PixelTrait traits=GetPixelChannelTraits(image,channel);
3250           PixelTrait scale_traits=GetPixelChannelTraits(scale_image,channel);
3251           if ((traits == UndefinedPixelTrait) ||
3252               (scale_traits == UndefinedPixelTrait))
3253             continue;
3254           if ((traits & BlendPixelTrait) == 0)
3255             {
3256               SetPixelChannel(scale_image,channel,ClampToQuantum(
3257                 scale_scanline[x*MaxPixelChannels+channel]),q);
3258               continue;
3259             }
3260           SetPixelChannel(scale_image,channel,ClampToQuantum(
3261             scale_scanline[x*MaxPixelChannels+channel]),q);
3262         }
3263         q+=GetPixelChannels(scale_image);
3264       }
3265     }
3266     if( IfMagickFalse(SyncCacheViewAuthenticPixels(scale_view,exception)) )
3267       {
3268         status=MagickFalse;
3269         break;
3270       }
3271     proceed=SetImageProgress(image,ScaleImageTag,(MagickOffsetType) y,
3272       image->rows);
3273     if( IfMagickFalse(proceed) )
3274       {
3275         status=MagickFalse;
3276         break;
3277       }
3278   }
3279   scale_view=DestroyCacheView(scale_view);
3280   image_view=DestroyCacheView(image_view);
3281   /*
3282     Free allocated memory.
3283   */
3284   y_vector=(double *) RelinquishMagickMemory(y_vector);
3285   scale_scanline=(double *) RelinquishMagickMemory(scale_scanline);
3286   if (scale_image->rows != image->rows)
3287     scanline=(double *) RelinquishMagickMemory(scanline);
3288   x_vector=(double *) RelinquishMagickMemory(x_vector);
3289   scale_image->type=image->type;
3290   if (status == MagickFalse)
3291     scale_image=DestroyImage(scale_image);
3292   return(scale_image);
3293 }
3294 \f
3295 /*
3296 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3297 %                                                                             %
3298 %                                                                             %
3299 %                                                                             %
3300 %   T h u m b n a i l I m a g e                                               %
3301 %                                                                             %
3302 %                                                                             %
3303 %                                                                             %
3304 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3305 %
3306 %  ThumbnailImage() changes the size of an image to the given dimensions and
3307 %  removes any associated profiles.  The goal is to produce small low cost
3308 %  thumbnail images suited for display on the Web.
3309 %
3310 %  The format of the ThumbnailImage method is:
3311 %
3312 %      Image *ThumbnailImage(const Image *image,const size_t columns,
3313 %        const size_t rows,ExceptionInfo *exception)
3314 %
3315 %  A description of each parameter follows:
3316 %
3317 %    o image: the image.
3318 %
3319 %    o columns: the number of columns in the scaled image.
3320 %
3321 %    o rows: the number of rows in the scaled image.
3322 %
3323 %    o exception: return any errors or warnings in this structure.
3324 %
3325 */
3326 MagickExport Image *ThumbnailImage(const Image *image,const size_t columns,
3327   const size_t rows,ExceptionInfo *exception)
3328 {
3329 #define SampleFactor  5
3330
3331   char
3332     value[MaxTextExtent];
3333
3334   const char
3335     *name;
3336
3337   Image
3338     *thumbnail_image;
3339
3340   double
3341     x_factor,
3342     y_factor;
3343
3344   size_t
3345     version;
3346
3347   struct stat
3348     attributes;
3349
3350   assert(image != (Image *) NULL);
3351   assert(image->signature == MagickSignature);
3352   if( IfMagickTrue(image->debug) )
3353     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3354   assert(exception != (ExceptionInfo *) NULL);
3355   assert(exception->signature == MagickSignature);
3356   x_factor=(double) columns/(double) image->columns;
3357   y_factor=(double) rows/(double) image->rows;
3358   if ((x_factor*y_factor) > 0.1)
3359     thumbnail_image=ResizeImage(image,columns,rows,image->filter,exception);
3360   else
3361     if (((SampleFactor*columns) < 128) || ((SampleFactor*rows) < 128))
3362       thumbnail_image=ResizeImage(image,columns,rows,image->filter,exception);
3363     else
3364       {
3365         Image
3366           *sample_image;
3367
3368         sample_image=SampleImage(image,SampleFactor*columns,SampleFactor*rows,
3369           exception);
3370         if (sample_image == (Image *) NULL)
3371           return((Image *) NULL);
3372         thumbnail_image=ResizeImage(sample_image,columns,rows,image->filter,
3373           exception);
3374         sample_image=DestroyImage(sample_image);
3375       }
3376   if (thumbnail_image == (Image *) NULL)
3377     return(thumbnail_image);
3378   (void) ParseAbsoluteGeometry("0x0+0+0",&thumbnail_image->page);
3379   if(thumbnail_image->alpha_trait != BlendPixelTrait)
3380     (void) SetImageAlphaChannel(thumbnail_image,OpaqueAlphaChannel,exception);
3381   thumbnail_image->depth=8;
3382   thumbnail_image->interlace=NoInterlace;
3383   /*
3384     Strip all profiles except color profiles.
3385   */
3386   ResetImageProfileIterator(thumbnail_image);
3387   for (name=GetNextImageProfile(thumbnail_image); name != (const char *) NULL; )
3388   {
3389     if ((LocaleCompare(name,"icc") != 0) && (LocaleCompare(name,"icm") != 0))
3390      {
3391        (void) DeleteImageProfile(thumbnail_image,name);
3392        ResetImageProfileIterator(thumbnail_image);
3393      }
3394     name=GetNextImageProfile(thumbnail_image);
3395   }
3396   (void) DeleteImageProperty(thumbnail_image,"comment");
3397   (void) CopyMagickString(value,image->magick_filename,MaxTextExtent);
3398   if (strstr(image->magick_filename,"//") == (char *) NULL)
3399     (void) FormatLocaleString(value,MaxTextExtent,"file://%s",
3400       image->magick_filename);
3401   (void) SetImageProperty(thumbnail_image,"Thumb::URI",value,exception);
3402   (void) CopyMagickString(value,image->magick_filename,MaxTextExtent);
3403   if( IfMagickTrue(GetPathAttributes(image->filename,&attributes)) )
3404     {
3405       (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3406         attributes.st_mtime);
3407       (void) SetImageProperty(thumbnail_image,"Thumb::MTime",value,exception);
3408     }
3409   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3410     attributes.st_mtime);
3411   (void) FormatMagickSize(GetBlobSize(image),MagickFalse,value);
3412   (void) ConcatenateMagickString(value,"B",MaxTextExtent);
3413   (void) SetImageProperty(thumbnail_image,"Thumb::Size",value,exception);
3414   (void) FormatLocaleString(value,MaxTextExtent,"image/%s",image->magick);
3415   LocaleLower(value);
3416   (void) SetImageProperty(thumbnail_image,"Thumb::Mimetype",value,exception);
3417   (void) SetImageProperty(thumbnail_image,"software",GetMagickVersion(&version),
3418     exception);
3419   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3420     image->magick_columns);
3421   (void) SetImageProperty(thumbnail_image,"Thumb::Image::Width",value,
3422     exception);
3423   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3424     image->magick_rows);
3425   (void) SetImageProperty(thumbnail_image,"Thumb::Image::Height",value,
3426     exception);
3427   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3428     GetImageListLength(image));
3429   (void) SetImageProperty(thumbnail_image,"Thumb::Document::Pages",value,
3430     exception);
3431   return(thumbnail_image);
3432 }