Cheaper sinc computation for resize.c

[anthony]

In cylindrical filters the Bessel (Jinc()) function is generally used instead of Sinc.

Basically it is a single pass filter using a circular (or elliptical) weighted sampling of the pixel neighbourhood, instead of a 2 pass linear sampling. It is used for image distortions where rotations can be involved (as well as scale).

From IM Examples...

The Jinc() function (more commonly known as a 'Bessel' filter) is the equivalent filter for use in a two dimensional cylindrical, or radial filtering operation, though it is not a perfect fit in this regard. Though very similar and closely related to Sinc() (see the graph provided) it is designed to filter a rectangular array of values using a radial or cylindrical distance, rather than orthogonal (axis aligned) distance. As such the first 'zero-crossing' falls between the values of 1.0 (for orthogonal neighbours) and the square root of 2. That is it crosses zero at approximatally the value of 1.2196, so as to match as close as posible all its 'close' neighbours.

The Lanczos equivalent would thus be a Bessel windowed Bessel filter.

I have seen on paper that did use this specific variant, and the improvement in a rotation distorted image, exampled in the paper, was quite simply amazing. However I have not reached the point of actually duplicating those results.
Hmmm... Can't seem to find the reference to the paper in my bookmarks.
Perhaps I have a physical copy. I'll have a look.

Actually I am more interested in a Mitchell like cylindrical varient.

[anthony]

As such for 'typical' resizing and distorts, the caching vastly reduces the cost of trig functions.
Only in the extreme case where 'caching' provides little benefit, does the polynomial form show a major improvement. However that does not mean we would not want to make use of the improvements. Especially when the error difference is so low. We may even make the polynomial function the default!

May I suggest calling the polynomial approximations "Sinc" and their more accurate versions "SincTrigonometric"? (!)

My thoughts exactly. Though explaining this complication in IM examples will be tricky, especially now you can lookup exactly what functions were selected using the filter:verbose setting.

[NicolasRobidoux]

May I suggest calling the polynomial approximations "Sinc" and their more accurate versions "SincTrigonometric"? (!)

My thoughts exactly. Though explaining this complication in IM examples will be tricky, especially now you can lookup exactly what functions were selected using the filter:verbose setting.

I was actually joking!

If you're going to do something like that, I may boost the orders of the approximations just a little so that people can be assured that there is ABSOLUTELY NO SIGNIFICANT DIFFERENCE (except that the Q64 version is actually more precise).

A great test image is the large "rings" picture http://www.imagemagick.org/Usage/resize ... g_orig.png. It really reveals how close things are (the only thing it's lacking is "text like" features and straight monochrome lines, although for such features Sinc sucks so reproducing its behavior exactly is nonsense). I'll do some additional testing to figure out what "you can't tell the difference between genuine and the imitation" requires.

[NicolasRobidoux]

Of course I'd be delighted if the Poor Man's Sincs became the Sinc defaults, and "standard trig sincs" became SincTrigonometric.

[NicolasRobidoux]

In cylindrical filters the Bessel (Jinc()) function is generally used instead of Sinc.

Basically it is a single pass filter using a circular (or elliptical) weighted sampling of the pixel neighbourhood, instead of a 2 pass linear sampling. It is used for image distortions where rotations can be involved (as well as scale).

From IM Examples...

The Jinc() function (more commonly known as a 'Bessel' filter) is the equivalent filter for use in a two dimensional cylindrical, or radial filtering operation, though it is not a perfect fit in this regard. Though very similar and closely related to Sinc() (see the graph provided) it is designed to filter a rectangular array of values using a radial or cylindrical distance, rather than orthogonal (axis aligned) distance. As such the first 'zero-crossing' falls between the values of 1.0 (for orthogonal neighbours) and the square root of 2. That is it crosses zero at approximatally the value of 1.2196, so as to match as close as posible all its 'close' neighbours.

The Lanczos equivalent would thus be a Bessel windowed Bessel filter.

I have seen on paper that did use this specific variant, and the improvement in a rotation distorted image, exampled in the paper, was quite simply amazing. However I have not reached the point of actually duplicating those results.
Hmmm... Can't seem to find the reference to the paper in my bookmarks.
Perhaps I have a physical copy. I'll have a look.

Actually I am more interested in a Mitchell like cylindrical varient.

This is just an intuition, but I do believe that lanczos2(r^2) (where r is the distance to the point about which the filter is cylindrically symmetric, that is, r^2 = x^2+y^2) would give a high quality filter. Its support would be the disk of radius sqrt(2) (0 outside of that). Kind of a cylindrical Catmull-Rom (a bit less sharp), and consquently sharper than small radius Bessel (I think). This "cylindrical lanczos2" should be quite fast because of the tight support, and have a reasonable balance between blur and halo.

