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