BC-splines with 2C+B=1 are optimal for EWA resampling

[NicolasRobidoux]

BC-splines with 2C+B=1 are sometime called "Keys filters" (http://www.imagemagick.org/Usage/resize/#filter_cubics). We'll use this terminology in this thread.

Suppose that one decides to do EWA (Elliptical Weighted Averaging) resampling using a cubic Hermite kernel with breaks at the integers.

Suppose, in addition, that one requires the kernel to have the following properties:

(1) It is continuous.

(2) Its derivative is continuous. ((1) and (2) together mean that the kernel is C^1.)

(3) It is even.

An even piecewise cubic with support [-2,2] is defined by 8 parameters. Condition (1) gives two conditions (left = right at x=1 and at x+2). Condition (2) adds three conditions (left derivative = right derivative at x=0, x=1 and x=2). So, the family is defined by 8-5 = 3 parameters. If, in addition, one does not care what value is taken at x=0 (provided it's not equal to zero), this leaves 2 parameters, often called B and C.

Suppose now that one wants to minimize the resampling error when the data is affine.

It is "common knowledge" that if one is "interpolating" in the usual, tensor, way, the resampling error is 0 if and only if the piecewise cubics are actually the BC-splines with 2C+B = 1 (http://www.cg.tuwien.ac.at/~theussl/DA/node11.html).

What is not common knowledge is the following (I am actually guessing that it is new), which came as a surprise to me:

If you are using even C^1 piecewise cubics with breaks at the integers and support [-2,2] as kernels for EWA resampling, and you want to minimize the resampling error when the data is univariate and affine (that is, the data is of the form a + b*x or a + b*y), then the same Keys piecewise cubics are, once again, the best ones. In that case, the error is not zero (unless b=0) for any choice of parameters (it can't). But the error for BC-splines with 2C+B=1 is smaller than for the nearby piecewise cubic splines satisfying (1)-(3).

Specifically, if one normalizes a piecewise cubic kernel satisfying (1)-(3) so that it's value at the origin is 1 (we can always normalize to any non-zero value when a kernel is used to filter), then the optimal kernels are defined by the relation: (value at 1) + (value of the slope at 1) = -1/2. This last relationship is basically equivalent to 2C+B=1.

(If I have time, I'll explain how I got this result: It's through careful numerical experiments based on proper simplification of the problem.)


Moral of the story:

If you are going to upsample or resample with a piecewise cubic kernel with breaks at the integers and support [-2,2], used as an EWA (Elliptical Weighted Averaging) filter weight, use a BC-spline with 2C+B=1. (If downsampling, I don't know.)

The current ImageMagick EWA default, the Robidoux filter, is one such cubic kernel. When I formulated the Robidoux filter, I made a quick and dirty guess that using a BC-spline with 2C+B=1 would be good (because "best" as an orthogonal interpolator should give something at least OK as an EWA kernel), with the intent of changing the parameters when I have time to tune it better. I had no idea that actually the Keys filters were actually truly, not only approximately, optimal in terms of exactness on linears a.k.a. second order accuracy, and consequently I optimized the right family (with respect to preservation of constant on rows/columns data under no-op).

[NicolasRobidoux]

My last post begs the question: Which, of the BC-splines with 2C+B=1, minimizes the error when resampling affine data? That is, which one is the "best of the best" w.r.t. to this criterion?

I'm in the process of answering this question, and will report later. The short answer: a quite blurry BC-spline with 2C+B=1.

[NicolasRobidoux]

(Still crunching number for the Keys cubic that minimizes the error on
univariate affine data.)

Now that I am certain that if one uses cubic splines (with breaks at
the integers and support [-2,2]) as EWA filter kernels they should be
Keys cubics, I can confidently answer the following question:

Which one gives the sharpest results when used in a No-Op?

(I believe that Anthony cares.)

Specifically, which one leads to the least maximum possible deviation
from the original data when the data is bounded (as image data is
always) under No-Op?

I actually discussed the answer somewhere which I can't track down
just now. It turns out that I applied the same criterion to the tuning
of the Jinc lanczoses but decided at the time that vertical/horizontal
preservation under No-Op was more important, decision which led to the
EWA Jinc LanczosSharp and Lanczos2Sharp blur values as well as the
Robidoux cubic filter which is currently the ImageMagick EWA default.

In any case, here is the answer:

The sharpest Keys cubic under No-Op is the one with

B = (-42*sqrt(2)+78)/71 = 0.26201451239901422

and

C = (1-B)/2 = (42*sqrt(2)-7)/142 = 0.3689927438004929

Compare to the Mitchell filter's B=C=1/3=0.33333333333333333 and the
Robidoux filter's

B = (228-108*sqrt(2))/199 = 0.37821575509399867

and

C = (1-B)/2 = (42*sqrt(2)-7)/142 = 0.31089212245300067

So, the new cubic is closer to Catmull-Rom than both Robidoux and Mitchell
(and, of course, B-Spline smoothing).

This sharpest Keys cubic leads to a maximum possible error under No-Op
of 24 in 8 bit (data in [0,255]). (This error, actually, occurs with
the checkerboard mode.) P.S. Correction: It's actually 49. I was halving my error bounds for a while.

In addition, it does look quite good:

Code: Select all

convert rose: -define filter:c=.3689927438004929 -filter Cubic -distort Resize 3000% rose_sharpestaCubic.png

A bit on the aliased side, which is not surprising.

So, I think this new cubic deserves a name: RobidouxSharp.

-----

Now, you may think that Catmull-Rom will lead to less maximum error
when used with EWA No-Op, but it's not the case. I can do the
computation if someone asks, but the large negative weights at sqrt(2)
will lead to large deviations from the original data.

The new cubic filter is defined by having values at 1 and sqrt(2)
exactly equal in size. At 1, "RobidouxSharp" is equal to
4.7848046234275417E-2, and at sqrt(2), it is equal to
-4.7848046234275417E-2. That is, RobidouxSharp is as close to being a
cardinal piecewise cubic (the way one must understand being cardinal
when working radially) as can be without losing "optimal"
reconstruction of linear gradients (roughly speaking).

[NicolasRobidoux]

Not sure "RobidouxSharp" is worth adding as a default: Yes it's sharp. But it's a little too aliased for my taste.

[anthony]

Here is the 'Mitchell-Netravali' diagram for comparison. The dashed line are keys filters.


The proposed RobidouxSharp would be roughly on the other side of the Mitchell-Netravali filter to the current Robidoux filter.

As the Mitchell-Netravali filter was a 'good guess' based on a 'survey of image professionals', I am amazed just how close these mathematically determined results are to that filter!

[NicolasRobidoux]

As the Mitchell-Netravali filter was a 'good guess' based on a 'survey of image professionals', I am amazed just how close these mathematically determined results are to that filter!

... except that Mitchell-Netravali was determined to be a good compromise for "direct interpolatory" use, not EWA. The amazing thing to me is that it could clearly be argued to be a good (better than both Robidoux and RobidouxSharp?) compromise for EWA too.

It's the carry over from orthogonal to polar that blows my mind. (Not that the carry over is complete: Catmull-Rom, which is excellent in orthogonal use, is not very good for EWA.)

[NicolasRobidoux]

The above RobidouxSharp comes from minimizing the worst case deviation from the original value when the original value is one of the two extremes (0 and 255 in typical 8 bit).

An alternate way is to minimize the worst case deviation when the original value is the average of the largest and smallest possible values (for example, 127.5 in typical 8 bit). I will give this a try later.

It would not surprise me if this one lands almost on top of Mitchell-Netravali.

P.S. It's not. I should have remembered that minimizing the L^1 norm generally gives sparsity, and indeed the minimizer is the Keys cubic which is equal to 0 at sqrt(2):

B = (9-3*sqrt(2))/7 = 0.67962275898295921
and
C = 0.1601886205085204

which, in principle, should be quite blurry.

P.S. 2... and indeed it is pleasant, and very blurry:

Code: Select all

convert rose: -define filter:c=.1601886205085204 -filter Cubic -distort Resize 3000% rose_sharpestaCubic.png

[NicolasRobidoux]

Now that I've taken the time to compare with a number of other schemes, classical (orthogonal) Lanczos in particular, I have to say that IMHO RobidouxSharp is probably worth adding as a "named" scheme. It is a bit more aliased than Sinc-windowed Sinc 3, but there is no "bounce back" halo, and the first halo is both narrower and slightly less deep.

[NicolasRobidoux]

Anthony:

IMHO, RobidouxSharp (discussed above) is worth adding as a named built-in. It's a good sharp EWA scheme, with haloing and aliasing (jaggies) that are not overly offensive given the sharpness.

Code: Select all

convert rose: -define filter:c=.3689927438004929 -filter Cubic -distort Resize 3000% rose_sharpestaCubic.png

[anthony]

I need the values for B, C and the first zero crossing (for using it as a window!)

[NicolasRobidoux]

RobidouxSharp Keys cubic spline:

B = 1-2*C = (78-42*sqrt(2))/71 = 0.2620145123990142

C = (42*sqrt(2)-7)/142 = 0.3689927438004929

first crossing: x = (14*sqrt(2)-1)/17 = 1.105822933719019

Axiom code to get these values:

Code: Select all

)cl a

inner x == _
eval( ( (12 - 9 * B - 6 * C) * x * x * x _
        + (-18 + 12 * B + 6 * C) * x * x _
        + (6 - 2 * B) ) / 6, B = 1-2*C )

outer x == _
eval( ( (-B - 6 * C) * x * x * x _
      + (6 * B + 30 * C) * x * x _
      + (-12 * B - 48 * C) * x _
      + (8 * B + 24 * C) ) / 6, B = 1-2*C )

outer 1 = - outer sqrt 2

solve(%,C)

eval(1-2*C,%)

eval(outer x,%%(4)) :: POLY AlgebraicNumber

rhs solve(%,x).1

[NicolasRobidoux]

Anthony:

Let me quadruple check something before moving on this.

P.S. All good. Go!

[NicolasRobidoux]

Anthony:

Even though I tuned RobidouxSharp using a minimax principle, there is an interpretation of its defining property which I think you would appreciate:

EWA with RobidouxSharp is the only Keys cubic based EWA such that under no-op, the coefficients contributed by the black pixels of a hash pattern are exactly as large as those contributed by the white pixels, excluding, as should be, the central coefficient (which we often normalize to 1, although this makes no difference).

If I read between the lines, you (a long time ago?) figured that "something like that" would be good.

Well, here it is.

If you think about it for a minute, this is the same as saying that under no-op the normalizing EWA denominator is exactly equal to the central coefficient (equal to 1, if normalized), which happens to be equivalent to the minimax criterion I also use to derive the new, as yet unimplemented, Jinc-windowed Jinc EWA LanczosSharps.

[NicolasRobidoux]

Sorry Anthony:

Let me try to more things and get back to you. I'm pretty sure RobidouxSharp is a good as I can make it, but I want to be really sure before you take time to change the source/API.

[NicolasRobidoux]

Anthony:

Now, it's a go.

I double checked that RobidouxSharp is also the l1 minimizer. And it is. So, it does not only minimize the maximum possible difference from extreme data, it actually minimizes the maximum possible difference for any data.

The Axiom code that verifies this (once other things have been set up):

Code: Select all

[ eval( ( ( abs( inner 0 - ( inner 0 + 4 * ( outer 1 + outer sqrt 2 ) ) ) + 4* ( abs outer 1 + abs outer sqrt 2) ) / ( inner 0 + 4 * ( outer 1 + outer sqrt 2 ) ) ) :: EXPR SF, C=(42*sqrt 2 - 7)/142 + i*.001 ) for i in -25..25 ]

If you know how to read the formula, it's actually obvious that the "extreme" minimizer is also the "any" minimizer.