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 function used as a windowing function
 722                 Typically equal to the location of the first zero crossing. */
 723       B,C;   /* BC-spline coefficients, ignored if not a CubicBC filter. */
 724   } const filters[SentinelFilter] =
 725   {
 726     { Box,       0.5, 0.5,     0.0, 0.0 }, /* Undefined (default to Box)  */
 727     { Box,       0.0, 0.5,     0.0, 0.0 }, /* Point (special handling)    */
 728     { Box,       0.5, 0.5,     0.0, 0.0 }, /* Box                         */
 729     { Triangle,  1.0, 1.0,     0.0, 0.0 }, /* Triangle                    */
 730     { CubicBC,   1.0, 1.0,     0.0, 0.0 }, /* Hermite (cubic  B=C=0)      */
 731     { Hanning,   1.0, 1.0,     0.0, 0.0 }, /* Hanning, cosine window      */
 732     { Hamming,   1.0, 1.0,     0.0, 0.0 }, /* Hamming, '' variation       */
 733     { Blackman,  1.0, 1.0,     0.0, 0.0 }, /* Blackman, 2*cosine window   */
 734     { Gaussian,  2.0, 1.5,     0.0, 0.0 }, /* Gaussian                    */
 735     { Quadratic, 1.5, 1.5,     0.0, 0.0 }, /* Quadratic gaussian          */
 736     { CubicBC,   2.0, 2.0,     1.0, 0.0 }, /* Cubic B-Spline (B=1,C=0)    */
 737     { CubicBC,   2.0, 1.0,     0.0, 0.5 }, /* Catmull-Rom    (B=0,C=1/2)  */
 738     { CubicBC,   2.0, 1.0, 1./3., 1./3. }, /* Mitchell       (B=C=1/3)    */
 739     { SincFast,  3.0, 1.0,     0.0, 0.0 }, /* Lanczos, 3-lobed Sinc-Sinc  */
 740     { Jinc,      3.0, 1.2196698912665045, 0.0, 0.0 }, /* Raw 3-lobed Jinc */
 741     { Sinc,      4.0, 1.0,     0.0, 0.0 }, /* Raw 4-lobed Sinc            */
 742     { Kaiser,    1.0, 1.0,     0.0, 0.0 }, /* Kaiser (square root window) */
 743     { Welsh,     1.0, 1.0,     0.0, 0.0 }, /* Welsh (parabolic window)    */
 744     { CubicBC,   2.0, 2.0,     1.0, 0.0 }, /* Parzen (B-Spline window)    */
 745     { Lagrange,  2.0, 1.0,     0.0, 0.0 }, /* Lagrange sinc approximation */
 746     { Bohman,    1.0, 1.0,     0.0, 0.0 }, /* Bohman, 2*Cosine window     */
 747     { Triangle,  1.0, 1.0,     0.0, 0.0 }, /* Bartlett (triangle window)  */
 748     { SincFast,  4.0, 1.0,     0.0, 0.0 }, /* Raw fast sinc ("Pade"-type) */
 749     { Jinc,      2.0, 1.2196698912665045, 0.0, 0.0 },
 750                                                  /* Lanczos2D (Jinc-Jinc) */
 751     { Jinc,      2.0, 1.1684849904329952, 0.0, 0.0 },
 752                              /* Lanczos2D sharpened with blur=0.958033808 */
 753     { CubicBC,   2.0, 1.1685777620836932,
 754                           0.37821575509399867, 0.31089212245300067 }
 755                      /* Robidoux: Keys cubic close to Lanczos2D sharpened */
 756   };
 757   /*
 758     The known zero crossings of the Jinc() or more accurately the Jinc(x*PI)
 759     function being used as a filter. It is used by the "filter:lobes" expert
 760     setting and for 'lobes' for Jinc functions in the previous table. This
 761     way users do not have to deal with the highly irrational lobe sizes of the
 762     Jinc filter.
 763
 764     Values taken from
 765       http://cose.math.bas.bg/webMathematica/webComputing/BesselZeros.jsp
 766     using Jv-function with v=1, then dividing by PI.
 767   */
 768   static MagickRealType
 769     jinc_zeros[16] =
 770     {
 771       1.2196698912665045,
 772       2.2331305943815286,
 773       3.2383154841662362,
 774       4.2410628637960699,
 775       5.2427643768701817,
 776       6.2439216898644877,
 777       7.244759868719957,
 778       8.2453949139520427,
 779       9.2458926849494673,
 780       10.246293348754916,
 781       11.246622794877883,
 782       12.246898461138105,
 783       13.247132522181061,
 784       14.247333735806849,
 785       15.2475085630373,
 786       16.247661874700962
 787    };
 788
 789   /*
 790     Allocate resize filter.
 791   */
 792   assert(image != (const Image *) NULL);
 793   assert(image->signature == MagickSignature);
 794   if (image->debug != MagickFalse)
 795     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 796   assert(UndefinedFilter < filter && filter < SentinelFilter);
 797   assert(exception != (ExceptionInfo *) NULL);
 798   assert(exception->signature == MagickSignature);
 799   resize_filter=(ResizeFilter *) AcquireMagickMemory(sizeof(*resize_filter));
 800   if (resize_filter == (ResizeFilter *) NULL)
 801     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 802   /*
 803     Defaults for the requested filter.
 804   */
 805   filter_type=mapping[filter].filter;
 806   window_type=mapping[filter].window;
 807   resize_filter->blur = blur;
 808   sigma = 0.5;
 809   /* Cylindrical Filters should use Jinc instead of Sinc */
 810   if (cylindrical != MagickFalse)
 811     switch (filter_type)
 812     {
 813       case SincFilter:
 814         /* Promote 1D Sinc Filter to a 2D Jinc filter. */
 815         if ( filter != SincFilter )
 816           filter_type=JincFilter;
 817         break;
 818       case SincFastFilter:
 819         /* Ditto for SincFast variant */
 820         if ( filter != SincFastFilter )
 821           filter_type=JincFilter;
 822         break;
 823       case LanczosFilter:
 824         /* Promote Lanczos from a Sinc-Sinc to a Jinc-Jinc. */
 825         filter_type=JincFilter;
 826         window_type=JincFilter;
 827         break;
 828       case Lanczos2DSharpFilter:
 829         /* Sharpened by Nicholas Robidoux so as to optimize for
 830          * minimal blurring of orthogonal lines
 831          */
 832         resize_filter->blur *= 0.958033808;
 833         break;
 834       case GaussianFilter:
 835         sigma = (MagickRealType) (MagickSQ2/2.0);  /* Cylindrical Gaussian sigma is sqrt(2)/2 */
 836         break;
 837       default:
 838         break;
 839     }
 840   else
 841     switch (filter_type)
 842     {
 843       case Lanczos2DFilter:
 844       case Lanczos2DSharpFilter:
 845         /* Demote to a 2-lobe Sinc-Sinc for orthogonal use. */
 846         window_type=SincFastFilter;
 847         break;
 848       default:
 849         break;
 850     }
 851
 852   artifact=GetImageArtifact(image,"filter:filter");
 853   if (artifact != (const char *) NULL)
 854     {
 855       option=ParseMagickOption(MagickFilterOptions,MagickFalse,artifact);
 856       if ((UndefinedFilter < option) && (option < SentinelFilter))
 857         { /* Raw filter request - no window function. */
 858           filter_type=(FilterTypes) option;
 859           window_type=BoxFilter;
 860         }
 861       if (option == LanczosFilter)
 862         { /* Lanczos is not a real filter but a self windowing Sinc/Jinc. */
 863           filter_type=cylindrical != MagickFalse ? JincFilter : LanczosFilter;
 864           window_type=cylindrical != MagickFalse ? JincFilter : SincFastFilter;
 865         }
 866       /* Filter override with a specific window function. */
 867       artifact=GetImageArtifact(image,"filter:window");
 868       if (artifact != (const char *) NULL)
 869         {
 870           option=ParseMagickOption(MagickFilterOptions,MagickFalse,artifact);
 871           if ((UndefinedFilter < option) && (option < SentinelFilter))
 872             {
 873               if (option != LanczosFilter)
 874                 window_type=(FilterTypes) option;
 875               else
 876                 window_type=cylindrical != MagickFalse ? JincFilter :
 877                   SincFastFilter;
 878             }
 879         }
 880     }
 881   else
 882     {
 883       /* Window specified, but no filter function?  Assume Sinc/Jinc. */
 884       artifact=GetImageArtifact(image,"filter:window");
 885       if (artifact != (const char *) NULL)
 886         {
 887           option=ParseMagickOption(MagickFilterOptions,MagickFalse,
 888             artifact);
 889           if ((UndefinedFilter < option) && (option < SentinelFilter))
 890             {
 891               filter_type=cylindrical != MagickFalse ?
 892                          JincFilter : SincFastFilter;
 893               window_type=(FilterTypes) option;
 894             }
 895         }
 896     }
 897   /* Assign the real functions to use for the filters selected. */
 898   resize_filter->filter=filters[filter_type].function;
 899   resize_filter->support=filters[filter_type].lobes;
 900   resize_filter->window=filters[window_type].function;
 901   resize_filter->scale=filters[window_type].scale;
 902   resize_filter->signature=MagickSignature;
 903
 904   /* Filter Modifications for orthogonal/cylindrical usage */
 905   if (cylindrical != MagickFalse)
 906     switch (filter_type)
 907     {
 908       case PointFilter:
 909       case BoxFilter:
 910         /* Support for Cylindrical Box should be sqrt(2)/2 */
 911         resize_filter->support=(MagickRealType) MagickSQ1_2;
 912         break;
 913       default:
 914         break;
 915     }
 916   else
 917     switch (filter_type)
 918     {
 919       case Lanczos2DFilter:
 920       case Lanczos2DSharpFilter:
 921         /* Demote to a 2-lobe Lanczos (Sinc-Sinc) for orthogonal use. */
 922         resize_filter->filter=SincFast;
 923         break;
 924       default:
 925         break;
 926     }
 927
 928   /*
 929   ** More Expert Option Modifications
 930   */
 931
 932   /* User Sigma Override - no support change */
 933   artifact=GetImageArtifact(image,"filter:sigma");
 934   if (artifact != (const char *) NULL)
 935     sigma=StringToDouble(artifact);
 936   /* Define coefficents for Gaussian (assumes no cubic window) */
 937   if ( GaussianFilter ) {
 938     resize_filter->coeff[0] = 1.0/(2.0*sigma*sigma);
 939     resize_filter->coeff[1] = (MagickRealType) (1.0/(Magick2PI*sigma*sigma)); /* unused */
 940   }
 941
 942   /* Blur Override */
 943   artifact=GetImageArtifact(image,"filter:blur");
 944   if (artifact != (const char *) NULL)
 945     resize_filter->blur=StringToDouble(artifact);
 946   if (resize_filter->blur < MagickEpsilon)
 947     resize_filter->blur=(MagickRealType) MagickEpsilon;
 948
 949   /* Support Overrides */
 950   artifact=GetImageArtifact(image,"filter:lobes");
 951   if (artifact != (const char *) NULL)
 952     {
 953       ssize_t
 954         lobes;
 955
 956       lobes=(ssize_t) StringToLong(artifact);
 957       if (lobes < 1)
 958         lobes=1;
 959       resize_filter->support=(MagickRealType) lobes;
 960     }
 961   /* convert Jinc lobes to a real support value */
 962   if (resize_filter->filter == Jinc)
 963     {
 964       if (resize_filter->support > 16)
 965         resize_filter->support=jinc_zeros[15];  /* largest entry in table */
 966       else
 967         resize_filter->support = jinc_zeros[((long)resize_filter->support)-1];
 968     }
 969   /* expert override of the support setting */
 970   artifact=GetImageArtifact(image,"filter:support");
 971   if (artifact != (const char *) NULL)
 972     resize_filter->support=fabs(StringToDouble(artifact));
 973   /*
 974     Scale windowing function separatally to the support 'clipping'
 975     window that calling operator is planning to actually use. (Expert
 976     override)
 977   */
 978   resize_filter->window_support=resize_filter->support; /* default */
 979   artifact=GetImageArtifact(image,"filter:win-support");
 980   if (artifact != (const char *) NULL)
 981     resize_filter->window_support=fabs(StringToDouble(artifact));
 982   /*
 983     Adjust window function scaling to the windowing support for
 984     weighting function.  This avoids a division on every filter call.
 985   */
 986   resize_filter->scale /= resize_filter->window_support;
 987
 988   /*
 989    * Set Cubic Spline B,C values, calculate Cubic coefficients.
 990   */
 991   B=0.0;
 992   C=0.0;
 993   if ((filters[filter_type].function == CubicBC) ||
 994       (filters[window_type].function == CubicBC))
 995     {
 996       B=filters[filter_type].B;
 997       C=filters[filter_type].C;
 998       if (filters[window_type].function == CubicBC)
 999         {
1000           B=filters[window_type].B;
1001           C=filters[window_type].C;
1002         }
1003       artifact=GetImageArtifact(image,"filter:b");
1004       if (artifact != (const char *) NULL)
1005         {
1006           B=StringToDouble(artifact);
1007           C=(1.0-B)/2.0; /* Calculate C as if it is a Keys cubic filter. */
1008           artifact=GetImageArtifact(image,"filter:c"); /* user C override */
1009           if (artifact != (const char *) NULL)
1010             C=StringToDouble(artifact);
1011         }
1012       else
1013         {
1014           artifact=GetImageArtifact(image,"filter:c");
1015           if (artifact != (const char *) NULL)
1016             {
1017               C=StringToDouble(artifact);
1018               B=1.0-2.0*C;  /* Calculate B as if it is a Keys cubic filter. */
1019             }
1020         }
1021     /* Convert B,C values into Cubic Coefficents.  See CubicBC().  */
1022     resize_filter->coeff[0]=(6.0-2.0*B)/6.0;
1023     resize_filter->coeff[1]=0.0;
1024     resize_filter->coeff[2]=(-18.0+12.0*B+6.0*C)/6.0;
1025     resize_filter->coeff[3]=(12.0-9.0*B-6.0*C)/6.0;
1026     resize_filter->coeff[4]=(8.0*B+24.0*C)/6.0;
1027     resize_filter->coeff[5]=(-12.0*B-48.0*C)/6.0;
1028     resize_filter->coeff[6]=(6.0*B+30.0*C)/6.0;
1029     resize_filter->coeff[7]=(-B-6.0*C)/6.0;
1030   }
1031
1032   /*
1033     Expert Option Request for verbose details of the resulting filter.
1034   */
1035 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1036   #pragma omp master
1037   {
1038 #endif
1039     artifact=GetImageArtifact(image,"filter:verbose");
1040     if (artifact != (const char *) NULL)
1041       {
1042         double
1043       support,
1044           x;
1045
1046         /*
1047           Set the weighting function properly when the weighting
1048           function may not exactly match the filter of the same name.
1049           EG: a Point filter really uses a Box weighting function
1050           with a different support than is typically used.
1051
1052         */
1053         if (resize_filter->filter == Box)       filter_type=BoxFilter;
1054         if (resize_filter->filter == Sinc)      filter_type=SincFilter;
1055         if (resize_filter->filter == SincFast)  filter_type=SincFastFilter;
1056         if (resize_filter->filter == Jinc)      filter_type=JincFilter;
1057         if (resize_filter->filter == CubicBC)   filter_type=CubicFilter;
1058         /*
1059           Report Filter Details.
1060         */
1061         support=GetResizeFilterSupport(resize_filter); /* practical_support */
1062         (void) fprintf(stdout,"# Resize Filter (for graphing)\n#\n");
1063         (void) fprintf(stdout,"# filter = %s\n",MagickOptionToMnemonic(
1064            MagickFilterOptions,filter_type));
1065         (void) fprintf(stdout,"# window = %s\n",MagickOptionToMnemonic(
1066            MagickFilterOptions, window_type));
1067         (void) fprintf(stdout,"# support = %.*g\n",GetMagickPrecision(),
1068            (double) resize_filter->support);
1069         (void) fprintf(stdout,"# win-support = %.*g\n",GetMagickPrecision(),
1070            (double) resize_filter->window_support);
1071         (void) fprintf(stdout,"# scale_blur = %.*g\n",GetMagickPrecision(),
1072            (double) resize_filter->blur);
1073         if ( filter_type == GaussianFilter )
1074           (void) fprintf(stdout,"# gaussian_sigma = %.*g\n",GetMagickPrecision(),
1075              (double) sigma);
1076         (void) fprintf(stdout,"# practical_support = %.*g\n",GetMagickPrecision(),
1077            (double) support);
1078         if ( filter_type == CubicFilter || window_type == CubicFilter )
1079           (void) fprintf(stdout,"# B,C = %.*g,%.*g\n",GetMagickPrecision(),
1080              (double) B,GetMagickPrecision(),(double) C);
1081         (void) fprintf(stdout,"\n");
1082         /*
1083           Output values of resulting filter graph -- for graphing
1084           filter result.
1085         */
1086         for (x=0.0; x <= support; x+=0.01f)
1087           (void) fprintf(stdout,"%5.2lf\t%.*g\n",x,GetMagickPrecision(),
1088             (double) GetResizeFilterWeight(resize_filter,x));
1089         /* A final value so gnuplot can graph the 'stop' properly. */
1090         (void) fprintf(stdout,"%5.2lf\t%.*g\n",support,GetMagickPrecision(),
1091           0.0);
1092       }
1093       /* Output the above once only for each image - remove setting */
1094     (void) DeleteImageArtifact((Image *) image,"filter:verbose");
1095 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1096   }
1097 #endif
1098   return(resize_filter);
1099 }
1100 \f
1101 /*
1102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1103 %                                                                             %
1104 %                                                                             %
1105 %                                                                             %
1106 %   A d a p t i v e R e s i z e I m a g e                                     %
1107 %                                                                             %
1108 %                                                                             %
1109 %                                                                             %
1110 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1111 %
1112 %  AdaptiveResizeImage() adaptively resize image with pixel resampling.
1113 %
1114 %  The format of the AdaptiveResizeImage method is:
1115 %
1116 %      Image *AdaptiveResizeImage(const Image *image,const size_t columns,
1117 %        const size_t rows,ExceptionInfo *exception)
1118 %
1119 %  A description of each parameter follows:
1120 %
1121 %    o image: the image.
1122 %
1123 %    o columns: the number of columns in the resized image.
1124 %
1125 %    o rows: the number of rows in the resized image.
1126 %
1127 %    o exception: return any errors or warnings in this structure.
1128 %
1129 */
1130 MagickExport Image *AdaptiveResizeImage(const Image *image,
1131   const size_t columns,const size_t rows,ExceptionInfo *exception)
1132 {
1133 #define AdaptiveResizeImageTag  "Resize/Image"
1134
1135   CacheView
1136     *resize_view;
1137
1138   Image
1139     *resize_image;
1140
1141   MagickBooleanType
1142     proceed;
1143
1144   MagickPixelPacket
1145     pixel;
1146
1147   PointInfo
1148     offset;
1149
1150   ResampleFilter
1151     *resample_filter;
1152
1153   ssize_t
1154     y;
1155
1156   /*
1157     Adaptively resize image.
1158   */
1159   assert(image != (const Image *) NULL);
1160   assert(image->signature == MagickSignature);
1161   if (image->debug != MagickFalse)
1162     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1163   assert(exception != (ExceptionInfo *) NULL);
1164   assert(exception->signature == MagickSignature);
1165   if ((columns == 0) || (rows == 0))
1166     return((Image *) NULL);
1167   if ((columns == image->columns) && (rows == image->rows))
1168     return(CloneImage(image,0,0,MagickTrue,exception));
1169   resize_image=CloneImage(image,columns,rows,MagickTrue,exception);
1170   if (resize_image == (Image *) NULL)
1171     return((Image *) NULL);
1172   if (SetImageStorageClass(resize_image,DirectClass) == MagickFalse)
1173     {
1174       InheritException(exception,&resize_image->exception);
1175       resize_image=DestroyImage(resize_image);
1176       return((Image *) NULL);
1177     }
1178   GetMagickPixelPacket(image,&pixel);
1179   resample_filter=AcquireResampleFilter(image,exception);
1180   (void) SetResampleFilter(resample_filter,PointFilter,1.0);
1181   (void) SetResampleFilterInterpolateMethod(resample_filter,
1182     MeshInterpolatePixel);
1183   resize_view=AcquireCacheView(resize_image);
1184   for (y=0; y < (ssize_t) resize_image->rows; y++)
1185   {
1186     register IndexPacket
1187       *restrict resize_indexes;
1188
1189     register ssize_t
1190       x;
1191
1192     register PixelPacket
1193       *restrict q;
1194
1195     q=QueueCacheViewAuthenticPixels(resize_view,0,y,resize_image->columns,1,
1196       exception);
1197     if (q == (PixelPacket *) NULL)
1198       break;
1199     resize_indexes=GetCacheViewAuthenticIndexQueue(resize_view);
1200     offset.y=((MagickRealType) y*image->rows/resize_image->rows);
1201     for (x=0; x < (ssize_t) resize_image->columns; x++)
1202     {
1203       offset.x=((MagickRealType) x*image->columns/resize_image->columns);
1204       (void) ResamplePixelColor(resample_filter,offset.x-0.5,offset.y-0.5,
1205         &pixel);
1206       SetPixelPacket(resize_image,&pixel,q,resize_indexes+x);
1207       q++;
1208     }
1209     if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
1210       break;
1211     proceed=SetImageProgress(image,AdaptiveResizeImageTag,(MagickOffsetType) y,
1212       image->rows);
1213     if (proceed == MagickFalse)
1214       break;
1215   }
1216   resample_filter=DestroyResampleFilter(resample_filter);
1217   resize_view=DestroyCacheView(resize_view);
1218   return(resize_image);
1219 }
1220 \f
1221 /*
1222 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1223 %                                                                             %
1224 %                                                                             %
1225 %                                                                             %
1226 +   B e s s e l O r d e r O n e                                               %
1227 %                                                                             %
1228 %                                                                             %
1229 %                                                                             %
1230 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1231 %
1232 %  BesselOrderOne() computes the Bessel function of x of the first kind of
1233 %  order 0.  This is used to create the Jinc() filter function below.
1234 %
1235 %    Reduce x to |x| since j1(x)= -j1(-x), and for x in (0,8]
1236 %
1237 %       j1(x) = x*j1(x);
1238 %
1239 %    For x in (8,inf)
1240 %
1241 %       j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
1242 %
1243 %    where x1 = x-3*pi/4. Compute sin(x1) and cos(x1) as follow:
1244 %
1245 %       cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
1246 %               =  1/sqrt(2) * (sin(x) - cos(x))
1247 %       sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
1248 %               = -1/sqrt(2) * (sin(x) + cos(x))
1249 %
1250 %  The format of the BesselOrderOne method is:
1251 %
1252 %      MagickRealType BesselOrderOne(MagickRealType x)
1253 %
1254 %  A description of each parameter follows:
1255 %
1256 %    o x: MagickRealType value.
1257 %
1258 */
1259
1260 #undef I0
1261 static MagickRealType I0(MagickRealType x)
1262 {
1263   MagickRealType
1264     sum,
1265     t,
1266     y;
1267
1268   register ssize_t
1269     i;
1270
1271   /*
1272     Zeroth order Bessel function of the first kind.
1273   */
1274   sum=1.0;
1275   y=x*x/4.0;
1276   t=y;
1277   for (i=2; t > MagickEpsilon; i++)
1278   {
1279     sum+=t;
1280     t*=y/((MagickRealType) i*i);
1281   }
1282   return(sum);
1283 }
1284
1285 #undef J1
1286 static MagickRealType J1(MagickRealType x)
1287 {
1288   MagickRealType
1289     p,
1290     q;
1291
1292   register ssize_t
1293     i;
1294
1295   static const double
1296     Pone[] =
1297     {
1298        0.581199354001606143928050809e+21,
1299       -0.6672106568924916298020941484e+20,
1300        0.2316433580634002297931815435e+19,
1301       -0.3588817569910106050743641413e+17,
1302        0.2908795263834775409737601689e+15,
1303       -0.1322983480332126453125473247e+13,
1304        0.3413234182301700539091292655e+10,
1305       -0.4695753530642995859767162166e+7,
1306        0.270112271089232341485679099e+4
1307     },
1308     Qone[] =
1309     {
1310       0.11623987080032122878585294e+22,
1311       0.1185770712190320999837113348e+20,
1312       0.6092061398917521746105196863e+17,
1313       0.2081661221307607351240184229e+15,
1314       0.5243710262167649715406728642e+12,
1315       0.1013863514358673989967045588e+10,
1316       0.1501793594998585505921097578e+7,
1317       0.1606931573481487801970916749e+4,
1318       0.1e+1
1319     };
1320
1321   p=Pone[8];
1322   q=Qone[8];
1323   for (i=7; i >= 0; i--)
1324   {
1325     p=p*x*x+Pone[i];
1326     q=q*x*x+Qone[i];
1327   }
1328   return(p/q);
1329 }
1330
1331 #undef P1
1332 static MagickRealType P1(MagickRealType x)
1333 {
1334   MagickRealType
1335     p,
1336     q;
1337
1338   register ssize_t
1339     i;
1340
1341   static const double
1342     Pone[] =
1343     {
1344       0.352246649133679798341724373e+5,
1345       0.62758845247161281269005675e+5,
1346       0.313539631109159574238669888e+5,
1347       0.49854832060594338434500455e+4,
1348       0.2111529182853962382105718e+3,
1349       0.12571716929145341558495e+1
1350     },
1351     Qone[] =
1352     {
1353       0.352246649133679798068390431e+5,
1354       0.626943469593560511888833731e+5,
1355       0.312404063819041039923015703e+5,
1356       0.4930396490181088979386097e+4,
1357       0.2030775189134759322293574e+3,
1358       0.1e+1
1359     };
1360
1361   p=Pone[5];
1362   q=Qone[5];
1363   for (i=4; i >= 0; i--)
1364   {
1365     p=p*(8.0/x)*(8.0/x)+Pone[i];
1366     q=q*(8.0/x)*(8.0/x)+Qone[i];
1367   }
1368   return(p/q);
1369 }
1370
1371 #undef Q1
1372 static MagickRealType Q1(MagickRealType x)
1373 {
1374   MagickRealType
1375     p,
1376     q;
1377
1378   register ssize_t
1379     i;
1380
1381   static const double
1382     Pone[] =
1383     {
1384       0.3511751914303552822533318e+3,
1385       0.7210391804904475039280863e+3,
1386       0.4259873011654442389886993e+3,
1387       0.831898957673850827325226e+2,
1388       0.45681716295512267064405e+1,
1389       0.3532840052740123642735e-1
1390     },
1391     Qone[] =
1392     {
1393       0.74917374171809127714519505e+4,
1394       0.154141773392650970499848051e+5,
1395       0.91522317015169922705904727e+4,
1396       0.18111867005523513506724158e+4,
1397       0.1038187585462133728776636e+3,
1398       0.1e+1
1399     };
1400
1401   p=Pone[5];
1402   q=Qone[5];
1403   for (i=4; i >= 0; i--)
1404   {
1405     p=p*(8.0/x)*(8.0/x)+Pone[i];
1406     q=q*(8.0/x)*(8.0/x)+Qone[i];
1407   }
1408   return(p/q);
1409 }
1410
1411 static MagickRealType BesselOrderOne(MagickRealType x)
1412 {
1413   MagickRealType
1414     p,
1415     q;
1416
1417   if (x == 0.0)
1418     return(0.0);
1419   p=x;
1420   if (x < 0.0)
1421     x=(-x);
1422   if (x < 8.0)
1423     return(p*J1(x));
1424   q=sqrt((double) (2.0/(MagickPI*x)))*(P1(x)*(1.0/sqrt(2.0)*(sin((double) x)-
1425     cos((double) x)))-8.0/x*Q1(x)*(-1.0/sqrt(2.0)*(sin((double) x)+
1426     cos((double) x))));
1427   if (p < 0.0)
1428     q=(-q);
1429   return(q);
1430 }
1431 \f
1432 /*
1433 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1434 %                                                                             %
1435 %                                                                             %
1436 %                                                                             %
1437 +   D e s t r o y R e s i z e F i l t e r                                     %
1438 %                                                                             %
1439 %                                                                             %
1440 %                                                                             %
1441 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1442 %
1443 %  DestroyResizeFilter() destroy the resize filter.
1444 %
1445 %  The format of the DestroyResizeFilter method is:
1446 %
1447 %      ResizeFilter *DestroyResizeFilter(ResizeFilter *resize_filter)
1448 %
1449 %  A description of each parameter follows:
1450 %
1451 %    o resize_filter: the resize filter.
1452 %
1453 */
1454 MagickExport ResizeFilter *DestroyResizeFilter(ResizeFilter *resize_filter)
1455 {
1456   assert(resize_filter != (ResizeFilter *) NULL);
1457   assert(resize_filter->signature == MagickSignature);
1458   resize_filter->signature=(~MagickSignature);
1459   resize_filter=(ResizeFilter *) RelinquishMagickMemory(resize_filter);
1460   return(resize_filter);
1461 }
1462 \f
1463 /*
1464 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1465 %                                                                             %
1466 %                                                                             %
1467 %                                                                             %
1468 +   G e t R e s i z e F i l t e r S u p p o r t                               %
1469 %                                                                             %
1470 %                                                                             %
1471 %                                                                             %
1472 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1473 %
1474 %  GetResizeFilterSupport() return the current support window size for this
1475 %  filter.  Note that this may have been enlarged by filter:blur factor.
1476 %
1477 %  The format of the GetResizeFilterSupport method is:
1478 %
1479 %      MagickRealType GetResizeFilterSupport(const ResizeFilter *resize_filter)
1480 %
1481 %  A description of each parameter follows:
1482 %
1483 %    o filter: Image filter to use.
1484 %
1485 */
1486 MagickExport MagickRealType GetResizeFilterSupport(
1487   const ResizeFilter *resize_filter)
1488 {
1489   assert(resize_filter != (ResizeFilter *) NULL);
1490   assert(resize_filter->signature == MagickSignature);
1491   return(resize_filter->support*resize_filter->blur);
1492 }
1493 \f
1494 /*
1495 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1496 %                                                                             %
1497 %                                                                             %
1498 %                                                                             %
1499 +   G e t R e s i z e F i l t e r W e i g h t                                 %
1500 %                                                                             %
1501 %                                                                             %
1502 %                                                                             %
1503 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1504 %
1505 %  GetResizeFilterWeight evaluates the specified resize filter at the point x
1506 %  which usally lies between zero and the filters current 'support' and
1507 %  returns the weight of the filter function at that point.
1508 %
1509 %  The format of the GetResizeFilterWeight method is:
1510 %
1511 %      MagickRealType GetResizeFilterWeight(const ResizeFilter *resize_filter,
1512 %        const MagickRealType x)
1513 %
1514 %  A description of each parameter follows:
1515 %
1516 %    o filter: the filter type.
1517 %
1518 %    o x: the point.
1519 %
1520 */
1521 MagickExport MagickRealType GetResizeFilterWeight(
1522   const ResizeFilter *resize_filter,const MagickRealType x)
1523 {
1524   MagickRealType
1525     scale,
1526     x_blur;
1527
1528   /*
1529     Windowing function - scale the weighting filter by this amount.
1530   */
1531   assert(resize_filter != (ResizeFilter *) NULL);
1532   assert(resize_filter->signature == MagickSignature);
1533   x_blur=fabs((double) x)/resize_filter->blur;  /* X offset with blur scaling */
1534   if ((resize_filter->window_support < MagickEpsilon) ||
1535       (resize_filter->window == Box))
1536     scale=1.0;  /* Point or Box Filter -- avoid division by zero */
1537   else
1538     {
1539       scale=resize_filter->scale;
1540       scale=resize_filter->window(x_blur*scale,resize_filter);
1541     }
1542   return(scale*resize_filter->filter(x_blur,resize_filter));
1543 }
1544 \f
1545 /*
1546 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1547 %                                                                             %
1548 %                                                                             %
1549 %                                                                             %
1550 %   M a g n i f y I m a g e                                                   %
1551 %                                                                             %
1552 %                                                                             %
1553 %                                                                             %
1554 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1555 %
1556 %  MagnifyImage() is a convenience method that scales an image proportionally
1557 %  to twice its size.
1558 %
1559 %  The format of the MagnifyImage method is:
1560 %
1561 %      Image *MagnifyImage(const Image *image,ExceptionInfo *exception)
1562 %
1563 %  A description of each parameter follows:
1564 %
1565 %    o image: the image.
1566 %
1567 %    o exception: return any errors or warnings in this structure.
1568 %
1569 */
1570 MagickExport Image *MagnifyImage(const Image *image,ExceptionInfo *exception)
1571 {
1572   Image
1573     *magnify_image;
1574
1575   assert(image != (Image *) NULL);
1576   assert(image->signature == MagickSignature);
1577   if (image->debug != MagickFalse)
1578     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1579   assert(exception != (ExceptionInfo *) NULL);
1580   assert(exception->signature == MagickSignature);
1581   magnify_image=ResizeImage(image,2*image->columns,2*image->rows,CubicFilter,
1582     1.0,exception);
1583   return(magnify_image);
1584 }
1585 \f
1586 /*
1587 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1588 %                                                                             %
1589 %                                                                             %
1590 %                                                                             %
1591 %   M i n i f y I m a g e                                                     %
1592 %                                                                             %
1593 %                                                                             %
1594 %                                                                             %
1595 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1596 %
1597 %  MinifyImage() is a convenience method that scales an image proportionally
1598 %  to half its size.
1599 %
1600 %  The format of the MinifyImage method is:
1601 %
1602 %      Image *MinifyImage(const Image *image,ExceptionInfo *exception)
1603 %
1604 %  A description of each parameter follows:
1605 %
1606 %    o image: the image.
1607 %
1608 %    o exception: return any errors or warnings in this structure.
1609 %
1610 */
1611 MagickExport Image *MinifyImage(const Image *image,ExceptionInfo *exception)
1612 {
1613   Image
1614     *minify_image;
1615
1616   assert(image != (Image *) NULL);
1617   assert(image->signature == MagickSignature);
1618   if (image->debug != MagickFalse)
1619     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1620   assert(exception != (ExceptionInfo *) NULL);
1621   assert(exception->signature == MagickSignature);
1622   minify_image=ResizeImage(image,image->columns/2,image->rows/2,CubicFilter,
1623     1.0,exception);
1624   return(minify_image);
1625 }
1626 \f
1627 /*
1628 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1629 %                                                                             %
1630 %                                                                             %
1631 %                                                                             %
1632 %   R e s a m p l e I m a g e                                                 %
1633 %                                                                             %
1634 %                                                                             %
1635 %                                                                             %
1636 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1637 %
1638 %  ResampleImage() resize image in terms of its pixel size, so that when
1639 %  displayed at the given resolution it will be the same size in terms of
1640 %  real world units as the original image at the original resolution.
1641 %
1642 %  The format of the ResampleImage method is:
1643 %
1644 %      Image *ResampleImage(Image *image,const double x_resolution,
1645 %        const double y_resolution,const FilterTypes filter,const double blur,
1646 %        ExceptionInfo *exception)
1647 %
1648 %  A description of each parameter follows:
1649 %
1650 %    o image: the image to be resized to fit the given resolution.
1651 %
1652 %    o x_resolution: the new image x resolution.
1653 %
1654 %    o y_resolution: the new image y resolution.
1655 %
1656 %    o filter: Image filter to use.
1657 %
1658 %    o blur: the blur factor where > 1 is blurry, < 1 is sharp.
1659 %
1660 */
1661 MagickExport Image *ResampleImage(const Image *image,const double x_resolution,
1662   const double y_resolution,const FilterTypes filter,const double blur,
1663   ExceptionInfo *exception)
1664 {
1665 #define ResampleImageTag  "Resample/Image"
1666
1667   Image
1668     *resample_image;
1669
1670   size_t
1671     height,
1672     width;
1673
1674   /*
1675     Initialize sampled image attributes.
1676   */
1677   assert(image != (const Image *) NULL);
1678   assert(image->signature == MagickSignature);
1679   if (image->debug != MagickFalse)
1680     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1681   assert(exception != (ExceptionInfo *) NULL);
1682   assert(exception->signature == MagickSignature);
1683   width=(size_t) (x_resolution*image->columns/(image->x_resolution == 0.0 ?
1684     72.0 : image->x_resolution)+0.5);
1685   height=(size_t) (y_resolution*image->rows/(image->y_resolution == 0.0 ?
1686     72.0 : image->y_resolution)+0.5);
1687   resample_image=ResizeImage(image,width,height,filter,blur,exception);
1688   if (resample_image != (Image *) NULL)
1689     {
1690       resample_image->x_resolution=x_resolution;
1691       resample_image->y_resolution=y_resolution;
1692     }
1693   return(resample_image);
1694 }
1695 #if defined(MAGICKCORE_LQR_DELEGATE)
1696 \f
1697 /*
1698 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1699 %                                                                             %
1700 %                                                                             %
1701 %                                                                             %
1702 %   L i q u i d R e s c a l e I m a g e                                       %
1703 %                                                                             %
1704 %                                                                             %
1705 %                                                                             %
1706 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1707 %
1708 %  LiquidRescaleImage() rescales image with seam carving.
1709 %
1710 %  The format of the LiquidRescaleImage method is:
1711 %
1712 %      Image *LiquidRescaleImage(const Image *image,
1713 %        const size_t columns,const size_t rows,
1714 %        const double delta_x,const double rigidity,ExceptionInfo *exception)
1715 %
1716 %  A description of each parameter follows:
1717 %
1718 %    o image: the image.
1719 %
1720 %    o columns: the number of columns in the rescaled image.
1721 %
1722 %    o rows: the number of rows in the rescaled image.
1723 %
1724 %    o delta_x: maximum seam transversal step (0 means straight seams).
1725 %
1726 %    o rigidity: introduce a bias for non-straight seams (typically 0).
1727 %
1728 %    o exception: return any errors or warnings in this structure.
1729 %
1730 */
1731 MagickExport Image *LiquidRescaleImage(const Image *image,const size_t columns,
1732   const size_t rows,const double delta_x,const double rigidity,
1733   ExceptionInfo *exception)
1734 {
1735 #define LiquidRescaleImageTag  "Rescale/Image"
1736
1737   CacheView
1738     *rescale_view;
1739
1740   const char
1741     *map;
1742
1743   guchar
1744     *packet;
1745
1746   Image
1747     *rescale_image;
1748
1749   int
1750     x,
1751     y;
1752
1753   LqrCarver
1754     *carver;
1755
1756   LqrRetVal
1757     lqr_status;
1758
1759   MagickBooleanType
1760     status;
1761
1762   MagickPixelPacket
1763     pixel;
1764
1765   unsigned char
1766     *pixels;
1767
1768   /*
1769     Liquid rescale image.
1770   */
1771   assert(image != (const Image *) NULL);
1772   assert(image->signature == MagickSignature);
1773   if (image->debug != MagickFalse)
1774     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1775   assert(exception != (ExceptionInfo *) NULL);
1776   assert(exception->signature == MagickSignature);
1777   if ((columns == 0) || (rows == 0))
1778     return((Image *) NULL);
1779   if ((columns == image->columns) && (rows == image->rows))
1780     return(CloneImage(image,0,0,MagickTrue,exception));
1781   if ((columns <= 2) || (rows <= 2))
1782     return(ResizeImage(image,columns,rows,image->filter,image->blur,exception));
1783   if ((columns >= (2*image->columns)) || (rows >= (2*image->rows)))
1784     {
1785       Image
1786         *resize_image;
1787
1788       size_t
1789         height,
1790         width;
1791
1792       /*
1793         Honor liquid resize size limitations.
1794       */
1795       for (width=image->columns; columns >= (2*width-1); width*=2);
1796       for (height=image->rows; rows >= (2*height-1); height*=2);
1797       resize_image=ResizeImage(image,width,height,image->filter,image->blur,
1798         exception);
1799       if (resize_image == (Image *) NULL)
1800         return((Image *) NULL);
1801       rescale_image=LiquidRescaleImage(resize_image,columns,rows,delta_x,
1802         rigidity,exception);
1803       resize_image=DestroyImage(resize_image);
1804       return(rescale_image);
1805     }
1806   map="RGB";
1807   if (image->matte == MagickFalse)
1808     map="RGBA";
1809   if (image->colorspace == CMYKColorspace)
1810     {
1811       map="CMYK";
1812       if (image->matte == MagickFalse)
1813         map="CMYKA";
1814     }
1815   pixels=(unsigned char *) AcquireQuantumMemory(image->columns,image->rows*
1816     strlen(map)*sizeof(*pixels));
1817   if (pixels == (unsigned char *) NULL)
1818     return((Image *) NULL);
1819   status=ExportImagePixels(image,0,0,image->columns,image->rows,map,CharPixel,
1820     pixels,exception);
1821   if (status == MagickFalse)
1822     {
1823       pixels=(unsigned char *) RelinquishMagickMemory(pixels);
1824       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1825     }
1826   carver=lqr_carver_new(pixels,image->columns,image->rows,strlen(map));
1827   if (carver == (LqrCarver *) NULL)
1828     {
1829       pixels=(unsigned char *) RelinquishMagickMemory(pixels);
1830       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1831     }
1832   lqr_status=lqr_carver_init(carver,(int) delta_x,rigidity);
1833   lqr_status=lqr_carver_resize(carver,columns,rows);
1834   rescale_image=CloneImage(image,lqr_carver_get_width(carver),
1835     lqr_carver_get_height(carver),MagickTrue,exception);
1836   if (rescale_image == (Image *) NULL)
1837     {
1838       pixels=(unsigned char *) RelinquishMagickMemory(pixels);
1839       return((Image *) NULL);
1840     }
1841   if (SetImageStorageClass(rescale_image,DirectClass) == MagickFalse)
1842     {
1843       InheritException(exception,&rescale_image->exception);
1844       rescale_image=DestroyImage(rescale_image);
1845       return((Image *) NULL);
1846     }
1847   GetMagickPixelPacket(rescale_image,&pixel);
1848   (void) lqr_carver_scan_reset(carver);
1849   rescale_view=AcquireCacheView(rescale_image);
1850   while (lqr_carver_scan(carver,&x,&y,&packet) != 0)
1851   {
1852     register IndexPacket
1853       *restrict rescale_indexes;
1854
1855     register PixelPacket
1856       *restrict q;
1857
1858     q=QueueCacheViewAuthenticPixels(rescale_view,x,y,1,1,exception);
1859     if (q == (PixelPacket *) NULL)
1860       break;
1861     rescale_indexes=GetCacheViewAuthenticIndexQueue(rescale_view);
1862     pixel.red=QuantumRange*(packet[0]/255.0);
1863     pixel.green=QuantumRange*(packet[1]/255.0);
1864     pixel.blue=QuantumRange*(packet[2]/255.0);
1865     if (image->colorspace != CMYKColorspace)
1866       {
1867         if (image->matte == MagickFalse)
1868           pixel.opacity=QuantumRange*(packet[3]/255.0);
1869       }
1870     else
1871       {
1872         pixel.index=QuantumRange*(packet[3]/255.0);
1873         if (image->matte == MagickFalse)
1874           pixel.opacity=QuantumRange*(packet[4]/255.0);
1875       }
1876     SetPixelPacket(rescale_image,&pixel,q,rescale_indexes);
1877     if (SyncCacheViewAuthenticPixels(rescale_view,exception) == MagickFalse)
1878       break;
1879   }
1880   rescale_view=DestroyCacheView(rescale_view);
1881   /*
1882     Relinquish resources.
1883   */
1884   lqr_carver_destroy(carver);
1885   return(rescale_image);
1886 }
1887 #else
1888 MagickExport Image *LiquidRescaleImage(const Image *image,
1889   const size_t magick_unused(columns),const size_t magick_unused(rows),
1890   const double magick_unused(delta_x),const double magick_unused(rigidity),
1891   ExceptionInfo *exception)
1892 {
1893   assert(image != (const Image *) NULL);
1894   assert(image->signature == MagickSignature);
1895   if (image->debug != MagickFalse)
1896     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1897   assert(exception != (ExceptionInfo *) NULL);
1898   assert(exception->signature == MagickSignature);
1899   (void) ThrowMagickException(exception,GetMagickModule(),MissingDelegateError,
1900     "DelegateLibrarySupportNotBuiltIn","`%s' (LQR)",image->filename);
1901   return((Image *) NULL);
1902 }
1903 #endif
1904 \f
1905 /*
1906 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1907 %                                                                             %
1908 %                                                                             %
1909 %                                                                             %
1910 %   R e s i z e I m a g e                                                     %
1911 %                                                                             %
1912 %                                                                             %
1913 %                                                                             %
1914 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1915 %
1916 %  ResizeImage() scales an image to the desired dimensions, using the given
1917 %  filter (see AcquireFilterInfo()).
1918 %
1919 %  If an undefined filter is given the filter defaults to Mitchell for a
1920 %  colormapped image, a image with a matte channel, or if the image is
1921 %  enlarged.  Otherwise the filter defaults to a Lanczos.
1922 %
1923 %  ResizeImage() was inspired by Paul Heckbert's "zoom" program.
1924 %
1925 %  The format of the ResizeImage method is:
1926 %
1927 %      Image *ResizeImage(Image *image,const size_t columns,
1928 %        const size_t rows,const FilterTypes filter,const double blur,
1929 %        ExceptionInfo *exception)
1930 %
1931 %  A description of each parameter follows:
1932 %
1933 %    o image: the image.
1934 %
1935 %    o columns: the number of columns in the scaled image.
1936 %
1937 %    o rows: the number of rows in the scaled image.
1938 %
1939 %    o filter: Image filter to use.
1940 %
1941 %    o blur: the blur factor where > 1 is blurry, < 1 is sharp.  Typically set
1942 %      this to 1.0.
1943 %
1944 %    o exception: return any errors or warnings in this structure.
1945 %
1946 */
1947
1948 typedef struct _ContributionInfo
1949 {
1950   MagickRealType
1951     weight;
1952
1953   ssize_t
1954     pixel;
1955 } ContributionInfo;
1956
1957 static ContributionInfo **DestroyContributionThreadSet(
1958   ContributionInfo **contribution)
1959 {
1960   register ssize_t
1961     i;
1962
1963   assert(contribution != (ContributionInfo **) NULL);
1964   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
1965     if (contribution[i] != (ContributionInfo *) NULL)
1966       contribution[i]=(ContributionInfo *) RelinquishMagickMemory(
1967         contribution[i]);
1968   contribution=(ContributionInfo **) RelinquishMagickMemory(contribution);
1969   return(contribution);
1970 }
1971
1972 static ContributionInfo **AcquireContributionThreadSet(const size_t count)
1973 {
1974   register ssize_t
1975     i;
1976
1977   ContributionInfo
1978     **contribution;
1979
1980   size_t
1981     number_threads;
1982
1983   number_threads=GetOpenMPMaximumThreads();
1984   contribution=(ContributionInfo **) AcquireQuantumMemory(number_threads,
1985     sizeof(*contribution));
1986   if (contribution == (ContributionInfo **) NULL)
1987     return((ContributionInfo **) NULL);
1988   (void) ResetMagickMemory(contribution,0,number_threads*sizeof(*contribution));
1989   for (i=0; i < (ssize_t) number_threads; i++)
1990   {
1991     contribution[i]=(ContributionInfo *) AcquireQuantumMemory(count,
1992       sizeof(**contribution));
1993     if (contribution[i] == (ContributionInfo *) NULL)
1994       return(DestroyContributionThreadSet(contribution));
1995   }
1996   return(contribution);
1997 }
1998
1999 static inline double MagickMax(const double x,const double y)
2000 {
2001   if (x > y)
2002     return(x);
2003   return(y);
2004 }
2005
2006 static inline double MagickMin(const double x,const double y)
2007 {
2008   if (x < y)
2009     return(x);
2010   return(y);
2011 }
2012
2013 static MagickBooleanType HorizontalFilter(const ResizeFilter *resize_filter,
2014   const Image *image,Image *resize_image,const MagickRealType x_factor,
2015   const MagickSizeType span,MagickOffsetType *offset,ExceptionInfo *exception)
2016 {
2017 #define ResizeImageTag  "Resize/Image"
2018
2019   CacheView
2020     *image_view,
2021     *resize_view;
2022
2023   ClassType
2024     storage_class;
2025
2026   ContributionInfo
2027     **restrict contributions;
2028
2029   MagickBooleanType
2030     status;
2031
2032   MagickPixelPacket
2033     zero;
2034
2035   MagickRealType
2036     scale,
2037     support;
2038
2039   ssize_t
2040     x;
2041
2042   /*
2043     Apply filter to resize horizontally from image to resize image.
2044   */
2045   scale=MagickMax(1.0/x_factor+MagickEpsilon,1.0);
2046   support=scale*GetResizeFilterSupport(resize_filter);
2047   storage_class=support > 0.5 ? DirectClass : image->storage_class;
2048   if (SetImageStorageClass(resize_image,storage_class) == MagickFalse)
2049     {
2050       InheritException(exception,&resize_image->exception);
2051       return(MagickFalse);
2052     }
2053   if (support < 0.5)
2054     {
2055       /*
2056         Support too small even for nearest neighbour: Reduce to point
2057         sampling.
2058       */
2059       support=(MagickRealType) 0.5;
2060       scale=1.0;
2061     }
2062   contributions=AcquireContributionThreadSet((size_t) (2.0*support+3.0));
2063   if (contributions == (ContributionInfo **) NULL)
2064     {
2065       (void) ThrowMagickException(exception,GetMagickModule(),
2066         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2067       return(MagickFalse);
2068     }
2069   status=MagickTrue;
2070   scale=1.0/scale;
2071   (void) ResetMagickMemory(&zero,0,sizeof(zero));
2072   image_view=AcquireCacheView(image);
2073   resize_view=AcquireCacheView(resize_image);
2074 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2075   #pragma omp parallel for shared(status)
2076 #endif
2077   for (x=0; x < (ssize_t) resize_image->columns; x++)
2078   {
2079     MagickRealType
2080       center,
2081       density;
2082
2083     register const IndexPacket
2084       *restrict indexes;
2085
2086     register const PixelPacket
2087       *restrict p;
2088
2089     register ContributionInfo
2090       *restrict contribution;
2091
2092     register IndexPacket
2093       *restrict resize_indexes;
2094
2095     register PixelPacket
2096       *restrict q;
2097
2098     register ssize_t
2099       y;
2100
2101     ssize_t
2102       n,
2103       start,
2104       stop;
2105
2106     if (status == MagickFalse)
2107       continue;
2108     center=(MagickRealType) (x+0.5)/x_factor;
2109     start=(ssize_t) MagickMax(center-support+0.5,0.0);
2110     stop=(ssize_t) MagickMin(center+support+0.5,(double) image->columns);
2111     density=0.0;
2112     contribution=contributions[GetOpenMPThreadId()];
2113     for (n=0; n < (stop-start); n++)
2114     {
2115       contribution[n].pixel=start+n;
2116       contribution[n].weight=GetResizeFilterWeight(resize_filter,scale*
2117         ((MagickRealType) (start+n)-center+0.5));
2118       density+=contribution[n].weight;
2119     }
2120     if ((density != 0.0) && (density != 1.0))
2121       {
2122         register ssize_t
2123           i;
2124
2125         /*
2126           Normalize.
2127         */
2128         density=1.0/density;
2129         for (i=0; i < n; i++)
2130           contribution[i].weight*=density;
2131       }
2132     p=GetCacheViewVirtualPixels(image_view,contribution[0].pixel,0,(size_t)
2133       (contribution[n-1].pixel-contribution[0].pixel+1),image->rows,exception);
2134     q=QueueCacheViewAuthenticPixels(resize_view,x,0,1,resize_image->rows,
2135       exception);
2136     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2137       {
2138         status=MagickFalse;
2139         continue;
2140       }
2141     indexes=GetCacheViewVirtualIndexQueue(image_view);
2142     resize_indexes=GetCacheViewAuthenticIndexQueue(resize_view);
2143     for (y=0; y < (ssize_t) resize_image->rows; y++)
2144     {
2145       MagickPixelPacket
2146         pixel;
2147
2148       MagickRealType
2149         alpha;
2150
2151       register ssize_t
2152         i;
2153
2154       ssize_t
2155         j;
2156
2157       pixel=zero;
2158       if (image->matte == MagickFalse)
2159         {
2160           for (i=0; i < n; i++)
2161           {
2162             j=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2163               (contribution[i].pixel-contribution[0].pixel);
2164             alpha=contribution[i].weight;
2165             pixel.red+=alpha*(p+j)->red;
2166             pixel.green+=alpha*(p+j)->green;
2167             pixel.blue+=alpha*(p+j)->blue;
2168             pixel.opacity+=alpha*(p+j)->opacity;
2169           }
2170           SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
2171           SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
2172           SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
2173           SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
2174           if ((image->colorspace == CMYKColorspace) &&
2175               (resize_image->colorspace == CMYKColorspace))
2176             {
2177               for (i=0; i < n; i++)
2178               {
2179                 j=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2180                   (contribution[i].pixel-contribution[0].pixel);
2181                 alpha=contribution[i].weight;
2182                 pixel.index+=alpha*indexes[j];
2183               }
2184               resize_indexes[y]=(IndexPacket) ClampToQuantum(pixel.index);
2185             }
2186         }
2187       else
2188         {
2189           MagickRealType
2190             gamma;
2191
2192           gamma=0.0;
2193           for (i=0; i < n; i++)
2194           {
2195             j=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2196               (contribution[i].pixel-contribution[0].pixel);
2197             alpha=contribution[i].weight*QuantumScale*
2198               GetAlphaPixelComponent(p+j);
2199             pixel.red+=alpha*(p+j)->red;
2200             pixel.green+=alpha*(p+j)->green;
2201             pixel.blue+=alpha*(p+j)->blue;
2202             pixel.opacity+=contribution[i].weight*(p+j)->opacity;
2203             gamma+=alpha;
2204           }
2205           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2206           q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
2207           q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
2208           q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
2209           SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
2210           if ((image->colorspace == CMYKColorspace) &&
2211               (resize_image->colorspace == CMYKColorspace))
2212             {
2213               for (i=0; i < n; i++)
2214               {
2215                 j=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2216                   (contribution[i].pixel-contribution[0].pixel);
2217                 alpha=contribution[i].weight*QuantumScale*
2218                   GetAlphaPixelComponent(p+j);
2219                 pixel.index+=alpha*indexes[j];
2220               }
2221               resize_indexes[y]=(IndexPacket) ClampToQuantum(gamma*
2222                GetIndexPixelComponent(&pixel));
2223             }
2224         }
2225       if ((resize_image->storage_class == PseudoClass) &&
2226           (image->storage_class == PseudoClass))
2227         {
2228           i=(ssize_t) (MagickMin(MagickMax(center,(double) start),(double) stop-
2229             1.0)+0.5);
2230           j=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2231             (contribution[i-start].pixel-contribution[0].pixel);
2232           resize_indexes[y]=indexes[j];
2233         }
2234       q++;
2235     }
2236     if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
2237       status=MagickFalse;
2238     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2239       {
2240         MagickBooleanType
2241           proceed;
2242
2243 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2244   #pragma omp critical (MagickCore_HorizontalFilter)
2245 #endif
2246         proceed=SetImageProgress(image,ResizeImageTag,(*offset)++,span);
2247         if (proceed == MagickFalse)
2248           status=MagickFalse;
2249       }
2250   }
2251   resize_view=DestroyCacheView(resize_view);
2252   image_view=DestroyCacheView(image_view);
2253   contributions=DestroyContributionThreadSet(contributions);
2254   return(status);
2255 }
2256
2257 static MagickBooleanType VerticalFilter(const ResizeFilter *resize_filter,
2258   const Image *image,Image *resize_image,const MagickRealType y_factor,
2259   const MagickSizeType span,MagickOffsetType *offset,ExceptionInfo *exception)
2260 {
2261   CacheView
2262     *image_view,
2263     *resize_view;
2264
2265   ClassType
2266     storage_class;
2267
2268   ContributionInfo
2269     **restrict contributions;
2270
2271   MagickBooleanType
2272     status;
2273
2274   MagickPixelPacket
2275     zero;
2276
2277   MagickRealType
2278     scale,
2279     support;
2280
2281   ssize_t
2282     y;
2283
2284   /*
2285     Apply filter to resize vertically from image to resize image.
2286   */
2287   scale=MagickMax(1.0/y_factor+MagickEpsilon,1.0);
2288   support=scale*GetResizeFilterSupport(resize_filter);
2289   storage_class=support > 0.5 ? DirectClass : image->storage_class;
2290   if (SetImageStorageClass(resize_image,storage_class) == MagickFalse)
2291     {
2292       InheritException(exception,&resize_image->exception);
2293       return(MagickFalse);
2294     }
2295   if (support < 0.5)
2296     {
2297       /*
2298         Support too small even for nearest neighbour: Reduce to point
2299         sampling.
2300       */
2301       support=(MagickRealType) 0.5;
2302       scale=1.0;
2303     }
2304   contributions=AcquireContributionThreadSet((size_t) (2.0*support+3.0));
2305   if (contributions == (ContributionInfo **) NULL)
2306     {
2307       (void) ThrowMagickException(exception,GetMagickModule(),
2308         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2309       return(MagickFalse);
2310     }
2311   status=MagickTrue;
2312   scale=1.0/scale;
2313   (void) ResetMagickMemory(&zero,0,sizeof(zero));
2314   image_view=AcquireCacheView(image);
2315   resize_view=AcquireCacheView(resize_image);
2316 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2317   #pragma omp parallel for shared(status)
2318 #endif
2319   for (y=0; y < (ssize_t) resize_image->rows; y++)
2320   {
2321     MagickRealType
2322       center,
2323       density;
2324
2325     register const IndexPacket
2326       *restrict indexes;
2327
2328     register const PixelPacket
2329       *restrict p;
2330
2331     register ContributionInfo
2332       *restrict contribution;
2333
2334     register IndexPacket
2335       *restrict resize_indexes;
2336
2337     register PixelPacket
2338       *restrict q;
2339
2340     register ssize_t
2341       x;
2342
2343     ssize_t
2344       n,
2345       start,
2346       stop;
2347
2348     if (status == MagickFalse)
2349       continue;
2350     center=(MagickRealType) (y+0.5)/y_factor;
2351     start=(ssize_t) MagickMax(center-support+0.5,0.0);
2352     stop=(ssize_t) MagickMin(center+support+0.5,(double) image->rows);
2353     density=0.0;
2354     contribution=contributions[GetOpenMPThreadId()];
2355     for (n=0; n < (stop-start); n++)
2356     {
2357       contribution[n].pixel=start+n;
2358       contribution[n].weight=GetResizeFilterWeight(resize_filter,scale*
2359         ((MagickRealType) (start+n)-center+0.5));
2360       density+=contribution[n].weight;
2361     }
2362     if ((density != 0.0) && (density != 1.0))
2363       {
2364         register ssize_t
2365           i;
2366
2367         /*
2368           Normalize.
2369         */
2370         density=1.0/density;
2371         for (i=0; i < n; i++)
2372           contribution[i].weight*=density;
2373       }
2374     p=GetCacheViewVirtualPixels(image_view,0,contribution[0].pixel,
2375       image->columns,(size_t) (contribution[n-1].pixel-contribution[0].pixel+1),
2376       exception);
2377     q=QueueCacheViewAuthenticPixels(resize_view,0,y,resize_image->columns,1,
2378       exception);
2379     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2380       {
2381         status=MagickFalse;
2382         continue;
2383       }
2384     indexes=GetCacheViewVirtualIndexQueue(image_view);
2385     resize_indexes=GetCacheViewAuthenticIndexQueue(resize_view);
2386     for (x=0; x < (ssize_t) resize_image->columns; x++)
2387     {
2388       MagickPixelPacket
2389         pixel;
2390
2391       MagickRealType
2392         alpha;
2393
2394       register ssize_t
2395         i;
2396
2397       ssize_t
2398         j;
2399
2400       pixel=zero;
2401       if (image->matte == MagickFalse)
2402         {
2403           for (i=0; i < n; i++)
2404           {
2405             j=(ssize_t) ((contribution[i].pixel-contribution[0].pixel)*
2406               image->columns+x);
2407             alpha=contribution[i].weight;
2408             pixel.red+=alpha*(p+j)->red;
2409             pixel.green+=alpha*(p+j)->green;
2410             pixel.blue+=alpha*(p+j)->blue;
2411             pixel.opacity+=alpha*(p+j)->opacity;
2412           }
2413           SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
2414           SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
2415           SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
2416           SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
2417           if ((image->colorspace == CMYKColorspace) &&
2418               (resize_image->colorspace == CMYKColorspace))
2419             {
2420               for (i=0; i < n; i++)
2421               {
2422                 j=(ssize_t) ((contribution[i].pixel-contribution[0].pixel)*
2423                   image->columns+x);
2424                 alpha=contribution[i].weight;
2425                 pixel.index+=alpha*indexes[j];
2426               }
2427               resize_indexes[x]=(IndexPacket) ClampToQuantum(pixel.index);
2428             }
2429         }
2430       else
2431         {
2432           MagickRealType
2433             gamma;
2434
2435           gamma=0.0;
2436           for (i=0; i < n; i++)
2437           {
2438             j=(ssize_t) ((contribution[i].pixel-contribution[0].pixel)*
2439               image->columns+x);
2440             alpha=contribution[i].weight*QuantumScale*
2441               GetAlphaPixelComponent(p+j);
2442             pixel.red+=alpha*(p+j)->red;
2443             pixel.green+=alpha*(p+j)->green;
2444             pixel.blue+=alpha*(p+j)->blue;
2445             pixel.opacity+=contribution[i].weight*(p+j)->opacity;
2446             gamma+=alpha;
2447           }
2448           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2449           q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
2450           q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
2451           q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
2452           SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
2453           if ((image->colorspace == CMYKColorspace) &&
2454               (resize_image->colorspace == CMYKColorspace))
2455             {
2456               for (i=0; i < n; i++)
2457               {
2458                 j=(ssize_t) ((contribution[i].pixel-contribution[0].pixel)*
2459                   image->columns+x);
2460                 alpha=contribution[i].weight*QuantumScale*
2461                   GetAlphaPixelComponent(p+j);
2462                 pixel.index+=alpha*indexes[j];
2463               }
2464               resize_indexes[x]=(IndexPacket) ClampToQuantum(gamma*
2465                 GetIndexPixelComponent(&pixel));
2466             }
2467         }
2468       if ((resize_image->storage_class == PseudoClass) &&
2469           (image->storage_class == PseudoClass))
2470         {
2471           i=(ssize_t) (MagickMin(MagickMax(center,(double) start),(double) stop-
2472             1.0)+0.5);
2473           j=(ssize_t) ((contribution[i-start].pixel-contribution[0].pixel)*
2474             image->columns+x);
2475           resize_indexes[x]=indexes[j];
2476         }
2477       q++;
2478     }
2479     if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
2480       status=MagickFalse;
2481     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2482       {
2483         MagickBooleanType
2484           proceed;
2485
2486 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2487   #pragma omp critical (MagickCore_VerticalFilter)
2488 #endif
2489         proceed=SetImageProgress(image,ResizeImageTag,(*offset)++,span);
2490         if (proceed == MagickFalse)
2491           status=MagickFalse;
2492       }
2493   }
2494   resize_view=DestroyCacheView(resize_view);
2495   image_view=DestroyCacheView(image_view);
2496   contributions=DestroyContributionThreadSet(contributions);
2497   return(status);
2498 }
2499
2500 MagickExport Image *ResizeImage(const Image *image,const size_t columns,
2501   const size_t rows,const FilterTypes filter,const double blur,
2502   ExceptionInfo *exception)
2503 {
2504 #define WorkLoadFactor  0.265
2505
2506   FilterTypes
2507     filter_type;
2508
2509   Image
2510     *filter_image,
2511     *resize_image;
2512
2513   MagickOffsetType
2514     offset;
2515
2516   MagickRealType
2517     x_factor,
2518     y_factor;
2519
2520   MagickSizeType
2521     span;
2522
2523   MagickStatusType
2524     status;
2525
2526   ResizeFilter
2527     *resize_filter;
2528
2529   /*
2530     Acquire resize image.
2531   */
2532   assert(image != (Image *) NULL);
2533   assert(image->signature == MagickSignature);
2534   if (image->debug != MagickFalse)
2535     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2536   assert(exception != (ExceptionInfo *) NULL);
2537   assert(exception->signature == MagickSignature);
2538   if ((columns == 0) || (rows == 0))
2539     ThrowImageException(ImageError,"NegativeOrZeroImageSize");
2540   if ((columns == image->columns) && (rows == image->rows) &&
2541       (filter == UndefinedFilter) && (blur == 1.0))
2542     return(CloneImage(image,0,0,MagickTrue,exception));
2543   resize_image=CloneImage(image,columns,rows,MagickTrue,exception);
2544   if (resize_image == (Image *) NULL)
2545     return(resize_image);
2546   /*
2547     Acquire resize filter.
2548   */
2549   x_factor=(MagickRealType) columns/(MagickRealType) image->columns;
2550   y_factor=(MagickRealType) rows/(MagickRealType) image->rows;
2551   if ((x_factor*y_factor) > WorkLoadFactor)
2552     filter_image=CloneImage(image,columns,image->rows,MagickTrue,exception);
2553   else
2554     filter_image=CloneImage(image,image->columns,rows,MagickTrue,exception);
2555   if (filter_image == (Image *) NULL)
2556     return(DestroyImage(resize_image));
2557   filter_type=LanczosFilter;
2558   if (filter != UndefinedFilter)
2559     filter_type=filter;
2560   else
2561     if ((x_factor == 1.0) && (y_factor == 1.0))
2562       filter_type=PointFilter;
2563     else
2564       if ((image->storage_class == PseudoClass) ||
2565           (image->matte != MagickFalse) || ((x_factor*y_factor) > 1.0))
2566         filter_type=MitchellFilter;
2567   resize_filter=AcquireResizeFilter(image,filter_type,blur,MagickFalse,
2568     exception);
2569   /*
2570     Resize image.
2571   */
2572   offset=0;
2573   if ((x_factor*y_factor) > WorkLoadFactor)
2574     {
2575       span=(MagickSizeType) (filter_image->columns+rows);
2576       status=HorizontalFilter(resize_filter,image,filter_image,x_factor,span,
2577         &offset,exception);
2578       status&=VerticalFilter(resize_filter,filter_image,resize_image,y_factor,
2579         span,&offset,exception);
2580     }
2581   else
2582     {
2583       span=(MagickSizeType) (filter_image->rows+columns);
2584       status=VerticalFilter(resize_filter,image,filter_image,y_factor,span,
2585         &offset,exception);
2586       status&=HorizontalFilter(resize_filter,filter_image,resize_image,x_factor,
2587         span,&offset,exception);
2588     }
2589   /*
2590     Free resources.
2591   */
2592   filter_image=DestroyImage(filter_image);
2593   resize_filter=DestroyResizeFilter(resize_filter);
2594   if ((status == MagickFalse) || (resize_image == (Image *) NULL))
2595     return((Image *) NULL);
2596   resize_image->type=image->type;
2597   return(resize_image);
2598 }
2599 \f
2600 /*
2601 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2602 %                                                                             %
2603 %                                                                             %
2604 %                                                                             %
2605 %   S a m p l e I m a g e                                                     %
2606 %                                                                             %
2607 %                                                                             %
2608 %                                                                             %
2609 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2610 %
2611 %  SampleImage() scales an image to the desired dimensions with pixel
2612 %  sampling.  Unlike other scaling methods, this method does not introduce
2613 %  any additional color into the scaled image.
2614 %
2615 %  The format of the SampleImage method is:
2616 %
2617 %      Image *SampleImage(const Image *image,const size_t columns,
2618 %        const size_t rows,ExceptionInfo *exception)
2619 %
2620 %  A description of each parameter follows:
2621 %
2622 %    o image: the image.
2623 %
2624 %    o columns: the number of columns in the sampled image.
2625 %
2626 %    o rows: the number of rows in the sampled image.
2627 %
2628 %    o exception: return any errors or warnings in this structure.
2629 %
2630 */
2631 MagickExport Image *SampleImage(const Image *image,const size_t columns,
2632   const size_t rows,ExceptionInfo *exception)
2633 {
2634 #define SampleImageTag  "Sample/Image"
2635
2636   CacheView
2637     *image_view,
2638     *sample_view;
2639
2640   Image
2641     *sample_image;
2642
2643   MagickBooleanType
2644     status;
2645
2646   MagickOffsetType
2647     progress;
2648
2649   register ssize_t
2650     x;
2651
2652   ssize_t
2653     *x_offset,
2654     y;
2655
2656   /*
2657     Initialize sampled image attributes.
2658   */
2659   assert(image != (const Image *) NULL);
2660   assert(image->signature == MagickSignature);
2661   if (image->debug != MagickFalse)
2662     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2663   assert(exception != (ExceptionInfo *) NULL);
2664   assert(exception->signature == MagickSignature);
2665   if ((columns == 0) || (rows == 0))
2666     ThrowImageException(ImageError,"NegativeOrZeroImageSize");
2667   if ((columns == image->columns) && (rows == image->rows))
2668     return(CloneImage(image,0,0,MagickTrue,exception));
2669   sample_image=CloneImage(image,columns,rows,MagickTrue,exception);
2670   if (sample_image == (Image *) NULL)
2671     return((Image *) NULL);
2672   /*
2673     Allocate scan line buffer and column offset buffers.
2674   */
2675   x_offset=(ssize_t *) AcquireQuantumMemory((size_t) sample_image->columns,
2676     sizeof(*x_offset));
2677   if (x_offset == (ssize_t *) NULL)
2678     {
2679       sample_image=DestroyImage(sample_image);
2680       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2681     }
2682   for (x=0; x < (ssize_t) sample_image->columns; x++)
2683     x_offset[x]=(ssize_t) (((MagickRealType) x+0.5)*image->columns/
2684       sample_image->columns);
2685   /*
2686     Sample each row.
2687   */
2688   status=MagickTrue;
2689   progress=0;
2690   image_view=AcquireCacheView(image);
2691   sample_view=AcquireCacheView(sample_image);
2692 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2693   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2694 #endif
2695   for (y=0; y < (ssize_t) sample_image->rows; y++)
2696   {
2697     register const IndexPacket
2698       *restrict indexes;
2699
2700     register const PixelPacket
2701       *restrict p;
2702
2703     register IndexPacket
2704       *restrict sample_indexes;
2705
2706     register PixelPacket
2707       *restrict q;
2708
2709     register ssize_t
2710       x;
2711
2712     ssize_t
2713       y_offset;
2714
2715     if (status == MagickFalse)
2716       continue;
2717     y_offset=(ssize_t) (((MagickRealType) y+0.5)*image->rows/
2718       sample_image->rows);
2719     p=GetCacheViewVirtualPixels(image_view,0,y_offset,image->columns,1,
2720       exception);
2721     q=QueueCacheViewAuthenticPixels(sample_view,0,y,sample_image->columns,1,
2722       exception);
2723     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2724       {
2725         status=MagickFalse;
2726         continue;
2727       }
2728     indexes=GetCacheViewAuthenticIndexQueue(image_view);
2729     sample_indexes=GetCacheViewAuthenticIndexQueue(sample_view);
2730     /*
2731       Sample each column.
2732     */
2733     for (x=0; x < (ssize_t) sample_image->columns; x++)
2734       *q++=p[x_offset[x]];
2735     if ((image->storage_class == PseudoClass) ||
2736         (image->colorspace == CMYKColorspace))
2737       for (x=0; x < (ssize_t) sample_image->columns; x++)
2738         sample_indexes[x]=indexes[x_offset[x]];
2739     if (SyncCacheViewAuthenticPixels(sample_view,exception) == MagickFalse)
2740       status=MagickFalse;
2741     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2742       {
2743         MagickBooleanType
2744           proceed;
2745
2746 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2747         #pragma omp critical (MagickCore_SampleImage)
2748 #endif
2749         proceed=SetImageProgress(image,SampleImageTag,progress++,image->rows);
2750         if (proceed == MagickFalse)
2751           status=MagickFalse;
2752       }
2753   }
2754   image_view=DestroyCacheView(image_view);
2755   sample_view=DestroyCacheView(sample_view);
2756   x_offset=(ssize_t *) RelinquishMagickMemory(x_offset);
2757   sample_image->type=image->type;
2758   return(sample_image);
2759 }
2760 \f
2761 /*
2762 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2763 %                                                                             %
2764 %                                                                             %
2765 %                                                                             %
2766 %   S c a l e I m a g e                                                       %
2767 %                                                                             %
2768 %                                                                             %
2769 %                                                                             %
2770 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2771 %
2772 %  ScaleImage() changes the size of an image to the given dimensions.
2773 %
2774 %  The format of the ScaleImage method is:
2775 %
2776 %      Image *ScaleImage(const Image *image,const size_t columns,
2777 %        const size_t rows,ExceptionInfo *exception)
2778 %
2779 %  A description of each parameter follows:
2780 %
2781 %    o image: the image.
2782 %
2783 %    o columns: the number of columns in the scaled image.
2784 %
2785 %    o rows: the number of rows in the scaled image.
2786 %
2787 %    o exception: return any errors or warnings in this structure.
2788 %
2789 */
2790 MagickExport Image *ScaleImage(const Image *image,const size_t columns,
2791   const size_t rows,ExceptionInfo *exception)
2792 {
2793 #define ScaleImageTag  "Scale/Image"
2794
2795   CacheView
2796     *image_view,
2797     *scale_view;
2798
2799   Image
2800     *scale_image;
2801
2802   MagickBooleanType
2803     next_column,
2804     next_row,
2805     proceed;
2806
2807   MagickPixelPacket
2808     pixel,
2809     *scale_scanline,
2810     *scanline,
2811     *x_vector,
2812     *y_vector,
2813     zero;
2814
2815   PointInfo
2816     scale,
2817     span;
2818
2819   register ssize_t
2820     i;
2821
2822   ssize_t
2823     number_rows,
2824     y;
2825
2826   /*
2827     Initialize scaled image attributes.
2828   */
2829   assert(image != (const Image *) NULL);
2830   assert(image->signature == MagickSignature);
2831   if (image->debug != MagickFalse)
2832     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2833   assert(exception != (ExceptionInfo *) NULL);
2834   assert(exception->signature == MagickSignature);
2835   if ((columns == 0) || (rows == 0))
2836     return((Image *) NULL);
2837   if ((columns == image->columns) && (rows == image->rows))
2838     return(CloneImage(image,0,0,MagickTrue,exception));
2839   scale_image=CloneImage(image,columns,rows,MagickTrue,exception);
2840   if (scale_image == (Image *) NULL)
2841     return((Image *) NULL);
2842   if (SetImageStorageClass(scale_image,DirectClass) == MagickFalse)
2843     {
2844       InheritException(exception,&scale_image->exception);
2845       scale_image=DestroyImage(scale_image);
2846       return((Image *) NULL);
2847     }
2848   /*
2849     Allocate memory.
2850   */
2851   x_vector=(MagickPixelPacket *) AcquireQuantumMemory((size_t) image->columns,
2852     sizeof(*x_vector));
2853   scanline=x_vector;
2854   if (image->rows != scale_image->rows)
2855     scanline=(MagickPixelPacket *) AcquireQuantumMemory((size_t) image->columns,
2856       sizeof(*scanline));
2857   scale_scanline=(MagickPixelPacket *) AcquireQuantumMemory((size_t)
2858     scale_image->columns,sizeof(*scale_scanline));
2859   y_vector=(MagickPixelPacket *) AcquireQuantumMemory((size_t) image->columns,
2860     sizeof(*y_vector));
2861   if ((scanline == (MagickPixelPacket *) NULL) ||
2862       (scale_scanline == (MagickPixelPacket *) NULL) ||
2863       (x_vector == (MagickPixelPacket *) NULL) ||
2864       (y_vector == (MagickPixelPacket *) NULL))
2865     {
2866       scale_image=DestroyImage(scale_image);
2867       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2868     }
2869   /*
2870     Scale image.
2871   */
2872   number_rows=0;
2873   next_row=MagickTrue;
2874   span.y=1.0;
2875   scale.y=(double) scale_image->rows/(double) image->rows;
2876   (void) ResetMagickMemory(y_vector,0,(size_t) image->columns*
2877     sizeof(*y_vector));
2878   GetMagickPixelPacket(image,&pixel);
2879   (void) ResetMagickMemory(&zero,0,sizeof(zero));
2880   i=0;
2881   image_view=AcquireCacheView(image);
2882   scale_view=AcquireCacheView(scale_image);
2883   for (y=0; y < (ssize_t) scale_image->rows; y++)
2884   {
2885     register const IndexPacket
2886       *restrict indexes;
2887
2888     register const PixelPacket
2889       *restrict p;
2890
2891     register IndexPacket
2892       *restrict scale_indexes;
2893
2894     register MagickPixelPacket
2895       *restrict s,
2896       *restrict t;
2897
2898     register PixelPacket
2899       *restrict q;
2900
2901     register ssize_t
2902       x;
2903
2904     q=QueueCacheViewAuthenticPixels(scale_view,0,y,scale_image->columns,1,
2905       exception);
2906     if (q == (PixelPacket *) NULL)
2907       break;
2908     scale_indexes=GetAuthenticIndexQueue(scale_image);
2909     if (scale_image->rows == image->rows)
2910       {
2911         /*
2912           Read a new scanline.
2913         */
2914         p=GetCacheViewVirtualPixels(image_view,0,i++,image->columns,1,
2915           exception);
2916         if (p == (const PixelPacket *) NULL)
2917           break;
2918         indexes=GetCacheViewVirtualIndexQueue(image_view);
2919         for (x=0; x < (ssize_t) image->columns; x++)
2920         {
2921           x_vector[x].red=(MagickRealType) GetRedPixelComponent(p);
2922           x_vector[x].green=(MagickRealType) GetGreenPixelComponent(p);
2923           x_vector[x].blue=(MagickRealType) GetBluePixelComponent(p);
2924           if (image->matte != MagickFalse)
2925             x_vector[x].opacity=(MagickRealType) GetOpacityPixelComponent(p);
2926           if (indexes != (IndexPacket *) NULL)
2927             x_vector[x].index=(MagickRealType) indexes[x];
2928           p++;
2929         }
2930       }
2931     else
2932       {
2933         /*
2934           Scale Y direction.
2935         */
2936         while (scale.y < span.y)
2937         {
2938           if ((next_row != MagickFalse) &&
2939               (number_rows < (ssize_t) image->rows))
2940             {
2941               /*
2942                 Read a new scanline.
2943               */
2944               p=GetCacheViewVirtualPixels(image_view,0,i++,image->columns,1,
2945                 exception);
2946               if (p == (const PixelPacket *) NULL)
2947                 break;
2948               indexes=GetCacheViewVirtualIndexQueue(image_view);
2949               for (x=0; x < (ssize_t) image->columns; x++)
2950               {
2951                 x_vector[x].red=(MagickRealType) GetRedPixelComponent(p);
2952                 x_vector[x].green=(MagickRealType) GetGreenPixelComponent(p);
2953                 x_vector[x].blue=(MagickRealType) GetBluePixelComponent(p);
2954                 if (image->matte != MagickFalse)
2955                   x_vector[x].opacity=(MagickRealType)
2956                     GetOpacityPixelComponent(p);
2957                 if (indexes != (IndexPacket *) NULL)
2958                   x_vector[x].index=(MagickRealType) indexes[x];
2959                 p++;
2960               }
2961               number_rows++;
2962             }
2963           for (x=0; x < (ssize_t) image->columns; x++)
2964           {
2965             y_vector[x].red+=scale.y*x_vector[x].red;
2966             y_vector[x].green+=scale.y*x_vector[x].green;
2967             y_vector[x].blue+=scale.y*x_vector[x].blue;
2968             if (scale_image->matte != MagickFalse)
2969               y_vector[x].opacity+=scale.y*x_vector[x].opacity;
2970             if (scale_indexes != (IndexPacket *) NULL)
2971               y_vector[x].index+=scale.y*x_vector[x].index;
2972           }
2973           span.y-=scale.y;
2974           scale.y=(double) scale_image->rows/(double) image->rows;
2975           next_row=MagickTrue;
2976         }
2977         if ((next_row != MagickFalse) && (number_rows < (ssize_t) image->rows))
2978           {
2979             /*
2980               Read a new scanline.
2981             */
2982             p=GetCacheViewVirtualPixels(image_view,0,i++,image->columns,1,
2983               exception);
2984             if (p == (const PixelPacket *) NULL)
2985               break;
2986             indexes=GetCacheViewVirtualIndexQueue(image_view);
2987             for (x=0; x < (ssize_t) image->columns; x++)
2988             {
2989               x_vector[x].red=(MagickRealType) GetRedPixelComponent(p);
2990               x_vector[x].green=(MagickRealType) GetGreenPixelComponent(p);
2991               x_vector[x].blue=(MagickRealType) GetBluePixelComponent(p);
2992               if (image->matte != MagickFalse)
2993                 x_vector[x].opacity=(MagickRealType)
2994                   GetOpacityPixelComponent(p);
2995               if (indexes != (IndexPacket *) NULL)
2996                 x_vector[x].index=(MagickRealType) indexes[x];
2997               p++;
2998             }
2999             number_rows++;
3000             next_row=MagickFalse;
3001           }
3002         s=scanline;
3003         for (x=0; x < (ssize_t) image->columns; x++)
3004         {
3005           pixel.red=y_vector[x].red+span.y*x_vector[x].red;
3006           pixel.green=y_vector[x].green+span.y*x_vector[x].green;
3007           pixel.blue=y_vector[x].blue+span.y*x_vector[x].blue;
3008           if (image->matte != MagickFalse)
3009             pixel.opacity=y_vector[x].opacity+span.y*x_vector[x].opacity;
3010           if (scale_indexes != (IndexPacket *) NULL)
3011             pixel.index=y_vector[x].index+span.y*x_vector[x].index;
3012           s->red=pixel.red;
3013           s->green=pixel.green;
3014           s->blue=pixel.blue;
3015           if (scale_image->matte != MagickFalse)
3016             s->opacity=pixel.opacity;
3017           if (scale_indexes != (IndexPacket *) NULL)
3018             s->index=pixel.index;
3019           s++;
3020           y_vector[x]=zero;
3021         }
3022         scale.y-=span.y;
3023         if (scale.y <= 0)
3024           {
3025             scale.y=(double) scale_image->rows/(double) image->rows;
3026             next_row=MagickTrue;
3027           }
3028         span.y=1.0;
3029       }
3030     if (scale_image->columns == image->columns)
3031       {
3032         /*
3033           Transfer scanline to scaled image.
3034         */
3035         s=scanline;
3036         for (x=0; x < (ssize_t) scale_image->columns; x++)
3037         {
3038           q->red=ClampToQuantum(s->red);
3039           q->green=ClampToQuantum(s->green);
3040           q->blue=ClampToQuantum(s->blue);
3041           if (scale_image->matte != MagickFalse)
3042             q->opacity=ClampToQuantum(s->opacity);
3043           if (scale_indexes != (IndexPacket *) NULL)
3044             scale_indexes[x]=(IndexPacket) ClampToQuantum(s->index);
3045           q++;
3046           s++;
3047         }
3048       }
3049     else
3050       {
3051         /*
3052           Scale X direction.
3053         */
3054         pixel=zero;
3055         next_column=MagickFalse;
3056         span.x=1.0;
3057         s=scanline;
3058         t=scale_scanline;
3059         for (x=0; x < (ssize_t) image->columns; x++)
3060         {
3061           scale.x=(double) scale_image->columns/(double) image->columns;
3062           while (scale.x >= span.x)
3063           {
3064             if (next_column != MagickFalse)
3065               {
3066                 pixel=zero;
3067                 t++;
3068               }
3069             pixel.red+=span.x*s->red;
3070             pixel.green+=span.x*s->green;
3071             pixel.blue+=span.x*s->blue;
3072             if (image->matte != MagickFalse)
3073               pixel.opacity+=span.x*s->opacity;
3074             if (scale_indexes != (IndexPacket *) NULL)
3075               pixel.index+=span.x*s->index;
3076             t->red=pixel.red;
3077             t->green=pixel.green;
3078             t->blue=pixel.blue;
3079             if (scale_image->matte != MagickFalse)
3080               t->opacity=pixel.opacity;
3081             if (scale_indexes != (IndexPacket *) NULL)
3082               t->index=pixel.index;
3083             scale.x-=span.x;
3084             span.x=1.0;
3085             next_column=MagickTrue;
3086           }
3087         if (scale.x > 0)
3088           {
3089             if (next_column != MagickFalse)
3090               {
3091                 pixel=zero;
3092                 next_column=MagickFalse;
3093                 t++;
3094               }
3095             pixel.red+=scale.x*s->red;
3096             pixel.green+=scale.x*s->green;
3097             pixel.blue+=scale.x*s->blue;
3098             if (scale_image->matte != MagickFalse)
3099               pixel.opacity+=scale.x*s->opacity;
3100             if (scale_indexes != (IndexPacket *) NULL)
3101               pixel.index+=scale.x*s->index;
3102             span.x-=scale.x;
3103           }
3104         s++;
3105       }
3106       if (span.x > 0)
3107         {
3108           s--;
3109           pixel.red+=span.x*s->red;
3110           pixel.green+=span.x*s->green;
3111           pixel.blue+=span.x*s->blue;
3112           if (scale_image->matte != MagickFalse)
3113             pixel.opacity+=span.x*s->opacity;
3114           if (scale_indexes != (IndexPacket *) NULL)
3115             pixel.index+=span.x*s->index;
3116         }
3117       if ((next_column == MagickFalse) &&
3118           ((ssize_t) (t-scale_scanline) < (ssize_t) scale_image->columns))
3119         {
3120           t->red=pixel.red;
3121           t->green=pixel.green;
3122           t->blue=pixel.blue;
3123           if (scale_image->matte != MagickFalse)
3124             t->opacity=pixel.opacity;
3125           if (scale_indexes != (IndexPacket *) NULL)
3126             t->index=pixel.index;
3127         }
3128       /*
3129         Transfer scanline to scaled image.
3130       */
3131       t=scale_scanline;
3132       for (x=0; x < (ssize_t) scale_image->columns; x++)
3133       {
3134         q->red=ClampToQuantum(t->red);
3135         q->green=ClampToQuantum(t->green);
3136         q->blue=ClampToQuantum(t->blue);
3137         if (scale_image->matte != MagickFalse)
3138           q->opacity=ClampToQuantum(t->opacity);
3139         if (scale_indexes != (IndexPacket *) NULL)
3140           scale_indexes[x]=(IndexPacket) ClampToQuantum(t->index);
3141         t++;
3142         q++;
3143       }
3144     }
3145     if (SyncCacheViewAuthenticPixels(scale_view,exception) == MagickFalse)
3146       break;
3147     proceed=SetImageProgress(image,ScaleImageTag,(MagickOffsetType) y,
3148       image->rows);
3149     if (proceed == MagickFalse)
3150       break;
3151   }
3152   scale_view=DestroyCacheView(scale_view);
3153   image_view=DestroyCacheView(image_view);
3154   /*
3155     Free allocated memory.
3156   */
3157   y_vector=(MagickPixelPacket *) RelinquishMagickMemory(y_vector);
3158   scale_scanline=(MagickPixelPacket *) RelinquishMagickMemory(scale_scanline);
3159   if (scale_image->rows != image->rows)
3160     scanline=(MagickPixelPacket *) RelinquishMagickMemory(scanline);
3161   x_vector=(MagickPixelPacket *) RelinquishMagickMemory(x_vector);
3162   scale_image->type=image->type;
3163   return(scale_image);
3164 }
3165 \f
3166 #if 0
3167       THIS IS NOT USED  --  to be removed
3168 /*
3169 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3170 %                                                                             %
3171 %                                                                             %
3172 %                                                                             %
3173 +   S e t R e s i z e F i l t e r S u p p o r t                               %
3174 %                                                                             %
3175 %                                                                             %
3176 %                                                                             %
3177 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3178 %
3179 %  SetResizeFilterSupport() specifies which IR filter to use to window
3180 %
3181 %  The format of the SetResizeFilterSupport method is:
3182 %
3183 %      void SetResizeFilterSupport(ResizeFilter *resize_filter,
3184 %        const MagickRealType support)
3185 %
3186 %  A description of each parameter follows:
3187 %
3188 %    o resize_filter: the resize filter.
3189 %
3190 %    o support: the filter spport radius.
3191 %
3192 */
3193 MagickExport void SetResizeFilterSupport(ResizeFilter *resize_filter,
3194   const MagickRealType support)
3195 {
3196   assert(resize_filter != (ResizeFilter *) NULL);
3197   assert(resize_filter->signature == MagickSignature);
3198   resize_filter->support=support;
3199 }
3200 #endif
3201 \f
3202 /*
3203 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3204 %                                                                             %
3205 %                                                                             %
3206 %                                                                             %
3207 %   T h u m b n a i l I m a g e                                               %
3208 %                                                                             %
3209 %                                                                             %
3210 %                                                                             %
3211 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3212 %
3213 %  ThumbnailImage() changes the size of an image to the given dimensions and
3214 %  removes any associated profiles.  The goal is to produce small low cost
3215 %  thumbnail images suited for display on the Web.
3216 %
3217 %  The format of the ThumbnailImage method is:
3218 %
3219 %      Image *ThumbnailImage(const Image *image,const size_t columns,
3220 %        const size_t rows,ExceptionInfo *exception)
3221 %
3222 %  A description of each parameter follows:
3223 %
3224 %    o image: the image.
3225 %
3226 %    o columns: the number of columns in the scaled image.
3227 %
3228 %    o rows: the number of rows in the scaled image.
3229 %
3230 %    o exception: return any errors or warnings in this structure.
3231 %
3232 */
3233 MagickExport Image *ThumbnailImage(const Image *image,const size_t columns,
3234   const size_t rows,ExceptionInfo *exception)
3235 {
3236 #define SampleFactor  5
3237
3238   char
3239     value[MaxTextExtent];
3240
3241   const char
3242     *name;
3243
3244   Image
3245     *thumbnail_image;
3246
3247   MagickRealType
3248     x_factor,
3249     y_factor;
3250
3251   size_t
3252     version;
3253
3254   struct stat
3255     attributes;
3256
3257   assert(image != (Image *) NULL);
3258   assert(image->signature == MagickSignature);
3259   if (image->debug != MagickFalse)
3260     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3261   assert(exception != (ExceptionInfo *) NULL);
3262   assert(exception->signature == MagickSignature);
3263   x_factor=(MagickRealType) columns/(MagickRealType) image->columns;
3264   y_factor=(MagickRealType) rows/(MagickRealType) image->rows;
3265   if ((x_factor*y_factor) > 0.1)
3266     thumbnail_image=ResizeImage(image,columns,rows,image->filter,image->blur,
3267       exception);
3268   else
3269     if (((SampleFactor*columns) < 128) || ((SampleFactor*rows) < 128))
3270       thumbnail_image=ResizeImage(image,columns,rows,image->filter,
3271         image->blur,exception);
3272     else
3273       {
3274         Image
3275           *sample_image;
3276
3277         sample_image=SampleImage(image,SampleFactor*columns,SampleFactor*rows,
3278           exception);
3279         if (sample_image == (Image *) NULL)
3280           return((Image *) NULL);
3281         thumbnail_image=ResizeImage(sample_image,columns,rows,image->filter,
3282           image->blur,exception);
3283         sample_image=DestroyImage(sample_image);
3284       }
3285   if (thumbnail_image == (Image *) NULL)
3286     return(thumbnail_image);
3287   (void) ParseAbsoluteGeometry("0x0+0+0",&thumbnail_image->page);
3288   if (thumbnail_image->matte == MagickFalse)
3289     (void) SetImageAlphaChannel(thumbnail_image,OpaqueAlphaChannel);
3290   thumbnail_image->depth=8;
3291   thumbnail_image->interlace=NoInterlace;
3292   /*
3293     Strip all profiles except color profiles.
3294   */
3295   ResetImageProfileIterator(thumbnail_image);
3296   for (name=GetNextImageProfile(thumbnail_image); name != (const char *) NULL; )
3297   {
3298     if ((LocaleCompare(name,"icc") != 0) && (LocaleCompare(name,"icm") != 0))
3299      {
3300        (void) DeleteImageProfile(thumbnail_image,name);
3301        ResetImageProfileIterator(thumbnail_image);
3302      }
3303     name=GetNextImageProfile(thumbnail_image);
3304   }
3305   (void) DeleteImageProperty(thumbnail_image,"comment");
3306   (void) CopyMagickString(value,image->magick_filename,MaxTextExtent);
3307   if (strstr(image->magick_filename,"//") == (char *) NULL)
3308     (void) FormatMagickString(value,MaxTextExtent,"file://%s",
3309       image->magick_filename);
3310   (void) SetImageProperty(thumbnail_image,"Thumb::URI",value);
3311   (void) CopyMagickString(value,image->magick_filename,MaxTextExtent);
3312   if (GetPathAttributes(image->filename,&attributes) != MagickFalse)
3313     {
3314       (void) FormatMagickString(value,MaxTextExtent,"%.20g",(double)
3315         attributes.st_mtime);
3316       (void) SetImageProperty(thumbnail_image,"Thumb::MTime",value);
3317     }
3318   (void) FormatMagickString(value,MaxTextExtent,"%.20g",(double)
3319     attributes.st_mtime);
3320   (void) FormatMagickSize(GetBlobSize(image),MagickFalse,value);
3321   (void) ConcatenateMagickString(value,"B",MaxTextExtent);
3322   (void) SetImageProperty(thumbnail_image,"Thumb::Size",value);
3323   (void) FormatMagickString(value,MaxTextExtent,"image/%s",image->magick);
3324   LocaleLower(value);
3325   (void) SetImageProperty(thumbnail_image,"Thumb::Mimetype",value);
3326   (void) SetImageProperty(thumbnail_image,"software",
3327     GetMagickVersion(&version));
3328   (void) FormatMagickString(value,MaxTextExtent,"%.20g",(double)
3329     image->magick_columns);
3330   (void) SetImageProperty(thumbnail_image,"Thumb::Image::Width",value);
3331   (void) FormatMagickString(value,MaxTextExtent,"%.20g",(double)
3332     image->magick_rows);
3333   (void) SetImageProperty(thumbnail_image,"Thumb::Image::height",value);
3334   (void) FormatMagickString(value,MaxTextExtent,"%.20g",(double)
3335     GetImageListLength(image));
3336   (void) SetImageProperty(thumbnail_image,"Thumb::Document::Pages",value);
3337   return(thumbnail_image);
3338 }