on the theory of the Spline16/36/64 resizers

[Wilbert]

I contacted the author of Panorama Tools (this is where the Spline16/36/64 appear to come from) and asked him about the coefficient of the Spline polynomials. He sent me a short reply which i have worked out in detail:

Let me first say that the polynomials are not the interpolating polynomials, but they are certain blending polynomials as will be explained below.

Let's consider Spline4 for simplicity. Consider one interpolating (cubic) spline for the interval [0,1]:

S(x) = a[3]*x^3 + a[2]*x^2 + a[1]*x + a[0] for x in [0,1]

Let's define: y0 = S(0), y1 = S(1), z0 = S''(0), z1 = S''(1). The coefficients (a[j])__j depend on (y0,y1,z0,z1). In fact, the values (a[j])_j are fixed by (y0,y1,z0,z1). Once you know (y0,y1,z0,z1), you can calculate the coefficients (a[j])_j, and thus the polynomial S.

Ok, let's set z0=z1=0. In other words, the second derivatives at the end points are zero. Under this condition (a[j])_j are fixed by (y0,y1).

What happens if you replace y0 by y1 and y1 by y0? In that case the polynomial S will be mirrored at x=1/2, since (1) S is fixed by (y0,y1) and (2) z0=z1. A consequence of this is that:

S[y0;y1](x) = S[y1;y0](1-x)

since the mapping x |-> 1-x results in the same mirroring of S at x=1/2.

It turns out that S can be expressed as a linear combination of the control points (y0,y1):

S(x) = w0(x) * y0 + w1(x) * y1

Since w0,w1 blend the points y0,y1 together, they are called blending polynomials (or basis functions). As a result of the boundary conditions we have:

y0 = S(0) = w0(0) * y0 + w1(0) * y1 => w0(0)=1, w1(0)=0
y1 = S(1) = w0(1) * y0 + w1(1) * y1 => w0(1)=0, w1(1)=1

It's also true that the sum of the blending polynomials is equal to one, that is

w0(x)+w1(x)=1

but i'm unable to prove that for the general case.

Because of the mirroring property we also have that

S(x) = w0(1-x) * y1 + w1(1-x) * y0

It follows that

w1(x) = w0(1-x)
w0(x) = w1(1-x)

Subtracting the segment bias results in

w1(x) = w0(1-x) -> w0[-x]

getting a symmetrical kernel.

Solving (w0,w1) using Maple:
The boundary conditions result in:

> eqn:={a0=y0, a3+a2+a1+a0=y1, 2*a2=0,6*a3+2*a2=0};

Solving for (a[j])_j:

> solve(eqn,{a0,a1,a2,a3});
{a0 = y0, a2 = 0, a3 = 0, a1 = -y0 + y1}

Calculating the coefficients of (y0,y1), that is (w0,w1):

> a0s:=y0; a1s:=-y0+y1; a2s:=0; a3s:=0;
> S:=a3s*x^3 + a2s*x^2 + a1s*x + a0s;
S:= (-y0 + y1) x + y0

> w0:=coeff(S,y0); w1:=coeff(S,y1);
w0 := -x + 1
w1 := x

Yes, that's the bilinear resizer (because a3=a2=0). I never realized that.

Two additional remarks:

1) In general for Spline k^2, we start with k-1 interpolating polynomials. They are located at (-k/2+1,-k/2+2), ..., (k/2-1,k/2). The bending polynomials are calculated from the polynomial located at [0,1], which is expressed as a linear combination of (y0, ..., y(k-1)). The boundary conditions are:

y[1](-k/2+1) = y0
y[1](-k/2+2) = y1
y[2](-k/2+2) = y1
...
y[k-2](k/2-1) = y(k-2)
y[k-2](k/2-1) = y(k-2)
y[k-1](k/2) = y(k-1)

and also continuity of first and second derivatives, where the second derivatives of the two endpoints are zero.

I have no idea how to derive a general formula. I can recalculate the known coefficients for Spline16 and Spline36. However for Spline64 I get slightly different ones. So, i assume that the known ones for Spline64 are simply wrong.

