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=StringToDouble(artifact,(char **) NULL);  /* override sigma */
 936   if (GaussianFilter)
 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*=StringToDouble(artifact,(char **) NULL);  /* override blur */
 948   if (resize_filter->blur < MagickEpsilon)
 949     resize_filter->blur=(MagickRealType) MagickEpsilon;
 950   artifact=GetImageArtifact(image,"filter:lobes");
 951   if (artifact != (const char *) NULL)
 952     {
 953       ssize_t
 954         lobes;
 955
 956       /*
 957         Override lobes.
 958       */
 959       lobes=(ssize_t) StringToLong(artifact);
 960       if (lobes < 1)
 961         lobes=1;
 962       resize_filter->support=(MagickRealType) lobes;
 963     }
 964   if (resize_filter->filter == Jinc)
 965     {
 966       /*
 967         Convert a Jinc function lobes value to a real support value.
 968       */
 969       if (resize_filter->support > 16)
 970         resize_filter->support=jinc_zeros[15];  /* largest entry in table */
 971       else
 972         resize_filter->support=jinc_zeros[((long)resize_filter->support)-1];
 973     }
 974   artifact=GetImageArtifact(image,"filter:support");
 975   if (artifact != (const char *) NULL)
 976     resize_filter->support=fabs(StringToDouble(artifact,(char **) NULL)); /* override support */
 977   /*
 978     Scale windowing function separately to the support 'clipping' window that
 979     calling operator is planning to actually use (expert override).
 980   */
 981   resize_filter->window_support=resize_filter->support; /* default */
 982   artifact=GetImageArtifact(image,"filter:win-support");
 983   if (artifact != (const char *) NULL)
 984     resize_filter->window_support=fabs(StringToDouble(artifact,(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=StringToDouble(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=StringToDouble(artifact,(char **) NULL);
1013         }
1014       else
1015         {
1016           artifact=GetImageArtifact(image,"filter:c");
1017           if (artifact != (const char *) NULL)
1018             {
1019               C=StringToDouble(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,const PixelInterpolateMethod method,
1144 %        ExceptionInfo *exception)
1145 %
1146 %  A description of each parameter follows:
1147 %
1148 %    o image: the image.
1149 %
1150 %    o columns: the number of columns in the resized image.
1151 %
1152 %    o rows: the number of rows in the resized image.
1153 %
1154 %    o method: the pixel interpolation method.
1155 %
1156 %    o exception: return any errors or warnings in this structure.
1157 %
1158 */
1159 MagickExport Image *AdaptiveResizeImage(const Image *image,
1160   const size_t columns,const size_t rows,const PixelInterpolateMethod method,
1161   ExceptionInfo *exception)
1162 {
1163 #define AdaptiveResizeImageTag  "Resize/Image"
1164
1165   CacheView
1166     *image_view,
1167     *resize_view;
1168
1169   Image
1170     *resize_image;
1171
1172   MagickBooleanType
1173     status;
1174
1175   MagickOffsetType
1176     progress;
1177
1178   ssize_t
1179     y;
1180
1181   /*
1182     Adaptively resize image.
1183   */
1184   assert(image != (const Image *) NULL);
1185   assert(image->signature == MagickSignature);
1186   if (image->debug != MagickFalse)
1187     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1188   assert(exception != (ExceptionInfo *) NULL);
1189   assert(exception->signature == MagickSignature);
1190   if ((columns == 0) || (rows == 0))
1191     return((Image *) NULL);
1192   if ((columns == image->columns) && (rows == image->rows))
1193     return(CloneImage(image,0,0,MagickTrue,exception));
1194   resize_image=CloneImage(image,columns,rows,MagickTrue,exception);
1195   if (resize_image == (Image *) NULL)
1196     return((Image *) NULL);
1197   if (SetImageStorageClass(resize_image,DirectClass,exception) == MagickFalse)
1198     {
1199       resize_image=DestroyImage(resize_image);
1200       return((Image *) NULL);
1201     }
1202   status=MagickTrue;
1203   progress=0;
1204   image_view=AcquireCacheView(image);
1205   resize_view=AcquireCacheView(resize_image);
1206 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1207   #pragma omp parallel for schedule(dynamic,4) shared(progress,status) omp_throttle(1)
1208 #endif
1209   for (y=0; y < (ssize_t) resize_image->rows; y++)
1210   {
1211     PointInfo
1212       offset;
1213
1214     register Quantum
1215       *restrict q;
1216
1217     register ssize_t
1218       x;
1219
1220     if (status == MagickFalse)
1221       continue;
1222     q=QueueCacheViewAuthenticPixels(resize_view,0,y,resize_image->columns,1,
1223       exception);
1224     if (q == (Quantum *) NULL)
1225       continue;
1226     offset.y=((MagickRealType) (y+0.5)*image->rows/resize_image->rows);
1227     for (x=0; x < (ssize_t) resize_image->columns; x++)
1228     {
1229       offset.x=((MagickRealType) (x+0.5)*image->columns/resize_image->columns);
1230       status=InterpolatePixelChannels(image,image_view,resize_image,
1231         method,offset.x-0.5,offset.y-0.5,q,exception);
1232       q+=GetPixelChannels(resize_image);
1233     }
1234     if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
1235       continue;
1236     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1237       {
1238         MagickBooleanType
1239           proceed;
1240
1241 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1242   #pragma omp critical (MagickCore_AdaptiveResizeImage)
1243 #endif
1244         proceed=SetImageProgress(image,AdaptiveResizeImageTag,progress++,
1245           image->rows);
1246         if (proceed == MagickFalse)
1247           status=MagickFalse;
1248       }
1249   }
1250   resize_view=DestroyCacheView(resize_view);
1251   image_view=DestroyCacheView(image_view);
1252   if (status == MagickFalse)
1253     resize_image=DestroyImage(resize_image);
1254   return(resize_image);
1255 }
1256 \f
1257 /*
1258 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1259 %                                                                             %
1260 %                                                                             %
1261 %                                                                             %
1262 +   B e s s e l O r d e r O n e                                               %
1263 %                                                                             %
1264 %                                                                             %
1265 %                                                                             %
1266 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1267 %
1268 %  BesselOrderOne() computes the Bessel function of x of the first kind of
1269 %  order 0.  This is used to create the Jinc() filter function below.
1270 %
1271 %    Reduce x to |x| since j1(x)= -j1(-x), and for x in (0,8]
1272 %
1273 %       j1(x) = x*j1(x);
1274 %
1275 %    For x in (8,inf)
1276 %
1277 %       j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
1278 %
1279 %    where x1 = x-3*pi/4. Compute sin(x1) and cos(x1) as follow:
1280 %
1281 %       cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
1282 %               =  1/sqrt(2) * (sin(x) - cos(x))
1283 %       sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
1284 %               = -1/sqrt(2) * (sin(x) + cos(x))
1285 %
1286 %  The format of the BesselOrderOne method is:
1287 %
1288 %      MagickRealType BesselOrderOne(MagickRealType x)
1289 %
1290 %  A description of each parameter follows:
1291 %
1292 %    o x: MagickRealType value.
1293 %
1294 */
1295
1296 #undef I0
1297 static MagickRealType I0(MagickRealType x)
1298 {
1299   MagickRealType
1300     sum,
1301     t,
1302     y;
1303
1304   register ssize_t
1305     i;
1306
1307   /*
1308     Zeroth order Bessel function of the first kind.
1309   */
1310   sum=1.0;
1311   y=x*x/4.0;
1312   t=y;
1313   for (i=2; t > MagickEpsilon; i++)
1314   {
1315     sum+=t;
1316     t*=y/((MagickRealType) i*i);
1317   }
1318   return(sum);
1319 }
1320
1321 #undef J1
1322 static MagickRealType J1(MagickRealType x)
1323 {
1324   MagickRealType
1325     p,
1326     q;
1327
1328   register ssize_t
1329     i;
1330
1331   static const double
1332     Pone[] =
1333     {
1334        0.581199354001606143928050809e+21,
1335       -0.6672106568924916298020941484e+20,
1336        0.2316433580634002297931815435e+19,
1337       -0.3588817569910106050743641413e+17,
1338        0.2908795263834775409737601689e+15,
1339       -0.1322983480332126453125473247e+13,
1340        0.3413234182301700539091292655e+10,
1341       -0.4695753530642995859767162166e+7,
1342        0.270112271089232341485679099e+4
1343     },
1344     Qone[] =
1345     {
1346       0.11623987080032122878585294e+22,
1347       0.1185770712190320999837113348e+20,
1348       0.6092061398917521746105196863e+17,
1349       0.2081661221307607351240184229e+15,
1350       0.5243710262167649715406728642e+12,
1351       0.1013863514358673989967045588e+10,
1352       0.1501793594998585505921097578e+7,
1353       0.1606931573481487801970916749e+4,
1354       0.1e+1
1355     };
1356
1357   p=Pone[8];
1358   q=Qone[8];
1359   for (i=7; i >= 0; i--)
1360   {
1361     p=p*x*x+Pone[i];
1362     q=q*x*x+Qone[i];
1363   }
1364   return(p/q);
1365 }
1366
1367 #undef P1
1368 static MagickRealType P1(MagickRealType x)
1369 {
1370   MagickRealType
1371     p,
1372     q;
1373
1374   register ssize_t
1375     i;
1376
1377   static const double
1378     Pone[] =
1379     {
1380       0.352246649133679798341724373e+5,
1381       0.62758845247161281269005675e+5,
1382       0.313539631109159574238669888e+5,
1383       0.49854832060594338434500455e+4,
1384       0.2111529182853962382105718e+3,
1385       0.12571716929145341558495e+1
1386     },
1387     Qone[] =
1388     {
1389       0.352246649133679798068390431e+5,
1390       0.626943469593560511888833731e+5,
1391       0.312404063819041039923015703e+5,
1392       0.4930396490181088979386097e+4,
1393       0.2030775189134759322293574e+3,
1394       0.1e+1
1395     };
1396
1397   p=Pone[5];
1398   q=Qone[5];
1399   for (i=4; i >= 0; i--)
1400   {
1401     p=p*(8.0/x)*(8.0/x)+Pone[i];
1402     q=q*(8.0/x)*(8.0/x)+Qone[i];
1403   }
1404   return(p/q);
1405 }
1406
1407 #undef Q1
1408 static MagickRealType Q1(MagickRealType x)
1409 {
1410   MagickRealType
1411     p,
1412     q;
1413
1414   register ssize_t
1415     i;
1416
1417   static const double
1418     Pone[] =
1419     {
1420       0.3511751914303552822533318e+3,
1421       0.7210391804904475039280863e+3,
1422       0.4259873011654442389886993e+3,
1423       0.831898957673850827325226e+2,
1424       0.45681716295512267064405e+1,
1425       0.3532840052740123642735e-1
1426     },
1427     Qone[] =
1428     {
1429       0.74917374171809127714519505e+4,
1430       0.154141773392650970499848051e+5,
1431       0.91522317015169922705904727e+4,
1432       0.18111867005523513506724158e+4,
1433       0.1038187585462133728776636e+3,
1434       0.1e+1
1435     };
1436
1437   p=Pone[5];
1438   q=Qone[5];
1439   for (i=4; i >= 0; i--)
1440   {
1441     p=p*(8.0/x)*(8.0/x)+Pone[i];
1442     q=q*(8.0/x)*(8.0/x)+Qone[i];
1443   }
1444   return(p/q);
1445 }
1446
1447 static MagickRealType BesselOrderOne(MagickRealType x)
1448 {
1449   MagickRealType
1450     p,
1451     q;
1452
1453   if (x == 0.0)
1454     return(0.0);
1455   p=x;
1456   if (x < 0.0)
1457     x=(-x);
1458   if (x < 8.0)
1459     return(p*J1(x));
1460   q=sqrt((double) (2.0/(MagickPI*x)))*(P1(x)*(1.0/sqrt(2.0)*(sin((double) x)-
1461     cos((double) x)))-8.0/x*Q1(x)*(-1.0/sqrt(2.0)*(sin((double) x)+
1462     cos((double) x))));
1463   if (p < 0.0)
1464     q=(-q);
1465   return(q);
1466 }
1467 \f
1468 /*
1469 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1470 %                                                                             %
1471 %                                                                             %
1472 %                                                                             %
1473 +   D e s t r o y R e s i z e F i l t e r                                     %
1474 %                                                                             %
1475 %                                                                             %
1476 %                                                                             %
1477 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1478 %
1479 %  DestroyResizeFilter() destroy the resize filter.
1480 %
1481 %  The format of the DestroyResizeFilter method is:
1482 %
1483 %      ResizeFilter *DestroyResizeFilter(ResizeFilter *resize_filter)
1484 %
1485 %  A description of each parameter follows:
1486 %
1487 %    o resize_filter: the resize filter.
1488 %
1489 */
1490 MagickPrivate ResizeFilter *DestroyResizeFilter(ResizeFilter *resize_filter)
1491 {
1492   assert(resize_filter != (ResizeFilter *) NULL);
1493   assert(resize_filter->signature == MagickSignature);
1494   resize_filter->signature=(~MagickSignature);
1495   resize_filter=(ResizeFilter *) RelinquishMagickMemory(resize_filter);
1496   return(resize_filter);
1497 }
1498 \f
1499 /*
1500 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1501 %                                                                             %
1502 %                                                                             %
1503 %                                                                             %
1504 +   G e t R e s i z e F i l t e r S u p p o r t                               %
1505 %                                                                             %
1506 %                                                                             %
1507 %                                                                             %
1508 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1509 %
1510 %  GetResizeFilterSupport() return the current support window size for this
1511 %  filter.  Note that this may have been enlarged by filter:blur factor.
1512 %
1513 %  The format of the GetResizeFilterSupport method is:
1514 %
1515 %      MagickRealType GetResizeFilterSupport(const ResizeFilter *resize_filter)
1516 %
1517 %  A description of each parameter follows:
1518 %
1519 %    o filter: Image filter to use.
1520 %
1521 */
1522 MagickPrivate MagickRealType GetResizeFilterSupport(
1523   const ResizeFilter *resize_filter)
1524 {
1525   assert(resize_filter != (ResizeFilter *) NULL);
1526   assert(resize_filter->signature == MagickSignature);
1527   return(resize_filter->support*resize_filter->blur);
1528 }
1529 \f
1530 /*
1531 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1532 %                                                                             %
1533 %                                                                             %
1534 %                                                                             %
1535 +   G e t R e s i z e F i l t e r W e i g h t                                 %
1536 %                                                                             %
1537 %                                                                             %
1538 %                                                                             %
1539 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1540 %
1541 %  GetResizeFilterWeight evaluates the specified resize filter at the point x
1542 %  which usally lies between zero and the filters current 'support' and
1543 %  returns the weight of the filter function at that point.
1544 %
1545 %  The format of the GetResizeFilterWeight method is:
1546 %
1547 %      MagickRealType GetResizeFilterWeight(const ResizeFilter *resize_filter,
1548 %        const MagickRealType x)
1549 %
1550 %  A description of each parameter follows:
1551 %
1552 %    o filter: the filter type.
1553 %
1554 %    o x: the point.
1555 %
1556 */
1557 MagickPrivate MagickRealType GetResizeFilterWeight(
1558   const ResizeFilter *resize_filter,const MagickRealType x)
1559 {
1560   MagickRealType
1561     scale,
1562     weight,
1563     x_blur;
1564
1565   /*
1566     Windowing function - scale the weighting filter by this amount.
1567   */
1568   assert(resize_filter != (ResizeFilter *) NULL);
1569   assert(resize_filter->signature == MagickSignature);
1570   x_blur=fabs((double) x)/resize_filter->blur;  /* X offset with blur scaling */
1571   if ((resize_filter->window_support < MagickEpsilon) ||
1572       (resize_filter->window == Box))
1573     scale=1.0;  /* Point or Box Filter -- avoid division by zero */
1574   else
1575     {
1576       scale=resize_filter->scale;
1577       scale=resize_filter->window(x_blur*scale,resize_filter);
1578     }
1579   weight=scale*resize_filter->filter(x_blur,resize_filter);
1580   return(weight);
1581 }
1582 #if defined(MAGICKCORE_LQR_DELEGATE)
1583 \f
1584 /*
1585 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1586 %                                                                             %
1587 %                                                                             %
1588 %                                                                             %
1589 %   L i q u i d R e s c a l e I m a g e                                       %
1590 %                                                                             %
1591 %                                                                             %
1592 %                                                                             %
1593 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1594 %
1595 %  LiquidRescaleImage() rescales image with seam carving.
1596 %
1597 %  The format of the LiquidRescaleImage method is:
1598 %
1599 %      Image *LiquidRescaleImage(const Image *image,const size_t columns,
1600 %        const size_t rows,const double delta_x,const double rigidity,
1601 %        ExceptionInfo *exception)
1602 %
1603 %  A description of each parameter follows:
1604 %
1605 %    o image: the image.
1606 %
1607 %    o columns: the number of columns in the rescaled image.
1608 %
1609 %    o rows: the number of rows in the rescaled image.
1610 %
1611 %    o delta_x: maximum seam transversal step (0 means straight seams).
1612 %
1613 %    o rigidity: introduce a bias for non-straight seams (typically 0).
1614 %
1615 %    o exception: return any errors or warnings in this structure.
1616 %
1617 */
1618 MagickExport 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 #define LiquidRescaleImageTag  "Rescale/Image"
1623
1624   CacheView
1625     *image_view,
1626     *rescale_view;
1627
1628   gfloat
1629     *packet,
1630     *pixels;
1631
1632   Image
1633     *rescale_image;
1634
1635   int
1636     x_offset,
1637     y_offset;
1638
1639   LqrCarver
1640     *carver;
1641
1642   LqrRetVal
1643     lqr_status;
1644
1645   MagickBooleanType
1646     status;
1647
1648   register gfloat
1649     *q;
1650
1651   ssize_t
1652     y;
1653
1654   /*
1655     Liquid rescale image.
1656   */
1657   assert(image != (const Image *) NULL);
1658   assert(image->signature == MagickSignature);
1659   if (image->debug != MagickFalse)
1660     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1661   assert(exception != (ExceptionInfo *) NULL);
1662   assert(exception->signature == MagickSignature);
1663   if ((columns == 0) || (rows == 0))
1664     return((Image *) NULL);
1665   if ((columns == image->columns) && (rows == image->rows))
1666     return(CloneImage(image,0,0,MagickTrue,exception));
1667   if ((columns <= 2) || (rows <= 2))
1668     return(ResizeImage(image,columns,rows,image->filter,image->blur,exception));
1669   if ((columns >= (2*image->columns)) || (rows >= (2*image->rows)))
1670     {
1671       Image
1672         *resize_image;
1673
1674       size_t
1675         height,
1676         width;
1677
1678       /*
1679         Honor liquid resize size limitations.
1680       */
1681       for (width=image->columns; columns >= (2*width-1); width*=2);
1682       for (height=image->rows; rows >= (2*height-1); height*=2);
1683       resize_image=ResizeImage(image,width,height,image->filter,image->blur,
1684         exception);
1685       if (resize_image == (Image *) NULL)
1686         return((Image *) NULL);
1687       rescale_image=LiquidRescaleImage(resize_image,columns,rows,delta_x,
1688         rigidity,exception);
1689       resize_image=DestroyImage(resize_image);
1690       return(rescale_image);
1691     }
1692   pixels=(gfloat *) AcquireQuantumMemory(image->columns,image->rows*
1693     GetPixelChannels(image)*sizeof(*pixels));
1694   if (pixels == (gfloat *) NULL)
1695     return((Image *) NULL);
1696   status=MagickTrue;
1697   q=pixels;
1698   image_view=AcquireCacheView(image);
1699   for (y=0; y < (ssize_t) image->rows; y++)
1700   {
1701     register const Quantum
1702       *restrict p;
1703
1704     register ssize_t
1705       x;
1706
1707     if (status == MagickFalse)
1708       continue;
1709     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1710     if (p == (const Quantum *) NULL)
1711       {
1712         status=MagickFalse;
1713         continue;
1714       }
1715     for (x=0; x < (ssize_t) image->columns; x++)
1716     {
1717       register ssize_t
1718         i;
1719
1720       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1721         *q++=QuantumScale*p[i];
1722       p+=GetPixelChannels(image);
1723     }
1724   }
1725   image_view=DestroyCacheView(image_view);
1726   carver=lqr_carver_new_ext(pixels,image->columns,image->rows,
1727     GetPixelChannels(image),LQR_COLDEPTH_32F);
1728   if (carver == (LqrCarver *) NULL)
1729     {
1730       pixels=(gfloat *) RelinquishMagickMemory(pixels);
1731       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1732     }
1733   lqr_status=lqr_carver_init(carver,(int) delta_x,rigidity);
1734   lqr_status=lqr_carver_resize(carver,columns,rows);
1735   (void) lqr_status;
1736   rescale_image=CloneImage(image,lqr_carver_get_width(carver),
1737     lqr_carver_get_height(carver),MagickTrue,exception);
1738   if (rescale_image == (Image *) NULL)
1739     {
1740       pixels=(gfloat *) RelinquishMagickMemory(pixels);
1741       return((Image *) NULL);
1742     }
1743   if (SetImageStorageClass(rescale_image,DirectClass,exception) == MagickFalse)
1744     {
1745       pixels=(gfloat *) RelinquishMagickMemory(pixels);
1746       rescale_image=DestroyImage(rescale_image);
1747       return((Image *) NULL);
1748     }
1749   rescale_view=AcquireCacheView(rescale_image);
1750   (void) lqr_carver_scan_reset(carver);
1751   while (lqr_carver_scan_ext(carver,&x_offset,&y_offset,(void **) &packet) != 0)
1752   {
1753     register Quantum
1754       *restrict q;
1755
1756     register ssize_t
1757       i;
1758
1759     q=QueueCacheViewAuthenticPixels(rescale_view,x_offset,y_offset,1,1,
1760       exception);
1761     if (q == (Quantum *) NULL)
1762       break;
1763     for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1764     {
1765       PixelChannel
1766         channel;
1767
1768       PixelTrait
1769         rescale_traits,
1770         traits;
1771
1772       traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
1773       channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
1774       rescale_traits=GetPixelChannelMapTraits(rescale_image,channel);
1775       if ((traits == UndefinedPixelTrait) ||
1776           (rescale_traits == UndefinedPixelTrait))
1777         continue;
1778       SetPixelChannel(rescale_image,channel,ClampToQuantum(QuantumRange*
1779         packet[i]),q);
1780     }
1781     if (SyncCacheViewAuthenticPixels(rescale_view,exception) == MagickFalse)
1782       break;
1783   }
1784   rescale_view=DestroyCacheView(rescale_view);
1785   lqr_carver_destroy(carver);
1786   return(rescale_image);
1787 }
1788 #else
1789 MagickExport Image *LiquidRescaleImage(const Image *image,
1790   const size_t magick_unused(columns),const size_t magick_unused(rows),
1791   const double magick_unused(delta_x),const double magick_unused(rigidity),
1792   ExceptionInfo *exception)
1793 {
1794   assert(image != (const Image *) NULL);
1795   assert(image->signature == MagickSignature);
1796   if (image->debug != MagickFalse)
1797     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1798   assert(exception != (ExceptionInfo *) NULL);
1799   assert(exception->signature == MagickSignature);
1800   (void) ThrowMagickException(exception,GetMagickModule(),MissingDelegateError,
1801     "DelegateLibrarySupportNotBuiltIn","`%s' (LQR)",image->filename);
1802   return((Image *) NULL);
1803 }
1804 #endif
1805 \f
1806 /*
1807 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1808 %                                                                             %
1809 %                                                                             %
1810 %                                                                             %
1811 %   M a g n i f y I m a g e                                                   %
1812 %                                                                             %
1813 %                                                                             %
1814 %                                                                             %
1815 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1816 %
1817 %  MagnifyImage() is a convenience method that scales an image proportionally
1818 %  to twice its size.
1819 %
1820 %  The format of the MagnifyImage method is:
1821 %
1822 %      Image *MagnifyImage(const Image *image,ExceptionInfo *exception)
1823 %
1824 %  A description of each parameter follows:
1825 %
1826 %    o image: the image.
1827 %
1828 %    o exception: return any errors or warnings in this structure.
1829 %
1830 */
1831 MagickExport Image *MagnifyImage(const Image *image,ExceptionInfo *exception)
1832 {
1833   Image
1834     *magnify_image;
1835
1836   assert(image != (Image *) NULL);
1837   assert(image->signature == MagickSignature);
1838   if (image->debug != MagickFalse)
1839     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1840   assert(exception != (ExceptionInfo *) NULL);
1841   assert(exception->signature == MagickSignature);
1842   magnify_image=ResizeImage(image,2*image->columns,2*image->rows,CubicFilter,
1843     1.0,exception);
1844   return(magnify_image);
1845 }
1846 \f
1847 /*
1848 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1849 %                                                                             %
1850 %                                                                             %
1851 %                                                                             %
1852 %   M i n i f y I m a g e                                                     %
1853 %                                                                             %
1854 %                                                                             %
1855 %                                                                             %
1856 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1857 %
1858 %  MinifyImage() is a convenience method that scales an image proportionally to
1859 %  half its size.
1860 %
1861 %  The format of the MinifyImage method is:
1862 %
1863 %      Image *MinifyImage(const Image *image,ExceptionInfo *exception)
1864 %
1865 %  A description of each parameter follows:
1866 %
1867 %    o image: the image.
1868 %
1869 %    o exception: return any errors or warnings in this structure.
1870 %
1871 */
1872 MagickExport Image *MinifyImage(const Image *image,ExceptionInfo *exception)
1873 {
1874   Image
1875     *minify_image;
1876
1877   assert(image != (Image *) NULL);
1878   assert(image->signature == MagickSignature);
1879   if (image->debug != MagickFalse)
1880     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1881   assert(exception != (ExceptionInfo *) NULL);
1882   assert(exception->signature == MagickSignature);
1883   minify_image=ResizeImage(image,image->columns/2,image->rows/2,CubicFilter,1.0,
1884     exception);
1885   return(minify_image);
1886 }
1887 \f
1888 /*
1889 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1890 %                                                                             %
1891 %                                                                             %
1892 %                                                                             %
1893 %   R e s a m p l e I m a g e                                                 %
1894 %                                                                             %
1895 %                                                                             %
1896 %                                                                             %
1897 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1898 %
1899 %  ResampleImage() resize image in terms of its pixel size, so that when
1900 %  displayed at the given resolution it will be the same size in terms of
1901 %  real world units as the original image at the original resolution.
1902 %
1903 %  The format of the ResampleImage method is:
1904 %
1905 %      Image *ResampleImage(Image *image,const double x_resolution,
1906 %        const double y_resolution,const FilterTypes filter,const double blur,
1907 %        ExceptionInfo *exception)
1908 %
1909 %  A description of each parameter follows:
1910 %
1911 %    o image: the image to be resized to fit the given resolution.
1912 %
1913 %    o x_resolution: the new image x resolution.
1914 %
1915 %    o y_resolution: the new image y resolution.
1916 %
1917 %    o filter: Image filter to use.
1918 %
1919 %    o blur: the blur factor where > 1 is blurry, < 1 is sharp.
1920 %
1921 */
1922 MagickExport Image *ResampleImage(const Image *image,const double x_resolution,
1923   const double y_resolution,const FilterTypes filter,const double blur,
1924   ExceptionInfo *exception)
1925 {
1926 #define ResampleImageTag  "Resample/Image"
1927
1928   Image
1929     *resample_image;
1930
1931   size_t
1932     height,
1933     width;
1934
1935   /*
1936     Initialize sampled image attributes.
1937   */
1938   assert(image != (const Image *) NULL);
1939   assert(image->signature == MagickSignature);
1940   if (image->debug != MagickFalse)
1941     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1942   assert(exception != (ExceptionInfo *) NULL);
1943   assert(exception->signature == MagickSignature);
1944   width=(size_t) (x_resolution*image->columns/(image->resolution.x == 0.0 ?
1945     72.0 : image->resolution.x)+0.5);
1946   height=(size_t) (y_resolution*image->rows/(image->resolution.y == 0.0 ?
1947     72.0 : image->resolution.y)+0.5);
1948   resample_image=ResizeImage(image,width,height,filter,blur,exception);
1949   if (resample_image != (Image *) NULL)
1950     {
1951       resample_image->resolution.x=x_resolution;
1952       resample_image->resolution.y=y_resolution;
1953     }
1954   return(resample_image);
1955 }
1956 \f
1957 /*
1958 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1959 %                                                                             %
1960 %                                                                             %
1961 %                                                                             %
1962 %   R e s i z e I m a g e                                                     %
1963 %                                                                             %
1964 %                                                                             %
1965 %                                                                             %
1966 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1967 %
1968 %  ResizeImage() scales an image to the desired dimensions, using the given
1969 %  filter (see AcquireFilterInfo()).
1970 %
1971 %  If an undefined filter is given the filter defaults to Mitchell for a
1972 %  colormapped image, a image with a matte channel, or if the image is
1973 %  enlarged.  Otherwise the filter defaults to a Lanczos.
1974 %
1975 %  ResizeImage() was inspired by Paul Heckbert's "zoom" program.
1976 %
1977 %  The format of the ResizeImage method is:
1978 %
1979 %      Image *ResizeImage(Image *image,const size_t columns,
1980 %        const size_t rows,const FilterTypes filter,const double blur,
1981 %        ExceptionInfo *exception)
1982 %
1983 %  A description of each parameter follows:
1984 %
1985 %    o image: the image.
1986 %
1987 %    o columns: the number of columns in the scaled image.
1988 %
1989 %    o rows: the number of rows in the scaled image.
1990 %
1991 %    o filter: Image filter to use.
1992 %
1993 %    o blur: the blur factor where > 1 is blurry, < 1 is sharp.  Typically set
1994 %      this to 1.0.
1995 %
1996 %    o exception: return any errors or warnings in this structure.
1997 %
1998 */
1999
2000 typedef struct _ContributionInfo
2001 {
2002   MagickRealType
2003     weight;
2004
2005   ssize_t
2006     pixel;
2007 } ContributionInfo;
2008
2009 static ContributionInfo **DestroyContributionThreadSet(
2010   ContributionInfo **contribution)
2011 {
2012   register ssize_t
2013     i;
2014
2015   assert(contribution != (ContributionInfo **) NULL);
2016   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2017     if (contribution[i] != (ContributionInfo *) NULL)
2018       contribution[i]=(ContributionInfo *) RelinquishMagickMemory(
2019         contribution[i]);
2020   contribution=(ContributionInfo **) RelinquishMagickMemory(contribution);
2021   return(contribution);
2022 }
2023
2024 static ContributionInfo **AcquireContributionThreadSet(const size_t count)
2025 {
2026   register ssize_t
2027     i;
2028
2029   ContributionInfo
2030     **contribution;
2031
2032   size_t
2033     number_threads;
2034
2035   number_threads=GetOpenMPMaximumThreads();
2036   contribution=(ContributionInfo **) AcquireQuantumMemory(number_threads,
2037     sizeof(*contribution));
2038   if (contribution == (ContributionInfo **) NULL)
2039     return((ContributionInfo **) NULL);
2040   (void) ResetMagickMemory(contribution,0,number_threads*sizeof(*contribution));
2041   for (i=0; i < (ssize_t) number_threads; i++)
2042   {
2043     contribution[i]=(ContributionInfo *) AcquireQuantumMemory(count,
2044       sizeof(**contribution));
2045     if (contribution[i] == (ContributionInfo *) NULL)
2046       return(DestroyContributionThreadSet(contribution));
2047   }
2048   return(contribution);
2049 }
2050
2051 static inline double MagickMax(const double x,const double y)
2052 {
2053   if (x > y)
2054     return(x);
2055   return(y);
2056 }
2057
2058 static inline double MagickMin(const double x,const double y)
2059 {
2060   if (x < y)
2061     return(x);
2062   return(y);
2063 }
2064
2065 static MagickBooleanType HorizontalFilter(const ResizeFilter *resize_filter,
2066   const Image *image,Image *resize_image,const MagickRealType x_factor,
2067   const MagickSizeType span,MagickOffsetType *offset,ExceptionInfo *exception)
2068 {
2069 #define ResizeImageTag  "Resize/Image"
2070
2071   CacheView
2072     *image_view,
2073     *resize_view;
2074
2075   ClassType
2076     storage_class;
2077
2078   ContributionInfo
2079     **restrict contributions;
2080
2081   MagickBooleanType
2082     status;
2083
2084   MagickRealType
2085     scale,
2086     support;
2087
2088   ssize_t
2089     x;
2090
2091   /*
2092     Apply filter to resize horizontally from image to resize image.
2093   */
2094   scale=MagickMax(1.0/x_factor+MagickEpsilon,1.0);
2095   support=scale*GetResizeFilterSupport(resize_filter);
2096   storage_class=support > 0.5 ? DirectClass : image->storage_class;
2097   if (SetImageStorageClass(resize_image,storage_class,exception) == MagickFalse)
2098     return(MagickFalse);
2099   if (support < 0.5)
2100     {
2101       /*
2102         Support too small even for nearest neighbour: Reduce to point sampling.
2103       */
2104       support=(MagickRealType) 0.5;
2105       scale=1.0;
2106     }
2107   contributions=AcquireContributionThreadSet((size_t) (2.0*support+3.0));
2108   if (contributions == (ContributionInfo **) NULL)
2109     {
2110       (void) ThrowMagickException(exception,GetMagickModule(),
2111         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2112       return(MagickFalse);
2113     }
2114   status=MagickTrue;
2115   scale=1.0/scale;
2116   image_view=AcquireCacheView(image);
2117   resize_view=AcquireCacheView(resize_image);
2118 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2119   #pragma omp parallel for shared(status)
2120 #endif
2121   for (x=0; x < (ssize_t) resize_image->columns; x++)
2122   {
2123     MagickRealType
2124       bisect,
2125       density;
2126
2127     register const Quantum
2128       *restrict p;
2129
2130     register ContributionInfo
2131       *restrict contribution;
2132
2133     register Quantum
2134       *restrict q;
2135
2136     register ssize_t
2137       y;
2138
2139     ssize_t
2140       n,
2141       start,
2142       stop;
2143
2144     if (status == MagickFalse)
2145       continue;
2146     bisect=(MagickRealType) (x+0.5)/x_factor;
2147     start=(ssize_t) MagickMax(bisect-support+0.5,0.0);
2148     stop=(ssize_t) MagickMin(bisect+support+0.5,(double) image->columns);
2149     density=0.0;
2150     contribution=contributions[GetOpenMPThreadId()];
2151     for (n=0; n < (stop-start); n++)
2152     {
2153       contribution[n].pixel=start+n;
2154       contribution[n].weight=GetResizeFilterWeight(resize_filter,scale*
2155         ((MagickRealType) (start+n)-bisect+0.5));
2156       density+=contribution[n].weight;
2157     }
2158     if ((density != 0.0) && (density != 1.0))
2159       {
2160         register ssize_t
2161           i;
2162
2163         /*
2164           Normalize.
2165         */
2166         density=1.0/density;
2167         for (i=0; i < n; i++)
2168           contribution[i].weight*=density;
2169       }
2170     p=GetCacheViewVirtualPixels(image_view,contribution[0].pixel,0,(size_t)
2171       (contribution[n-1].pixel-contribution[0].pixel+1),image->rows,exception);
2172     q=QueueCacheViewAuthenticPixels(resize_view,x,0,1,resize_image->rows,
2173       exception);
2174     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2175       {
2176         status=MagickFalse;
2177         continue;
2178       }
2179     for (y=0; y < (ssize_t) resize_image->rows; y++)
2180     {
2181       register ssize_t
2182         i;
2183
2184       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2185       {
2186         MagickRealType
2187           alpha,
2188           gamma,
2189           pixel;
2190
2191         PixelChannel
2192           channel;
2193
2194         PixelTrait
2195           resize_traits,
2196           traits;
2197
2198         register ssize_t
2199           j;
2200
2201         ssize_t
2202           k;
2203
2204         traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
2205         channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
2206         resize_traits=GetPixelChannelMapTraits(resize_image,channel);
2207         if ((traits == UndefinedPixelTrait) ||
2208             (resize_traits == UndefinedPixelTrait))
2209           continue;
2210         if ((resize_traits & CopyPixelTrait) != 0)
2211           {
2212             j=(ssize_t) (MagickMin(MagickMax(bisect,(double) start),(double)
2213               stop-1.0)+0.5);
2214             k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2215               (contribution[j-start].pixel-contribution[0].pixel);
2216             SetPixelChannel(resize_image,channel,p[k*GetPixelChannels(image)+i],
2217               q);
2218             continue;
2219           }
2220         pixel=0.0;
2221         if ((resize_traits & BlendPixelTrait) == 0)
2222           {
2223             /*
2224               No alpha blending.
2225             */
2226             for (j=0; j < n; j++)
2227             {
2228               k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2229                 (contribution[j].pixel-contribution[0].pixel);
2230               alpha=contribution[j].weight;
2231               pixel+=alpha*p[k*GetPixelChannels(image)+i];
2232             }
2233             SetPixelChannel(resize_image,channel,ClampToQuantum(pixel),q);
2234             continue;
2235           }
2236         /*
2237           Alpha blending.
2238         */
2239         gamma=0.0;
2240         for (j=0; j < n; j++)
2241         {
2242           k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2243             (contribution[j].pixel-contribution[0].pixel);
2244           alpha=contribution[j].weight*QuantumScale*
2245             GetPixelAlpha(image,p+k*GetPixelChannels(image));
2246           pixel+=alpha*p[k*GetPixelChannels(image)+i];
2247           gamma+=alpha;
2248         }
2249         gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2250         SetPixelChannel(resize_image,channel,ClampToQuantum(gamma*pixel),q);
2251       }
2252       q+=GetPixelChannels(resize_image);
2253     }
2254     if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
2255       status=MagickFalse;
2256     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2257       {
2258         MagickBooleanType
2259           proceed;
2260
2261 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2262   #pragma omp critical (MagickCore_HorizontalFilter)
2263 #endif
2264         proceed=SetImageProgress(image,ResizeImageTag,(*offset)++,span);
2265         if (proceed == MagickFalse)
2266           status=MagickFalse;
2267       }
2268   }
2269   resize_view=DestroyCacheView(resize_view);
2270   image_view=DestroyCacheView(image_view);
2271   contributions=DestroyContributionThreadSet(contributions);
2272   return(status);
2273 }
2274
2275 static MagickBooleanType VerticalFilter(const ResizeFilter *resize_filter,
2276   const Image *image,Image *resize_image,const MagickRealType y_factor,
2277   const MagickSizeType span,MagickOffsetType *offset,ExceptionInfo *exception)
2278 {
2279   CacheView
2280     *image_view,
2281     *resize_view;
2282
2283   ClassType
2284     storage_class;
2285
2286   ContributionInfo
2287     **restrict contributions;
2288
2289   MagickBooleanType
2290     status;
2291
2292   PixelInfo
2293     zero;
2294
2295   MagickRealType
2296     scale,
2297     support;
2298
2299   ssize_t
2300     y;
2301
2302   /*
2303     Apply filter to resize vertically from image to resize image.
2304   */
2305   scale=MagickMax(1.0/y_factor+MagickEpsilon,1.0);
2306   support=scale*GetResizeFilterSupport(resize_filter);
2307   storage_class=support > 0.5 ? DirectClass : image->storage_class;
2308   if (SetImageStorageClass(resize_image,storage_class,exception) == MagickFalse)
2309     return(MagickFalse);
2310   if (support < 0.5)
2311     {
2312       /*
2313         Support too small even for nearest neighbour: Reduce to point sampling.
2314       */
2315       support=(MagickRealType) 0.5;
2316       scale=1.0;
2317     }
2318   contributions=AcquireContributionThreadSet((size_t) (2.0*support+3.0));
2319   if (contributions == (ContributionInfo **) NULL)
2320     {
2321       (void) ThrowMagickException(exception,GetMagickModule(),
2322         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2323       return(MagickFalse);
2324     }
2325   status=MagickTrue;
2326   scale=1.0/scale;
2327   (void) ResetMagickMemory(&zero,0,sizeof(zero));
2328   image_view=AcquireCacheView(image);
2329   resize_view=AcquireCacheView(resize_image);
2330 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2331   #pragma omp parallel for shared(status)
2332 #endif
2333   for (y=0; y < (ssize_t) resize_image->rows; y++)
2334   {
2335     MagickRealType
2336       bisect,
2337       density;
2338
2339     register const Quantum
2340       *restrict p;
2341
2342     register ContributionInfo
2343       *restrict contribution;
2344
2345     register Quantum
2346       *restrict q;
2347
2348     register ssize_t
2349       x;
2350
2351     ssize_t
2352       n,
2353       start,
2354       stop;
2355
2356     if (status == MagickFalse)
2357       continue;
2358     bisect=(MagickRealType) (y+0.5)/y_factor;
2359     start=(ssize_t) MagickMax(bisect-support+0.5,0.0);
2360     stop=(ssize_t) MagickMin(bisect+support+0.5,(double) image->rows);
2361     density=0.0;
2362     contribution=contributions[GetOpenMPThreadId()];
2363     for (n=0; n < (stop-start); n++)
2364     {
2365       contribution[n].pixel=start+n;
2366       contribution[n].weight=GetResizeFilterWeight(resize_filter,scale*
2367         ((MagickRealType) (start+n)-bisect+0.5));
2368       density+=contribution[n].weight;
2369     }
2370     if ((density != 0.0) && (density != 1.0))
2371       {
2372         register ssize_t
2373           i;
2374
2375         /*
2376           Normalize.
2377         */
2378         density=1.0/density;
2379         for (i=0; i < n; i++)
2380           contribution[i].weight*=density;
2381       }
2382     p=GetCacheViewVirtualPixels(image_view,0,contribution[0].pixel,
2383       image->columns,(size_t) (contribution[n-1].pixel-contribution[0].pixel+1),
2384       exception);
2385     q=QueueCacheViewAuthenticPixels(resize_view,0,y,resize_image->columns,1,
2386       exception);
2387     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2388       {
2389         status=MagickFalse;
2390         continue;
2391       }
2392     for (x=0; x < (ssize_t) resize_image->columns; x++)
2393     {
2394       register ssize_t
2395         i;
2396
2397       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2398       {
2399         MagickRealType
2400           alpha,
2401           gamma,
2402           pixel;
2403
2404         PixelChannel
2405           channel;
2406
2407         PixelTrait
2408           resize_traits,
2409           traits;
2410
2411         register ssize_t
2412           j;
2413
2414         ssize_t
2415           k;
2416
2417         traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
2418         channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
2419         resize_traits=GetPixelChannelMapTraits(resize_image,channel);
2420         if ((traits == UndefinedPixelTrait) ||
2421             (resize_traits == UndefinedPixelTrait))
2422           continue;
2423         if ((resize_traits & CopyPixelTrait) != 0)
2424           {
2425             j=(ssize_t) (MagickMin(MagickMax(bisect,(double) start),(double)
2426               stop-1.0)+0.5);
2427             k=(ssize_t) ((contribution[j-start].pixel-contribution[0].pixel)*
2428               image->columns+x);
2429             SetPixelChannel(resize_image,channel,p[k*GetPixelChannels(image)+i],
2430               q);
2431             continue;
2432           }
2433         pixel=0.0;
2434         if ((resize_traits & BlendPixelTrait) == 0)
2435           {
2436             /*
2437               No alpha blending.
2438             */
2439             for (j=0; j < n; j++)
2440             {
2441               k=(ssize_t) ((contribution[j].pixel-contribution[0].pixel)*
2442                 image->columns+x);
2443               alpha=contribution[j].weight;
2444               pixel+=alpha*p[k*GetPixelChannels(image)+i];
2445             }
2446             SetPixelChannel(resize_image,channel,ClampToQuantum(pixel),q);
2447             continue;
2448           }
2449         gamma=0.0;
2450         for (j=0; j < n; j++)
2451         {
2452           k=(ssize_t) ((contribution[j].pixel-contribution[0].pixel)*
2453             image->columns+x);
2454           alpha=contribution[j].weight*QuantumScale*
2455             GetPixelAlpha(image,p+k*GetPixelChannels(image));
2456           pixel+=alpha*p[k*GetPixelChannels(image)+i];
2457           gamma+=alpha;
2458         }
2459         gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2460         SetPixelChannel(resize_image,channel,ClampToQuantum(gamma*pixel),q);
2461       }
2462       q+=GetPixelChannels(resize_image);
2463     }
2464     if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
2465       status=MagickFalse;
2466     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2467       {
2468         MagickBooleanType
2469           proceed;
2470
2471 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2472   #pragma omp critical (MagickCore_VerticalFilter)
2473 #endif
2474         proceed=SetImageProgress(image,ResizeImageTag,(*offset)++,span);
2475         if (proceed == MagickFalse)
2476           status=MagickFalse;
2477       }
2478   }
2479   resize_view=DestroyCacheView(resize_view);
2480   image_view=DestroyCacheView(image_view);
2481   contributions=DestroyContributionThreadSet(contributions);
2482   return(status);
2483 }
2484
2485 MagickExport Image *ResizeImage(const Image *image,const size_t columns,
2486   const size_t rows,const FilterTypes filter,const double blur,
2487   ExceptionInfo *exception)
2488 {
2489 #define WorkLoadFactor  0.265
2490
2491   FilterTypes
2492     filter_type;
2493
2494   Image
2495     *filter_image,
2496     *resize_image;
2497
2498   MagickOffsetType
2499     offset;
2500
2501   MagickRealType
2502     x_factor,
2503     y_factor;
2504
2505   MagickSizeType
2506     span;
2507
2508   MagickStatusType
2509     status;
2510
2511   ResizeFilter
2512     *resize_filter;
2513
2514   /*
2515     Acquire resize image.
2516   */
2517   assert(image != (Image *) NULL);
2518   assert(image->signature == MagickSignature);
2519   if (image->debug != MagickFalse)
2520     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2521   assert(exception != (ExceptionInfo *) NULL);
2522   assert(exception->signature == MagickSignature);
2523   if ((columns == 0) || (rows == 0))
2524     ThrowImageException(ImageError,"NegativeOrZeroImageSize");
2525   if ((columns == image->columns) && (rows == image->rows) &&
2526       (filter == UndefinedFilter) && (blur == 1.0))
2527     return(CloneImage(image,0,0,MagickTrue,exception));
2528   resize_image=CloneImage(image,columns,rows,MagickTrue,exception);
2529   if (resize_image == (Image *) NULL)
2530     return(resize_image);
2531   /*
2532     Acquire resize filter.
2533   */
2534   x_factor=(MagickRealType) columns/(MagickRealType) image->columns;
2535   y_factor=(MagickRealType) rows/(MagickRealType) image->rows;
2536   if ((x_factor*y_factor) > WorkLoadFactor)
2537     filter_image=CloneImage(image,columns,image->rows,MagickTrue,exception);
2538   else
2539     filter_image=CloneImage(image,image->columns,rows,MagickTrue,exception);
2540   if (filter_image == (Image *) NULL)
2541     return(DestroyImage(resize_image));
2542   filter_type=LanczosFilter;
2543   if (filter != UndefinedFilter)
2544     filter_type=filter;
2545   else
2546     if ((x_factor == 1.0) && (y_factor == 1.0))
2547       filter_type=PointFilter;
2548     else
2549       if ((image->storage_class == PseudoClass) ||
2550           (image->matte != MagickFalse) || ((x_factor*y_factor) > 1.0))
2551         filter_type=MitchellFilter;
2552   resize_filter=AcquireResizeFilter(image,filter_type,blur,MagickFalse,
2553     exception);
2554   /*
2555     Resize image.
2556   */
2557   offset=0;
2558   if ((x_factor*y_factor) > WorkLoadFactor)
2559     {
2560       span=(MagickSizeType) (filter_image->columns+rows);
2561       status=HorizontalFilter(resize_filter,image,filter_image,x_factor,span,
2562         &offset,exception);
2563       status&=VerticalFilter(resize_filter,filter_image,resize_image,y_factor,
2564         span,&offset,exception);
2565     }
2566   else
2567     {
2568       span=(MagickSizeType) (filter_image->rows+columns);
2569       status=VerticalFilter(resize_filter,image,filter_image,y_factor,span,
2570         &offset,exception);
2571       status&=HorizontalFilter(resize_filter,filter_image,resize_image,x_factor,
2572         span,&offset,exception);
2573     }
2574   /*
2575     Free resources.
2576   */
2577   filter_image=DestroyImage(filter_image);
2578   resize_filter=DestroyResizeFilter(resize_filter);
2579   if (status == MagickFalse)
2580     {
2581       resize_image=DestroyImage(resize_image);
2582       return((Image *) NULL);
2583     }
2584   resize_image->type=image->type;
2585   return(resize_image);
2586 }
2587 \f
2588 /*
2589 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2590 %                                                                             %
2591 %                                                                             %
2592 %                                                                             %
2593 %   S a m p l e I m a g e                                                     %
2594 %                                                                             %
2595 %                                                                             %
2596 %                                                                             %
2597 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2598 %
2599 %  SampleImage() scales an image to the desired dimensions with pixel
2600 %  sampling.  Unlike other scaling methods, this method does not introduce
2601 %  any additional color into the scaled image.
2602 %
2603 %  The format of the SampleImage method is:
2604 %
2605 %      Image *SampleImage(const Image *image,const size_t columns,
2606 %        const size_t rows,ExceptionInfo *exception)
2607 %
2608 %  A description of each parameter follows:
2609 %
2610 %    o image: the image.
2611 %
2612 %    o columns: the number of columns in the sampled image.
2613 %
2614 %    o rows: the number of rows in the sampled image.
2615 %
2616 %    o exception: return any errors or warnings in this structure.
2617 %
2618 */
2619 MagickExport Image *SampleImage(const Image *image,const size_t columns,
2620   const size_t rows,ExceptionInfo *exception)
2621 {
2622 #define SampleImageTag  "Sample/Image"
2623
2624   CacheView
2625     *image_view,
2626     *sample_view;
2627
2628   Image
2629     *sample_image;
2630
2631   MagickBooleanType
2632     status;
2633
2634   MagickOffsetType
2635     progress;
2636
2637   register ssize_t
2638     x;
2639
2640   ssize_t
2641     *x_offset,
2642     y;
2643
2644   /*
2645     Initialize sampled image attributes.
2646   */
2647   assert(image != (const Image *) NULL);
2648   assert(image->signature == MagickSignature);
2649   if (image->debug != MagickFalse)
2650     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2651   assert(exception != (ExceptionInfo *) NULL);
2652   assert(exception->signature == MagickSignature);
2653   if ((columns == 0) || (rows == 0))
2654     ThrowImageException(ImageError,"NegativeOrZeroImageSize");
2655   if ((columns == image->columns) && (rows == image->rows))
2656     return(CloneImage(image,0,0,MagickTrue,exception));
2657   sample_image=CloneImage(image,columns,rows,MagickTrue,exception);
2658   if (sample_image == (Image *) NULL)
2659     return((Image *) NULL);
2660   /*
2661     Allocate scan line buffer and column offset buffers.
2662   */
2663   x_offset=(ssize_t *) AcquireQuantumMemory((size_t) sample_image->columns,
2664     sizeof(*x_offset));
2665   if (x_offset == (ssize_t *) NULL)
2666     {
2667       sample_image=DestroyImage(sample_image);
2668       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2669     }
2670   for (x=0; x < (ssize_t) sample_image->columns; x++)
2671     x_offset[x]=(ssize_t) (((MagickRealType) x+0.5)*image->columns/
2672       sample_image->columns);
2673   /*
2674     Sample each row.
2675   */
2676   status=MagickTrue;
2677   progress=0;
2678   image_view=AcquireCacheView(image);
2679   sample_view=AcquireCacheView(sample_image);
2680 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2681   #pragma omp parallel for schedule(dynamic,4) shared(progress,status) omp_throttle(1)
2682 #endif
2683   for (y=0; y < (ssize_t) sample_image->rows; y++)
2684   {
2685     register const Quantum
2686       *restrict p;
2687
2688     register Quantum
2689       *restrict q;
2690
2691     register ssize_t
2692       x;
2693
2694     ssize_t
2695       y_offset;
2696
2697     if (status == MagickFalse)
2698       continue;
2699     y_offset=(ssize_t) (((MagickRealType) y+0.5)*image->rows/
2700       sample_image->rows);
2701     p=GetCacheViewVirtualPixels(image_view,0,y_offset,image->columns,1,
2702       exception);
2703     q=QueueCacheViewAuthenticPixels(sample_view,0,y,sample_image->columns,1,
2704       exception);
2705     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2706       {
2707         status=MagickFalse;
2708         continue;
2709       }
2710     /*
2711       Sample each column.
2712     */
2713     for (x=0; x < (ssize_t) sample_image->columns; x++)
2714     {
2715       register ssize_t
2716         i;
2717
2718       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2719       {
2720         PixelChannel
2721           channel;
2722
2723         PixelTrait
2724           sample_traits,
2725           traits;
2726
2727         traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
2728         channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
2729         sample_traits=GetPixelChannelMapTraits(sample_image,channel);
2730         if ((traits == UndefinedPixelTrait) ||
2731             (sample_traits == UndefinedPixelTrait))
2732           continue;
2733         SetPixelChannel(sample_image,channel,p[x_offset[x]*GetPixelChannels(
2734           image)+i],q);
2735       }
2736       q+=GetPixelChannels(sample_image);
2737     }
2738     if (SyncCacheViewAuthenticPixels(sample_view,exception) == MagickFalse)
2739       status=MagickFalse;
2740     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2741       {
2742         MagickBooleanType
2743           proceed;
2744
2745 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2746         #pragma omp critical (MagickCore_SampleImage)
2747 #endif
2748         proceed=SetImageProgress(image,SampleImageTag,progress++,image->rows);
2749         if (proceed == MagickFalse)
2750           status=MagickFalse;
2751       }
2752   }
2753   image_view=DestroyCacheView(image_view);
2754   sample_view=DestroyCacheView(sample_view);
2755   x_offset=(ssize_t *) RelinquishMagickMemory(x_offset);
2756   sample_image->type=image->type;
2757   return(sample_image);
2758 }
2759 \f
2760 /*
2761 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2762 %                                                                             %
2763 %                                                                             %
2764 %                                                                             %
2765 %   S c a l e I m a g e                                                       %
2766 %                                                                             %
2767 %                                                                             %
2768 %                                                                             %
2769 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2770 %
2771 %  ScaleImage() changes the size of an image to the given dimensions.
2772 %
2773 %  The format of the ScaleImage method is:
2774 %
2775 %      Image *ScaleImage(const Image *image,const size_t columns,
2776 %        const size_t rows,ExceptionInfo *exception)
2777 %
2778 %  A description of each parameter follows:
2779 %
2780 %    o image: the image.
2781 %
2782 %    o columns: the number of columns in the scaled image.
2783 %
2784 %    o rows: the number of rows in the scaled image.
2785 %
2786 %    o exception: return any errors or warnings in this structure.
2787 %
2788 */
2789 MagickExport Image *ScaleImage(const Image *image,const size_t columns,
2790   const size_t rows,ExceptionInfo *exception)
2791 {
2792 #define ScaleImageTag  "Scale/Image"
2793
2794   CacheView
2795     *image_view,
2796     *scale_view;
2797
2798   Image
2799     *scale_image;
2800
2801   MagickBooleanType
2802     next_column,
2803     next_row,
2804     proceed;
2805
2806   MagickRealType
2807     alpha,
2808     gamma,
2809     pixel[CompositePixelChannel],
2810     *scale_scanline,
2811     *scanline,
2812     *x_vector,
2813     *y_vector;
2814
2815   PixelChannel
2816     channel;
2817
2818   PixelTrait
2819     scale_traits,
2820     traits;
2821
2822   PointInfo
2823     scale,
2824     span;
2825
2826   register ssize_t
2827     i;
2828
2829   ssize_t
2830     n,
2831     number_rows,
2832     y;
2833
2834   /*
2835     Initialize scaled image attributes.
2836   */
2837   assert(image != (const Image *) NULL);
2838   assert(image->signature == MagickSignature);
2839   if (image->debug != MagickFalse)
2840     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2841   assert(exception != (ExceptionInfo *) NULL);
2842   assert(exception->signature == MagickSignature);
2843   if ((columns == 0) || (rows == 0))
2844     return((Image *) NULL);
2845   if ((columns == image->columns) && (rows == image->rows))
2846     return(CloneImage(image,0,0,MagickTrue,exception));
2847   scale_image=CloneImage(image,columns,rows,MagickTrue,exception);
2848   if (scale_image == (Image *) NULL)
2849     return((Image *) NULL);
2850   if (SetImageStorageClass(scale_image,DirectClass,exception) == MagickFalse)
2851     {
2852       scale_image=DestroyImage(scale_image);
2853       return((Image *) NULL);
2854     }
2855   /*
2856     Allocate memory.
2857   */
2858   x_vector=(MagickRealType *) AcquireQuantumMemory((size_t) image->columns,
2859     GetPixelChannels(image)*sizeof(*x_vector));
2860   scanline=x_vector;
2861   if (image->rows != scale_image->rows)
2862     scanline=(MagickRealType *) AcquireQuantumMemory((size_t) image->columns,
2863       GetPixelChannels(image)*sizeof(*scanline));
2864   scale_scanline=(MagickRealType *) AcquireQuantumMemory((size_t)
2865     scale_image->columns,GetPixelChannels(scale_image)*sizeof(*scale_scanline));
2866   y_vector=(MagickRealType *) AcquireQuantumMemory((size_t) image->columns,
2867     GetPixelChannels(image)*sizeof(*y_vector));
2868   if ((scanline == (MagickRealType *) NULL) ||
2869       (scale_scanline == (MagickRealType *) NULL) ||
2870       (x_vector == (MagickRealType *) NULL) ||
2871       (y_vector == (MagickRealType *) NULL))
2872     {
2873       scale_image=DestroyImage(scale_image);
2874       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2875     }
2876   /*
2877     Scale image.
2878   */
2879   number_rows=0;
2880   next_row=MagickTrue;
2881   span.y=1.0;
2882   scale.y=(double) scale_image->rows/(double) image->rows;
2883   for (i=0; i < (ssize_t) (GetPixelChannels(image)*image->columns); i++)
2884     y_vector[i]=0.0;
2885   n=0;
2886   image_view=AcquireCacheView(image);
2887   scale_view=AcquireCacheView(scale_image);
2888   for (y=0; y < (ssize_t) scale_image->rows; y++)
2889   {
2890     register const Quantum
2891       *restrict p;
2892
2893     register Quantum
2894       *restrict q;
2895
2896     register ssize_t
2897       x;
2898
2899     q=QueueCacheViewAuthenticPixels(scale_view,0,y,scale_image->columns,1,
2900       exception);
2901     if (q == (Quantum *) NULL)
2902       break;
2903     alpha=1.0;
2904     if (scale_image->rows == image->rows)
2905       {
2906         /*
2907           Read a new scanline.
2908         */
2909         p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
2910           exception);
2911         if (p == (const Quantum *) NULL)
2912           break;
2913         for (x=0; x < (ssize_t) image->columns; x++)
2914         {
2915           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2916           {
2917             traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
2918             if ((traits & BlendPixelTrait) == 0)
2919               {
2920                 x_vector[x*GetPixelChannels(image)+i]=(MagickRealType) p[i];
2921                 continue;
2922               }
2923             alpha=QuantumScale*GetPixelAlpha(image,p);
2924             x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
2925           }
2926           p+=GetPixelChannels(image);
2927         }
2928       }
2929     else
2930       {
2931         /*
2932           Scale Y direction.
2933         */
2934         while (scale.y < span.y)
2935         {
2936           if ((next_row != MagickFalse) &&
2937               (number_rows < (ssize_t) image->rows))
2938             {
2939               /*
2940                 Read a new scanline.
2941               */
2942               p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
2943                 exception);
2944               if (p == (const Quantum *) NULL)
2945                 break;
2946               for (x=0; x < (ssize_t) image->columns; x++)
2947               {
2948                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2949                 {
2950                   traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
2951                   if ((traits & BlendPixelTrait) == 0)
2952                     {
2953                       x_vector[x*GetPixelChannels(image)+i]=(MagickRealType)
2954                         p[i];
2955                       continue;
2956                     }
2957                   alpha=QuantumScale*GetPixelAlpha(image,p);
2958                   x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
2959                 }
2960                 p+=GetPixelChannels(image);
2961               }
2962               number_rows++;
2963             }
2964           for (x=0; x < (ssize_t) image->columns; x++)
2965             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2966               y_vector[x*GetPixelChannels(image)+i]+=scale.y*
2967                 x_vector[x*GetPixelChannels(image)+i];
2968           span.y-=scale.y;
2969           scale.y=(double) scale_image->rows/(double) image->rows;
2970           next_row=MagickTrue;
2971         }
2972         if ((next_row != MagickFalse) && (number_rows < (ssize_t) image->rows))
2973           {
2974             /*
2975               Read a new scanline.
2976             */
2977             p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
2978               exception);
2979             if (p == (const Quantum *) NULL)
2980               break;
2981             for (x=0; x < (ssize_t) image->columns; x++)
2982             {
2983               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2984               {
2985                 traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
2986                 if ((traits & BlendPixelTrait) == 0)
2987                   {
2988                     x_vector[x*GetPixelChannels(image)+i]=(MagickRealType) p[i];
2989                     continue;
2990                   }
2991                 alpha=QuantumScale*GetPixelAlpha(image,p);
2992                 x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
2993               }
2994               p+=GetPixelChannels(image);
2995             }
2996             number_rows++;
2997             next_row=MagickFalse;
2998           }
2999         for (x=0; x < (ssize_t) image->columns; x++)
3000         {
3001           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3002           {
3003             pixel[i]=y_vector[x*GetPixelChannels(image)+i]+span.y*
3004               x_vector[x*GetPixelChannels(image)+i];
3005             scanline[x*GetPixelChannels(image)+i]=pixel[i];
3006             y_vector[x*GetPixelChannels(image)+i]=0.0;
3007           }
3008         }
3009         scale.y-=span.y;
3010         if (scale.y <= 0)
3011           {
3012             scale.y=(double) scale_image->rows/(double) image->rows;
3013             next_row=MagickTrue;
3014           }
3015         span.y=1.0;
3016       }
3017     if (scale_image->columns == image->columns)
3018       {
3019         /*
3020           Transfer scanline to scaled image.
3021         */
3022         for (x=0; x < (ssize_t) scale_image->columns; x++)
3023         {
3024           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3025           {
3026             traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
3027             channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
3028             scale_traits=GetPixelChannelMapTraits(scale_image,channel);
3029             if ((traits == UndefinedPixelTrait) ||
3030                 (scale_traits == UndefinedPixelTrait))
3031               continue;
3032             if ((scale_traits & BlendPixelTrait) == 0)
3033               {
3034                 SetPixelChannel(scale_image,channel,ClampToQuantum(scanline[
3035                   x*GetPixelChannels(image)+i]),q);
3036                 continue;
3037               }
3038             alpha=QuantumScale*scanline[x*GetPixelChannels(image)+
3039               GetPixelChannelMapChannel(image,AlphaPixelChannel)];
3040             gamma=1.0/(fabs((double) alpha) <= MagickEpsilon ? 1.0 : alpha);
3041             SetPixelChannel(scale_image,channel,ClampToQuantum(gamma*scanline[
3042               x*GetPixelChannels(image)+i]),q);
3043           }
3044           q+=GetPixelChannels(scale_image);
3045         }
3046       }
3047     else
3048       {
3049         ssize_t
3050           n;
3051
3052         /*
3053           Scale X direction.
3054         */
3055         next_column=MagickFalse;
3056         n=0;
3057         span.x=1.0;
3058         for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3059           pixel[i]=0.0;
3060         for (x=0; x < (ssize_t) image->columns; x++)
3061         {
3062           scale.x=(double) scale_image->columns/(double) image->columns;
3063           while (scale.x >= span.x)
3064           {
3065             if (next_column != MagickFalse)
3066               {
3067                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3068                   pixel[i]=0.0;
3069                 n++;
3070               }
3071             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3072             {
3073               traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
3074               if (traits == UndefinedPixelTrait)
3075                 continue;
3076               channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
3077               pixel[i]+=span.x*scanline[x*GetPixelChannels(image)+i];
3078               scale_scanline[n*GetPixelChannels(scale_image)+channel]=pixel[i];
3079             }
3080             scale.x-=span.x;
3081             span.x=1.0;
3082             next_column=MagickTrue;
3083           }
3084         if (scale.x > 0)
3085           {
3086             if (next_column != MagickFalse)
3087               {
3088                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3089                   pixel[i]=0.0;
3090                 n++;
3091                 next_column=MagickFalse;
3092               }
3093             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3094               pixel[i]+=scale.x*scanline[x*GetPixelChannels(image)+i];
3095             span.x-=scale.x;
3096           }
3097       }
3098       if (span.x > 0)
3099         {
3100           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3101             pixel[i]+=span.x*scanline[(x-1)*GetPixelChannels(image)+i];
3102         }
3103       if ((next_column == MagickFalse) &&
3104           ((ssize_t) n < (ssize_t) scale_image->columns))
3105         for (i=0; i < (ssize_t) GetPixelChannels(scale_image); i++)
3106           scale_scanline[n*GetPixelChannels(scale_image)+i]=pixel[i];
3107       /*
3108         Transfer scanline to scaled image.
3109       */
3110       for (x=0; x < (ssize_t) scale_image->columns; x++)
3111       {
3112         for (i=0; i < (ssize_t) GetPixelChannels(scale_image); i++)
3113         {
3114           traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
3115           channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
3116           scale_traits=GetPixelChannelMapTraits(scale_image,channel);
3117           if ((traits == UndefinedPixelTrait) ||
3118               (scale_traits == UndefinedPixelTrait))
3119             continue;
3120           if ((scale_traits & BlendPixelTrait) == 0)
3121             {
3122               SetPixelChannel(scale_image,channel,ClampToQuantum(
3123                 scale_scanline[x*GetPixelChannels(scale_image)+channel]),q);
3124               continue;
3125             }
3126           alpha=QuantumScale*scanline[x*GetPixelChannels(image)+
3127             GetPixelChannelMapChannel(image,AlphaPixelChannel)];
3128           gamma=1.0/(fabs((double) alpha) <= MagickEpsilon ? 1.0 : alpha);
3129           SetPixelChannel(scale_image,channel,ClampToQuantum(gamma*
3130             scale_scanline[x*GetPixelChannels(scale_image)+channel]),q);
3131         }
3132         q+=GetPixelChannels(scale_image);
3133       }
3134     }
3135     if (SyncCacheViewAuthenticPixels(scale_view,exception) == MagickFalse)
3136       break;
3137     proceed=SetImageProgress(image,ScaleImageTag,(MagickOffsetType) y,
3138       image->rows);
3139     if (proceed == MagickFalse)
3140       break;
3141   }
3142   scale_view=DestroyCacheView(scale_view);
3143   image_view=DestroyCacheView(image_view);
3144   /*
3145     Free allocated memory.
3146   */
3147   y_vector=(MagickRealType *) RelinquishMagickMemory(y_vector);
3148   scale_scanline=(MagickRealType *) RelinquishMagickMemory(scale_scanline);
3149   if (scale_image->rows != image->rows)
3150     scanline=(MagickRealType *) RelinquishMagickMemory(scanline);
3151   x_vector=(MagickRealType *) RelinquishMagickMemory(x_vector);
3152   scale_image->type=image->type;
3153   return(scale_image);
3154 }
3155 \f
3156 /*
3157 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3158 %                                                                             %
3159 %                                                                             %
3160 %                                                                             %
3161 %   T h u m b n a i l I m a g e                                               %
3162 %                                                                             %
3163 %                                                                             %
3164 %                                                                             %
3165 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3166 %
3167 %  ThumbnailImage() changes the size of an image to the given dimensions and
3168 %  removes any associated profiles.  The goal is to produce small low cost
3169 %  thumbnail images suited for display on the Web.
3170 %
3171 %  The format of the ThumbnailImage method is:
3172 %
3173 %      Image *ThumbnailImage(const Image *image,const size_t columns,
3174 %        const size_t rows,ExceptionInfo *exception)
3175 %
3176 %  A description of each parameter follows:
3177 %
3178 %    o image: the image.
3179 %
3180 %    o columns: the number of columns in the scaled image.
3181 %
3182 %    o rows: the number of rows in the scaled image.
3183 %
3184 %    o exception: return any errors or warnings in this structure.
3185 %
3186 */
3187 MagickExport Image *ThumbnailImage(const Image *image,const size_t columns,
3188   const size_t rows,ExceptionInfo *exception)
3189 {
3190 #define SampleFactor  5
3191
3192   char
3193     value[MaxTextExtent];
3194
3195   const char
3196     *name;
3197
3198   Image
3199     *thumbnail_image;
3200
3201   MagickRealType
3202     x_factor,
3203     y_factor;
3204
3205   size_t
3206     version;
3207
3208   struct stat
3209     attributes;
3210
3211   assert(image != (Image *) NULL);
3212   assert(image->signature == MagickSignature);
3213   if (image->debug != MagickFalse)
3214     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3215   assert(exception != (ExceptionInfo *) NULL);
3216   assert(exception->signature == MagickSignature);
3217   x_factor=(MagickRealType) columns/(MagickRealType) image->columns;
3218   y_factor=(MagickRealType) rows/(MagickRealType) image->rows;
3219   if ((x_factor*y_factor) > 0.1)
3220     thumbnail_image=ResizeImage(image,columns,rows,image->filter,image->blur,
3221       exception);
3222   else
3223     if (((SampleFactor*columns) < 128) || ((SampleFactor*rows) < 128))
3224       thumbnail_image=ResizeImage(image,columns,rows,image->filter,
3225         image->blur,exception);
3226     else
3227       {
3228         Image
3229           *sample_image;
3230
3231         sample_image=SampleImage(image,SampleFactor*columns,SampleFactor*rows,
3232           exception);
3233         if (sample_image == (Image *) NULL)
3234           return((Image *) NULL);
3235         thumbnail_image=ResizeImage(sample_image,columns,rows,image->filter,
3236           image->blur,exception);
3237         sample_image=DestroyImage(sample_image);
3238       }
3239   if (thumbnail_image == (Image *) NULL)
3240     return(thumbnail_image);
3241   (void) ParseAbsoluteGeometry("0x0+0+0",&thumbnail_image->page);
3242   if (thumbnail_image->matte == MagickFalse)
3243     (void) SetImageAlphaChannel(thumbnail_image,OpaqueAlphaChannel,exception);
3244   thumbnail_image->depth=8;
3245   thumbnail_image->interlace=NoInterlace;
3246   /*
3247     Strip all profiles except color profiles.
3248   */
3249   ResetImageProfileIterator(thumbnail_image);
3250   for (name=GetNextImageProfile(thumbnail_image); name != (const char *) NULL; )
3251   {
3252     if ((LocaleCompare(name,"icc") != 0) && (LocaleCompare(name,"icm") != 0))
3253      {
3254        (void) DeleteImageProfile(thumbnail_image,name);
3255        ResetImageProfileIterator(thumbnail_image);
3256      }
3257     name=GetNextImageProfile(thumbnail_image);
3258   }
3259   (void) DeleteImageProperty(thumbnail_image,"comment");
3260   (void) CopyMagickString(value,image->magick_filename,MaxTextExtent);
3261   if (strstr(image->magick_filename,"//") == (char *) NULL)
3262     (void) FormatLocaleString(value,MaxTextExtent,"file://%s",
3263       image->magick_filename);
3264   (void) SetImageProperty(thumbnail_image,"Thumb::URI",value,exception);
3265   (void) CopyMagickString(value,image->magick_filename,MaxTextExtent);
3266   if (GetPathAttributes(image->filename,&attributes) != MagickFalse)
3267     {
3268       (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3269         attributes.st_mtime);
3270       (void) SetImageProperty(thumbnail_image,"Thumb::MTime",value,exception);
3271     }
3272   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3273     attributes.st_mtime);
3274   (void) FormatMagickSize(GetBlobSize(image),MagickFalse,value);
3275   (void) ConcatenateMagickString(value,"B",MaxTextExtent);
3276   (void) SetImageProperty(thumbnail_image,"Thumb::Size",value,exception);
3277   (void) FormatLocaleString(value,MaxTextExtent,"image/%s",image->magick);
3278   LocaleLower(value);
3279   (void) SetImageProperty(thumbnail_image,"Thumb::Mimetype",value,exception);
3280   (void) SetImageProperty(thumbnail_image,"software",GetMagickVersion(&version),
3281     exception);
3282   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3283     image->magick_columns);
3284   (void) SetImageProperty(thumbnail_image,"Thumb::Image::Width",value,
3285     exception);
3286   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3287     image->magick_rows);
3288   (void) SetImageProperty(thumbnail_image,"Thumb::Image::height",value,
3289     exception);
3290   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3291     GetImageListLength(image));
3292   (void) SetImageProperty(thumbnail_image,"Thumb::Document::Pages",value,
3293     exception);
3294   return(thumbnail_image);
3295 }