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

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