edit Wilbert: I made a mistake with Spline64. Its coefficients are correct now.

[Egh]

Great job!

I will forward this explanation to a contact of mine. I wondered if we can derive higher order of splines, for instance, spline144 etc, and if they are on any use in practice, of course.

[MfA]

Originally Posted by Wilbert

and also continuity of first and second derivatives, where the second derivatives of the two endpoints are zero.

Derivative of the interpolation kernel of spline 16 goes from -4/5 to -7/15 at the crossover ... not very continuous.

[madshi]

@Wilbert, good job! Although I barely understand most of what you wrote...

Originally Posted by Wilbert

However for Spline64 I get slightly different ones. So, i assume that the known ones for Spline64 are simply wrong.

What Spline64 coefficients do you get? Maybe they look better than those used now by everyone? Also can you calculate coefficients for Spline256?

Thanks!

[Wilbert]

@MfA,

Derivative of the interpolation kernel of spline 16 goes from -4/5 to -7/15 at the crossover ... not very continuous.

No, you are referring to the blending polynomials here, their derivatives are indeed not continuous.

The derivatives of the interpolating functions (S(x) and the other ones) are continuous at the crossovers. Well, in case you have a multiple of them (for k>2). Note that only one of them is used (namely the one for [0,1]) for interpolating pixels. The source pixels are set to (y0, ..., y(k-1)).

[MfA]

The polynomials in the code aren't used as blending functions, it's a piecewise cubic interpolator ... they are used as a convolution kernel (a non C1 continuous convolution kernel).

This whole methodology seems a weird mishmash of disparate theories to me, there is a rhyme to it ... but very little reason.

[Wilbert]

@madshi,

Start with the polynomials

Y_I(x) = g3*x^3 + g2*x^2 + g1*x + g0, for x in [-3,-2]
Y_I(x) = f3*x^3 + f2*x^2 + f1*x + f0, for x in [-2,-1]
Y_I(x) = e3*x^3 + e2*x^2 + e1*x + e0, for x in [-1,0]
Y_I(x) = a3*x^3 + a2*x^2 + a1*x + a0, for x in [0,1]
Y_I(x) = b3*x^3 + b2*x^2 + b1*x + b0, for x in [1,2]
Y_I(x) = c3*x^3 + c2*x^2 + c1*x + c0, for x in [2,3]
Y_I(x) = d3*x^3 + d2*x^2 + d1*x + d0, for x in [3,4]

The boundary conditions should result in:

> eqn:={
> y0 = -27*g3 + 9*g2 - 3*g1 + g0,
> y1 = -8*g3 + 4*g2 - 2*g1 + g0,
> y1 = -8*f3 + 4*f2 - 2*f1 + f0,
> y2 = -f3 + f2 - f1 + f0,
> y2 = -e3 + e2 - e1 + e0,
> y3 = e0,
> y3 = a0,
> y4 = a3 + a2 + a1 + a0,
> y4 = b3 + b2 + b1 + b0,
> y5 = 8*b3 + 4*b2 + 2*b1 + b0,
> y5 = 8*c3 + 4*c2 + 2*c1 + c0,
> y6 = 27*c3 + 9*c2 + 3*c1 + c0,
> y6 = 27*d3 + 9*d2 + 3*d1 + d0,
> y7 = 64*d3 + 16*d2 + 4*d1 + d0,
> 12*g3 - 4*g2 + g1 = 12*f3 - 4*f2 + f1,
> 3*f3 - 2*f2 + f1 = 3*e3 - 2*e2 + e1,
> e1 = a1,
> 3*a3 + 2*a2 + a1 = 3*b3 + 2*b2 + b1,
> 12*b3 + 4*b2 + b1 = 12*c3 + 4*c2 + c1,
> 27*c3 + 6*c2 + c1 = 27*d3 + 6*d2 + d1,
> -12*g3 + 2*g2 = -12*f3 + 2*f2,
> -6*f3 + 2*f2 = -6*e3 + 2*e2, 
> 2*e2 = 2*a2,
> 6*a3 + 2*a2 = 6*b3 + 2*b2,
> 12*b3 + 2*b2 = 12*c3 + 2*c2,
> 18*c3 + 2*c2 = 18*d3 + 2*d2,
> -18*g3 + 2*g2 = 0,
> 24*d3 + 2*d2 = 0};

