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