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

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