The remaining Maple calculations for Spline64:

> solution:=solve(eqn,{a0,a1,a2,a3,b0,b1,b2,b3,c0,c1,c2,c3,d0,d1,d2,d3,e0,e1,e2,e3,f0,f1,f2,f3,g0,g1,g2,g3});

> a0s:=subs(solution, a0); a1s:=subs(solution, a1); a2s:=subs(solution, a2); a3s:=subs(solution, a3);
                              a0s := y3


         2340      156       624        26        97       582
  a1s := ---- y4 + ---- y6 - ---- y5 - ---- y7 - ---- y0 + ---- y1
         2911      2911      2911      2911      2911      2911

           2328
         - ---- y2 - 3/2911 y3
           2911


         4050      270       1080       45       168       1008
  a2s := ---- y4 + ---- y6 - ---- y5 - ---- y7 + ---- y0 - ---- y1
         2911      2911      2911      2911      2911      2911

           4032      6387
         + ---- y2 - ---- y3
           2911      2911


           49                24
  a3s := - -- y4 - 6/41 y6 + -- y5 + 1/41 y7 - 1/41 y0 + 6/41 y1
           41                41

           24      49
         - -- y2 + -- y3
           41      41

> w:= a3s*x^3 + a2s*x^2 + a1s*x + a0s;

       /  49                24
  w := |- -- y4 - 6/41 y6 + -- y5 + 1/41 y7 - 1/41 y0 + 6/41 y1
       \  41                41

           24      49   \  3   /4050      270       1080       45
         - -- y2 + -- y3| x  + |---- y4 + ---- y6 - ---- y5 - ---- y7
           41      41   /      \2911      2911      2911      2911

           168       1008      4032      6387   \  2   /2340
         + ---- y0 - ---- y1 + ---- y2 - ---- y3| x  + |---- y4
           2911      2911      2911      2911   /      \2911

           156       624        26        97       582       2328
         + ---- y6 - ---- y5 - ---- y7 - ---- y0 + ---- y1 - ---- y2
           2911      2911      2911      2911      2911      2911

                    \
         - 3/2911 y3| x + y3
                    /

> w0:=coeff(w,y0); w1:=coeff(w,y1); w2:=coeff(w,y2); w3:=coeff(w,y3);

                                3   168   2    97
                  w0 := - 1/41 x  + ---- x  - ---- x
                                    2911      2911


                               3   1008  2   582
                   w1 := 6/41 x  - ---- x  + ---- x
                                   2911      2911


                           24  3   4032  2   2328
                   w2 := - -- x  + ---- x  - ---- x
                           41      2911      2911


                       49  3   6387  2
                 w3 := -- x  - ---- x  - 3/2911 x + 1
                       41      2911

> w4:=coeff(w,y4); w5:=coeff(w,y5); w6:=coeff(w,y6); w7:=coeff(w,y7);

                           49  3   4050  2   2340
                   w4 := - -- x  + ---- x  + ---- x
                           41      2911      2911


                          24  3   1080  2   624
                    w5 := -- x  - ---- x  - ---- x
                          41      2911      2911


                                3   270   2   156
                  w6 := - 6/41 x  + ---- x  + ---- x
                                    2911      2911


                               3    45   2    26
                   w7 := 1/41 x  - ---- x  - ---- x
                                   2911      2911

> w0+w1+w2+w3+w4+w5+w6+w7;

                                  1

Note that (w4,w5,w6,w7) are mirrored versions of (w3,w2,w1,w0). You can see that by replacing (y0,y1,y2,y3) by (y7,y6,y5,y4) and x by (1-x).

[Wilbert]

The polynomials in the code aren't used as blending functions, it's a piecewise cubic interpolator ... they are used as a convolution kernel.

