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]

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... :o

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,

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.
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/triggs/pubs/Triggs-iccv01-subpix.pdf

[lych_necross]

:eek: 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]

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/resources/technical_notes/time_domain_em/pdfs/Akima_tension_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 :D Now we finally have our definition.

[Egh]

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? :confused:

[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/graphics/asst5/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.

[MfA]

Hmm? The first tangent for the Akima spline constraints on page 4 of that pdf you linked was :

1/2*(z2-z1/y2-y1+z3-z2/y3-y2), so essentially the tangent between P1 and P3 with uniform control points (y2-y1=y3-y2).

[MfA]

BTW, I should add that as far as linear interpolators go this is all really old hat ... stuff like this is a little more fresh :
http://www-ist.massey.ac.nz/elai/Publications/ISPACS_beilei_2008.pdf

This method is pretty cute and easy to implement (basically doing a little unsharp before interpolation) :
http://www.ing.unibs.it/~marco.dalai/pub/DLM_VLBV05.pdf

[mikenadia]

Trying to demonstrate that some of the weights equal 1.
I assume that all the weights are polynomial order n-1 and do not depends on the yi values.
Will take then all the yi=1 and try to prove that the polynomial G(x)= (sum weights) -1 has n zeros
and therefore is null.
Based on this paper [http://online.redwoods.cc.ca.us/instruct/darnold/laproj/Fall98/SkyMeg/Proj.PDF (equation 23),
for i=1to n-1: ai+bi+ci=0 di=1 (h=1 and all the yi=1)
Si(xi+h)=ai(h)^3 +bi(h)^2+ci(h)+di=1 because h=1 and Si(xi+1)=sum of weigths(xi) because all the yi=1
so all the xi(for 1 to n-1) are zero of G(x).
Also S1(xn)=dn-1=1 so xn is also a zero of G(x).
so G(x) has n zeros (polynomial order n-1) , and G(x)=0.

Edit:different algorithm for Akima ( that is not a Catmull-Rom ) from the one from "fugroairbone": http://www.iue.tuwien.ac.at/phd/rottinger/node60.html
They compute different derivatives at all points and interpolation is a parabole in the "end-intervals.
They use five points to compute derivatives. It is probably too hard compared to Catmull-Rom (only two points)).
We probably could modify Catmull-Rom to make it harder . Tangent in 2 = (s3-s1)/2.We could take (s3-s1)/k or k*ln(s3/s1) .
We couls also compute first and second derivative at all points and based on a "magic formula", decide to increase radius, to favor natural cubic spline or Catmull-Rom.

Edit2: I read that pixels that are in a low-gradient environment should be given more weight.
May be when we interpolate along the axis X, we should keep an "indicator" of the gradient along X of the interpolated "pixel" and before interpolating along Y, take into account those "gradient along X" and normalize.
And if there is value to it, is there a difference (human perception) by starting the interpolation along the X or Y axis.

[Gavino]

Trying to demonstrate that some [sum!] of the weights equal 1.
Here's what I think is a simpler proof:

If all the yi are equal (to k, say), then it is clear that the solution Si(x)=k (for all i,x) satisfies the initial constraints.
(intuitively, interpolating a set of constant values gives that same value).
Since S1(x) (= sum(Wn(x)*yn)) in this case reduces to k times the sum of the weights, it follows that the sum is 1.

[madshi]

BTW, I should add that as far as linear interpolators go this is all really old hat ... stuff like this is a little more fresh :
http://www-ist.massey.ac.nz/elai/Publications/ISPACS_beilei_2008.pdf

This method is pretty cute and easy to implement (basically doing a little unsharp before interpolation) :
http://www.ing.unibs.it/~marco.dalai/pub/DLM_VLBV05.pdf
Thanks, these papers do look interesting. Unfortunately "easy to implement" probably applies to math professors. At least me, being a programmer with only school math knowledge, I don't understand most of those papers... :( I guess there's no source code available for any of these techniques? :o

[MfA]

Well I didn't say the first one was simple :) For the second one you simply apply a FIR filter with taps 7/720, -11/90, 49/40, -11/90, 7/720 in both dimensions and then you do bi-linear interpolation.

[madshi]

Well I didn't say the first one was simple :)
:D

For the second one you simply apply a FIR filter with taps 7/720, -11/90, 49/40, -11/90, 7/720 in both dimensions and then you do bi-linear interpolation.
That sounds easy enough!! But are those taps identical for any up- and downscaling ratios? Normally when doing (big) downscaling a lot more taps are used than when doing upscaling...

Thanks!

[MfA]

No, but as I've said before interpolation and resizing are two fundamentally different tasks. For downsizing you can simply add an "appropriate" gaussian blur before doing interpolation.

[madshi]

No, but as I've said before interpolation and resizing are two fundamentally different tasks. For downsizing you can simply add an "appropriate" gaussian blur before doing interpolation.
Ok, but the FIR configuration you posted would be alright for any upscaling factor, correct?

For downscaling should I also use Bilinear as the interpolation filter (after having done Gaussian blur)? To be honest, I'm not sure about doing Gaussian blur + Bilinear. I have already tried using Gaussian resampling for downscaling and the results were much too soft.

[MfA]

As I said, appropriate amount ... you can make it as strong or weak as you want.

Direct Gaussian resampling only makes sense for very large ratios.

[madshi]

BTW, I should add that as far as linear interpolators go this is all really old hat ...
Maybe it's old hat. But I've just tried the 2nd "fresh" paper which you kindly explained to me. And I don't like the results much. Old hat Spline36 looks worlds better to me. For those interested, here is a photo upsampled 425% by various resampling filters:

http://madshi.net/prefiltering.rar (25MB)

- Bilinear.bmp
- Prefiltering + Bilinear.bmp
- Bicubic50.bmp
- Bicubic75.bmp
- Spline36.bmp
- Lanczos4 - madshi tweaked.bmp

The "Prefiltering + Bilinear" is what MfA's 2nd "fresh" paper produces. It's a lot less blurry compared to simple Bilinear resampling, but it has quite visible ringing and isn't any better in terms of aliasing compared to simple Bilinear filtering. Spline36 has a similar amount of ringing, but has *MUCH* reduced aliasing.

---------------------

Sooooo, back to Spline16/36/64/256? :)

