resize.c

Go to the documentation of this file.

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%                 RRRR   EEEEE  SSSSS  IIIII  ZZZZZ  EEEEE                    %
%                 R   R  E      SS       I       ZZ  E                        %
%                 RRRR   EEE     SSS     I     ZZZ   EEE                      %
%                 R R    E         SS    I    ZZ     E                        %
%                 R  R   EEEEE  SSSSS  IIIII  ZZZZZ  EEEEE                    %
%                                                                             %
%                                                                             %
%                      MagickCore Image Resize Methods                        %
%                                                                             %
%                              Software Design                                %
%                                John Cristy                                  %
%                                 July 1992                                   %
%                                                                             %
%                                                                             %
%  Copyright 1999-2008 ImageMagick Studio LLC, a non-profit organization      %
%  dedicated to making software imaging solutions freely available.           %
%                                                                             %
%  You may not use this file except in compliance with the License.  You may  %
%  obtain a copy of the License at                                            %
%                                                                             %
%    http://www.imagemagick.org/script/license.php                            %
%                                                                             %
%  Unless required by applicable law or agreed to in writing, software        %
%  distributed under the License is distributed on an "AS IS" BASIS,          %
%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
%  See the License for the specific language governing permissions and        %
%  limitations under the License.                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
*/

/*
  Include declarations.
*/
#include "magick/studio.h"
#include "magick/artifact.h"
#include "magick/blob.h"
#include "magick/cache.h"
#include "magick/cache-view.h"
#include "magick/color.h"
#include "magick/color-private.h"
#include "magick/draw.h"
#include "magick/exception.h"
#include "magick/exception-private.h"
#include "magick/gem.h"
#include "magick/image.h"
#include "magick/image-private.h"
#include "magick/list.h"
#include "magick/memory_.h"
#include "magick/pixel-private.h"
#include "magick/property.h"
#include "magick/monitor.h"
#include "magick/monitor-private.h"
#include "magick/pixel.h"
#include "magick/option.h"
#include "magick/resample.h"
#include "magick/resize.h"
#include "magick/resize-private.h"
#include "magick/string_.h"
#include "magick/utility.h"
#include "magick/version.h"
#if defined(MAGICKCORE_LQR_DELEGATE)
#include <lqr.h>
#endif

/*
  Typedef declarations.
*/
struct _ResizeFilter
{
  MagickRealType
    (*filter)(const MagickRealType,const ResizeFilter *),
    (*window)(const MagickRealType,const ResizeFilter *),
    support,   /* filter region of support - the filter support limit */
    window_support,  /* window support, usally equal to support (expert only) */
    scale,     /* dimension to scale to fit window support (usally 1.0) */
    blur,      /* x-scale (blur-sharpen) */
    cubic[8];  /* cubic coefficents for smooth Cubic filters */

  unsigned long
    signature;
};

/*
  Forward declaractions.
*/
static MagickRealType
  I0(MagickRealType x),
  BesselOrderOne(MagickRealType);

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   F i l t e r F u n c t i o n s                                             %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  These are the various filter and windowing functions that are provided,
%
%  They are all internal to this module only.  See AcquireResizeFilterInfo()
%  for details of the access to these functions, via the
%  GetResizeFilterSupport() and GetResizeFilterWeight() API interface.
%
%  The individual filter functions have this format...
%
%     static MagickRealtype *FilterName(const MagickRealType x,
%        const MagickRealType support)
%
%    o x: the distance from the sampling point
%         generally in the range of  0 to support
%         The GetResizeFilterWeight() ensures this a positive value.
%
%    o resize_filter: Current Filter Information
%        This allows function to access support, and posibly other
%        pre-calculated information defineding the functions.
%
*/

static MagickRealType Bessel(const MagickRealType x,
  const ResizeFilter *magick_unused(resize_filter))
{
  /*
    See Pratt "Digital Image Processing" p.97 for Bessel functions

    This function actually a X-scaled Jinc(x) function.
      http://mathworld.wolfram.com/JincFunction.html
    And on page 11 of...
      http://www.ph.ed.ac.uk/%7ewjh/teaching/mo/slides/lens/lens.pdf
  */
  if (x == 0.0)
    return((MagickRealType) (MagickPI/4.0));
  return(BesselOrderOne(MagickPI*x)/(2.0*x));
}

static MagickRealType Blackman(const MagickRealType x,
     const ResizeFilter *magick_unused(resize_filter))
{
  /*
    Blackman: 2rd Order cosine windowing function.
  */
  return(0.42+0.5*cos(MagickPI*(double) x)+0.08*cos(2.0*MagickPI*(double) x));
}

