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