Cheaper sinc computation for resize.c

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Cheaper sinc computation for resize.c
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 1e8 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.2777777749828003702607348742700210983971e1
#define FAST_LANCZOS_COEF_2 0.7883967237692122941462552028914990810270e2
#define FAST_LANCZOS_COEF_3 0.1014962147766290313113433635286152701178e3
#define FAST_LANCZOS_COEF_4 0.7957122890011771235521289628748676216784e4
#define FAST_LANCZOS_COEF_5 0.4297822393998422516448009437099172798829e5
#define FAST_LANCZOS_COEF_6 0.1707846962680966498461068335001167780751e6
#define FAST_LANCZOS_COEF_7 0.5110955098479687596630062395922169657449e8
#define FAST_LANCZOS_COEF_8 0.1091121789464303759146861882233623802483e9
#define FAST_LANCZOS_COEF_9 0.1270023023877167259514316150814763206939e11
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 onepiece 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
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 1e8 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.2777777749828003702607348742700210983971e1
#define FAST_LANCZOS_COEF_2 0.7883967237692122941462552028914990810270e2
#define FAST_LANCZOS_COEF_3 0.1014962147766290313113433635286152701178e3
#define FAST_LANCZOS_COEF_4 0.7957122890011771235521289628748676216784e4
#define FAST_LANCZOS_COEF_5 0.4297822393998422516448009437099172798829e5
#define FAST_LANCZOS_COEF_6 0.1707846962680966498461068335001167780751e6
#define FAST_LANCZOS_COEF_7 0.5110955098479687596630062395922169657449e8
#define FAST_LANCZOS_COEF_8 0.1091121789464303759146861882233623802483e9
#define FAST_LANCZOS_COEF_9 0.1270023023877167259514316150814763206939e11
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 onepiece 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
Re: Cheaper sinc computation for resize.c
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.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'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?

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: Cheaper sinc computation for resize.c
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 3lobes).
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.8e14 over that interval.
const double x_squared = (double) x * (double) x;
if (x_squared <= 1.0)
{
const double p =
0.2777777777777777484967514793229935886858e1 + x_squared * (
0.7883970992697542830560986463962874704581e2 + x_squared * (
0.1014971048673801299627391800253629784307e2 + x_squared * (
0.7957975170651347522690517580617391446624e4 + x_squared * (
0.4302064159535280660121265849050342513500e5 + x_squared * (
0.1720088580385734039018550565889345187441e6 + x_squared * (
0.5325247340694045282279767977535166005127e8 + x_squared * (
0.1318969739543994035763170253658588795558e9 + x_squared * (
0.2678779601269262783351382208935247318217e11 + x_squared * (
0.4547427567791873812472026076072662961706e13 + x_squared * (
0.6542016827740848562863869296875889077432e15 + x_squared * (
0.7978989414064956224188997380470875190815e17 + x_squared * (
0.7824271041263904097239661475065843020761e19 + x_squared * (
0.4764330182417253304763220033014732190061e21 )))))))))))));
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.0e8 over that interval.
const double x_squared = (double) x * (double) x;
if (x_squared <= 1.0)
{
const double p =
0.2777777749828003702607348742700210983971e1 + x_squared * (
0.7883967237692122941462552028914990810270e2 + x_squared * (
0.1014962147766290313113433635286152701178e3 + x_squared * (
0.7957122890011771235521289628748676216784e4 + x_squared * (
0.4297822393998422516448009437099172798829e5 + x_squared * (
0.1707846962680966498461068335001167780751e6 + x_squared * (
0.5110955098479687596630062395922169657449e8 + x_squared * (
0.1091121789464303759146861882233623802483e9 + x_squared * (
0.1270023023877167259514316150814763206939e11 ))))))));
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.8e5 over that interval.
const double x_squared = (double) x * (double) x;
if (x_squared <= 1.0)
{
const double p =
0.2777533762119262916850997993215487249694e1 + x_squared * (
0.7871060719417389265828115362113334868483e2 + x_squared * (
0.1002039946063031209329660819970837885402e2 + x_squared * (
0.7437810130466799848673421417105135076745e4 + x_squared * (
0.3275486331327512805059915658357413415728e5 + x_squared * (
0.6742010644073333825302432891245471531150e7 )))));
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 Xscaled Sinc(x) function.
*/
if (x == 0.0)
return(1.0);
return(sin(MagickPI*(double) x)/(MagickPI*(double) x));
}
#endif
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 3lobes).
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.8e14 over that interval.
const double x_squared = (double) x * (double) x;
if (x_squared <= 1.0)
{
const double p =
0.2777777777777777484967514793229935886858e1 + x_squared * (
0.7883970992697542830560986463962874704581e2 + x_squared * (
0.1014971048673801299627391800253629784307e2 + x_squared * (
0.7957975170651347522690517580617391446624e4 + x_squared * (
0.4302064159535280660121265849050342513500e5 + x_squared * (
0.1720088580385734039018550565889345187441e6 + x_squared * (
0.5325247340694045282279767977535166005127e8 + x_squared * (
0.1318969739543994035763170253658588795558e9 + x_squared * (
0.2678779601269262783351382208935247318217e11 + x_squared * (
0.4547427567791873812472026076072662961706e13 + x_squared * (
0.6542016827740848562863869296875889077432e15 + x_squared * (
0.7978989414064956224188997380470875190815e17 + x_squared * (
0.7824271041263904097239661475065843020761e19 + x_squared * (
0.4764330182417253304763220033014732190061e21 )))))))))))));
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.0e8 over that interval.
const double x_squared = (double) x * (double) x;
if (x_squared <= 1.0)
{
const double p =
0.2777777749828003702607348742700210983971e1 + x_squared * (
0.7883967237692122941462552028914990810270e2 + x_squared * (
0.1014962147766290313113433635286152701178e3 + x_squared * (
0.7957122890011771235521289628748676216784e4 + x_squared * (
0.4297822393998422516448009437099172798829e5 + x_squared * (
0.1707846962680966498461068335001167780751e6 + x_squared * (
0.5110955098479687596630062395922169657449e8 + x_squared * (
0.1091121789464303759146861882233623802483e9 + x_squared * (
0.1270023023877167259514316150814763206939e11 ))))))));
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.8e5 over that interval.
const double x_squared = (double) x * (double) x;
if (x_squared <= 1.0)
{
const double p =
0.2777533762119262916850997993215487249694e1 + x_squared * (
0.7871060719417389265828115362113334868483e2 + x_squared * (
0.1002039946063031209329660819970837885402e2 + x_squared * (
0.7437810130466799848673421417105135076745e4 + x_squared * (
0.3275486331327512805059915658357413415728e5 + x_squared * (
0.6742010644073333825302432891245471531150e7 )))));
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 Xscaled Sinc(x) function.
*/
if (x == 0.0)
return(1.0);
return(sin(MagickPI*(double) x)/(MagickPI*(double) x));
}
#endif

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: Cheaper sinc computation for resize.c
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
I fixed the version posted at
http://web.cs.laurentian.ca/nrobidoux/misc/resize.c
Re: Cheaper sinc computation for resize.c
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.
On behalf of the ImageMagick user community, thank you for your contribution.
Re: Cheaper sinc computation for resize.c
Grab ImageMagick6.6.310 Beta which includes your patch. We still need to do more testing before we officially accept the patch. In the meantime if you have any other filter approximations, feel free to post them here.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: Cheaper sinc computation for resize.c
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.)

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.)
Re: Cheaper sinc computation for resize.c
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.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: Cheaper sinc computation for resize.c
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?
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?
Re: Cheaper sinc computation for resize.c
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
 Posts: 25764
 Joined: 20070702T17:14:5107:00
 Authentication code: 1152
 Location: Sunnyvale, California, USA
Re: Cheaper sinc computation for resize.c
Magick,magick wrote: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.
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
 Posts: 25764
 Joined: 20070702T17:14:5107:00
 Authentication code: 1152
 Location: Sunnyvale, California, USA
Re: Cheaper sinc computation for resize.c
fmw42 wrote:Magick,magick wrote: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.
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?
Re: Cheaper sinc computation for resize.c
MagickRealType is not affected by HDRI.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: Cheaper sinc computation for resize.c
ETA: Thursday at the latest.
Nicolas Robidoux
Nicolas Robidoux
 fmw42
 Posts: 25764
 Joined: 20070702T17:14:5107:00
 Authentication code: 1152
 Location: Sunnyvale, California, USA
Re: Cheaper sinc computation for resize.c
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 sinclike 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.
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 sinclike 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.