static MagickRealType Bohman(const MagickRealType x,
     const ResizeFilter *magick_unused(resize_filter))
{
  /*
    Bohman: 2rd Order cosine windowing function.
  */
  return((1-x)*cos(MagickPI*(double) x)+sin(MagickPI*(double) x)/MagickPI);
}

static MagickRealType Box(const MagickRealType magick_unused(x),
  const ResizeFilter *magick_unused(resize_filter))
{
  /*
    Just return 1.0, filter will still be clipped by its support window.
  */
  return(1.0);
}

static MagickRealType CubicBC(const MagickRealType x,
  const ResizeFilter *resize_filter)
{
  /*
    Cubic Filters using B,C determined values:

    Mitchell-Netravali  B=1/3 C=1/3   Qualitively ideal Cubic Filter
    Catmull-Rom         B= 0  C=1/2   Cublic Interpolation Function
    Cubic B-Spline      B= 1  C= 0    Spline Approximation of Gaussian
    Hermite             B= 0  C= 0    Quadratic Spline (support = 1)

    See paper by Mitchell and Netravali,
      Reconstruction Filters in Computer Graphics
      Computer Graphics, Volume 22, Number 4, August 1988
        http://www.cs.utexas.edu/users/fussell/courses/cs384g/
                 lectures/mitchell/Mitchell.pdf

    Coefficents are determined from B,C values
       P0 = (  6 - 2*B       )/6
       P1 =         0
       P2 = (-18 +12*B + 6*C )/6
       P3 = ( 12 - 9*B - 6*C )/6
       Q0 = (      8*B +24*C )/6
       Q1 = (    -12*B -48*C )/6
       Q2 = (      6*B +30*C )/6
       Q3 = (    - 1*B - 6*C )/6

    Which is used to define the filter...
       P0 + P1*x + P2*x^2 + P3*x^3      0 <= x < 1
       Q0 + Q1*x + Q2*x^2 + Q3*x^3      1 <= x <= 2

    Which ensures function is continuous in value and derivative (slope).
  */
  if (x < 1.0)
    return(resize_filter->cubic[0]+x*(resize_filter->cubic[1]+x*
      (resize_filter->cubic[2]+x*resize_filter->cubic[3])));
  if (x < 2.0)
    return(resize_filter->cubic[4] +x*(resize_filter->cubic[5]+x*
      (resize_filter->cubic[6] +x*resize_filter->cubic[7])));
  return(0.0);
}

static MagickRealType Gaussian(const MagickRealType x,
  const ResizeFilter *magick_unused(resize_filter))
{
  return(exp((double) (-2.0*x*x))*sqrt(2.0/MagickPI));
}

static MagickRealType Hanning(const MagickRealType x,
  const ResizeFilter *magick_unused(resize_filter))
{
  /*
    A Cosine windowing function.
  */
  return(0.5+0.5*cos(MagickPI*(double) x));
}

static MagickRealType Hamming(const MagickRealType x,
  const ResizeFilter *magick_unused(resize_filter))
{
  /*
    A offset Cosine windowing function.
  */
  return(0.54+0.46*cos(MagickPI*(double) x));
}

static MagickRealType Kaiser(const MagickRealType x,
  const ResizeFilter *magick_unused(resize_filter))
{
#define Alpha  6.5
#define I0A  (1.0/I0(Alpha))

  /*
    Kaiser Windowing Function (bessel windowing):
      Alpha is a free value  from 5 to 8 (currently hardcoded to 6.5)
      Future: make alphand the IOA pre-calculation, a 'expert' setting.
  */
  return(I0A*I0(Alpha*sqrt((double) (1.0-x*x))));
}

static MagickRealType Lagrange(const MagickRealType x,
  const ResizeFilter *resize_filter)
{
  long
    n,
    order;

  MagickRealType
    value;

  register long
    i;

  /*
    Lagrange Piece-Wise polynomial fit of Sinc:
      N is the 'order' of the lagrange function and depends on
      the overall support window size of the filter. That is for
      a support of 2, gives a lagrange-4 or piece-wise cubic functions

      Note that n is the specific piece of the piece-wise function to calculate.

      See Survey: Interpolation Methods, IEEE Transactions on Medical Imaging,
      Vol 18, No 11, November 1999, p1049-1075, -- Equation 27 on p1064
  */
  if (x > resize_filter->support)
    return(0.0);
  order=(long) (2.0*resize_filter->window_support);  /* number of pieces */
  n=(long) ((1.0*order)/2.0+x);  /* which piece does x belong to */
  value=1.0f;
  for (i=0; i < order; i++)
    if (i != n)
      value*=(n-i-x)/(n-i);
  return(value);
}

