Cheaper sinc computation for resize.c

[NicolasRobidoux]

I am an applied mathematician working on image filtering. For example, I am a regular contributor to the free open source VIPS/NIP2 system, and me and my students have contributed to the free open source GEGL library.

My masters student Chantal Racette and I have just discovered an highly accurate approximation method for the sinc (and related) functions on finite ranges.

I believe that this could be useful for your recize.c code, and would be delighted to see it integrated in IM right away. (We'll write the article later.)

Are you interested?

To give you an idea of what I'm talking about, here is an single piece polynomial approximation with 1e-8 maximum relative error over the interval [-3,3] (note that getting a low relative error is not trivial because the sinc function has zeros in the interval), which is a direct replacement for the corresponding Sinc function in resize.c:

static MagickRealType Sinc(const MagickRealType x,
const ResizeFilter *magick_unused(resize_filter))
{
#define FAST_LANCZOS_COEF_1 -0.2777777749828003702607348742700210983971e-1
#define FAST_LANCZOS_COEF_2 0.7883967237692122941462552028914990810270e-2
#define FAST_LANCZOS_COEF_3 -0.1014962147766290313113433635286152701178e-3
#define FAST_LANCZOS_COEF_4 0.7957122890011771235521289628748676216784e-4
#define FAST_LANCZOS_COEF_5 -0.4297822393998422516448009437099172798829e-5
#define FAST_LANCZOS_COEF_6 0.1707846962680966498461068335001167780751e-6
#define FAST_LANCZOS_COEF_7 -0.5110955098479687596630062395922169657449e-8
#define FAST_LANCZOS_COEF_8 0.1091121789464303759146861882233623802483e-9
#define FAST_LANCZOS_COEF_9 -0.1270023023877167259514316150814763206939e-11
const double x_squared = x*x;
if (x_squared <= 1.0)
{
const double p =
FAST_LANCZOS_COEF_1 + x_squared * (
FAST_LANCZOS_COEF_2 + x_squared * (
FAST_LANCZOS_COEF_3 + x_squared * (
FAST_LANCZOS_COEF_4 + x_squared * (
FAST_LANCZOS_COEF_5 + x_squared * (
FAST_LANCZOS_COEF_6 + x_squared * (
FAST_LANCZOS_COEF_7 + x_squared * (
FAST_LANCZOS_COEF_8 + x_squared * (
FAST_LANCZOS_COEF_9 ))))))));
const double s =
( (x_squared + -1.0) * (x_squared + -4.0) )
*
( (x_squared + -9.0) * p );
return(s);
}
else
return(sin(MagickPI*x)/(MagickPI*x));
}

(I have not taken the time to triple check this code. Apologies if it is not quite right.)
We can produce one-piece or piecewise polynomial approximations of arbitrary accuracy, and that we minimize the relative error because we have found that, in general, for an equal computation cost, the minimizer of relative error has a frequency response closer to the original than the minimizer of absolute error when used as a filter kernel.)

Basically, you tell us what you'd like, we compute for a bit, and you get the code.

We also have highly accurate approximations of the lanczos functions. Generally, 30 or so flops gives machine precision.

Nicolas Robidoux
Department of Mathematics and Computer Science
Universite Laurentienne/Laurentian University
Sudbuy Ontario Canada

[magick]

I believe that this could be useful for your resize.c code, and would be delighted to see it integrated in IM right away.

We're always interested in contributions to ImageMagick, particularly if they speed up computations. You must accept the ImageMagick license, particularly the Submission of Contribution clause @ http://www.imagemagick.org/script/licen ... tributions in order for us to to include your work in the main ImageMagick distribution.

We'll take a look the your Sinc function over the next few days. You can submit the image filters individually or you could submit your entire contribution as a patch. Just post a URL to the patch that we can download and apply against the Subversion trunk.

Anthony, do you have time to review the proposed Sinc function and ensure its a suitable replacement inside resize.c?

[NicolasRobidoux]