[Egh]

Sooooo, back to Spline16/36/64/256?

has somebody calculated spline256 coefficients by now? I was also thinking about less computationally expensive option of spline144....

[Gavino]

I was also thinking about less computationally expensive option of spline144....
Well, Spline100 (5 taps) would actually be the next one in the sequence. ;)

[madshi]

has somebody calculated spline256 coefficients by now? I was also thinking about less computationally expensive option of spline144....
In my experience the more taps you use, the less obvious the difference is. For my eyes going from 3 to 4 taps is a bigger difference than going from 4 taps to 8 taps. Oh well, I'm talking about positive differences here. Negative differences (ringing) are quite obvious when going from 4 taps to 8 taps. Personally, I don't think more than 4 taps is really useful, because the positive aspects are extremely small while the negative aspects are quite noticeable. But that's just my personal opinion...

[MfA]

For those interested, here is a photo upsampled 425% by various resampling filters
That's not pure interpolation though ... with that amount of resizing you are testing the filters on how believably they can make up extra detail.

[Gavino]

That's not pure interpolation though ...
Why not? What do you mean by 'pure interpolation' then?

[MfA]

Small ratios, shifts, rotations.

For large ratio upsizing something like Tritical's filters seem more appropriate.

[madshi]

That's not pure interpolation though ... with that amount of resizing you are testing the filters on how believably they can make up extra detail.
Even if I upscale by "only" 33% (anamorphic stretch) I still get similar results, the differences are just smaller. I think the main point of Spline16/36/64 (which is the topic of this thread) is doing DVD -> HD upscaling. That's still over 200% enlargement. So I think my comparison screenshots are valid.

I see in the paper that they tested by rotating images. Maybe the results would be different in that case, I don't know. But I do find it kind of funny that both papers seem to use some simple cubic interpolation (probably bicubic 0.5) for the "conventional" interpolation comparison screenshots. And then they wonder why the results are so blurry?! Are they not aware that e.g. Lanczos exists? Which usually produces much better results when doing multiple operations after each other (like rotating multiple times or upscaling/downscaling multiple times)...

[MfA]

Edit:different algorithm for Akima ( that is not a Catmull-Rom ) from the one from "fugroairbone"
Unfortunately this one no longer simplifies to a convolution kernel ... so it's going to be quite slow.

[tetsuo55]

- Lanczos4 - madshi tweaked.bmp

Are you using that one in MadVR?

[madshi]

Are you using that one in MadVR?
Nope.

[tetsuo55]

Nope.

too bad, i hope you will

[IanB]

@Wilbert,

From your post 7 you missed out on the polynomial for the [-4,-3]segment. Also some of the constraints were not present. I have marked these in Red. Also some of the constraints you have expressed in negated form, these I have adjusted these and marked them in Blue, the original expression is above in Light Gray.
Y_I(x) = h3*x3 + h2*x2 + h1*x + h0, for x in [-4,-3]
Y_I(x) = g3*x3 + g2*x2 + g1*x + g0, for x in [-3,-2]
Y_I(x) = f3*x3 + f2*x2 + f1*x + f0, for x in [-2,-1]
Y_I(x) = e3*x3 + e2*x2 + e1*x + e0, for x in [-1,0]
Y_I(x) = a3*x3 + a2*x2 + a1*x + a0, for x in [0,1]
Y_I(x) = b3*x3 + b2*x2 + b1*x + b0, for x in [1,2]
Y_I(x) = c3*x3 + c2*x2 + c1*x + c0, for x in [2,3]
Y_I(x) = d3*x3 + d2*x2 + d1*x + d0, for x in [3,4]
Analysing the set of constraints it appears that half are redundant, this make sense as the overall function is symmetrical about zero.

