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

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