proper scaling of the Jinc filter for EWA use

[NicolasRobidoux]

Thank you Anthony!

-----

Argh! I_0 is the 0th modified Bessel function of the first kind. (Spent the night going how in the world does one get a decent window how of Kaiser. Mistery solved.)

My choice of alpha is then fairly arbitrary. (Still works well, though.)

Back to the saddle.

[NicolasRobidoux]

If one modifies the formula for the Kaiser window so that it is continuous, that is, if one uses

Code: Select all

  return(I0(Alpha*sqrt((double) (1.0-x*x)))-1.);

instead of

Code: Select all

  return(I0A*I0(Alpha*sqrt((double) (1.0-x*x))));

one gets a pleasant family of "modified Kaiser"-windowed Jinc 3 EWA filters for which the Alpha value (in (0,infinity)) gives low aliasing for small values (but larger haloing) and less haloing (but more, but still mild, jaggies) for large values.

Note: I have not bothered normalizing the window formula so as to get 1 at zero, but this is, of course, easy.

P.S. These windowed Jincs (with the modified Kaiser window fixed so that it is continuous) work really well, and they allow one to minimize the second halo without too much increase in aliasing/blur---Alpha values between about 1.5 and 8 work well---but I still like the current distort Lanczos and LanczosSharp EWA, Jinc-windowed Jinc 3 (the second one unblurred a touch). They just antialias beautifully.

[NicolasRobidoux]

I'm getting closer to some conclusions for this round of attempts at improving things.

[NicolasRobidoux]

Before moving on, I computed the blur for Jinc-windowed Jinc 3 which would make the hash pattern stay the same under no-op (the "Anthony criterion", more or less).

It is 0.79841305642749343307. It gives horridly aliased results:

Code: Select all

~/src/ImageMagick-6.7.0/utilities/convert rose: -define filter:filter=Jinc -define filter:window=Jinc -define filter:blur=0.79841305642749343307 -define filter:lobes=3 -distort Resize 3000% rose_JincJinc3blur0p79841305642749343307.png

The axiom code to find it (its reciprocal, actually) is here:

Code: Select all

)cl a

digits 100

R1 := 3.831705970207512315614435886308160766564545274287801928762298989918839309519011470214112874757423127

r1 := R1 / %pi

R3 := 10.17346813506272207718571177677584406981951250019168555611465006811578704378288387382891893264510929

r3 := R3 / %pi

jinc x == besselJ(1,%pi*x)/x

sinc x == besselJ(1,x*(R1/r3))/x

l x == if (x<r3) then ( sinc(x) * jinc(x) ) else 0.

f r == [ r, sqrt(1.+1)*r, 2*r, sqrt(1.+4)*r, sqrt(4.+4)*r, 3*r, sqrt(9.+1.)*r, sqrt(9.+4)*r, 4*r, sqrt(16.+1)*r, sqrt(9.+9)*r, sqrt(16.+4)*r, 5*r ]

z := l 1.e-128

gplus(r) == _
  ( _
  z + _
  4 * _
  ( _
  - l((f(r)).1) _
  + l((f(r)).2) _
  + l((f(r)).3) _
  + l((f(r)).5) _
  - l((f(r)).6) _
  + l((f(r)).9) _
  + l((f(r)).11) _
  - l((f(r)).13) _
  + 2 * _
  ( _
  - l((f(r)).4) _
  + l((f(r)).7) _
  - l((f(r)).8) _
  - l((f(r)).10) _
  + l((f(r)).12) ) ) _
  ) / _
  ( _
  z + _
  4 * _
  ( l((f(r)).1) _
  + l((f(r)).2) _
  + l((f(r)).3) _
  + l((f(r)).5) _
  + l((f(r)).6) _
  + l((f(r)).9) _
  + l((f(r)).11) _
  + l((f(r)).13) _
  + 2 * _
  ( l((f(r)).4) _
  + l((f(r)).7) _
  + l((f(r)).8) _
  + l((f(r)).10) _
  + l((f(r)).12) ) ) )