I have put a modified resize.c with three alternate sinc( functions commented out with #if 0 ... #endif here:

http://web.cs.laurentian.ca/nrobidoux/misc/resize.c

The key modified code follows at the end of this post.

All three sinc approximations approximate sinc over [-3,3] because it appears to me that this is interval most used (because of lanczos 3-lobes).

As mentioned earlier, I can approximate over other intervals, I can approximate using piecewise polynomials instead of one piece over the interval of interest, and I can approximate other functions.

Nicolas Robidoux

Here is the code:

#if 0

static MagickRealType Sinc(const MagickRealType x,
const ResizeFilter *magick_unused(resize_filter))
{
// The following approximation of the sinc function over the
// interval [-3,3], constructed by Nicolas Robidoux with the
// assistance of Chantal Racette, has maximum relative error of
// 1.8e-14 over that interval.
const double x_squared = (double) x * (double) x;
if (x_squared <= 1.0)
{
const double p =
-0.2777777777777777484967514793229935886858e-1 + x_squared * (
0.7883970992697542830560986463962874704581e-2 + x_squared * (
-0.1014971048673801299627391800253629784307e-2 + x_squared * (
0.7957975170651347522690517580617391446624e-4 + x_squared * (
-0.4302064159535280660121265849050342513500e-5 + x_squared * (
0.1720088580385734039018550565889345187441e-6 + x_squared * (
-0.5325247340694045282279767977535166005127e-8 + x_squared * (
0.1318969739543994035763170253658588795558e-9 + x_squared * (
-0.2678779601269262783351382208935247318217e-11 + x_squared * (
0.4547427567791873812472026076072662961706e-13 + x_squared * (
-0.6542016827740848562863869296875889077432e-15 + x_squared * (
0.7978989414064956224188997380470875190815e-17 + x_squared * (
-0.7824271041263904097239661475065843020761e-19 + x_squared * (
0.4764330182417253304763220033014732190061e-21 )))))))))))));
const double s =
( ( x_squared + -1.0 ) * ( x_squared + -4.0 ) )
*
( ( x_squared + -9.0 ) * p );
return(s);
}
else
return(sin(MagickPI*x)/(MagickPI*x));
}

#endif

#if 0
static MagickRealType Sinc(const MagickRealType x,
const ResizeFilter *magick_unused(resize_filter))
{
// The following approximation of the sinc function over the
// interval [-3,3], constructed by Nicolas Robidoux with the
// assistance of Chantal Racette, has maximum relative error of
// 1.0e-8 over that interval.
const double x_squared = (double) x * (double) x;
if (x_squared <= 1.0)
{
const double p =
-0.2777777749828003702607348742700210983971e-1 + x_squared * (
0.7883967237692122941462552028914990810270e-2 + x_squared * (
-0.1014962147766290313113433635286152701178e-3 + x_squared * (
0.7957122890011771235521289628748676216784e-4 + x_squared * (
-0.4297822393998422516448009437099172798829e-5 + x_squared * (
0.1707846962680966498461068335001167780751e-6 + x_squared * (
-0.5110955098479687596630062395922169657449e-8 + x_squared * (
0.1091121789464303759146861882233623802483e-9 + x_squared * (
-0.1270023023877167259514316150814763206939e-11 ))))))));
const double s =
( ( x_squared + -1.0 ) * ( x_squared + -4.0 ) )
*
( ( x_squared + -9.0 ) * p );
return(s);
}
else
return(sin(MagickPI*x)/(MagickPI*x));
}
#endif

#if 0
static MagickRealType Sinc(const MagickRealType x,
const ResizeFilter *magick_unused(resize_filter))
{
// The following approximation of the sinc function over the
// interval [-3,3], constructed by Nicolas Robidoux with the
// assistance of Chantal Racette, has maximum relative error of
// 8.8e-5 over that interval.
const double x_squared = (double) x * (double) x;
if (x_squared <= 1.0)
{
const double p =
-0.2777533762119262916850997993215487249694e-1 + x_squared * (
0.7871060719417389265828115362113334868483e-2 + x_squared * (
-0.1002039946063031209329660819970837885402e-2 + x_squared * (
0.7437810130466799848673421417105135076745e-4 + x_squared * (
-0.3275486331327512805059915658357413415728e-5 + x_squared * (
0.6742010644073333825302432891245471531150e-7 )))));
const double s =
( ( x_squared + -1.0 ) * ( x_squared + -4.0 ) )
*
( ( x_squared + -9.0 ) * p );
return(s);
}
else
return(sin(MagickPI*x)/(MagickPI*x));
}
#endif

#if 1
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));
}
#endif