Extracting the sets of constraint equations to do with each polynomial and pairing it with it's reflection about zero gives :-

Equations for ai & ei y2 = -e3 + e2 - e1 + e0
y4 = a3 + a2 + a1 + a0

y3 = e0
y3 = a0

3*f3 - 2*f2 + f1 = 3*e3 - 2*e2 + e1
-3*e3 + 2*e2 - e1 = -3*f3 + 2*f2 - f1
3*a3 + 2*a2 + a1 = 3*b3 + 2*b2 + b1

-6*f3 + 2*f2 = -6*e3 + 2*e2
6*a3 + 2*a2 = 6*b3 + 2*b2

e1 = a1

2*e2 = 2*a2
Equations for bi & fi y1 = -8*f3 + 4*f2 - 2*f1 + f0
y5 = 8*b3 + 4*b2 + 2*b1 + b0

y2 = -f3 + f2 - f1 + f0
y4 = b3 + b2 + b1 + b0

12*g3 - 4*g2 + g1 = 12*f3 - 4*f2 + f1
-12*f3 + 4*f2 - f1 = -12*g3 + 4*g2 - g1
12*b3 + 4*b2 + b1 = 12*c3 + 4*c2 + c1

-12*g3 + 2*g2 = -12*f3 + 2*f2
12*b3 + 2*b2 = 12*c3 + 2*c2

3*f3 - 2*f2 + f1 = 3*e3 - 2*e2 + e1
-3*e3 + 2*e2 - e1 = -3*f3 + 2*f2 - f1
3*a3 + 2*a2 + a1 = 3*b3 + 2*b2 + b1

-6*f3 + 2*f2 = -6*e3 + 2*e2
6*a3 + 2*a2 = 6*b3 + 2*b2
Equations for ci & gi y0 = -27*g3 + 9*g2 - 3*g1 + g0
y6 = 27*c3 + 9*c2 + 3*c1 + c0

y1 = -8*g3 + 4*g2 - 2*g1 + g0
y5 = 8*c3 + 4*c2 + 2*c1 + c0

-27*g3 + 6*g2 + g1 = -27*h3 + 6*h2 + h1
27*c3 + 6*c2 + c1 = 27*d3 + 6*d2 + d1

-18*g3 + 2*g2 = -18*h3 + 2*h2
18*c3 + 2*c2 = 18*d3 + 2*d2

-18*g3 + 2*g2 = 0
18*c3 + 2*c2 = 0

12*g3 - 4*g2 + g1 = 12*f3 - 4*f2 + f1
-12*f3 + 4*f2 - f1 = -12*g3 + 4*g2 - g1
12*b3 + 4*b2 + b1 = 12*c3 + 4*c2 + c1

-12*g3 + 2*g2 = -12*f3 + 2*f2
12*b3 + 2*b2 = 12*c3 + 2*c2
Equations for di & hi y-1= -64*h3 + 16*h2 - 4*h1 + h0
y7 = 64*d3 + 16*d2 + 4*d1 + d0

y0 = -27*h3 + 9*h2 - 3*h1 + h0
y6 = 27*d3 + 9*d2 + 3*d1 + d0

-24*h3 + 2*h2 = 0
24*d3 + 2*d2 = 0

-27*g3 + 6*g2 + g1 = -27*h3 + 6*h2 + h1
27*c3 + 6*c2 + c1 = 27*d3 + 6*d2 + d1

-18*g3 + 2*g2 = -18*h3 + 2*h2
18*c3 + 2*c2 = 18*d3 + 2*d2

Thus the system should be solvable with just this half of the equations.> eqn:={
> 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,
> 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,
> 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*c3 + 2*c2 = 0,
> 24*d3 + 2*d2 = 0};Can you please confirm this with your Maple Wilbert.

Also I think this constraint is also available 12*b3 + 2*b2 = 0, not sure about the 12 maybe it should be 9, my brain is quite rusty on this stuff.

Also we know a0=1 and b0,c0,d0=0 so we may be able to reduce the equation clutter even more, making it feasible to deduce spline100, spline144, spline196 and spline256 coefficients

[Wilbert]

