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