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