Cheaper sinc computation for resize.c

[anthony]

Filters generally have a "normalized scaling most appropriate for translations/rotations" which basically is the "support" variable.

If you use the "convolution" as opposed to the "interpolation" approach when enlarging (more specifically, when locally the Jacobian matrix of the geometrical transformation has at least one singular value larger than 1) you should clamp UP the width/height of the "convolution mask" so it is AT LEAST twice the support. If you do that, you do not need to "switch" between EWA and Interpolation (unless, of course, you find that interpolation is better).

This being said, for many methods (at least those for which the weights are normalized, either before or during the "interpolation"), this clamping is exactly equivalent to using the kernel as if it was the cardinal basis function for an interpolatory method. So, maybe you are already doing the above.

yes I noticed this with Resize. It 'clamps up' as you put it when magnification becomes involved. As such the same filter is used for both magnification and minification. I think the terms used in research papers was decimation, and reconstruction filters.

In EWA filtering Distort it is a little harder as a direct implantation could generate very long thin ellipses that just 'miss' every source pixel! The solution to that was to ensure ellipses have a minimum thickness.

That same technique however was expanded by Paul Heckbert in a 'High Quality EWA'.
What he did was modify the ellipse so it has a minimal size, with a supposedly smooth transition between the two 'modes' .

The problem is I think I may have implemented it wrong, resulting in a bad ellipse size when magnification becomes involved. At the time I did not understand filtering properly enough to fix it (it was before my research for resize filters that started soon after), and the problem was being hidden by a coping of code from another badly implemented EWA, and hidden by that same codes over blurry Gaussian filter.

Basically the result was the two wrongs cancelled each other out to an extent, and it has only appeared not that I attempted to fix the over-blurriness of distorts, and provided a method to turn off the switch to interpolation for areas of image enlargements. Untill I get this fix I am holding off updating distorts.

Enlarging an image is a good way of checking how well the filter is working as you can directly see the convolution effects of the filter in extreme enlargements. (Ringing, blocking, filter support clipping, etc).
As is demonstrated in the the IM Examples, Resize page.


For example...

Code: Select all

  convert -size 1x1 xc: -bordercolor '#444' -border 4x4 \
          -set option:filter:filter Sinc -resize 100x100\! \
          dot_sinc.gif

Clearly shows the 'ringing effects of a raw Sinc on a enlargement of a single pixel.

Repeating this with distort (before my changes)

Code: Select all

  convert -size 1x1 xc: -bordercolor '#444' -border 4x4 \
          -set option:filter:filter Sinc +distort SRT 12,0 \
          dot_distort.gif

shows the result of a bilinear repeating this however with interpolation turned off to force EWA filtering
however results in...

which is clearly WRONG! I Tracking down the problem.

UPDATE: The problem was figured out and fixed. See
viewtopic.php?f=2&t=17061

[NicolasRobidoux]

I'm an idiot:

Of course a two-term recursion should be extended with a loop. What am I doing doing a dozen nested ifs?

Fixing LanczosChebyshev when I have a minute.

nicolas

[anthony]

On my problems with distort. I was right the ellipse support radius was wrong.
I have now fixed this and suddenly most of the original blurriess, and later aliasing has disappeared.
Looks like two wrongs can sometimes make a "mostly right".

However any filter involving negative weighting is still not working as expected. Working on it.
UPDATE: Also now fixed

[anthony]

I'm an idiot:

Of course a two-term recursion should be extended with a loop. What am I doing doing a dozen nested ifs?

Fixing LanczosChebyshev when I have a minute.

nicolas

I sort of wondered about that. Especially a loop being feed with the array of constants you are using.

[anthony]

Note for people looking at this thread in the future

The Polynomial Sinc function is now accessable as SincFast, and is also now the internal default used for Sinc Windowing functions, including lanczos. However as discussed above, this only produces a very minor speed up for the same result due to filter caching by the Resize and Distort image processing operations.

The LancosChebyshev function is now known as LanczosFast, and at the time of writing (IM v6.6.4-5) is still classed as experimental, and as such is not the default for the "Lanczos" filter.
Lanczos is thus currently definined as a "SincFast windowed SincFast" filter for image resize, and a "Bessel windowed Bessel" filter for image distort.

You can use a -set option:filter:verbose 1 expert setting to see the actual functions internally defined according to the user supplied filter and expert settings given.
See Resize Filters, Expert Settings...
http://www.imagemagick.org/Usage/resize/#filter_options

For example...

Code: Select all

  convert xc: -filter Lanczos -set option:filter:verbose 1 -resize 200% null:  | grep '^#'

#
# Resize Filter (for graphing)
#
# filter = SincFast
# window = SincFast
# support = 3
# win-support = 3
# blur = 1
# blurred_support = 3
# B,C = 0,0
#

[NicolasRobidoux]

As it turns out, computing lanczos(n) = n sin(pi x) sin(pi x/n) / (pi x)^2 using the Chebyshev recursive formula for sines of multiple of an angle and/or high accuracy polynomial approximations of cos(pi x) looks good on paper, but actually SincFast is so, well, fast, that the recursion loses.

The final version of LanczosFast is found in ImageMagick svn version 2683. It runs just fine; it's just a little slower than stock Lanczos based on SincFast.

LanczosFast is gone as of version 2684.

I'm not totally clear on why LanczosFast runs more slowly (!). I'm guessing it's because of the while loop (and additional branching).

For reference, here is the code:

Code: Select all

static MagickRealType LanczosFast(const MagickRealType x,
  const ResizeFilter *resize_filter)
{
  /*
    WARNING:
    LanczosFast only outputs correct values if support is a
    POSITIVE INTEGER.
  */
  /*
    Computing the Lanczos kernel directly from its definition

      sinc(x)*sinc(x/n) = sin(pi*x)*sin(pi*x/n)/((pi*x)*(pi*x/n))

    requires two calls to a trigonometric function (slow).

    LanczosFast uses a recursive formula for sin(x)*sin(x/n) (n a
    positive integer) based on the Chebyshev method for the
    computation of sines and cosines of multiples of an angle:
    http://en.wikipedia.org/wiki/List_of_trigonometric_identities...
    ...#Chebyshev_method

    The two-term recursion allows computing the needed product of
    sines with only one trig call (cos(pi*x/n)). It is, however, only
    applicable when n=support is a positive integer (the only truly
    useful case).  If not, the code reverts to the usual
    product-of-(fast)sincs formula. It also reverts to the usual if
    x>support (an expert-only possibility).

    In order to eliminate the remaining trig call, high precision
    "minimax" polynomial approximants are used. See the SincFast
    comments for details.

    Recursive formula for the product of the sine of an angle and the
    sine of a multiple of the angle discovered by Nicolas Robidoux
    (pending the discovery of an earlier discoverer) with the
    assistance of Chantal Racette.  Approximations of cos(pi x) over
    the interval [-1,1] constructed by Nicolas Robidoux and Chantal
    Racette with funding from the Natural Sciences and Engineering
    Research Council of Canada.
  */
  /*
    We assume that support > 0 and x >= 0.
  */
  const MagickRealType support = resize_filter->support;
  const MagickRealType xos = x/support;
  const MagickRealType xx = xos*xos;
  const MagickRealType pi2xx = (MagickPIL*MagickPIL)*x*x;
  if (pi2xx == (MagickRealType) 0.0)
    return(1.0);
  {
    MagickRealType c;
    if (xos>1.0)
      {
	c = (MagickRealType) cos((double) (MagickPIL*xos));
      }
    else
      {
#if MAGICKCORE_QUANTUM_DEPTH <= 8
	/*
	  Maximum absolute relative error 1.8e-8 < 1/2^25.
	*/
	const MagickRealType c0 =
	  3.99999992643142920211338040253631444680L;
	const MagickRealType c1 =
	  -3.73920425616710679261196824501740620846L;
	const MagickRealType c2 =
	  1.27796557877584024350399820493902002676L;
	const MagickRealType c3 =
	  -0.228809086406287387417934333009611810036L;
	const MagickRealType c4 =
	  0.249852427433082398345689482146854999198e-1L;
	const MagickRealType c5 =
	  -0.160409656670710471758484348423047271587e-2L;
	const MagickRealType p =
	  c0+xx*(c1+xx*(c2+xx*(c3+xx*(c4+xx*c5))));
#elif MAGICKCORE_QUANTUM_DEPTH <= 16
	/*
	  Max. abs. rel. error 1.6e-12 < 1/2^39.
	*/
	const MagickRealType c0 =
	  3.99999999999369507875837503656312381210L;
	const MagickRealType c1 =
	  -3.73920880146193419336832119898956139180L;
	const MagickRealType c2 =
	  1.27801328321259865639714208686127144980L;
	const MagickRealType c3 =
	  -0.228997785334013943051673846919195886430L;
	const MagickRealType c4 =
	  0.253305916511203326575532808976653777888e-1L;
	const MagickRealType c5 =
	  -0.190290817480650124504465021065476218977e-2L;
	const MagickRealType c6 =
	  0.102686888764524439708563743432769693793e-3L;
	const MagickRealType c7 =
	  -0.373344419226166828092642504009826630593e-5L;
	const MagickRealType p =
	  c0+xx*(c1+xx*(c2+xx*(c3+xx*(c4+xx*(c5+xx*(c6+xx*c7))))));
#else
	/*
	  Max. abs. rel. error 6.1e-16 < 1/2^50 if computed with
	  "true" long doubles, 5.6e-17 < 1/2^53 if long doubles are
	  IEEE doubles.
	*/
	const MagickRealType c0 =
	  3.99999999999999977653487536495855668234L;
	const MagickRealType c1 =
	  -3.73920880217867665744995903521564136498L;
	const MagickRealType c2 =
	  1.27801329695096620491074078493372289128L;
	const MagickRealType c3 =
	  -0.228997887594702717198738793082652705850L;
	const MagickRealType c4 =
	  0.253309709061066296103141786511835130143e-1L;
	const MagickRealType c5 =
	  -0.190368124507833605338686581554879104404e-2L;
	const MagickRealType c6 =
	  0.103570244775199129899717783834074546882e-3L;
	const MagickRealType c7 =
	  -0.426762781982225449936672742185198500277e-5L;
	const MagickRealType c8 =
	  0.137089756109639545966503908032203389025e-6L;
	const MagickRealType c9 =
	  -0.321199312802383349700375491874286104564e-8L;
	const MagickRealType p =
	  c0+xx*(c1+xx*(c2+xx*(c3+xx*(c4+xx*(c5+xx*(c6+xx*(c7+xx*(c8+xx*c9
	  ))))))));
#endif    
	c = (0.25-xx)*p;
      }
    {
      const MagickRealType cpc = c+c;
      MagickRealType ss = 1.0-c*c;
      MagickRealType n = support-1.0;
      MagickRealType ss1 = 0.0;
      MagickRealType temp;
      while (n>0.0)
	{
	  temp = ss;
	  ss = cpc*temp-ss1;
	  ss1 = temp;
	  n = n-1.0;
	}
      return(support/pi2xx*ss);
  } }
}