[NicolasRobidoux]

Apologies for the multiple posts: the x_squared <= 1.0 tests were wrong. (Still gave correct results, so was hard to detect.)

I fixed the version posted at

http://web.cs.laurentian.ca/nrobidoux/misc/resize.c

[magick]

Got it. We'll review the changes over the next few days. You of course are not restricted in anyway with your code. You can offer it to any other project you want. The license gives us rights to utilize your code to use within ImageMagick.

On behalf of the ImageMagick user community, thank you for your contribution.

[magick]

Grab ImageMagick-6.6.3-10 Beta which includes your patch. We still need to do more testing before we officially accept the patch. In the mean-time if you have any other filter approximations, feel free to post them here.

[NicolasRobidoux]

Thank for you the (tentative) inclusion.

------

There is something I oversaw in the code I gave, which is a result of having never worked on filtering where the widest sinc extent is, at least in principle, arbitrarily large (as is the case in IM).

The long and short of it is that because the approximation will sometimes be used in IM with values of x outside of [-3,3], I need to match the approximation inside the interval (where the polynomial is used) with the outside values (where sin(pi x)/(pi x) is used) more carefully.

I'll update the code when me and Chantal have completed the necessary computation.

Note that you may not find a measurable difference, at least with the most precise polynomial approximation (the one which is not #if 0'd out). But I know it's there, so I'll fix it.

(Note: If you want me to explain the math, let me know. For now, I'll assume that you only care about the code and the results it gives. One side effect of the change will be that p(0) will be exactly 1.)

[magick]

We'd like to see an short explanation of the math posted here. Let us know when your corrections are complete for the Sinc approximation. We'll hold off on the next point release of ImageMagick until Sinc() is complete and tested.

[NicolasRobidoux]

Two questions:

1) Is MagickRealType ever larger (in terms of digits in the significand a.k.a. mantissa) than a standard double?

2) Is MagickRealType ever smaller than a standard double?

Are there guidelines somewhere RE: whether computations should be done with double or MagickRealType?

[magick]

MagickRealType is double unless the QuantumDepth is set to 32 (rarely) at build time. In that case its long double. The default QuantumDepth is 16 bits per pixel per pixel component.

[fmw42]

MagickRealType is double unless the QuantumDepth is set to 32 (rarely) at build time. In that case its long double. The default QuantumDepth is 16 bits per pixel per pixel component.

Magick,

Sorry to stick my nose in, but what about under HDRI mode? Any potential issues and does that have any bearing on the MagickRealType?

[fmw42]

MagickRealType is double unless the QuantumDepth is set to 32 (rarely) at build time. In that case its long double. The default QuantumDepth is 16 bits per pixel per pixel component.

Magick,

Sorry to stick my nose in. I am not sure this is relevant, but what about under HDRI mode? Any potential issues and does that have any bearing on the MagickRealType?

[magick]

MagickRealType is not affected by HDRI.

[NicolasRobidoux]

ETA: Thursday at the latest.

Nicolas Robidoux

[fmw42]

Forgive me again for sticking my nose in.

I would suggest that it might be useful for evaluation purposes to see a graph showing the difference in the "fast" sinc from the more formal "exact" one and a difference graph and error measure. Thus, if differences are significant enough (or results differ when using in -resize), one could offer it as an alternate sinc filter for -resize so that users can choose speed over quality.

Anthony has graphed the various sinc-like filters that are used with -resize. See http://www.imagemagick.org/Usage/resize ... r_windowed as well as the windowing functions slightly below that section and other filters above that.