static MagickRealType Quadratic(const MagickRealType x,
  const ResizeFilter *magick_unused(resize_filter))
{
  /*
    2rd order (quadratic) B-Spline approximation of Gaussian.
  */
  if (x < 0.5)
    return(0.75-x*x);
  if (x < 1.5)
    return(0.5*(x-1.5)*(x-1.5));
  return(0.0);
}

static MagickRealType Sinc(const MagickRealType x,
  const ResizeFilter *magick_unused(resize_filter))
{
  /*
    This function actually a X-scaled Sinc(x) function.
  */
  if (x == 0.0)
    return(1.0);
  return(sin(MagickPI*(double) x)/(MagickPI*(double) x));
}

static MagickRealType Triangle(const MagickRealType x,
  const ResizeFilter *magick_unused(resize_filter))
{
  /*
    1rd order (linear) B-Spline,  bilinear interpolation,
    Tent 1D filter, or a Bartlett 2D Cone filter
  */
  if (x < 1.0)
    return(1.0-x);
  return(0.0);
}

static MagickRealType Welsh(const MagickRealType x,
  const ResizeFilter *magick_unused(resize_filter))
{
  /*
    Welsh parabolic windowing filter.
  */
  if (x <  1.0)
    return(1.0-x*x);
  return(0.0);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   A c q u i r e R e s i z e F i l t e r                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AcquireResizeFilter() allocates the ResizeFilter structure.  Choose from
%  these filters:
%
%  FIR (Finite impulse Response) Filters
%      Box         Triangle   Quadratic
%      Cubic       Hermite    Catrom
%      Mitchell
%
%  IIR (Infinite impulse Response) Filters
%      Gaussian     Sinc        Bessel
%
%  Windowed Sinc/Bessel Method
%      Blackman     Hanning     Hamming
%      Kaiser       Lancos (Sinc)
%
%  FIR filters are used as is, and are limited by that filters support window
%  (unless over-ridden).  'Gaussian' while classed as an IIR filter, is also
%  simply clipped by its support size (1.5).
%
%  Requesting a windowed filter will return either a windowed Sinc, for a one
%  dimentional orthogonal filtering method, such as ResizeImage(), or a
%  windowed Bessel for image operations requiring a two dimentional
%  cylindrical filtering method, such a DistortImage().  Which function is
%  is used set by the "cylindrical" boolean argument.
%
%  Directly requesting 'Sinc' or 'Bessel' will force the use of that filter
%  function, with a default 'Blackman' windowing method.  This not however
%  recommended as it removes the correct filter selection for different
%  filtering image operations.  Selecting a window filtering method is better.
%
%  Lanczos is purely special case of a Sinc windowed Sinc, but defulting to
%  a 3 lobe support, rather that the default 4 lobe support.
%
%  Special options can be used to override specific, or all the filter
%  settings.   However doing so is not advisible unless you have expert
%  knowledge of the use of resampling filtered techniques. Extreme caution is
%  advised.
%
%    "filter:filter"    Select this function as the filter.
%        If a "filter:window" operation is not provided, then no windowing
%        will be performed on the selected filter, (support clipped)
%
%        This can be used to force the use of a windowing method as filter,
%        request a 'Sinc' filter in a radially filtered operation, or the
%        'Bessel' filter for a othogonal filtered operation.
%
%    "filter:window"   Select this windowing function for the filter.
%        While any filter could be used as a windowing function,
%        using that filters first lobe over the whole support window,
%        using a non-windowing method is not advisible.
%
%    "filter:lobes"    Number of lobes to use for the Sinc/Bessel filter.
%        This a simper method of setting filter support size that will
%        correctly handle the Sinc/Bessel switch for an operators filtering
%        requirements.
%
%    "filter:support"  Set the support size for filtering to the size given
%        This not recomented for Sinc/Bessel windowed filters, but is
%        used for simple filters like FIR filters, and the Gaussian Filter.
%        This will override any 'filter:lobes' option.
%
%    "filter:blur"     Scale the filter and support window by this amount.
%        A value >1 will generally result in a more burred image with
%        more ringing effects, while a value <1 will sharpen the
%        resulting image with more aliasing and Morie effects.
%
%    "filter:win-support"  Scale windowing function to this size instead.
%        This causes the windowing (or self-windowing Lagrange filter)
%        to act is if the support winodw it much much larger than what
%        is actually supplied to the calling operator.  The filter however
%        is still clipped to the real support size given.  If unset this
%        will equal the normal filter support size.
%
%    "filter:b"
%    "filter:c"    Override the preset B,C values for a Cubic type of filter
%         If only one of these are given it is assumes to be a 'Keys'
%         type of filter such that B+2C=1, where Keys 'alpha' value = C
%
%    "filter:verbose"   Output verbose plotting data for graphing the
%         resulting filter over the whole support range (with blur effect).
%
%  Set a true un-windowed Sinc filter with 10 lobes (very slow)
%     -set option:filter:filter  Sinc
%     -set option:filter:lobes   8
%
%  For example force an 8 lobe Lanczos (Sinc or Bessel) filter...
%     -filter Lanczos
%     -set option:filter:lobes   8
%
%  The format of the AcquireResizeFilter method is:
%
%      ResizeFilter *AcquireResizeFilter(const Image *image,
%        const FilterTypes filter_type, const MagickBooleanType radial,
%        ExceptionInfo *exception)
%
%    o image: the image.
%
%    o filter: the filter type, defining a preset filter, window and support.
%
%    o blur: blur the filter by this amount, use 1.0 if unknown.
%            Image artifact "filter:blur"  will override this old usage
%
%    o radial: 1D orthogonal filter (Sinc) or 2D radial filter (Bessel)
%
%    o exception: return any errors or warnings in this structure.
%
*/
MagickExport ResizeFilter *AcquireResizeFilter(const Image *image,
  const FilterTypes filter, const MagickRealType blur,
  const MagickBooleanType cylindrical,ExceptionInfo *exception)
{
  const char
    *artifact;

  FilterTypes
    filter_type,
    window_type;

  long
    filter_artifact;

  MagickRealType
    B,
    C;

  register ResizeFilter
    *resize_filter;

  /*
    Table Mapping given Filter, into  Weighting and Windowing functions.
    A 'Box' windowing function means its a simble non-windowed filter.
    A 'Sinc' filter function (must be windowed) could be upgraded to a
    'Bessel' filter if a "cylindrical" filter is requested, unless a "Sinc"
    filter specifically request.
  */
  static struct
  {
    FilterTypes
      filter,
      window;
  } const mapping[SentinelFilter] =
  {
    { UndefinedFilter, BoxFilter },  /* undefined */
    { PointFilter,     BoxFilter },  /* special, nearest-neighbour filter */
    { BoxFilter,       BoxFilter },  /* Box averaging Filter */
    { TriangleFilter,  BoxFilter },  /* Linear Interpolation Filter */
    { HermiteFilter,   BoxFilter },      /* Hermite interpolation filter */
    { SincFilter,      HanningFilter },  /* Hanning -- Cosine-Sinc */
    { SincFilter,      HammingFilter },  /* Hamming --  '' variation */
    { SincFilter,      BlackmanFilter }, /* Blackman -- 2*Cosine-Sinc */
    { GaussianFilter,  BoxFilter },      /* Gaussain Blurring filter */
    { QuadraticFilter, BoxFilter },      /* Quadratic Gaussian approximation */
    { CubicFilter,     BoxFilter },      /* Cubic Gaussian approximation */
    { CatromFilter,    BoxFilter },      /* Cubic Interpolator */
    { MitchellFilter,  BoxFilter },      /* 'ideal' Cubic Filter */
    { LanczosFilter,   SincFilter },     /* Special, 3 lobed Sinc-Sinc */
    { BesselFilter,    BlackmanFilter }, /* 3 lobed bessel -specific request */
    { SincFilter,      BlackmanFilter }, /* 4 lobed sinc - specific request */
    { SincFilter,      KaiserFilter },   /* Kaiser --  SqRoot-Sinc */
    { SincFilter,      WelshFilter },    /* Welsh -- Parabolic-Sinc */
    { SincFilter,      CubicFilter },    /* Parzen -- Cubic-Sinc */
    { LagrangeFilter,  BoxFilter },      /* Lagrange self-windowing filter */
    { SincFilter,      BohmanFilter },   /* Bohman -- 2*Cosine-Sinc */
    { SincFilter,      TriangleFilter }  /* Bartlett -- Triangle-Sinc */
  };
  /*
    Table maping the filter/window function from the above table to the actual
    filter/window function call to use.  The default support size for that
    filter as a weighting function, and the point to scale when that function
    is used as a windowing function (typ 1.0).
  */
  static struct
  {
    MagickRealType
      (*function)(const MagickRealType, const ResizeFilter*),
      support,  /* default support size for function as a filter */
      scale,    /* size windowing function, for scaling windowing function */
      B,
      C;        /* Cubic Filter factors for a CubicBC function, else ignored */
  } const filters[SentinelFilter] =
  {
    { Box,       0.0f,  0.5f, 0.0f, 0.0f }, /* Undefined */
    { Box,       0.0f,  0.5f, 0.0f, 0.0f }, /* Point */
    { Box,       0.5f,  0.5f, 0.0f, 0.0f }, /* Box */
    { Triangle,  1.0f,  1.0f, 0.0f, 0.0f }, /* Triangle */
    { CubicBC,   1.0f,  1.0f, 0.0f, 0.0f }, /* Hermite, Cubic B=C=0 */
    { Hanning,   1.0f,  1.0f, 0.0f, 0.0f }, /* Hanning, Cosine window */
    { Hamming,   1.0f,  1.0f, 0.0f, 0.0f }, /* Hamming, '' variation */
    { Blackman,  1.0f,  1.0f, 0.0f, 0.0f }, /* Blackman, 2*cos window */
    { Gaussian,  1.5f,  1.5f, 0.0f, 0.0f }, /* Gaussian */
    { Quadratic, 1.5f,  1.5f, 0.0f, 0.0f }, /* Quadratic Gaussian */
    { CubicBC,   2.0f,  2.0f, 1.0f, 0.0f }, /* B-Spline of Gaussian B=1 C=0 */
    { CubicBC,   2.0f,  1.0f, 0.0f, 0.5f }, /* Catmull-Rom  B=0 C=1/2 */
    { CubicBC,   2.0f,  1.0f, 1.0f/3.0f, 1.0f/3.0f }, /* Mitchel B=C=1/3 */
    { Sinc,      3.0f,  1.0f, 0.0f, 0.0f }, /* Lanczos, 3 lobed Sinc-Sinc */
    { Bessel,    3.2383f,1.2197f,.0f,.0f }, /* 3 lobed Blackman-Bessel */
    { Sinc,      4.0f,  1.0f, 0.0f, 0.0f }, /* 4 lobed Blackman-Sinc   */
    { Kaiser,    1.0f,  1.0f, 0.0f, 0.0f }, /* Kaiser, sq-root windowing */
    { Welsh,     1.0f,  1.0f, 0.0f, 0.0f }, /* Welsh, Parabolic windowing */
    { CubicBC,   2.0f,  2.0f, 1.0f, 0.0f }, /* Parzen, B-Spline windowing */
    { Lagrange,  2.0f,  1.0f, 0.0f, 0.0f }, /* Lagrangian Filter */
    { Bohman,    1.0f,  1.0f, 0.0f, 0.0f }, /* Bohman, 2*Cosine windowing */
    { Triangle,  1.0f,  1.0f, 0.0f, 0.0f }  /* Bartlett, Triangle windowing */
  };
  /*
    The known zero crossings of the Bessel() or the Jinc(x*PI) function
    Found by using
      http://cose.math.bas.bg/webMathematica/webComputing/BesselZeros.jsp
    for Jv-function with v=1,  then dividing X-roots by PI (tabled below)
  */
  static MagickRealType
    bessel_zeros[16] =
    {
      1.21966989126651f,
      2.23313059438153f,
      3.23831548416624f,
      4.24106286379607f,
      5.24276437687019f,
      6.24392168986449f,
      7.24475986871996f,
      8.24539491395205f,
      9.24589268494948f,
      10.2462933487549f,
      11.2466227948779f,
      12.2468984611381f,
      13.2471325221811f,
      14.2473337358069f,
      15.2475085630373f,
      16.247661874701f
   };

  assert(image != (const Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  assert(UndefinedFilter < filter && filter < SentinelFilter);
  assert(exception != (ExceptionInfo *) NULL);
  assert(exception->signature == MagickSignature);

  resize_filter=(ResizeFilter *) AcquireMagickMemory(sizeof(*resize_filter));
  if (resize_filter == (ResizeFilter *) NULL)
    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");

  /* defaults for the requested filter */
  filter_type = mapping[filter].filter;
  window_type = mapping[filter].window;


  /* Filter blur -- scaling both filter and support window */
  resize_filter->blur = blur;
  artifact=GetImageArtifact(image,"filter:blur");
  if (artifact != (const char *) NULL)
    resize_filter->blur = atof(artifact);
  if ( resize_filter->blur < MagickEpsilon )
    resize_filter->blur = (MagickRealType) MagickEpsilon;

  /* Modifications for Cylindrical filter use */
  if ( cylindrical != MagickFalse && filter != SincFilter ) {
    /* promote 1D Sinc Filter to a 2D Bessel filter */
    if ( filter_type == SincFilter )
      filter_type = BesselFilter;
    /* Prompote Lanczos (Sinc-Sinc) to Lanczos (Bessel-Bessel) */
    else if ( filter_type == LanczosFilter ) {
      filter_type = BesselFilter;
      window_type = BesselFilter;
    }
   /* Blur other filters appropriatally correct cylindrical usage */
    else if ( filter_type == GaussianFilter )
      /* Gaussian is scaled by  4*ln(2) and not 4*sqrt(2/MagickPI)
         - according to Paul Heckbert's paper on EWA resampling */
      resize_filter->blur *= 2.0*log(2.0)/sqrt(2.0/MagickPI);
    else if ( filter_type != BesselFilter )
      /* filters with a 1.0 zero root crossing by the first bessel_zero */
      resize_filter->blur *= bessel_zeros[0];
  }

  /* Override Filter Selection */
  artifact=GetImageArtifact(image,"filter:filter");
  if (artifact != (const char *) NULL) {
    /* raw filter request - no window function */
    filter_artifact=ParseMagickOption(MagickFilterOptions,
         MagickFalse,artifact);
    if ( UndefinedFilter < filter_artifact &&
             filter_artifact < SentinelFilter ) {
      filter_type = (FilterTypes) filter_artifact;
      window_type = BoxFilter;
    }
    /* Lanczos is nor a real filter but a self windowing Sinc/Bessel */
    if ( filter_artifact == LanczosFilter ) {
      filter_type = (cylindrical!=MagickFalse) ? BesselFilter : LanczosFilter;
      window_type = (cylindrical!=MagickFalse) ? BesselFilter : SincFilter;
    }
    /* Filter overwide with a specific window function? */
    artifact=GetImageArtifact(image,"filter:window");
    if (artifact != (const char *) NULL) {
      filter_artifact=ParseMagickOption(MagickFilterOptions,
            MagickFalse,artifact);
      if ( UndefinedFilter < filter_artifact &&
               filter_artifact < SentinelFilter ) {
        if ( filter_artifact != LanczosFilter )
          window_type = (FilterTypes) filter_artifact;
        else
          window_type = (cylindrical!=MagickFalse) ? BesselFilter : SincFilter;
      }
    }
  }
  else {
    /* window specified, but no filter function?  Assume Sinc/Bessel */
    artifact=GetImageArtifact(image,"filter:window");
    if (artifact != (const char *) NULL) {
      filter_artifact=ParseMagickOption(MagickFilterOptions,MagickFalse,
        artifact);
      if ( UndefinedFilter < filter_artifact &&
               filter_artifact < SentinelFilter ) {
        filter_type = (cylindrical!=MagickFalse) ? BesselFilter : SincFilter;
        if ( filter_artifact != LanczosFilter )
          window_type = (FilterTypes) filter_artifact;
        else
          window_type = filter_type;
      }
    }
  }

  resize_filter->filter   = filters[filter_type].function;
  resize_filter->support  = filters[filter_type].support;
  resize_filter->window   = filters[window_type].function;
  resize_filter->scale    = filters[window_type].scale;
  resize_filter->signature=MagickSignature;

  /* Filter support overrides */
  artifact=GetImageArtifact(image,"filter:lobes");
  if (artifact != (const char *) NULL) {
    long lobes = atol(artifact);
    if ( lobes < 1  ) lobes = 1;
    resize_filter->support = (MagickRealType) lobes;
    if ( filter_type == BesselFilter ) {
      if ( lobes > 16 ) lobes = 16;
      resize_filter->support = bessel_zeros[lobes-1];
    }
  }
  artifact=GetImageArtifact(image,"filter:support");
  if (artifact != (const char *) NULL)
    resize_filter->support = fabs(atof(artifact));

  /* Scale windowing function separatally to the support 'clipping' window
     that calling operator is planning to actually use. - Expert Use Only
  */
  resize_filter->window_support = resize_filter->support;
  artifact=GetImageArtifact(image,"filter:win-support");
  if (artifact != (const char *) NULL)
    resize_filter->window_support = fabs(atof(artifact));

  /* Set Cubic Spline B,C values, calculate Cubic coefficents */
  B=0.0;
  C=0.0;
  if ( filters[filter_type].function == CubicBC
       || filters[window_type].function == CubicBC ) {
    if ( filters[filter_type].function == CubicBC ) {
      B=filters[filter_type].B;
      C=filters[filter_type].C;
    }
    else if ( filters[window_type].function == CubicBC ) {
      B=filters[window_type].B;
      C=filters[window_type].C;
    }
    artifact=GetImageArtifact(image,"filter:b");
    if (artifact != (const char *) NULL) {
      B=atof(artifact);
      C=1.0-2.0*B;  /* Calculate C as if it is a Keys cubic filter */
      artifact=GetImageArtifact(image,"filter:c");
      if (artifact != (const char *) NULL)
        C=atof(artifact);
    }
    else {
      artifact=GetImageArtifact(image,"filter:c");
      if (artifact != (const char *) NULL) {
        C=atof(artifact);
        B=(1.0-C)/2.0; /* Calculate B as if it is a Keys cubic filter */
      }
    }
    /* Convert B,C values into Cubic Coefficents - See CubicBC() */
    resize_filter->cubic[0]=(  6.0 -2.0*B       )/6.0;
    resize_filter->cubic[1]=0.0;
    resize_filter->cubic[2]=(-18.0+12.0*B+ 6.0*C)/6.0;
    resize_filter->cubic[3]=( 12.0- 9.0*B- 6.0*C)/6.0;
    resize_filter->cubic[4]=(       8.0*B+24.0*C)/6.0;
    resize_filter->cubic[5]=(     -12.0*B-48.0*C)/6.0;
    resize_filter->cubic[6]=(       6.0*B+30.0*C)/6.0;
    resize_filter->cubic[7]=(     - 1.0*B- 6.0*C)/6.0;
  }
  artifact=GetImageArtifact(image,"filter:verbose");
  if (artifact != (const char *) NULL)
    {
      double
        support,
        x;

      /*
        Output filter graph -- for graphing filter result.
      */
      support=GetResizeFilterSupport(resize_filter);
      (void) printf("# support = %lg\n",support);
      for (x=0.0; x <= support; x+=0.01f)
        (void) printf("%5.2lf\t%lf\n",x,GetResizeFilterWeight(resize_filter,x));
      (void) printf("%5.2lf\t%lf\n",support,0.0);
    }
  return(resize_filter);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   A d a p t i v e R e s i z e I m a g e                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AdaptiveResizeImage() adaptively resize image with pixel resampling.
%
%  The format of the AdaptiveResizeImage method is:
%
%      Image *AdaptiveResizeImage(const Image *image,
%        const unsigned long columns,const unsigned long rows,
%        ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o columns: the number of columns in the resized image.
%
%    o rows: the number of rows in the resized image.
%
%    o exception: return any errors or warnings in this structure.
%
*/
MagickExport Image *AdaptiveResizeImage(const Image *image,
  const unsigned long columns,const unsigned long rows,ExceptionInfo *exception)
{
#define AdaptiveResizeImageTag  "Resize/Image"

  Image
    *resize_image;

  long
    y;

  MagickBooleanType
    proceed;

  MagickPixelPacket
    pixel;

  PointInfo
    offset;

  register IndexPacket
    *resize_indexes;

  register long
    x;

  register PixelPacket
    *q;

  ResampleFilter
    *resample_filter;

  ViewInfo
    *resize_view;

  /*
    Adaptively resize image.
  */
  assert(image != (const Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  assert(exception != (ExceptionInfo *) NULL);
  assert(exception->signature == MagickSignature);
  if ((columns == 0) || (rows == 0))
    return((Image *) NULL);
  if ((columns == image->columns) && (rows == image->rows))
    return(CloneImage(image,0,0,MagickTrue,exception));
  resize_image=CloneImage(image,columns,rows,MagickTrue,exception);
  if (resize_image == (Image *) NULL)
    return((Image *) NULL);
  if (SetImageStorageClass(resize_image,DirectClass) == MagickFalse)
    {
      InheritException(exception,&resize_image->exception);
      resize_image=DestroyImage(resize_image);
      return((Image *) NULL);
    }
  GetMagickPixelPacket(image,&pixel);
  resample_filter=AcquireResampleFilter(image,exception);
  if (image->interpolate == UndefinedInterpolatePixel)
    (void) SetResampleFilterInterpolateMethod(resample_filter,
      MeshInterpolatePixel);
  resize_view=AcquireCacheView(resize_image);
  for (y=0; y < (long) resize_image->rows; y++)
  {
    q=QueueCacheViewAuthenticPixels(resize_view,0,y,resize_image->columns,1,
      exception);
    if (q == (PixelPacket *) NULL)
      break;
    resize_indexes=GetCacheViewAuthenticIndexQueue(resize_view);
    offset.y=((MagickRealType) y*image->rows/resize_image->rows);
    for (x=0; x < (long) resize_image->columns; x++)
    {
      offset.x=((MagickRealType) x*image->columns/resize_image->columns);
      (void) ResamplePixelColor(resample_filter,offset.x-0.5,offset.y-0.5,
        &pixel);
      SetPixelPacket(resize_image,&pixel,q,resize_indexes+x);
      q++;
    }
    if (SyncCacheViewAuthenticPixels(resize_view,exception) == MagickFalse)
      break;
    proceed=SetImageProgress(image,AdaptiveResizeImageTag,y,image->rows);
    if (proceed == MagickFalse)
      break;
  }
  resample_filter=DestroyResampleFilter(resample_filter);
  resize_view=DestroyCacheView(resize_view);
  return(resize_image);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   B e s s e l O r d e r O n e                                               %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  BesselOrderOne() computes the Bessel function of x of the first kind of
%  order 0:
%
%    Reduce x to |x| since j1(x)= -j1(-x), and for x in (0,8]
%
%       j1(x) = x*j1(x);
%
%    For x in (8,inf)
%
%       j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
%
%    where x1 = x-3*pi/4. Compute sin(x1) and cos(x1) as follow:
%
%       cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
%               =  1/sqrt(2) * (sin(x) - cos(x))
%       sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
%               = -1/sqrt(2) * (sin(x) + cos(x))
%
%  The format of the BesselOrderOne method is:
%
%      MagickRealType BesselOrderOne(MagickRealType x)
%
%  A description of each parameter follows:
%
%    o x: MagickRealType value.
%
*/

#undef I0
static MagickRealType I0(MagickRealType x)
{
  MagickRealType
    sum,
    t,
    y;

  register long
    i;

  /*
    Zeroth order Bessel function of the first kind.
  */
  sum=1.0;
  y=x*x/4.0;
  t=y;
  for (i=2; t > MagickEpsilon; i++)
  {
    sum+=t;
    t*=y/((MagickRealType) i*i);
  }
  return(sum);
}

#undef J1
static MagickRealType J1(MagickRealType x)
{
  MagickRealType
    p,
    q;

  register long
    i;

  static const double
    Pone[] =
    {
       0.581199354001606143928050809e+21,
      -0.6672106568924916298020941484e+20,
       0.2316433580634002297931815435e+19,
      -0.3588817569910106050743641413e+17,
       0.2908795263834775409737601689e+15,
      -0.1322983480332126453125473247e+13,
       0.3413234182301700539091292655e+10,
      -0.4695753530642995859767162166e+7,
       0.270112271089232341485679099e+4
    },
    Qone[] =
    {
      0.11623987080032122878585294e+22,
      0.1185770712190320999837113348e+20,
      0.6092061398917521746105196863e+17,
      0.2081661221307607351240184229e+15,
      0.5243710262167649715406728642e+12,
      0.1013863514358673989967045588e+10,
      0.1501793594998585505921097578e+7,
      0.1606931573481487801970916749e+4,
      0.1e+1
    };

  p=Pone[8];
  q=Qone[8];
  for (i=7; i >= 0; i--)
  {
    p=p*x*x+Pone[i];
    q=q*x*x+Qone[i];
  }
  return(p/q);
}

#undef P1
static MagickRealType P1(MagickRealType x)
{
  MagickRealType
    p,
    q;

  register long
    i;

  static const double
    Pone[] =
    {
      0.352246649133679798341724373e+5,
      0.62758845247161281269005675e+5,
      0.313539631109159574238669888e+5,
      0.49854832060594338434500455e+4,
      0.2111529182853962382105718e+3,
      0.12571716929145341558495e+1
    },
    Qone[] =
    {
      0.352246649133679798068390431e+5,
      0.626943469593560511888833731e+5,
      0.312404063819041039923015703e+5,
      0.4930396490181088979386097e+4,
      0.2030775189134759322293574e+3,
      0.1e+1
    };

  p=Pone[5];
  q=Qone[5];
  for (i=4; i >= 0; i--)
  {
    p=p*(8.0/x)*(8.0/x)+Pone[i];
    q=q*(8.0/x)*(8.0/x)+Qone[i];
  }
  return(p/q);
}

#undef Q1
static MagickRealType Q1(MagickRealType x)
{
  MagickRealType
    p,
    q;

  register long
    i;

  static const double
    Pone[] =
    {
      0.3511751914303552822533318e+3,
      0.7210391804904475039280863e+3,
      0.4259873011654442389886993e+3,
      0.831898957673850827325226e+2,
      0.45681716295512267064405e+1,
      0.3532840052740123642735e-1
    },
    Qone[] =
    {
      0.74917374171809127714519505e+4,
      0.154141773392650970499848051e+5,
      0.91522317015169922705904727e+4,
      0.18111867005523513506724158e+4,
      0.1038187585462133728776636e+3,
      0.1e+1
    };

  p=Pone[5];
  q=Qone[5];
  for (i=4; i >= 0; i--)
  {
    p=p*(8.0/x)*(8.0/x)+Pone[i];
    q=q*(8.0/x)*(8.0/x)+Qone[i];
  }
  return(p/q);
}

static MagickRealType BesselOrderOne(MagickRealType x)
{
  MagickRealType
    p,
    q;

  if (x == 0.0)
    return(0.0);
  p=x;
  if (x < 0.0)
    x=(-x);
  if (x < 8.0)
    return(p*J1(x));
  q=sqrt((double) (2.0/(MagickPI*x)))*(P1(x)*(1.0/sqrt(2.0)*(sin((double) x)-
    cos((double) x)))-8.0/x*Q1(x)*(-1.0/sqrt(2.0)*(sin((double) x)+
    cos((double) x))));
  if (p < 0.0)
    q=(-q);
  return(q);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   D e s t r o y R e s i z e F i l t e r                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  DestroyResizeFilter() destroy the resize filter.
%
%  The format of the AcquireResizeFilter method is: