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

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