Welcome to Doom9's Forum, THE in-place to be for everyone interested in DVD conversion.

Before you start posting please read the forum rules. By posting to this forum you agree to abide by the rules.

 

Go Back   Doom9's Forum > Capturing and Editing Video > Avisynth Development

Reply
 
Thread Tools Search this Thread Display Modes
Old 15th May 2009, 23:54   #1  |  Link
Wilbert
Moderator
 
Join Date: Nov 2001
Location: Netherlands
Posts: 6,293
on the theory of the Spline16/36/64 resizers

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]:
Code:
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:
Code:
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):
Code:
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:
Code:
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
Code:
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
Code:
S(x) = w0(1-x) * y1 + w1(1-x) * y0
It follows that
Code:
w1(x) = w0(1-x)
w0(x) = w1(1-x)
Subtracting the segment bias results in
Code:
w1(x) = w0(1-x) -> w0[-x]
getting a symmetrical kernel.

Solving (w0,w1) using Maple:
The boundary conditions result in:
Code:
> eqn:={a0=y0, a3+a2+a1+a0=y1, 2*a2=0,6*a3+2*a2=0};
Solving for (a[j])_j:
Code:
> 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):
Code:
> 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:
Code:
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.

Last edited by Wilbert; 16th May 2009 at 22:52.
Wilbert is offline   Reply With Quote
Old 16th May 2009, 00:43   #2  |  Link
Egh
Registered User
 
Join Date: Jun 2005
Posts: 630
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.
Egh is offline   Reply With Quote
Old 16th May 2009, 02:56   #3  |  Link
MfA
Registered User
 
Join Date: Mar 2002
Posts: 1,075
Quote:
Originally Posted by Wilbert View Post
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.
MfA is offline   Reply With Quote
Old 16th May 2009, 07:31   #4  |  Link
madshi
Registered Developer
 
Join Date: Sep 2006
Posts: 8,797
@Wilbert, good job! Although I barely understand most of what you wrote...

Quote:
Originally Posted by Wilbert View Post
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!
madshi is offline   Reply With Quote
Old 16th May 2009, 16:05   #5  |  Link
Wilbert
Moderator
 
Join Date: Nov 2001
Location: Netherlands
Posts: 6,293
@MfA,
Quote:
Quote:
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)).

Last edited by Wilbert; 16th May 2009 at 16:08.
Wilbert is offline   Reply With Quote
Old 16th May 2009, 16:17   #6  |  Link
MfA
Registered User
 
Join Date: Mar 2002
Posts: 1,075
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.

Last edited by MfA; 16th May 2009 at 16:36.
MfA is offline   Reply With Quote
Old 16th May 2009, 16:35   #7  |  Link
Wilbert
Moderator
 
Join Date: Nov 2001
Location: Netherlands
Posts: 6,293
@madshi,

Start with the polynomials
Code:
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:
Code:
> 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:

Code:
> 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).

Last edited by Wilbert; 16th May 2009 at 23:02.
Wilbert is offline   Reply With Quote
Old 16th May 2009, 16:44   #8  |  Link
Wilbert
Moderator
 
Join Date: Nov 2001
Location: Netherlands
Posts: 6,293
Quote:
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?
Attached Files
File Type: zip panosrc.zip (157.0 KB, 235 views)
Wilbert is offline   Reply With Quote
Old 16th May 2009, 17:11   #9  |  Link
MfA
Registered User
 
Join Date: Mar 2002
Posts: 1,075
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.

Last edited by MfA; 16th May 2009 at 17:15.
MfA is offline   Reply With Quote
Old 16th May 2009, 20:58   #10  |  Link
Egh
Registered User
 
Join Date: Jun 2005
Posts: 630
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.

Last edited by Egh; 16th May 2009 at 21:14.
Egh is offline   Reply With Quote
Old 16th May 2009, 21:29   #11  |  Link
MfA
Registered User
 
Join Date: Mar 2002
Posts: 1,075
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
MfA is offline   Reply With Quote
Old 16th May 2009, 22:27   #12  |  Link
lych_necross
ZZZzzzz...
 
lych_necross's Avatar
 
Join Date: Jan 2007
Location: USA
Posts: 302
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?
lych_necross is offline   Reply With Quote
Old 16th May 2009, 22:47   #13  |  Link
Egh
Registered User
 
Join Date: Jun 2005
Posts: 630
Quote:
Originally Posted by MfA View Post
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.
Egh is offline   Reply With Quote
Old 16th May 2009, 22:53   #14  |  Link
Wilbert
Moderator
 
Join Date: Nov 2001
Location: Netherlands
Posts: 6,293
@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).
Wilbert is offline   Reply With Quote
Old 16th May 2009, 23:16   #15  |  Link
mikenadia
Registered User
 
Join Date: Nov 2007
Posts: 246
Smelling Catmull-Rom36 or Bicubicresize36
Code:
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

Last edited by mikenadia; 20th May 2009 at 01:54.
mikenadia is offline   Reply With Quote
Old 17th May 2009, 00:29   #16  |  Link
MfA
Registered User
 
Join Date: Mar 2002
Posts: 1,075
Mikenadia, I fail to see the difference between Akima and Catmull-Rom spline interpolation.

Last edited by MfA; 17th May 2009 at 00:51.
MfA is offline   Reply With Quote
Old 17th May 2009, 00:37   #17  |  Link
IanB
Avisynth Developer
 
Join Date: Jan 2003
Location: Melbourne, Australia
Posts: 3,168
@Wilbert,

Pretty cool stuff Now we finally have our definition.
IanB is offline   Reply With Quote
Old 17th May 2009, 01:32   #18  |  Link
Egh
Registered User
 
Join Date: Jun 2005
Posts: 630
Quote:
Originally Posted by MfA View Post
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?
Egh is offline   Reply With Quote
Old 17th May 2009, 02:08   #19  |  Link
MfA
Registered User
 
Join Date: Mar 2002
Posts: 1,075
The funny thing about the Catmull-Rom spline is the amount of different methods which lead to it.
MfA is offline   Reply With Quote
Old 17th May 2009, 02:33   #20  |  Link
mikenadia
Registered User
 
Join Date: Nov 2007
Posts: 246
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.

Last edited by mikenadia; 17th May 2009 at 13:58. Reason: MfA: you are right all along. My mistake.
mikenadia is offline   Reply With Quote
Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off

Forum Jump


All times are GMT +1. The time now is 18:14.


Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2017, vBulletin Solutions Inc.