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