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

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