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