The part of the code which uses a (new?) recursive formula for Lanczos n is the very last. (The part that starts right after the #endif.)
c is the cosine of pi x / n. From it, the product of sines ss is computed by looping. Finally, it's combined with the n/(pi x)^2 factor.

I am GUESSING that for a FIXED n, the recursive formula combined with a high efficiency cosine approximation may win. (Although, if n=3 for example, you are actually better off using

sin(t)*sin(3t) = sin^2(t)*(3-4*sin^2(t)) = -1 + c^2 ( 5 - 4 c^2 ) with c = cos(t)

than using the recursion directly.)

But I don't think this is something really useful for ImageMagick.

[anthony]

LanczosFast is gone as of version 2684.

Update to SVN version 2685...

Removed LanczosFast from filter types and command line option parsing
"resample.h" and "option.c" respecivally

Added comment in ChangeLog for svn 2686

[NicolasRobidoux]

LanczosFast is gone as of version 2684.

Update to SVN version 2685...

Removed LanczosFast from filter types and command line option parsing
"resample.h" and "option.c" respecivally

Added comment in ChangeLog for svn 2686

Apologies (I know better): I forgot to do a "find" on every occurence of LanczosFast in the library so as to remove it everywhere.

[anthony]

That is fine. I myself often do not know about things needed by other API's (like Perl, Ruby, etc). Usually they are picked up later by others more intimately involved with those API's.

This is part of the beauty of Open Source.