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