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