Just a quick response,
From your post 7 you missed out on the polynomial for the [-4,-3]segment.
The author of Panorama Tools doesn't define a polynomial on [-4,-3]. If you do it correctly (using 8 splines, setting yj=0 except at the center), you will get this one: http://forums.virtualdub.org/index.php?act=ST&f=11&t=17212& (see post: May 20 2009, 02:44 AM). I don't see any difference with the existing Spline64 though.

Getting the spline100, spline144, spline196 and spline256 coefficients (using the correct number of splines) is trivial in Maple. I will post them tomorrow. It's also possible to calculate an arbitrary one. If i've time i will try to implement that. I will comment on your boundary conditions later.

[Gavino]

@Wilbert,
From your post 7 you missed out on the polynomial for the [-4,-3]segment.
No, Wilbert is correct.
There is no [-4,-3] segment. For Spline 64, we have eight points, hence seven intervals/polynomials.
Also I think this constraint is also available 12*b3 + 2*b2 = 0, not sure about the 12 maybe it should be 9, my brain is quite rusty on this stuff.
Can you say where this constraint comes from?
Also we know a0=1 and b0,c0,d0=0 so we may be able to reduce the equation clutter even more, making it feasible to deduce spline100, spline144, spline196 and spline256 coefficients
If you use a tool like Maple, does the equation clutter matter?
Unless the tool has a limit on the size of equations it can solve, isn't it just a matter of following Wilbert's method with the appropriate number of polynomials?

[Gavino]

The author of Panorama Tools doesn't define a polynomial on [-4,-3]. If you do it correctly (using 8 splines, setting yj=0 except at the center), you will get this one: http://forums.virtualdub.org/index.php?act=ST&f=11&t=17212& (see post: May 20 2009, 02:44 AM). I don't see any difference with the existing Spline64 though.
Wilbert, I hadn't seen your reply when I made the previous post.
After looking at the vdub forum thread, I think there is some confusion over the underlying process. Your method (from Panorama tools) fits a spline through the sample points and then derives the filter kernel from the resulting blending polynomials. The approach suggested by phaeron derives the kernel itself as a spline. It's not clear to me which is more 'correct', or if there is an equivalence between the two approaches at some level.

[Wilbert]

After looking at the vdub forum thread, I think there is some confusion over the underlying process. Your method (from Panorama tools) fits a spline through the sample points and then derives the filter kernel from the resulting blending polynomials. The approach suggested by phaeron derives the kernel itself as a spline.
Yes, i know :) I should have made that more clear in my post above.

[Wilbert]

Spline100:


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

3 209 2 362
w0 := 1/153 x - ----- x + ----- x
13515 40545


3 418 2 724
w1 := - 2/51 x + ---- x - ----- x
4505 13515


3 1672 2 2896
w2 := 8/51 x - ---- x + ----- x
4505 13515


10 3 1254 2 724
w3 := - -- x + ---- x - --- x
17 901 901


61 3 9893 2
w4 := -- x - ---- x - 1/13515 x + 1
51 4505


Spline144:

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

3 2340 2 1351
w0 := - 1/571 x + ------ x - ------ x
564719 564719


3 14040 2 8106
w1 := 6/571 x - ------ x + ------ x
564719 564719


24 3 56160 2 32424
w2 := - --- x + ------ x - ------ x
571 564719 564719


90 3 210600 2 121590
w3 := --- x - ------ x + ------ x
571 564719 564719


336 3 786240 2 453936
w4 := - --- x + ------ x - ------ x
571 564719 564719


683 3 1240203 2
w5 := --- x - ------- x - 3/564719 x + 1
571 564719


The weights sum up to one for both splines (including their symmetric counterparts), so i suspect that they are correct.

[Gavino]

Spline100: ...
Spline144: ...
From the way these are expressed, I assume that they are calculated using the 'Panorama Tools' method (as for existing Spline16, etc), rather than the alternative method of fitting a spline through the kernel crossing points.

Any thoughts on the relative merits of the two approaches? It is clear that they give different coefficients and I suspect they are fundamentally different in nature, although they may both be equally valid as approximations to a sinc filter.

[Wilbert]

Sorry for the late response.

From the way these are expressed, I assume that they are calculated using the 'Panorama Tools' method (as for existing Spline16, etc), rather than the alternative method of fitting a spline through the kernel crossing points.
Yes, indeed.

Any thoughts on the relative merits of the two approaches? It is clear that they give different coefficients and I suspect they are fundamentally different in nature, although they may both be equally valid as approximations to a sinc filter.
The coefficients seem to be close to each other (well i only looked at Spline64, dunno about the other ones). I implemented both resizers here: http://www.geocities.com/wilbertdijkhof/SplineResize_v01.zip

For the alternate method you can choose the number of "taps". I leave it up to you guys to find and report about the differences :)