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

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