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-2013 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,image->columns,image->rows,
1865     GetPixelChannels(image),LQR_COLDEPTH_32F);
1866   if (carver == (LqrCarver *) NULL)
1867     {
1868       pixel_info=RelinquishVirtualMemory(pixel_info);
1869       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1870     }
1871   lqr_status=lqr_carver_init(carver,(int) delta_x,rigidity);
1872   lqr_status=lqr_carver_resize(carver,columns,rows);
1873   (void) lqr_status;
1874   rescale_image=CloneImage(image,lqr_carver_get_width(carver),
1875     lqr_carver_get_height(carver),MagickTrue,exception);
1876   if (rescale_image == (Image *) NULL)
1877     {
1878       pixel_info=RelinquishVirtualMemory(pixel_info);
1879       return((Image *) NULL);
1880     }
1881   if (SetImageStorageClass(rescale_image,DirectClass,exception) == MagickFalse)
1882     {
1883       pixel_info=RelinquishVirtualMemory(pixel_info);
1884       rescale_image=DestroyImage(rescale_image);
1885       return((Image *) NULL);
1886     }
1887   rescale_view=AcquireAuthenticCacheView(rescale_image,exception);
1888   (void) lqr_carver_scan_reset(carver);
1889   while (lqr_carver_scan_ext(carver,&x_offset,&y_offset,(void **) &packet) != 0)
1890   {
1891     register Quantum
1892       *restrict q;
1893
1894     register ssize_t
1895       i;
1896
1897     q=QueueCacheViewAuthenticPixels(rescale_view,x_offset,y_offset,1,1,
1898       exception);
1899     if (q == (Quantum *) NULL)
1900       break;
1901     for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1902     {
1903       PixelChannel
1904         channel;
1905
1906       PixelTrait
1907         rescale_traits,
1908         traits;
1909
1910       channel=GetPixelChannelChannel(image,i);
1911       traits=GetPixelChannelTraits(image,channel);
1912       rescale_traits=GetPixelChannelTraits(rescale_image,channel);
1913       if ((traits == UndefinedPixelTrait) ||
1914           (rescale_traits == UndefinedPixelTrait))
1915         continue;
1916       SetPixelChannel(rescale_image,channel,ClampToQuantum(QuantumRange*
1917         packet[i]),q);
1918     }
1919     if (SyncCacheViewAuthenticPixels(rescale_view,exception) == MagickFalse)
1920       break;
1921   }
1922   rescale_view=DestroyCacheView(rescale_view);
1923   pixel_info=RelinquishVirtualMemory(pixel_info);
1924   lqr_carver_destroy(carver);
1925   return(rescale_image);
1926 }
1927 #else
1928 MagickExport Image *LiquidRescaleImage(const Image *image,
1929   const size_t magick_unused(columns),const size_t magick_unused(rows),
1930   const double magick_unused(delta_x),const double magick_unused(rigidity),
1931   ExceptionInfo *exception)
1932 {
1933   assert(image != (const Image *) NULL);
1934   assert(image->signature == MagickSignature);
1935   if (image->debug != MagickFalse)
1936     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1937   assert(exception != (ExceptionInfo *) NULL);
1938   assert(exception->signature == MagickSignature);
1939   (void) ThrowMagickException(exception,GetMagickModule(),MissingDelegateError,
1940     "DelegateLibrarySupportNotBuiltIn","'%s' (LQR)",image->filename);
1941   return((Image *) NULL);
1942 }
1943 #endif
1944 \f
1945 /*
1946 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1947 %                                                                             %
1948 %                                                                             %
1949 %                                                                             %
1950 %   M a g n i f y I m a g e                                                   %
1951 %                                                                             %
1952 %                                                                             %
1953 %                                                                             %
1954 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1955 %
1956 %  MagnifyImage() doubles the size of the image with a pixel art scaling
1957 %  algorithm.
1958 %
1959 %  The format of the MagnifyImage method is:
1960 %
1961 %      Image *MagnifyImage(const Image *image,ExceptionInfo *exception)
1962 %
1963 %  A description of each parameter follows:
1964 %
1965 %    o image: the image.
1966 %
1967 %    o exception: return any errors or warnings in this structure.
1968 %
1969 */
1970 MagickExport Image *MagnifyImage(const Image *image,ExceptionInfo *exception)
1971 {
1972 #define MagnifyImageTag  "Magnify/Image"
1973
1974   CacheView
1975     *image_view,
1976     *magnify_view;
1977
1978   Image
1979     *magnify_image;
1980
1981   MagickBooleanType
1982     status;
1983
1984   MagickOffsetType
1985     progress;
1986
1987   ssize_t
1988     y;
1989
1990   /*
1991     Initialize magnified image attributes.
1992   */
1993   assert(image != (const Image *) NULL);
1994   assert(image->signature == MagickSignature);
1995   if (image->debug != MagickFalse)
1996     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1997   assert(exception != (ExceptionInfo *) NULL);
1998   assert(exception->signature == MagickSignature);
1999   magnify_image=CloneImage(image,2*image->columns,2*image->rows,MagickTrue,
2000     exception);
2001   if (magnify_image == (Image *) NULL)
2002     return((Image *) NULL);
2003   /*
2004     Magnify image.
2005   */
2006   status=MagickTrue;
2007   progress=0;
2008   image_view=AcquireVirtualCacheView(image,exception);
2009   magnify_view=AcquireAuthenticCacheView(magnify_image,exception);
2010 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2011   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2012     magick_threads(image,magnify_image,image->rows,1)
2013 #endif
2014   for (y=0; y < (ssize_t) image->rows; y++)
2015   {
2016     register Quantum
2017       *restrict q;
2018
2019     register ssize_t
2020       x;
2021
2022     if (status == MagickFalse)
2023       continue;
2024     q=QueueCacheViewAuthenticPixels(magnify_view,0,2*y,magnify_image->columns,2,
2025       exception);
2026     if (q == (Quantum *) NULL)
2027       {
2028         status=MagickFalse;
2029         continue;
2030       }
2031     /*
2032       Magnify this row of pixels.
2033     */
2034     for (x=0; x < (ssize_t) image->columns; x++)
2035     {
2036       MagickRealType
2037         intensity[9];
2038
2039       register const Quantum
2040         *restrict p;
2041
2042       register Quantum
2043         *restrict r;
2044
2045       register ssize_t
2046         i;
2047
2048       size_t
2049         channels;
2050
2051       p=GetCacheViewVirtualPixels(image_view,x-1,y-1,3,3,exception);
2052       if (p == (const Quantum *) NULL)
2053         {
2054           status=MagickFalse;
2055           continue;
2056         }
2057       channels=GetPixelChannels(image);
2058       for (i=0; i < 9; i++)
2059         intensity[i]=GetPixelIntensity(image,p+i*channels);
2060       r=q;
2061       if ((fabs(intensity[1]-intensity[7]) < MagickEpsilon) ||
2062           (fabs(intensity[3]-intensity[5]) < MagickEpsilon))
2063         {
2064           /*
2065             Clone center pixel.
2066           */
2067           for (i=0; i < (ssize_t) channels; i++)
2068             r[i]=p[4*channels+i];
2069           r+=GetPixelChannels(magnify_image);
2070           for (i=0; i < (ssize_t) channels; i++)
2071             r[i]=p[4*channels+i];
2072           r+=(magnify_image->columns-1)*GetPixelChannels(magnify_image);
2073           for (i=0; i < (ssize_t) channels; i++)
2074             r[i]=p[4*channels+i];
2075           r+=GetPixelChannels(magnify_image);
2076           for (i=0; i < (ssize_t) channels; i++)
2077             r[i]=p[4*channels+i];
2078         }
2079       else
2080         {
2081           /*
2082             Selectively clone pixel.
2083           */
2084           if (fabs(intensity[1]-intensity[3]) < MagickEpsilon)
2085             for (i=0; i < (ssize_t) channels; i++)
2086               r[i]=p[3*channels+i];
2087           else
2088             for (i=0; i < (ssize_t) channels; i++)
2089               r[i]=p[4*channels+i];
2090           r+=GetPixelChannels(magnify_image);
2091           if (fabs(intensity[1]-intensity[5]) < MagickEpsilon)
2092             for (i=0; i < (ssize_t) channels; i++)
2093               r[i]=p[5*channels+i];
2094           else
2095             for (i=0; i < (ssize_t) channels; i++)
2096               r[i]=p[4*channels+i];
2097           r+=(magnify_image->columns-1)*GetPixelChannels(magnify_image);
2098           if (fabs(intensity[3]-intensity[7]) < MagickEpsilon)
2099             for (i=0; i < (ssize_t) channels; i++)
2100               r[i]=p[3*channels+i];
2101           else
2102             for (i=0; i < (ssize_t) channels; i++)
2103               r[i]=p[4*channels+i];
2104           r+=GetPixelChannels(magnify_image);
2105           if (fabs(intensity[5]-intensity[7]) < MagickEpsilon)
2106             for (i=0; i < (ssize_t) channels; i++)
2107               r[i]=p[5*channels+i];
2108           else
2109             for (i=0; i < (ssize_t) channels; i++)
2110               r[i]=p[4*channels+i];
2111         }
2112       q+=2*GetPixelChannels(magnify_image);
2113     }
2114     if (SyncCacheViewAuthenticPixels(magnify_view,exception) == MagickFalse)
2115       status=MagickFalse;
2116     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2117       {
2118         MagickBooleanType
2119           proceed;
2120
2121 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2122         #pragma omp critical (MagickCore_MagnifyImage)
2123 #endif
2124         proceed=SetImageProgress(image,MagnifyImageTag,progress++,image->rows);
2125         if (proceed == MagickFalse)
2126           status=MagickFalse;
2127       }
2128   }
2129   magnify_view=DestroyCacheView(magnify_view);
2130   image_view=DestroyCacheView(image_view);
2131   if (status == MagickFalse)
2132     magnify_image=DestroyImage(magnify_image);
2133   return(magnify_image);
2134 }
2135 \f
2136 /*
2137 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2138 %                                                                             %
2139 %                                                                             %
2140 %                                                                             %
2141 %   M i n i f y I m a g e                                                     %
2142 %                                                                             %
2143 %                                                                             %
2144 %                                                                             %
2145 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2146 %
2147 %  MinifyImage() is a convenience method that scales an image proportionally to
2148 %  half its size.
2149 %
2150 %  The format of the MinifyImage method is:
2151 %
2152 %      Image *MinifyImage(const Image *image,ExceptionInfo *exception)
2153 %
2154 %  A description of each parameter follows:
2155 %
2156 %    o image: the image.
2157 %
2158 %    o exception: return any errors or warnings in this structure.
2159 %
2160 */
2161 MagickExport Image *MinifyImage(const Image *image,ExceptionInfo *exception)
2162 {
2163   Image
2164     *minify_image;
2165
2166   assert(image != (Image *) NULL);
2167   assert(image->signature == MagickSignature);
2168   if (image->debug != MagickFalse)
2169     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2170   assert(exception != (ExceptionInfo *) NULL);
2171   assert(exception->signature == MagickSignature);
2172   minify_image=ResizeImage(image,image->columns/2,image->rows/2,SplineFilter,
2173     exception);
2174   return(minify_image);
2175 }
2176 \f
2177 /*
2178 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2179 %                                                                             %
2180 %                                                                             %
2181 %                                                                             %
2182 %   R e s a m p l e I m a g e                                                 %
2183 %                                                                             %
2184 %                                                                             %
2185 %                                                                             %
2186 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2187 %
2188 %  ResampleImage() resize image in terms of its pixel size, so that when
2189 %  displayed at the given resolution it will be the same size in terms of
2190 %  real world units as the original image at the original resolution.
2191 %
2192 %  The format of the ResampleImage method is:
2193 %
2194 %      Image *ResampleImage(Image *image,const double x_resolution,
2195 %        const double y_resolution,const FilterTypes filter,
2196 %        ExceptionInfo *exception)
2197 %
2198 %  A description of each parameter follows:
2199 %
2200 %    o image: the image to be resized to fit the given resolution.
2201 %
2202 %    o x_resolution: the new image x resolution.
2203 %
2204 %    o y_resolution: the new image y resolution.
2205 %
2206 %    o filter: Image filter to use.
2207 %
2208 %    o exception: return any errors or warnings in this structure.
2209 %
2210 */
2211 MagickExport Image *ResampleImage(const Image *image,const double x_resolution,
2212   const double y_resolution,const FilterTypes filter,ExceptionInfo *exception)
2213 {
2214 #define ResampleImageTag  "Resample/Image"
2215
2216   Image
2217     *resample_image;
2218
2219   size_t
2220     height,
2221     width;
2222
2223   /*
2224     Initialize sampled image attributes.
2225   */
2226   assert(image != (const Image *) NULL);
2227   assert(image->signature == MagickSignature);
2228   if (image->debug != MagickFalse)
2229     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2230   assert(exception != (ExceptionInfo *) NULL);
2231   assert(exception->signature == MagickSignature);
2232   width=(size_t) (x_resolution*image->columns/(image->resolution.x == 0.0 ?
2233     72.0 : image->resolution.x)+0.5);
2234   height=(size_t) (y_resolution*image->rows/(image->resolution.y == 0.0 ?
2235     72.0 : image->resolution.y)+0.5);
2236   resample_image=ResizeImage(image,width,height,filter,exception);
2237   if (resample_image != (Image *) NULL)
2238     {
2239       resample_image->resolution.x=x_resolution;
2240       resample_image->resolution.y=y_resolution;
2241     }
2242   return(resample_image);
2243 }
2244 \f
2245 /*
2246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2247 %                                                                             %
2248 %                                                                             %
2249 %                                                                             %
2250 %   R e s i z e I m a g e                                                     %
2251 %                                                                             %
2252 %                                                                             %
2253 %                                                                             %
2254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2255 %
2256 %  ResizeImage() scales an image to the desired dimensions, using the given
2257 %  filter (see AcquireFilterInfo()).
2258 %
2259 %  If an undefined filter is given the filter defaults to Mitchell for a
2260 %  colormapped image, a image with a matte channel, or if the image is
2261 %  enlarged.  Otherwise the filter defaults to a Lanczos.
2262 %
2263 %  ResizeImage() was inspired by Paul Heckbert's "zoom" program.
2264 %
2265 %  The format of the ResizeImage method is:
2266 %
2267 %      Image *ResizeImage(Image *image,const size_t columns,
2268 %        const size_t rows,const FilterTypes filter,const double blur,
2269 %        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 blur: the blur factor where > 1 is blurry, < 1 is sharp.  Typically set
2282 %      this to 1.0.
2283 %
2284 %    o exception: return any errors or warnings in this structure.
2285 %
2286 */
2287
2288 typedef struct _ContributionInfo
2289 {
2290   double
2291     weight;
2292
2293   ssize_t
2294     pixel;
2295 } ContributionInfo;
2296
2297 static ContributionInfo **DestroyContributionThreadSet(
2298   ContributionInfo **contribution)
2299 {
2300   register ssize_t
2301     i;
2302
2303   assert(contribution != (ContributionInfo **) NULL);
2304   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2305     if (contribution[i] != (ContributionInfo *) NULL)
2306       contribution[i]=(ContributionInfo *) RelinquishAlignedMemory(
2307         contribution[i]);
2308   contribution=(ContributionInfo **) RelinquishMagickMemory(contribution);
2309   return(contribution);
2310 }
2311
2312 static ContributionInfo **AcquireContributionThreadSet(const size_t count)
2313 {
2314   register ssize_t
2315     i;
2316
2317   ContributionInfo
2318     **contribution;
2319
2320   size_t
2321     number_threads;
2322
2323   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2324   contribution=(ContributionInfo **) AcquireQuantumMemory(number_threads,
2325     sizeof(*contribution));
2326   if (contribution == (ContributionInfo **) NULL)
2327     return((ContributionInfo **) NULL);
2328   (void) ResetMagickMemory(contribution,0,number_threads*sizeof(*contribution));
2329   for (i=0; i < (ssize_t) number_threads; i++)
2330   {
2331     contribution[i]=(ContributionInfo *) MagickAssumeAligned(
2332       AcquireAlignedMemory(count,sizeof(**contribution)));
2333     if (contribution[i] == (ContributionInfo *) NULL)
2334       return(DestroyContributionThreadSet(contribution));
2335   }
2336   return(contribution);
2337 }
2338
2339 static inline double MagickMax(const double x,const double y)
2340 {
2341   if (x > y)
2342     return(x);
2343   return(y);
2344 }
2345
2346 static inline double MagickMin(const double x,const double y)
2347 {
2348   if (x < y)
2349     return(x);
2350   return(y);
2351 }
2352
2353 static MagickBooleanType HorizontalFilter(const ResizeFilter *resize_filter,
2354   const Image *image,Image *resize_image,const double x_factor,
2355   const MagickSizeType span,MagickOffsetType *offset,ExceptionInfo *exception)
2356 {
2357 #define ResizeImageTag  "Resize/Image"
2358
2359   CacheView
2360     *image_view,
2361     *resize_view;
2362
2363   ClassType
2364     storage_class;
2365
2366   ContributionInfo
2367     **restrict contributions;
2368
2369   MagickBooleanType
2370     status;
2371
2372   double
2373     scale,
2374     support;
2375
2376   ssize_t
2377     x;
2378
2379   /*
2380     Apply filter to resize horizontally from image to resize image.
2381   */
2382   scale=MagickMax(1.0/x_factor+MagickEpsilon,1.0);
2383   support=scale*GetResizeFilterSupport(resize_filter);
2384   storage_class=support > 0.5 ? DirectClass : image->storage_class;
2385   if (SetImageStorageClass(resize_image,storage_class,exception) == MagickFalse)
2386     return(MagickFalse);
2387   if (support < 0.5)
2388     {
2389       /*
2390         Support too small even for nearest neighbour: Reduce to point sampling.
2391       */
2392       support=(double) 0.5;
2393       scale=1.0;
2394     }
2395   contributions=AcquireContributionThreadSet((size_t) (2.0*support+3.0));
2396   if (contributions == (ContributionInfo **) NULL)
2397     {
2398       (void) ThrowMagickException(exception,GetMagickModule(),
2399         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2400       return(MagickFalse);
2401     }
2402   status=MagickTrue;
2403   scale=PerceptibleReciprocal(scale);
2404   image_view=AcquireVirtualCacheView(image,exception);
2405   resize_view=AcquireAuthenticCacheView(resize_image,exception);
2406 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2407   #pragma omp parallel for schedule(static,4) shared(status) \
2408     magick_threads(image,resize_image,resize_image->columns,1)
2409 #endif
2410   for (x=0; x < (ssize_t) resize_image->columns; x++)
2411   {
2412     double
2413       bisect,
2414       density;
2415
2416     register const Quantum
2417       *restrict p;
2418
2419     register ContributionInfo
2420       *restrict contribution;
2421
2422     register Quantum
2423       *restrict q;
2424
2425     register ssize_t
2426       y;
2427
2428     ssize_t
2429       n,
2430       start,
2431       stop;
2432
2433     if (status == MagickFalse)
2434       continue;
2435     bisect=(double) (x+0.5)/x_factor+MagickEpsilon;
2436     start=(ssize_t) MagickMax(bisect-support+0.5,0.0);
2437     stop=(ssize_t) MagickMin(bisect+support+0.5,(double) image->columns);
2438     density=0.0;
2439     contribution=contributions[GetOpenMPThreadId()];
2440     for (n=0; n < (stop-start); n++)
2441     {
2442       contribution[n].pixel=start+n;
2443       contribution[n].weight=GetResizeFilterWeight(resize_filter,scale*
2444         ((double) (start+n)-bisect+0.5));
2445       density+=contribution[n].weight;
2446     }
2447     if ((density != 0.0) && (density != 1.0))
2448       {
2449         register ssize_t
2450           i;
2451
2452         /*
2453           Normalize.
2454         */
2455         density=PerceptibleReciprocal(density);
2456         for (i=0; i < n; i++)
2457           contribution[i].weight*=density;
2458       }
2459     p=GetCacheViewVirtualPixels(image_view,contribution[0].pixel,0,(size_t)
2460       (contribution[n-1].pixel-contribution[0].pixel+1),image->rows,exception);
2461     q=QueueCacheViewAuthenticPixels(resize_view,x,0,1,resize_image->rows,
2462       exception);
2463     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2464       {
2465         status=MagickFalse;
2466         continue;
2467       }
2468     for (y=0; y < (ssize_t) resize_image->rows; y++)
2469     {
2470       register ssize_t
2471         i;
2472
2473       for (i=0; i < (ssize_t) GetPixelChannels(resize_image); i++)
2474       {
2475         double
2476           alpha,
2477           gamma,
2478           pixel;
2479
2480         PixelChannel
2481           channel;
2482
2483         PixelTrait
2484           resize_traits,
2485           traits;
2486
2487         register ssize_t
2488           j;
2489
2490         ssize_t
2491           k;
2492
2493         channel=GetPixelChannelChannel(image,i);
2494         traits=GetPixelChannelTraits(image,channel);
2495         resize_traits=GetPixelChannelTraits(resize_image,channel);
2496         if ((traits == UndefinedPixelTrait) ||
2497             (resize_traits == UndefinedPixelTrait))
2498           continue;
2499         if (((resize_traits & CopyPixelTrait) != 0) ||
2500             (GetPixelReadMask(resize_image,q) == 0))
2501           {
2502             j=(ssize_t) (MagickMin(MagickMax(bisect,(double) start),(double)
2503               stop-1.0)+0.5);
2504             k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2505               (contribution[j-start].pixel-contribution[0].pixel);
2506             SetPixelChannel(resize_image,channel,p[k*GetPixelChannels(image)+i],
2507               q);
2508             continue;
2509           }
2510         pixel=0.0;
2511         if ((resize_traits & BlendPixelTrait) == 0)
2512           {
2513             /*
2514               No alpha blending.
2515             */
2516             for (j=0; j < n; j++)
2517             {
2518               k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2519                 (contribution[j].pixel-contribution[0].pixel);
2520               alpha=contribution[j].weight;
2521               pixel+=alpha*p[k*GetPixelChannels(image)+i];
2522             }
2523             SetPixelChannel(resize_image,channel,ClampToQuantum(pixel),q);
2524             continue;
2525           }
2526         /*
2527           Alpha blending.
2528         */
2529         gamma=0.0;
2530         for (j=0; j < n; j++)
2531         {
2532           k=y*(contribution[n-1].pixel-contribution[0].pixel+1)+
2533             (contribution[j].pixel-contribution[0].pixel);
2534           alpha=contribution[j].weight*QuantumScale*
2535             GetPixelAlpha(image,p+k*GetPixelChannels(image));
2536           pixel+=alpha*p[k*GetPixelChannels(image)+i];
2537           gamma+=alpha;
2538         }
2539         gamma=PerceptibleReciprocal(gamma);
2540         SetPixelChannel(resize_image,channel,ClampToQuantum(gamma*pixel),q);
2541       }
2542       q+=GetPixelChannels(resize_image);
2543     }
2544     if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
2545       status=MagickFalse;
2546     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2547       {
2548         MagickBooleanType
2549           proceed;
2550
2551 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2552         #pragma omp critical (MagickCore_HorizontalFilter)
2553 #endif
2554         proceed=SetImageProgress(image,ResizeImageTag,(*offset)++,span);
2555         if (proceed == MagickFalse)
2556           status=MagickFalse;
2557       }
2558   }
2559   resize_view=DestroyCacheView(resize_view);
2560   image_view=DestroyCacheView(image_view);
2561   contributions=DestroyContributionThreadSet(contributions);
2562   return(status);
2563 }
2564
2565 static MagickBooleanType VerticalFilter(const ResizeFilter *resize_filter,
2566   const Image *image,Image *resize_image,const double y_factor,
2567   const MagickSizeType span,MagickOffsetType *offset,ExceptionInfo *exception)
2568 {
2569   CacheView
2570     *image_view,
2571     *resize_view;
2572
2573   ClassType
2574     storage_class;
2575
2576   ContributionInfo
2577     **restrict contributions;
2578
2579   double
2580     scale,
2581     support;
2582
2583   MagickBooleanType
2584     status;
2585
2586   ssize_t
2587     y;
2588
2589   /*
2590     Apply filter to resize vertically from image to resize image.
2591   */
2592   scale=MagickMax(1.0/y_factor+MagickEpsilon,1.0);
2593   support=scale*GetResizeFilterSupport(resize_filter);
2594   storage_class=support > 0.5 ? DirectClass : image->storage_class;
2595   if (SetImageStorageClass(resize_image,storage_class,exception) == MagickFalse)
2596     return(MagickFalse);
2597   if (support < 0.5)
2598     {
2599       /*
2600         Support too small even for nearest neighbour: Reduce to point sampling.
2601       */
2602       support=(double) 0.5;
2603       scale=1.0;
2604     }
2605   contributions=AcquireContributionThreadSet((size_t) (2.0*support+3.0));
2606   if (contributions == (ContributionInfo **) NULL)
2607     {
2608       (void) ThrowMagickException(exception,GetMagickModule(),
2609         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2610       return(MagickFalse);
2611     }
2612   status=MagickTrue;
2613   scale=PerceptibleReciprocal(scale);
2614   image_view=AcquireVirtualCacheView(image,exception);
2615   resize_view=AcquireAuthenticCacheView(resize_image,exception);
2616 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2617   #pragma omp parallel for schedule(static,4) shared(status) \
2618     magick_threads(image,resize_image,resize_image->rows,1)
2619 #endif
2620   for (y=0; y < (ssize_t) resize_image->rows; y++)
2621   {
2622     double
2623       bisect,
2624       density;
2625
2626     register const Quantum
2627       *restrict p;
2628
2629     register ContributionInfo
2630       *restrict contribution;
2631
2632     register Quantum
2633       *restrict q;
2634
2635     register ssize_t
2636       x;
2637
2638     ssize_t
2639       n,
2640       start,
2641       stop;
2642
2643     if (status == MagickFalse)
2644       continue;
2645     bisect=(double) (y+0.5)/y_factor+MagickEpsilon;
2646     start=(ssize_t) MagickMax(bisect-support+0.5,0.0);
2647     stop=(ssize_t) MagickMin(bisect+support+0.5,(double) image->rows);
2648     density=0.0;
2649     contribution=contributions[GetOpenMPThreadId()];
2650     for (n=0; n < (stop-start); n++)
2651     {
2652       contribution[n].pixel=start+n;
2653       contribution[n].weight=GetResizeFilterWeight(resize_filter,scale*
2654         ((double) (start+n)-bisect+0.5));
2655       density+=contribution[n].weight;
2656     }
2657     if ((density != 0.0) && (density != 1.0))
2658       {
2659         register ssize_t
2660           i;
2661
2662         /*
2663           Normalize.
2664         */
2665         density=PerceptibleReciprocal(density);
2666         for (i=0; i < n; i++)
2667           contribution[i].weight*=density;
2668       }
2669     p=GetCacheViewVirtualPixels(image_view,0,contribution[0].pixel,
2670       image->columns,(size_t) (contribution[n-1].pixel-contribution[0].pixel+1),
2671       exception);
2672     q=QueueCacheViewAuthenticPixels(resize_view,0,y,resize_image->columns,1,
2673       exception);
2674     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2675       {
2676         status=MagickFalse;
2677         continue;
2678       }
2679     for (x=0; x < (ssize_t) resize_image->columns; x++)
2680     {
2681       register ssize_t
2682         i;
2683
2684       for (i=0; i < (ssize_t) GetPixelChannels(resize_image); i++)
2685       {
2686         double
2687           alpha,
2688           gamma,
2689           pixel;
2690
2691         PixelChannel
2692           channel;
2693
2694         PixelTrait
2695           resize_traits,
2696           traits;
2697
2698         register ssize_t
2699           j;
2700
2701         ssize_t
2702           k;
2703
2704         channel=GetPixelChannelChannel(image,i);
2705         traits=GetPixelChannelTraits(image,channel);
2706         resize_traits=GetPixelChannelTraits(resize_image,channel);
2707         if ((traits == UndefinedPixelTrait) ||
2708             (resize_traits == UndefinedPixelTrait))
2709           continue;
2710         if (((resize_traits & CopyPixelTrait) != 0) ||
2711             (GetPixelReadMask(resize_image,q) == 0))
2712           {
2713             j=(ssize_t) (MagickMin(MagickMax(bisect,(double) start),(double)
2714               stop-1.0)+0.5);
2715             k=(ssize_t) ((contribution[j-start].pixel-contribution[0].pixel)*
2716               image->columns+x);
2717             SetPixelChannel(resize_image,channel,p[k*GetPixelChannels(image)+i],
2718               q);
2719             continue;
2720           }
2721         pixel=0.0;
2722         if ((resize_traits & BlendPixelTrait) == 0)
2723           {
2724             /*
2725               No alpha blending.
2726             */
2727             for (j=0; j < n; j++)
2728             {
2729               k=(ssize_t) ((contribution[j].pixel-contribution[0].pixel)*
2730                 image->columns+x);
2731               alpha=contribution[j].weight;
2732               pixel+=alpha*p[k*GetPixelChannels(image)+i];
2733             }
2734             SetPixelChannel(resize_image,channel,ClampToQuantum(pixel),q);
2735             continue;
2736           }
2737         gamma=0.0;
2738         for (j=0; j < n; j++)
2739         {
2740           k=(ssize_t) ((contribution[j].pixel-contribution[0].pixel)*
2741             image->columns+x);
2742           alpha=contribution[j].weight*QuantumScale*GetPixelAlpha(image,p+k*
2743             GetPixelChannels(image));
2744           pixel+=alpha*p[k*GetPixelChannels(image)+i];
2745           gamma+=alpha;
2746         }
2747         gamma=PerceptibleReciprocal(gamma);
2748         SetPixelChannel(resize_image,channel,ClampToQuantum(gamma*pixel),q);
2749       }
2750       q+=GetPixelChannels(resize_image);
2751     }
2752     if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
2753       status=MagickFalse;
2754     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2755       {
2756         MagickBooleanType
2757           proceed;
2758
2759 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2760         #pragma omp critical (MagickCore_VerticalFilter)
2761 #endif
2762         proceed=SetImageProgress(image,ResizeImageTag,(*offset)++,span);
2763         if (proceed == MagickFalse)
2764           status=MagickFalse;
2765       }
2766   }
2767   resize_view=DestroyCacheView(resize_view);
2768   image_view=DestroyCacheView(image_view);
2769   contributions=DestroyContributionThreadSet(contributions);
2770   return(status);
2771 }
2772
2773 MagickExport Image *ResizeImage(const Image *image,const size_t columns,
2774   const size_t rows,const FilterTypes filter,ExceptionInfo *exception)
2775 {
2776   double
2777     x_factor,
2778     y_factor;
2779
2780   FilterTypes
2781     filter_type;
2782
2783   Image
2784     *filter_image,
2785     *resize_image;
2786
2787   MagickOffsetType
2788     offset;
2789
2790   MagickSizeType
2791     span;
2792
2793   MagickStatusType
2794     status;
2795
2796   ResizeFilter
2797     *resize_filter;
2798
2799   /*
2800     Acquire resize image.
2801   */
2802   assert(image != (Image *) NULL);
2803   assert(image->signature == MagickSignature);
2804   if (image->debug != MagickFalse)
2805     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2806   assert(exception != (ExceptionInfo *) NULL);
2807   assert(exception->signature == MagickSignature);
2808   if ((columns == 0) || (rows == 0))
2809     ThrowImageException(ImageError,"NegativeOrZeroImageSize");
2810   if ((columns == image->columns) && (rows == image->rows) &&
2811       (filter == UndefinedFilter) && (blur == 1.0))
2812     return(CloneImage(image,0,0,MagickTrue,exception));
2813   resize_image=CloneImage(image,columns,rows,MagickTrue,exception);
2814   if (resize_image == (Image *) NULL)
2815     return(resize_image);
2816   /*
2817     Acquire resize filter.
2818   */
2819   x_factor=(double) columns/(double) image->columns;
2820   y_factor=(double) rows/(double) image->rows;
2821   if (x_factor > y_factor)
2822     filter_image=CloneImage(image,columns,image->rows,MagickTrue,exception);
2823   else
2824     filter_image=CloneImage(image,image->columns,rows,MagickTrue,exception);
2825   if (filter_image == (Image *) NULL)
2826     return(DestroyImage(resize_image));
2827   filter_type=LanczosFilter;
2828   if (filter != UndefinedFilter)
2829     filter_type=filter;
2830   else
2831     if ((x_factor == 1.0) && (y_factor == 1.0))
2832       filter_type=PointFilter;
2833     else
2834       if ((image->storage_class == PseudoClass) ||
2835           (image->alpha_trait == BlendPixelTrait) ||
2836           ((x_factor*y_factor) > 1.0))
2837         filter_type=MitchellFilter;
2838   resize_filter=AcquireResizeFilter(image,filter_type,MagickFalse,exception);
2839   /*
2840     Resize image.
2841   */
2842   offset=0;
2843   if (x_factor > y_factor)
2844     {
2845       span=(MagickSizeType) (filter_image->columns+rows);
2846       status=HorizontalFilter(resize_filter,image,filter_image,x_factor,span,
2847         &offset,exception);
2848       status&=VerticalFilter(resize_filter,filter_image,resize_image,y_factor,
2849         span,&offset,exception);
2850     }
2851   else
2852     {
2853       span=(MagickSizeType) (filter_image->rows+columns);
2854       status=VerticalFilter(resize_filter,image,filter_image,y_factor,span,
2855         &offset,exception);
2856       status&=HorizontalFilter(resize_filter,filter_image,resize_image,x_factor,
2857         span,&offset,exception);
2858     }
2859   /*
2860     Free resources.
2861   */
2862   filter_image=DestroyImage(filter_image);
2863   resize_filter=DestroyResizeFilter(resize_filter);
2864   if (status == MagickFalse)
2865     {
2866       resize_image=DestroyImage(resize_image);
2867       return((Image *) NULL);
2868     }
2869   resize_image->type=image->type;
2870   return(resize_image);
2871 }
2872 \f
2873 /*
2874 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2875 %                                                                             %
2876 %                                                                             %
2877 %                                                                             %
2878 %   S a m p l e I m a g e                                                     %
2879 %                                                                             %
2880 %                                                                             %
2881 %                                                                             %
2882 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2883 %
2884 %  SampleImage() scales an image to the desired dimensions with pixel
2885 %  sampling.  Unlike other scaling methods, this method does not introduce
2886 %  any additional color into the scaled image.
2887 %
2888 %  The format of the SampleImage method is:
2889 %
2890 %      Image *SampleImage(const Image *image,const size_t columns,
2891 %        const size_t rows,ExceptionInfo *exception)
2892 %
2893 %  A description of each parameter follows:
2894 %
2895 %    o image: the image.
2896 %
2897 %    o columns: the number of columns in the sampled image.
2898 %
2899 %    o rows: the number of rows in the sampled image.
2900 %
2901 %    o exception: return any errors or warnings in this structure.
2902 %
2903 */
2904 MagickExport Image *SampleImage(const Image *image,const size_t columns,
2905   const size_t rows,ExceptionInfo *exception)
2906 {
2907 #define SampleImageTag  "Sample/Image"
2908
2909   CacheView
2910     *image_view,
2911     *sample_view;
2912
2913   Image
2914     *sample_image;
2915
2916   MagickBooleanType
2917     status;
2918
2919   MagickOffsetType
2920     progress;
2921
2922   register ssize_t
2923     x;
2924
2925   ssize_t
2926     *x_offset,
2927     y;
2928
2929   PointInfo
2930     sample_offset;
2931
2932   /*
2933     Initialize sampled image attributes.
2934   */
2935   assert(image != (const Image *) NULL);
2936   assert(image->signature == MagickSignature);
2937   if (image->debug != MagickFalse)
2938     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2939   assert(exception != (ExceptionInfo *) NULL);
2940   assert(exception->signature == MagickSignature);
2941   if ((columns == 0) || (rows == 0))
2942     ThrowImageException(ImageError,"NegativeOrZeroImageSize");
2943   if ((columns == image->columns) && (rows == image->rows))
2944     return(CloneImage(image,0,0,MagickTrue,exception));
2945   sample_image=CloneImage(image,columns,rows,MagickTrue,exception);
2946   if (sample_image == (Image *) NULL)
2947     return((Image *) NULL);
2948   /*
2949     Set the sampling offset, default is in the mid-point of sample regions.
2950   */
2951   sample_offset.x=sample_offset.y=0.5-MagickEpsilon;
2952   {
2953     const char
2954       *value;
2955
2956     value=GetImageArtifact(image,"sample:offset");
2957     if (value != (char *) NULL)
2958       {
2959         GeometryInfo
2960           geometry_info;
2961
2962         MagickStatusType
2963           flags;
2964
2965         (void) ParseGeometry(value,&geometry_info);
2966         flags=ParseGeometry(value,&geometry_info);
2967         sample_offset.x=sample_offset.y=geometry_info.rho/100.0-MagickEpsilon;
2968         if ((flags & SigmaValue) != 0)
2969           sample_offset.y=geometry_info.sigma/100.0-MagickEpsilon;
2970       }
2971   }
2972   /*
2973     Allocate scan line buffer and column offset buffers.
2974   */
2975   x_offset=(ssize_t *) AcquireQuantumMemory((size_t) sample_image->columns,
2976     sizeof(*x_offset));
2977   if (x_offset == (ssize_t *) NULL)
2978     {
2979       sample_image=DestroyImage(sample_image);
2980       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2981     }
2982   for (x=0; x < (ssize_t) sample_image->columns; x++)
2983     x_offset[x]=(ssize_t) ((((double) x+sample_offset.x)*image->columns)/
2984       sample_image->columns);
2985   /*
2986     Sample each row.
2987   */
2988   status=MagickTrue;
2989   progress=0;
2990   image_view=AcquireVirtualCacheView(image,exception);
2991   sample_view=AcquireAuthenticCacheView(sample_image,exception);
2992 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2993   #pragma omp parallel for schedule(static,4) shared(status) \
2994     magick_threads(image,sample_image,1,1)
2995 #endif
2996   for (y=0; y < (ssize_t) sample_image->rows; y++)
2997   {
2998     register const Quantum
2999       *restrict p;
3000
3001     register Quantum
3002       *restrict q;
3003
3004     register ssize_t
3005       x;
3006
3007     ssize_t
3008       y_offset;
3009
3010     if (status == MagickFalse)
3011       continue;
3012     y_offset=(ssize_t) ((((double) y+sample_offset.y)*image->rows)/
3013       sample_image->rows);
3014     p=GetCacheViewVirtualPixels(image_view,0,y_offset,image->columns,1,
3015       exception);
3016     q=QueueCacheViewAuthenticPixels(sample_view,0,y,sample_image->columns,1,
3017       exception);
3018     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3019       {
3020         status=MagickFalse;
3021         continue;
3022       }
3023     /*
3024       Sample each column.
3025     */
3026     for (x=0; x < (ssize_t) sample_image->columns; x++)
3027     {
3028       register ssize_t
3029         i;
3030
3031       if (GetPixelReadMask(sample_image,q) == 0)
3032         {
3033           q+=GetPixelChannels(sample_image);
3034           continue;
3035         }
3036       for (i=0; i < (ssize_t) GetPixelChannels(sample_image); i++)
3037       {
3038         PixelChannel
3039           channel;
3040
3041         PixelTrait
3042           sample_traits,
3043           traits;
3044
3045         channel=GetPixelChannelChannel(image,i);
3046         traits=GetPixelChannelTraits(image,channel);
3047         sample_traits=GetPixelChannelTraits(sample_image,channel);
3048         if ((traits == UndefinedPixelTrait) ||
3049             (sample_traits == UndefinedPixelTrait))
3050           continue;
3051         SetPixelChannel(sample_image,channel,p[x_offset[x]*GetPixelChannels(
3052           image)+i],q);
3053       }
3054       q+=GetPixelChannels(sample_image);
3055     }
3056     if (SyncCacheViewAuthenticPixels(sample_view,exception) == MagickFalse)
3057       status=MagickFalse;
3058     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3059       {
3060         MagickBooleanType
3061           proceed;
3062
3063 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3064         #pragma omp critical (MagickCore_SampleImage)
3065 #endif
3066         proceed=SetImageProgress(image,SampleImageTag,progress++,image->rows);
3067         if (proceed == MagickFalse)
3068           status=MagickFalse;
3069       }
3070   }
3071   image_view=DestroyCacheView(image_view);
3072   sample_view=DestroyCacheView(sample_view);
3073   x_offset=(ssize_t *) RelinquishMagickMemory(x_offset);
3074   sample_image->type=image->type;
3075   if (status == MagickFalse)
3076     sample_image=DestroyImage(sample_image);
3077   return(sample_image);
3078 }
3079 \f
3080 /*
3081 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3082 %                                                                             %
3083 %                                                                             %
3084 %                                                                             %
3085 %   S c a l e I m a g e                                                       %
3086 %                                                                             %
3087 %                                                                             %
3088 %                                                                             %
3089 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3090 %
3091 %  ScaleImage() changes the size of an image to the given dimensions.
3092 %
3093 %  The format of the ScaleImage method is:
3094 %
3095 %      Image *ScaleImage(const Image *image,const size_t columns,
3096 %        const size_t rows,ExceptionInfo *exception)
3097 %
3098 %  A description of each parameter follows:
3099 %
3100 %    o image: the image.
3101 %
3102 %    o columns: the number of columns in the scaled image.
3103 %
3104 %    o rows: the number of rows in the scaled image.
3105 %
3106 %    o exception: return any errors or warnings in this structure.
3107 %
3108 */
3109 MagickExport Image *ScaleImage(const Image *image,const size_t columns,
3110   const size_t rows,ExceptionInfo *exception)
3111 {
3112 #define ScaleImageTag  "Scale/Image"
3113
3114   CacheView
3115     *image_view,
3116     *scale_view;
3117
3118   double
3119     alpha,
3120     gamma,
3121     pixel[CompositePixelChannel],
3122     *scale_scanline,
3123     *scanline,
3124     *x_vector,
3125     *y_vector;
3126
3127   Image
3128     *scale_image;
3129
3130   MagickBooleanType
3131     next_column,
3132     next_row,
3133     proceed,
3134     status;
3135
3136   PixelChannel
3137     channel;
3138
3139   PixelTrait
3140     scale_traits,
3141     traits;
3142
3143   PointInfo
3144     scale,
3145     span;
3146
3147   register ssize_t
3148     i;
3149
3150   ssize_t
3151     n,
3152     number_rows,
3153     y;
3154
3155   /*
3156     Initialize scaled image attributes.
3157   */
3158   assert(image != (const Image *) NULL);
3159   assert(image->signature == MagickSignature);
3160   if (image->debug != MagickFalse)
3161     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3162   assert(exception != (ExceptionInfo *) NULL);
3163   assert(exception->signature == MagickSignature);
3164   if ((columns == 0) || (rows == 0))
3165     return((Image *) NULL);
3166   if ((columns == image->columns) && (rows == image->rows))
3167     return(CloneImage(image,0,0,MagickTrue,exception));
3168   scale_image=CloneImage(image,columns,rows,MagickTrue,exception);
3169   if (scale_image == (Image *) NULL)
3170     return((Image *) NULL);
3171   if (SetImageStorageClass(scale_image,DirectClass,exception) == MagickFalse)
3172     {
3173       scale_image=DestroyImage(scale_image);
3174       return((Image *) NULL);
3175     }
3176   /*
3177     Allocate memory.
3178   */
3179   x_vector=(double *) AcquireQuantumMemory((size_t) image->columns,
3180     GetPixelChannels(image)*sizeof(*x_vector));
3181   scanline=x_vector;
3182   if (image->rows != scale_image->rows)
3183     scanline=(double *) AcquireQuantumMemory((size_t) image->columns,
3184       GetPixelChannels(image)*sizeof(*scanline));
3185   scale_scanline=(double *) AcquireQuantumMemory((size_t)
3186     scale_image->columns,MaxPixelChannels*sizeof(*scale_scanline));
3187   y_vector=(double *) AcquireQuantumMemory((size_t) image->columns,
3188     GetPixelChannels(image)*sizeof(*y_vector));
3189   if ((scanline == (double *) NULL) ||
3190       (scale_scanline == (double *) NULL) ||
3191       (x_vector == (double *) NULL) ||
3192       (y_vector == (double *) NULL))
3193     {
3194       scale_image=DestroyImage(scale_image);
3195       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3196     }
3197   /*
3198     Scale image.
3199   */
3200   number_rows=0;
3201   next_row=MagickTrue;
3202   span.y=1.0;
3203   scale.y=(double) scale_image->rows/(double) image->rows;
3204   for (i=0; i < (ssize_t) (GetPixelChannels(image)*image->columns); i++)
3205     y_vector[i]=0.0;
3206   n=0;
3207   status=MagickTrue;
3208   image_view=AcquireVirtualCacheView(image,exception);
3209   scale_view=AcquireAuthenticCacheView(scale_image,exception);
3210   for (y=0; y < (ssize_t) scale_image->rows; y++)
3211   {
3212     register const Quantum
3213       *restrict p;
3214
3215     register Quantum
3216       *restrict q;
3217
3218     register ssize_t
3219       x;
3220
3221     if (status == MagickFalse)
3222       break;
3223     q=QueueCacheViewAuthenticPixels(scale_view,0,y,scale_image->columns,1,
3224       exception);
3225     if (q == (Quantum *) NULL)
3226       {
3227         status=MagickFalse;
3228         break;
3229       }
3230     alpha=1.0;
3231     if (scale_image->rows == image->rows)
3232       {
3233         /*
3234           Read a new scanline.
3235         */
3236         p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
3237           exception);
3238         if (p == (const Quantum *) NULL)
3239           {
3240             status=MagickFalse;
3241             break;
3242           }
3243         for (x=0; x < (ssize_t) image->columns; x++)
3244         {
3245           if (GetPixelReadMask(image,p) == 0)
3246             {
3247               p+=GetPixelChannels(image);
3248               continue;
3249             }
3250           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3251           {
3252             PixelChannel
3253               channel;
3254
3255             PixelTrait
3256               traits;
3257
3258             channel=GetPixelChannelChannel(image,i);
3259             traits=GetPixelChannelTraits(image,channel);
3260             if ((traits & BlendPixelTrait) == 0)
3261               {
3262                 x_vector[x*GetPixelChannels(image)+i]=(double) p[i];
3263                 continue;
3264               }
3265             alpha=QuantumScale*GetPixelAlpha(image,p);
3266             x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
3267           }
3268           p+=GetPixelChannels(image);
3269         }
3270       }
3271     else
3272       {
3273         /*
3274           Scale Y direction.
3275         */
3276         while (scale.y < span.y)
3277         {
3278           if ((next_row != MagickFalse) &&
3279               (number_rows < (ssize_t) image->rows))
3280             {
3281               /*
3282                 Read a new scanline.
3283               */
3284               p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
3285                 exception);
3286               if (p == (const Quantum *) NULL)
3287                 {
3288                   status=MagickFalse;
3289                   break;
3290                 }
3291               for (x=0; x < (ssize_t) image->columns; x++)
3292               {
3293                 if (GetPixelReadMask(image,p) == 0)
3294                   {
3295                     p+=GetPixelChannels(image);
3296                     continue;
3297                   }
3298                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3299                 {
3300                   PixelChannel
3301                     channel;
3302
3303                   PixelTrait
3304                     traits;
3305
3306                   channel=GetPixelChannelChannel(image,i);
3307                   traits=GetPixelChannelTraits(image,channel);
3308                   if ((traits & BlendPixelTrait) == 0)
3309                     {
3310                       x_vector[x*GetPixelChannels(image)+i]=(double)
3311                         p[i];
3312                       continue;
3313                     }
3314                   alpha=QuantumScale*GetPixelAlpha(image,p);
3315                   x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
3316                 }
3317                 p+=GetPixelChannels(image);
3318               }
3319               number_rows++;
3320             }
3321           for (x=0; x < (ssize_t) image->columns; x++)
3322             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3323               y_vector[x*GetPixelChannels(image)+i]+=scale.y*
3324                 x_vector[x*GetPixelChannels(image)+i];
3325           span.y-=scale.y;
3326           scale.y=(double) scale_image->rows/(double) image->rows;
3327           next_row=MagickTrue;
3328         }
3329         if ((next_row != MagickFalse) && (number_rows < (ssize_t) image->rows))
3330           {
3331             /*
3332               Read a new scanline.
3333             */
3334             p=GetCacheViewVirtualPixels(image_view,0,n++,image->columns,1,
3335               exception);
3336             if (p == (const Quantum *) NULL)
3337               {
3338                 status=MagickFalse;
3339                 break;
3340               }
3341             for (x=0; x < (ssize_t) image->columns; x++)
3342             {
3343               if (GetPixelReadMask(image,p) == 0)
3344                 {
3345                   p+=GetPixelChannels(image);
3346                   continue;
3347                 }
3348               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3349               {
3350                 PixelChannel
3351                   channel;
3352
3353                 PixelTrait
3354                   traits;
3355
3356                 channel=GetPixelChannelChannel(image,i);
3357                 traits=GetPixelChannelTraits(image,channel);
3358                 if ((traits & BlendPixelTrait) == 0)
3359                   {
3360                     x_vector[x*GetPixelChannels(image)+i]=(double) p[i];
3361                     continue;
3362                   }
3363                 alpha=QuantumScale*GetPixelAlpha(image,p);
3364                 x_vector[x*GetPixelChannels(image)+i]=alpha*p[i];
3365               }
3366               p+=GetPixelChannels(image);
3367             }
3368             number_rows++;
3369             next_row=MagickFalse;
3370           }
3371         for (x=0; x < (ssize_t) image->columns; x++)
3372         {
3373           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3374           {
3375             pixel[i]=y_vector[x*GetPixelChannels(image)+i]+span.y*
3376               x_vector[x*GetPixelChannels(image)+i];
3377             scanline[x*GetPixelChannels(image)+i]=pixel[i];
3378             y_vector[x*GetPixelChannels(image)+i]=0.0;
3379           }
3380         }
3381         scale.y-=span.y;
3382         if (scale.y <= 0)
3383           {
3384             scale.y=(double) scale_image->rows/(double) image->rows;
3385             next_row=MagickTrue;
3386           }
3387         span.y=1.0;
3388       }
3389     if (scale_image->columns == image->columns)
3390       {
3391         /*
3392           Transfer scanline to scaled image.
3393         */
3394         for (x=0; x < (ssize_t) scale_image->columns; x++)
3395         {
3396           if (GetPixelReadMask(scale_image,q) == 0)
3397             {
3398               q+=GetPixelChannels(scale_image);
3399               continue;
3400             }
3401           for (i=0; i < (ssize_t) GetPixelChannels(scale_image); i++)
3402           {
3403             ssize_t
3404               offset;
3405
3406             channel=GetPixelChannelChannel(scale_image,i);
3407             traits=GetPixelChannelTraits(image,channel);
3408             scale_traits=GetPixelChannelTraits(scale_image,channel);
3409             if ((traits == UndefinedPixelTrait) ||
3410                 (scale_traits == UndefinedPixelTrait))
3411               continue;
3412             offset=GetPixelChannelOffset(image,channel);
3413             if ((traits & BlendPixelTrait) == 0)
3414               {
3415                 SetPixelChannel(scale_image,channel,ClampToQuantum(
3416                   scanline[x*GetPixelChannels(image)+offset]),q);
3417                 continue;
3418               }
3419             alpha=QuantumScale*scanline[x*GetPixelChannels(image)+
3420               GetPixelChannelChannel(image,AlphaPixelChannel)];
3421             gamma=PerceptibleReciprocal(alpha);
3422             SetPixelChannel(scale_image,channel,ClampToQuantum(gamma*scanline[
3423               x*GetPixelChannels(image)+offset]),q);
3424           }
3425           q+=GetPixelChannels(scale_image);
3426         }
3427       }
3428     else
3429       {
3430         ssize_t
3431           n;
3432
3433         /*
3434           Scale X direction.
3435         */
3436         next_column=MagickFalse;
3437         n=0;
3438         span.x=1.0;
3439         for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3440           pixel[i]=0.0;
3441         for (x=0; x < (ssize_t) image->columns; x++)
3442         {
3443           scale.x=(double) scale_image->columns/(double) image->columns;
3444           while (scale.x >= span.x)
3445           {
3446             if (next_column != MagickFalse)
3447               {
3448                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3449                   pixel[i]=0.0;
3450                 n++;
3451               }
3452             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3453             {
3454               PixelChannel
3455                 channel;
3456
3457               PixelTrait
3458                 traits;
3459
3460               channel=GetPixelChannelChannel(image,i);
3461               traits=GetPixelChannelTraits(image,channel);
3462               if (traits == UndefinedPixelTrait)
3463                 continue;
3464               pixel[i]+=span.x*scanline[x*GetPixelChannels(image)+i];
3465               scale_scanline[n*MaxPixelChannels+channel]=pixel[i];
3466             }
3467             scale.x-=span.x;
3468             span.x=1.0;
3469             next_column=MagickTrue;
3470           }
3471         if (scale.x > 0)
3472           {
3473             if (next_column != MagickFalse)
3474               {
3475                 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3476                   pixel[i]=0.0;
3477                 n++;
3478                 next_column=MagickFalse;
3479               }
3480             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3481               pixel[i]+=scale.x*scanline[x*GetPixelChannels(image)+i];
3482             span.x-=scale.x;
3483           }
3484       }
3485       if (span.x > 0)
3486         {
3487           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3488             pixel[i]+=span.x*scanline[(x-1)*GetPixelChannels(image)+i];
3489         }
3490       if ((next_column == MagickFalse) &&
3491           ((ssize_t) n < (ssize_t) scale_image->columns))
3492         for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3493         {
3494           channel=GetPixelChannelChannel(image,i);
3495           scale_scanline[n*MaxPixelChannels+channel]=pixel[i];
3496         }
3497       /*
3498         Transfer scanline to scaled image.
3499       */
3500       for (x=0; x < (ssize_t) scale_image->columns; x++)
3501       {
3502         if (GetPixelReadMask(scale_image,q) == 0)
3503           {
3504             q+=GetPixelChannels(scale_image);
3505             continue;
3506           }
3507         for (i=0; i < (ssize_t) GetPixelChannels(scale_image); i++)
3508         {
3509           channel=GetPixelChannelChannel(scale_image,i);
3510           traits=GetPixelChannelTraits(image,channel);
3511           scale_traits=GetPixelChannelTraits(scale_image,channel);
3512           if ((traits == UndefinedPixelTrait) ||
3513               (scale_traits == UndefinedPixelTrait))
3514             continue;
3515           if ((traits & BlendPixelTrait) == 0)
3516             {
3517               SetPixelChannel(scale_image,channel,ClampToQuantum(
3518                 scale_scanline[x*MaxPixelChannels+channel]),q);
3519               continue;
3520             }
3521           SetPixelChannel(scale_image,channel,ClampToQuantum(scale_scanline[
3522             x*MaxPixelChannels+channel]),q);
3523         }
3524         q+=GetPixelChannels(scale_image);
3525       }
3526     }
3527     if (SyncCacheViewAuthenticPixels(scale_view,exception) == MagickFalse)
3528       {
3529         status=MagickFalse;
3530         break;
3531       }
3532     proceed=SetImageProgress(image,ScaleImageTag,(MagickOffsetType) y,
3533       image->rows);
3534     if (proceed == MagickFalse)
3535       {
3536         status=MagickFalse;
3537         break;
3538       }
3539   }
3540   scale_view=DestroyCacheView(scale_view);
3541   image_view=DestroyCacheView(image_view);
3542   /*
3543     Free allocated memory.
3544   */
3545   y_vector=(double *) RelinquishMagickMemory(y_vector);
3546   scale_scanline=(double *) RelinquishMagickMemory(scale_scanline);
3547   if (scale_image->rows != image->rows)
3548     scanline=(double *) RelinquishMagickMemory(scanline);
3549   x_vector=(double *) RelinquishMagickMemory(x_vector);
3550   scale_image->type=image->type;
3551   if (status == MagickFalse)
3552     scale_image=DestroyImage(scale_image);
3553   return(scale_image);
3554 }
3555 \f
3556 /*
3557 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3558 %                                                                             %
3559 %                                                                             %
3560 %                                                                             %
3561 %   T h u m b n a i l I m a g e                                               %
3562 %                                                                             %
3563 %                                                                             %
3564 %                                                                             %
3565 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3566 %
3567 %  ThumbnailImage() changes the size of an image to the given dimensions and
3568 %  removes any associated profiles.  The goal is to produce small low cost
3569 %  thumbnail images suited for display on the Web.
3570 %
3571 %  The format of the ThumbnailImage method is:
3572 %
3573 %      Image *ThumbnailImage(const Image *image,const size_t columns,
3574 %        const size_t rows,ExceptionInfo *exception)
3575 %
3576 %  A description of each parameter follows:
3577 %
3578 %    o image: the image.
3579 %
3580 %    o columns: the number of columns in the scaled image.
3581 %
3582 %    o rows: the number of rows in the scaled image.
3583 %
3584 %    o exception: return any errors or warnings in this structure.
3585 %
3586 */
3587 MagickExport Image *ThumbnailImage(const Image *image,const size_t columns,
3588   const size_t rows,ExceptionInfo *exception)
3589 {
3590 #define SampleFactor  5
3591
3592   char
3593     value[MaxTextExtent];
3594
3595   const char
3596     *name;
3597
3598   Image
3599     *thumbnail_image;
3600
3601   double
3602     x_factor,
3603     y_factor;
3604
3605   size_t
3606     version;
3607
3608   struct stat
3609     attributes;
3610
3611   assert(image != (Image *) NULL);
3612   assert(image->signature == MagickSignature);
3613   if (image->debug != MagickFalse)
3614     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3615   assert(exception != (ExceptionInfo *) NULL);
3616   assert(exception->signature == MagickSignature);
3617   x_factor=(double) columns/(double) image->columns;
3618   y_factor=(double) rows/(double) image->rows;
3619   if ((x_factor*y_factor) > 0.1)
3620     thumbnail_image=ResizeImage(image,columns,rows,image->filter,exception);
3621   else
3622     if (((SampleFactor*columns) < 128) || ((SampleFactor*rows) < 128))
3623       thumbnail_image=ResizeImage(image,columns,rows,image->filter,exception);
3624     else
3625       {
3626         Image
3627           *sample_image;
3628
3629         sample_image=SampleImage(image,SampleFactor*columns,SampleFactor*rows,
3630           exception);
3631         if (sample_image == (Image *) NULL)
3632           return((Image *) NULL);
3633         thumbnail_image=ResizeImage(sample_image,columns,rows,image->filter,
3634           exception);
3635         sample_image=DestroyImage(sample_image);
3636       }
3637   if (thumbnail_image == (Image *) NULL)
3638     return(thumbnail_image);
3639   (void) ParseAbsoluteGeometry("0x0+0+0",&thumbnail_image->page);
3640   if (thumbnail_image->alpha_trait != BlendPixelTrait)
3641     (void) SetImageAlphaChannel(thumbnail_image,OpaqueAlphaChannel,exception);
3642   thumbnail_image->depth=8;
3643   thumbnail_image->interlace=NoInterlace;
3644   /*
3645     Strip all profiles except color profiles.
3646   */
3647   ResetImageProfileIterator(thumbnail_image);
3648   for (name=GetNextImageProfile(thumbnail_image); name != (const char *) NULL; )
3649   {
3650     if ((LocaleCompare(name,"icc") != 0) && (LocaleCompare(name,"icm") != 0))
3651      {
3652        (void) DeleteImageProfile(thumbnail_image,name);
3653        ResetImageProfileIterator(thumbnail_image);
3654      }
3655     name=GetNextImageProfile(thumbnail_image);
3656   }
3657   (void) DeleteImageProperty(thumbnail_image,"comment");
3658   (void) CopyMagickString(value,image->magick_filename,MaxTextExtent);
3659   if (strstr(image->magick_filename,"//") == (char *) NULL)
3660     (void) FormatLocaleString(value,MaxTextExtent,"file://%s",
3661       image->magick_filename);
3662   (void) SetImageProperty(thumbnail_image,"Thumb::URI",value,exception);
3663   (void) CopyMagickString(value,image->magick_filename,MaxTextExtent);
3664   if ( IfMagickTrue(GetPathAttributes(image->filename,&attributes)) )
3665     {
3666       (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3667         attributes.st_mtime);
3668       (void) SetImageProperty(thumbnail_image,"Thumb::MTime",value,exception);
3669     }
3670   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3671     attributes.st_mtime);
3672   (void) FormatMagickSize(GetBlobSize(image),MagickFalse,value);
3673   (void) ConcatenateMagickString(value,"B",MaxTextExtent);
3674   (void) SetImageProperty(thumbnail_image,"Thumb::Size",value,exception);
3675   (void) FormatLocaleString(value,MaxTextExtent,"image/%s",image->magick);
3676   LocaleLower(value);
3677   (void) SetImageProperty(thumbnail_image,"Thumb::Mimetype",value,exception);
3678   (void) SetImageProperty(thumbnail_image,"software",GetMagickVersion(&version),
3679     exception);
3680   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3681     image->magick_columns);
3682   (void) SetImageProperty(thumbnail_image,"Thumb::Image::Width",value,
3683     exception);
3684   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3685     image->magick_rows);
3686   (void) SetImageProperty(thumbnail_image,"Thumb::Image::Height",value,
3687     exception);
3688   (void) FormatLocaleString(value,MaxTextExtent,"%.20g",(double)
3689     GetImageListLength(image));
3690   (void) SetImageProperty(thumbnail_image,"Thumb::Document::Pages",value,
3691     exception);
3692   return(thumbnail_image);
3693 }