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