[ gplus(1+.000000000000000001*i) - 1 for i in 252484527838897316..252484527838897335 ]

[NicolasRobidoux]

My only solid recommendation at this point is this:

In resize.c, replace the blur for LanczosSharp from the current value (which is barely different from 1) to the following:

0.88549061701764

(I'll triple check the exactness of this value when I have more time.)

In 8 bit (data in [0,255]) the maximum change of a pixel value under no-op is 82. (Some of my ealier reported values may have been off by half. I'll have to go back and check all of these bounds.)

I have mentioned this value before. At the time, I was more concerned about preserving vertical/horizontal data under no-op than the worst case deviation under no-op.

This blur value minimizes the worst case deviation under no-op. It (or, rather, its reciprocal) was found using this axiom code:

Code: Select all

)cl a

digits 100

R1 := 3.831705970207512315614435886308160766564545274287801928762298989918839309519011470214112874757423127

r1 := R1 / %pi

R3 := 10.17346813506272207718571177677584406981951250019168555611465006811578704378288387382891893264510929

r3 := R3 / %pi

jinc x == besselJ(1,%pi*x)/x

sinc x == besselJ(1,x*(R1/r3))/x

l x == if (x<r3) then ( sinc(x) * jinc(x) ) else 0.

plus x == max(l x,0)

minus x == -min(l x, 0)

f r == [ r, sqrt(1.+1)*r, 2*r, sqrt(1.+4)*r, sqrt(4.+4)*r, 3*r, sqrt(9.+1.)*r, sqrt(9.+4)*r, 4*r, sqrt(16.+1)*r, sqrt(9.+9)*r, sqrt(16.+4)*r, 5*r ]

z := l 1.e-128

gplus(r) == _
  ( _
  4 * _
  ( plus((f(r)).1) _
  + plus((f(r)).2) _
  + plus((f(r)).3) _
  + plus((f(r)).5) _
  + plus((f(r)).6) _
  + plus((f(r)).9) _
  + plus((f(r)).11) _
  + plus((f(r)).13) _
  + 2 * _
  ( plus((f(r)).4) _
  + plus((f(r)).7) _
  + plus((f(r)).8) _
  + plus((f(r)).10) _
  + plus((f(r)).12) ) ) _
  ) / _
  ( _
  z + _
  4 * _
  ( l((f(r)).1) _
  + l((f(r)).2) _
  + l((f(r)).3) _
  + l((f(r)).5) _
  + l((f(r)).6) _
  + l((f(r)).9) _
  + l((f(r)).11) _
  + l((f(r)).13) _
  + 2 * _
  ( l((f(r)).4) _
  + l((f(r)).7) _
  + l((f(r)).8) _
  + l((f(r)).10) _
  + l((f(r)).12) ) ) )

gminus(r) == _
  ( _
  4 * _
  ( minus((f(r)).1) _
  + minus((f(r)).2) _
  + minus((f(r)).3) _
  + minus((f(r)).5) _
  + minus((f(r)).6) _
  + minus((f(r)).9) _
  + minus((f(r)).11) _
  + minus((f(r)).13) _
  + 2 * _
  ( minus((f(r)).4) _
  + minus((f(r)).7) _
  + minus((f(r)).8) _
  + minus((f(r)).10) _
  + minus((f(r)).12) ) ) _
  ) / _
  ( _
  z + _
  4 * _
  ( l((f(r)).1) _
  + l((f(r)).2) _
  + l((f(r)).3) _
  + l((f(r)).5) _
  + l((f(r)).6) _
  + l((f(r)).9) _
  + l((f(r)).11) _
  + l((f(r)).13) _
  + 2 * _
  ( l((f(r)).4) _
  + l((f(r)).7) _
  + l((f(r)).8) _
  + l((f(r)).10) _
  + l((f(r)).12) ) ) )

-- To find the one with least possible deviation (this is a converged result):

[ gplus(1+.00000000000000001*i) - gminus(1+.00000000000000001*i) for i in 12931744366533334..12931744366533349 ]

Now, with this blur, EWA LanczosSharp is a bit too aliased to my taste. It is also just a touch sharper than orthogonal Lanczos 3. On the other hand, "plain" EWA Lanczos is a bit too soft for my (and most people's) taste. I just don't know how to attach a "less than 1 but more than .88" blur value to EWA Lanczos without adding yet one more exception. So, I'll leave aside the issue of how to "officially" put a balanced value in IM for now.

With the usual test pictures, here is the code that generates the results:

Code: Select all

convert sl.png -define filter:filter=Jinc -define filter:window=Jinc -define filter:blur=0.88549061701764 -define filter:lobes=3 -distort Resize 3000% sl_JincJinc3blur0p885469.png

convert rose: -define filter:filter=Jinc -define filter:window=Jinc -define filter:blur=0.88549061701764 -define filter:lobes=3 -distort Resize 3000% rose_JincJinc3blur0p88549.png

Although I am not 100% satisfied with my attempt at thinking outside of the zero sum box, one thing is clear (at least to me): The above LanczosSharp (with blur=.88...) is unquestionably better than classical orthogonal Lanczos with pixel art (which is not saying much, I'll admit). In natural images, it appears to do better than classical Lanczos 3 with diagonal interfaces and lines (which is not much of an achievement) at the expense of straightening nearly vertical and horizontal interfaces and lines a little too much to my taste. This gives results a slightly "boxy" look compared to classical Lanczos (which, once again, "breaks" near diagonal lines and interfaces; the two methods are kind of complementary).

In any case, this whole thread is about exploratory work, so the door is wide open for criticism of the results and the methods.

[NicolasRobidoux]

For example, Jinc-windowed Jinc 3 with blur = 3/r3 = 0.9264075766146068 is totally deserving to be "special." This is the one blurred just enough so that the EWA disk has radius 3, no more and no less.

[anthony]

If one modifies the formula for the Kaiser window so that it is continuous, that is, if one uses

Code: Select all

  return(I0(Alpha*sqrt((double) (1.0-x*x)))-1.);

instead of

Code: Select all

  return(I0A*I0(Alpha*sqrt((double) (1.0-x*x))));

one gets a pleasant family of "modified Kaiser"-windowed Jinc 3 EWA filters for which the Alpha value (in (0,infinity)) gives low aliasing for small values (but larger haloing) and less haloing (but more, but still mild, jaggies) for large values.

Note: I have not bothered normalizing the window formula so as to get 1 at zero, but this is, of course, easy.

Actually the IOA is not even needed, as all filter function results are normalized before being used. Part of their defined use. In fact the normalization factor for Gaussian is not even applied so that you get a 1.0 at zero offset for comparision purposes.

The IOA is simply calculated to be 1/Kaiser() at x=0. For a short time during the above I even removed it, resulting in a graph that has a x peak of just over 200. It is easily calculated, but actually not needed.

Hmmm Check Wikipedia http://en.wikipedia.org/wiki/Kaiser_window
I believe it is defined right... It matches the original definition does match that of the original 'zoom' program from Paul Heckbert.

Code: Select all

typedef struct {	/* data for parameterized Kaiser window */
    double a;		/* = w*(N-1)/2 in Oppenheim&Schafer notation */
    double i0a;
    /*
     * typically 4<a<9
     * param a trades off main lobe width (sharpness)
     * for side lobe amplitude (ringing)
     */
} kaiser_data;

double filt_kaiser(double x, void *d)	/* parameterized Kaiser window */
{
    /* from Oppenheim & Schafer, Hamming */
    kaiser_data *k;

    k = (kaiser_data *)d;
    return bessel_i0(k->a*sqrt(1.-x*x))*k->i0a;
}

static void kaiser_init(double a, double b, void *d)
{
    kaiser_data *k;

    k = (kaiser_data *)d;
    k->a = a;
    k->i0a = 1./bessel_i0(a);
}

double bessel_i0(double x)
{
    /*
     * modified zeroth order Bessel function of the first kind.
     * Are there better ways to compute this than the power series?
     */
    int i;
    double sum, y, t;

    sum = 1.;
    y = x*x/4.;
    t = y;
    for (i=2; t>EPSILON; i++) {
	sum += t;
	t *= (double)y/(i*i);
    }
    return sum;
}

P.S. These windowed Jincs (with the modified Kaiser window fixed so that it is continuous) work really well, and they allow one to minimize the second halo without too much increase in aliasing/blur---Alpha values between about 1.5 and 8 work well---but I still like the current distort Lanczos and LanczosSharp EWA, Jinc-windowed Jinc 3 (the second one unblurred a touch). They just antialias beautifully.

I thought the Kaiser function was continuous, at least in the 0.0 to 1.0 range which is all that is required of a windowing function. The window function is scaled automatically to match the support (or windowed_support expert option) of the main filter function.

I am quite willing to 'unblur' the LanczosSharp 'a touch' if you want me to.

[NicolasRobidoux]

Anthony:

You programmed Kaiser just fine. (Although indeed I0A is superfluous for the reason you mentioned.)

It is just that Kaiser is nonzero at the endpoints of the window (at x=1) which makes it discontinuous when you look at the window in the context of the entire real line.

This nonzero value at the endpoints, when used with the Jinc filter (and, I imagine, Sinc as well) leads to "bounceback" haloing which has a "flat ring" feel to it (even though, generally, it is small). Which makes sense, because it's not tapered out as strongly near the edge of the support as if there is a double root (one that comes from the Jinc or Sinc filter itself, and one from most window functions).

The modified Kaiser is just that: a modified version to get rid of a specific artifact.

[NicolasRobidoux]

I am quite willing to 'unblur' the LanczosSharp 'a touch' if you want me to.

Please do: LanczosSharp should be defined by blur=0.88549061701764

(I currently don't have a properly set up svn to commit this change myself.)

I also suggest that in the IM documentation (... Examples) mention that blur values between the LanczosSharp value and 1 work well with Jinc-windowed Jinc EWA, and that, for example, blur = 3/r3 = 0.9264075766146068, a value set so that the no-op EWA disc has radius 3, strikes a good balance between blur and jaggies (r3 is the location of the third crossing of the Jinc function as documented in resize.c; I use a higher precision value in the axiom code above). Finally, mention that values closer to 1 (but slightly less: I'll post my favorite, which is a nudge bit larger than the current LanczosSharp value, when I will have recomputed it carefully) do an OK job at smoothly enlarging pixel art and text (if halo is not verboten, of course).

At this point, I want to express my thanks to Anthony, Magick and Fred (and Heckbert too) for having set up ImageMagick and this forum so that it is such an easy tool to use for investigation, and also direct special thanks to Anthony for the insights I gained from his feedback when trying to improve EWA resampling.

[NicolasRobidoux]

Anthony:

If you want LanczosSharp to be "good" as opposed to "extreme," just use the blur 0.9264075766146068 = 3/r3.

I don't seem to make up my mind on this, so probably you should test both this value and the "extreme" one and decide for yourself. I've not tested on a wide enough variety of images to know for sure. But on the rose, 3/r3 is definitely more "balanced." I do not like that the "extreme" .88... value emphasizes nearly vertical or horizontal features. .92... also does that, but less noticeably.

So, maybe, .92640... should be the LanczosSharp blur, with notes in the ImageMagick examples that if you want to minize worst case no-op variation .88... is the one, and that is is very sharp.

[NicolasRobidoux]

What this "emphasize vertical/horizontal" business suggests is that maybe I should try minimizing the error when the data comes from a linear gradient. However, doing this will likely lead to quite a blurry scheme.

The alternative is to play with windowing functions some more, or go to Jinc (or other)-windowed Jinc 4, which I don't really want to because it will widen haloing.

[NicolasRobidoux]

The blur to make Jinc-windowed Jinc 4 have least worst case deviation under no-op is 0.88451002338585141. It gives excellent results:

Code: Select all

convert rose:  -define filter:filter=Jinc -define filter:window=Jinc -define filter:blur=0.88451002338585141 -define filter:lobes=4 -distort Resize 3000% rose_JincJinc4blur0p88451002338585141.png

convert sl.png  -define filter:filter=Jinc -define filter:window=Jinc -define filter:blur=0.88451002338585141 -define filter:lobes=4 -distort Resize 3000% sl_JincJinc4blur0p88451002338585141.png

Like for Jinc-windowed Jinc 3, it looks just a bit too sharp to me.

The axiom code to find this value is at the end of this post.

-----

The blur to make Jinc-windowed Jinc 4-lobe have the same extent as Sinc-windowed Sinc 4-lobe, namely four, is 4/r4 = 0.943159799432847711 (r4 is the normalized Jinc root in resize.c). Just like for Jinc-windowed Jinc 3, this blur gives excellent results:

Code: Select all

convert rose:  -define filter:filter=Jinc -define filter:window=Jinc -define filter:blur=0.943159799432847711 -define filter:lobes=4 -distort Resize 3000% rose_JincJinc4blur0p943159799432847711.png

convert sl.png  -define filter:filter=Jinc -define filter:window=Jinc -define filter:blur=0.943159799432847711 -define filter:lobes=4 -distort Resize 3000% sl_JincJinc4blur0p943159799432847711.png

It appears that having the Jinc-windowed Jinc be blurred so that its extent matches the extent of the corresponding Sinc-windowed Sinc (that is: extent = number of lobes) is a good recipe to produce well rounded EWA filters. I am starting to think that these should be the ones that are highlighted by getting the LanczosSharp name, instead of the ones that are "maximally reasonably sharp."

This four lobe version, in particular, is just gorgeous. And the 3-lobe version (with extent exactly 3) is very nice if you want less haloing.

-----

Here is the axiom code that performs the "minimization of the worst possible deviation under no-op":

Code: Select all

)cl a

digits 100

R1 := 3.831705970207512315614435886308160766564545274287801928762298989918839309519011470214112874757423127

r1 := R1 / %pi

R4 := 13.32369193631422303239368412694787675121664473135786578547757152649656706334730478254711901794883995

r4 := R4 / %pi

jinc x == besselJ(1,%pi*x)/x

sinc x == besselJ(1,x*(R1/r4))/x

l x == if (x<r4) then ( sinc(x) * jinc(x) ) else 0.

plus x == max(l x,0)

minus x == -min(l x, 0)

f r == [ r, sqrt(1.+1)*r, 2*r, sqrt(1.+4)*r, sqrt(4.+4)*r, 3*r, sqrt(9.+1.)*r, sqrt(9.+4)*r, 4*r, sqrt(16.+1)*r, sqrt(9.+9)*r, sqrt(16.+4)*r, 5*r ]

z := l 1.e-128

gplus(r) == _
  ( _
  4 * _
  ( plus((f(r)).1) _
  + plus((f(r)).2) _
  + plus((f(r)).3) _
  + plus((f(r)).5) _
  + plus((f(r)).6) _
  + plus((f(r)).9) _
  + plus((f(r)).11) _
  + plus((f(r)).13) _
  + 2 * _
  ( plus((f(r)).4) _
  + plus((f(r)).7) _
  + plus((f(r)).8) _
  + plus((f(r)).10) _
  + plus((f(r)).12) ) ) _
  ) / _
  ( _
  z + _
  4 * _
  ( l((f(r)).1) _
  + l((f(r)).2) _
  + l((f(r)).3) _
  + l((f(r)).5) _
  + l((f(r)).6) _
  + l((f(r)).9) _
  + l((f(r)).11) _
  + l((f(r)).13) _
  + 2 * _
  ( l((f(r)).4) _
  + l((f(r)).7) _
  + l((f(r)).8) _
  + l((f(r)).10) _
  + l((f(r)).12) ) ) )

gminus(r) == _
  ( _
  4 * _
  ( minus((f(r)).1) _
  + minus((f(r)).2) _
  + minus((f(r)).3) _
  + minus((f(r)).5) _
  + minus((f(r)).6) _
  + minus((f(r)).9) _
  + minus((f(r)).11) _
  + minus((f(r)).13) _
  + 2 * _
  ( minus((f(r)).4) _
  + minus((f(r)).7) _
  + minus((f(r)).8) _
  + minus((f(r)).10) _
  + minus((f(r)).12) ) ) _
  ) / _
  ( _
  z + _
  4 * _
  ( l((f(r)).1) _
  + l((f(r)).2) _
  + l((f(r)).3) _
  + l((f(r)).5) _
  + l((f(r)).6) _
  + l((f(r)).9) _
  + l((f(r)).11) _
  + l((f(r)).13) _
  + 2 * _
  ( l((f(r)).4) _
  + l((f(r)).7) _
  + l((f(r)).8) _
  + l((f(r)).10) _
  + l((f(r)).12) ) ) )

g x == gplus x - gminus x

-- The one with least possible deviation:

g(1.13056943794945350)

gplus(1.13056943794945350)

gminus(1.13056943794945350)

-- Gives the reciprocal of .88451002338585141

1/1.13056943794945350

[NicolasRobidoux]

Anthony:

Here is a suggestion (warning: API change):

Whenever someone asks for Lanczos and it's an (orthogonal) -resize call, give them Sinc-windowed Sinc. If -define filter:lobes=n, given them Sinc-windowed Sinc n lobes. The default n is 3, so that without an explicit lobe specification, you get the current -filter Lanczos.

When it is an EWA (-distort) call, give them Jinc-windowed Jinc with blur = n/r_n, with r_n the nth root of Jinc. This gives a filter which is closer to what people expect (hope for, should I say) out of Lanczos, certainly less blurry than the current EWA Lanczos and Lanczos2. For example, the support of the EWA and orthogonal Lanczoses then have the same diameter. Again, the default number of lobes shoudl be 3.

If people want the unblurred Jinc-windowed Jinc, they can get it through Sinc/Jinc.

-----

I'll deal with the issue of what to do with LanczosSharp separately. What would be nice is to set the blurs for all the reasonable numbers of lobes (1 to what I believe is the current maximum eight), and then essentially the behavior would be defined very similarly to the above API change for Lanczos, just sharper. I can certainly do that, although of course it requires some programming and running (just bought two computers but had no time to install...).

-----

What do you think?

[NicolasRobidoux]

I am still stunned at the combination of naturality and sharpness of the result of enlarging the rose picture with Jinc-windowed Jinc 4-lobe EWA sharpened so as to minimize worst case deviation under no-op (blur=.88549...). It still looks like it may be oversharpening a touch. But the result fills me with amazement. And it is unquestionably much better on the sl.png test image than standard (orthogonal) Lanczos 3.

P.S. The "blurring" weakness of the above EWA is actually not surprising: Perfectly diagonal lines are blurred along their peak. I am starting to think that my method of determining "optimal" blur is actually the right one for "maximum bearable sharpness" for Jinc-windowed Jinc EWA.

P.S. 2 This "optimization" approach worked pretty well for the proposed "RobidouxSharp" cubic sampler. So, it may be reasonably general. Or else it could just be dumb luck.

[NicolasRobidoux]

I am happy with the results I got so far. When I have a minute, I may want to investigate this:

J_0 is approximately proportional to Jinc for small input values. This is a consequence of the continued fraction equation (77) of http://mathworld.wolfram.com/BesselFunc ... tKind.html.

What this means is that J_0-windowed Jinc may work well. Better than Jinc-windowed Jinc? Don't know.

P.S. I've not programmed J_0, but computations of the derivative at the windowing function's limit shows that haloing will be worse with J_0 as windowing function than with Sinc, which itself will be worse than with Jinc. At this point, I really doubt I can do much better than JInc-windowed EWA with a lot of work. So, I'm done. The only thing left to do in this round (I think) is to compute the "optimal" sharpening blurs for all numbers of lobes from 1 to 8.