Yes, i agree.

I have been looking at the code from Panorama_Tools in panosrc.zip (which i attached here). In particular at resample.c where the resampling is implemented. I still haven't figured it out how it is used there. Have you ever looked at that source code?

[MfA]

Not really, but then I have to say I don't really see any reason to look for some TRUE purely mathematical way of divining interpolation kernels ... it's an exercise in futility. The optimal kernel is almost always signal dependent and of course metric dependent (ie. a matter of taste). With the exception of some trivial examples such as the optimal interpolation of low order Taylor components (which can be done with 0 error).

The existing functions are given more worth than they are due really, you could play around with the coefficients and get something better or worse or the same looking ... in linear interpolators some things almost always make sense (symmetry, zero weights for integer positions not 0 and looking mostly like a low order Lanczos kernel) but everything else is pretty arbitrary.

[Egh]

However that is the point in analysis of the spline resampling coefficients. I.e. are the coefficients good enough for video resampling purposes or shall we tweak them to get better results for our case?

For instance, that 360degree rotation test has little relevance to madshi render, as we only do single resampling. for instance, in this case we can allow some oversharpening provided other artifacts are tolerable.

@willbert:

My contact implies you made a mistake when calculating spline64 coefficients.

His solution in Maple:
http://pastebin.ca/1424864

compare to the original equations in the code:
http://pastebin.ca/1424865

Only 4 equations have been calculated, as others are symmetrical.

[MfA]

Analysis of the coefficients gets you nowhere, the coefficients are a means to a cause and their effectiveness can only be derived from experiment.

You can find a optimal two lobe cubic interpolation function (like spline16) in a mean absolute error sense (for some images) in this paper : http://lear.inrialpes.fr/people/trig...v01-subpix.pdf

[lych_necross]

I was told there would be no math!
Given the complexity of Spline64's coefficients, isn't it reasonable to assume that this is where the law of diminishing returns kicks in?

[Egh]

Originally Posted by MfA

Analysis of the coefficients gets you nowhere, the coefficients are a means to a cause and their effectiveness can only be derived from experiment.

I agree. However we need these coefficients in the first place to test with. I.e. best thing now is to try several ways to derive these coefficients and then select best based on the visual tests.

[Wilbert]

@Egh,

Thanks! I made a small mistake in copying the boundary conditions into Maple. I should have checked them first. All is good now (i corrected my posts above).

[mikenadia]

Smelling Catmull-Rom36 or Bicubicresize36

From madshi:Great! yesgrey3 earlier reported that "Akima" splines may be interesting. 
I haven't been able to turn that documentation into a 
coefficient set that I could use for my resampling filters...

http://www.fugroairborne.com/resourc...ension_III.pdf

on page 4: equations 2 to 5 .
This will allow you to only interpolate between points 2 and 3.

Another one is the constrained spline (with different derivatives at endpoints) described in
http://www.korf.co.uk/spline.pdf (interesting because no need to solve a system of equations)

Other ideas (preserving convexity over C1 continuity , locally monotonous constraint...) in
http://math.lanl.gov/~mac/papers/numerics/H83.pdf

[MfA]

Mikenadia, I fail to see the difference between Akima and Catmull-Rom spline interpolation.

[IanB]

@Wilbert,

Pretty cool stuff Now we finally have our definition.

[Egh]

Originally Posted by MfA

Mikenadia, I fail to see the difference between Akima and Catmull-Rom spline interpolation.

i thought catmull-rom is just bicubic as written in the old thread (b=0, с=0.5). How come Akima splines are just bicubic then?

[MfA]

The funny thing about the Catmull-Rom spline is the amount of different methods which lead to it.

[mikenadia]

according to http://www.cs.cmu.edu/~fp/courses/gr...catmullRom.pdf

the difference between akima and catmull-rom is the way the tangents are computed at point 2 and 3 .
in akima . P1P2 in point 2
in catmull-rom. P1P3 in point 2 ( if tension=0.5) .

But obviously, they belong to the same family.