If you want something which is more Mitchell like, the first thing I'd try is actually mitchell(r) with support the disk of radius 2 (r = sqrt(x^2+y^2)).

If I figure out something better, I'll let you know.

You seem to have read more about this than I, so maybe these suggestions are "inferior." Take these as ignorant hunches, no more and no less.

[NicolasRobidoux]

Just to make sure there is no confusion:

All my suggestions so far where "using a kernel function as a weighing function" (Michell-Netravali is normally used "the other way"): Over the support, give weights to input pixels p(x_i,y_j) based on their relative position (x_i,y_j) w.r.t. the sampling location, normalize the weights, and convolve: sum weight(x_i,y_j)*p(x_i,y_j) / sum weight(x_i,y_j).

[anthony]

You seem to have read more about this than I, so maybe these suggestions are "inferior." Take these as ignorant hunches, no more and no less.

More that I have been plugging away at cylindrical filters for about 3 years, looking for information on something other than a blurry Gaussian. Wish I can find that article! It was the only one I ever found that dealt with a cylindrical filter (as they were looking at rotational distortion) and did not just simply follow Paul Heckbert's use of Gaussian, and they used a 8 lobe Lanczos basied on the Bessel function (which It thought was overkill).

I will need to get some time to go in an do some work on the EWA distortion filter setup, so I can properly do some experiments. Also to provide some means to turn off the use of an interpolation filter in areas the image is enlarged so it only relies on the cylindrical filter.

ASIDE: currently it switches automatically between EWA and Interpolation depending on if the pixel area is becoming enlarged or not. Paul Heckbert himself had a solution to that in a later paper. The only reason I did not use it was because Gaussian become very noticable when scaling is very low. That is images that were not distorted much turned out very blurry!


Sorry I am ranting. But then this is a problem I have been having for a long time.


Lanczos2(r^2) Nice idea.


Hmmm strange I never noticed that
1D Gaussian == exp( -(r^2)/(2*sigma^2) ) / (2*pi*sigma^2)
2D Gaussina == exp( -(r^2)/(2*sigma^2) ) / (sqrt(2*pi)*sigma)
That is only the normalization factor changed!! The filter is the same after normalizations!
I never notice that before! Which is strange as I have been programming 1 and 2 D gaussians quite heavily!

I was always under the impression that their was more of a difference as the Gaussian used by Paul Heckbert
did not seem to match up. But then they were creating a table of 1000 values which does a lookup with a r^2 value.
over a r=0 to 2 range. The exp() had weird things in it, like a 4*ln(2) constant that I could not understand and was not explained by Paul Heckbert as being already 'understood'. Of course I don't and didn't understand so had to use the formula blind (or semi-blind) which means I may have it all wrong!

Back to the books for me... I hate reading research papers -- they make too many assumptions about the background knowledge of the reader! At least I have kept hard copies of those papers!

[NicolasRobidoux]

ASIDE: currently it switches automatically between EWA and Interpolation depending on if the pixel area is becoming enlarged or not.

I would guess you know this already (I suppose I could figure it out by looking at the code):

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.

