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