(Forgive me if I'm preaching to the converted.)

[anthony]

Just to make sure there is no confusion:

All my suggestions so far where "using a kernel function as a weighing function" (Michell-Netravali is normally used "the other way"): Over the support, give weights to input pixels p(x_i,y_j) based on their relative position (x_i,y_j) w.r.t. the sampling location, normalize the weights, and convolve: sum weight(x_i,y_j)*p(x_i,y_j) / sum weight(x_i,y_j).

Yes I understand that, though it is more complex that than when transparency is involved!

i,j = destination x,y = source pixels being filtered
sum() = over all the x,y pixels in the scaled filter support area (be it linear or cylindrical)

p(i,j) = sum( weight(x,y)*a(x,j)*p(x,y) ) / sum( weight(x,y_j)*a(x,y) )
a(i.j) = sum( weight(x,y)*a(x,j) ) / sum( weight(x,y_j) )

where p() = color for that pixel
and a() = alpha or transparency for that pixel.

If you don't you end up with 'Transparency halo effects, such as was demonstrated in a OLD bug report for resize
http://www.imagemagick.org/Usage/bugs/resize_halo/
and again for convolution blurs
http://www.imagemagick.org/Usage/bugs/blur_trans/
Both problems were fixed in some way - though not fixed well for convolution blurs, though it is fixed properly in
Morphology Convolution (old convolution is not using the new version yet)

Sorry its late and I am ranting again.

[NicolasRobidoux]

(Forgive me if I'm preaching to the converted.)

Apologies: I now realize that I was preaching to the Pope

-----

This is interesting: What it says is that when resampling RGBA (say) must be treated kind of like CMYK.

[NicolasRobidoux]

May I suggest calling the polynomial approximations "Sinc" and their more accurate versions "SincTrigonometric"? (!)

My thoughts exactly. Though explaining this complication in IM examples will be tricky, especially now you can lookup exactly what functions were selected using the filter:verbose setting.

Actually, I don't think it's that tricky to explain:

Now that I just changed the code (it's not in svn yet: as soon as Anthony puts in the new lanczos n for integer n, I'll put it in) so that it's impossible to tell which one is computed with trig calls and which one is computed with direct polynomial approximations, maybe something like this would be enough:

A superfast method of computing the sinc function (based on polynomial approximations which are "optimal" with respect to relative error) is used by default. If you want to use the usual sin(pi*x)/(pi*x) approach (which relies on the sin function from math.c), use SincTrigonometric.

[NicolasRobidoux]

I double checked with Dr Ralf Meyer (who knows better about the standards and how modern compilers/chips behave under the hood than me) and my insistent explicit casting of MagickRealTypes to doubles prior to plugging into sin, cos etc was redundant:

Basically, math.h REQUIRES the input to be cast to double if it is not already, so the compiler MUST perform the cast even if it is not specified.

Likewise with assignments like:

const double out = in1 * in2;

If at least one of in1 and in2 is a long double, then the RHS is computed with long double arithmetic (which is very slow on current intel chips which farm out long double computations to the "x87" unit which makes long double arithmetic about as slow as if it was done in software). However, the result is automatically cast to double when it is "stored" in out.

I've started rewriting resize.c to remove the redundant casts (and relying more on automatic upcasting as required by the standard).

[NicolasRobidoux]

Suggestion:

It appears to me that there is no benefit to computing filter coefficients at long double accuracy, ever.

For one thing, most of them have been computed using double precision math.c routines, not their long (C99) versions, and nobody, apparently, has complained (unless there is some overloading I'm missing).

An obvious first suggestion is never to have the coefficients be larger than double. For one thing, this would make their cacheing more compact.

Another is to have the resize.c filters take in MagickRealTypes (which may be long doubles) and spit out MagickRealTypes, but actually use doubles to do the computation when it is nothing but trivial.

The motivation for this comment is the following:

In the not so old days, long doubles, with a good compiler, would turn now to be "extended" 80ish-bit floating point numbers, with which one could compute just as fast as with doubles.

With 64-bit chips/OS, however, extended 80-bit is not there any more (really), and computations with types larger than (standard) doubles are not offloaded to SSE2, but instead shipped off to the "x87" unit, which is, generally, brutally SLOW.

------

I have no idea whether long doubles should or should not be used at all in IM. I'm just saying that not something needs to be computed in long doubles just because it will "touch" something which is a long double.

[NicolasRobidoux]

I put a version NOT DOUBLE CHECKED of resize.c at the usual

http://web.cs.laurentian.ca/nrobidoux/misc/IM/

with the following features:

1)

I have eliminated SincPolynomial, basically by replacing Sinc with a version of it with tolerances tighter than they were before.

I have not triple checked the tolerances for Q32 and Q64, but for 8 and 16 there is no significant difference with using trigs, and the coefficient computation is much faster. (Granted, it is a small part of the whole computation...)

2)

I have rewritten LanczosChebyshev so that it always returns the right result. Unless the nested branching makes it slower than the "direct" approach on reasonably current chips (I doubt this very much, although I've not benchmarked enough), I think that it could be some sort of replacement for lanczos ("sinc-sinc").

3)

I have removed a lot of the explicit casts, given that "the compiler knows what to do."

--------

Given that the above changes could be argued for or against, I have NOT put them in svn. Once I know what is a go and what is a no, I'll triple check, and then commit.

[NicolasRobidoux]

SUMMARY: The current version of LanczosChebyshev is not efficient enough to be ready for prime time.

Details:

As you may know, my opinion is that Sinc should actually be "my" Sinc (a.k.a. SincPolynomial, but it's just Sinc in my home version of resize.c).

So, I compared LanczosChebyshev (which calls one trig then uses a short recursion to compute the needed product of sines) to "nicolas' plain vanilla lanczos" (which uses two calls to (polynomial) Sinc).

As it turns out (this is for lanczos 3), the one trig call still left in LanczosChebyshev makes it SLOWER than simply calling Sinc (the polynomial version) twice.

Conclusion:

IF, indeed, SincPolynomial becomes Sinc, there is NO REASON to use the recursion to compute the lanczoses UNLESS THE NUMBER OF LOBES IS LARGER THAN 4.

Or UNLESS I USE A POLYNOMIAL APPROXIMATION OF cos(pi t) for t in [-1,1] to start the recursion.

PS

It is indeed and issue of number of trig calls: If I compare lanczospolynomial to the lanczos computed with Sinc's which make trig calls, then lanczospolynomial comes out ahead.