Welcome to Doom9's Forum, THE inplace 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. 
15th May 2009, 23:54  #1  Link 
Moderator
Join Date: Nov 2001
Location: Netherlands
Posts: 6,312

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] 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](1x) 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 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 Code:
w0(x)+w1(x)=1 Because of the mirroring property we also have that Code:
S(x) = w0(1x) * y1 + w1(1x) * y0 Code:
w1(x) = w0(1x) w0(x) = w1(1x) Code:
w1(x) = w0(1x) > w0[x] 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}; Code:
> solve(eqn,{a0,a1,a2,a3}); {a0 = y0, a2 = 0, a3 = 0, a1 = y0 + y1} 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 Two additional remarks: 1) In general for Spline k^2, we start with k1 interpolating polynomials. They are located at (k/2+1,k/2+2), ..., (k/21,k/2). The bending polynomials are calculated from the polynomial located at [0,1], which is expressed as a linear combination of (y0, ..., y(k1)). The boundary conditions are: Code:
y[1](k/2+1) = y0 y[1](k/2+2) = y1 y[2](k/2+2) = y1 ... y[k2](k/21) = y(k2) y[k2](k/21) = y(k2) y[k1](k/2) = y(k1) 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. 
16th May 2009, 07:31  #4  Link  
Registered Developer
Join Date: Sep 2006
Posts: 9,126

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

16th May 2009, 16:05  #5  Link  
Moderator
Join Date: Nov 2001
Location: Netherlands
Posts: 6,312

@MfA,
Quote:
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(k1)). Last edited by Wilbert; 16th May 2009 at 16:08. 

16th May 2009, 16:17  #6  Link 
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. 
16th May 2009, 16:35  #7  Link 
Moderator
Join Date: Nov 2001
Location: Netherlands
Posts: 6,312

@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] 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}; 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 Last edited by Wilbert; 16th May 2009 at 23:02. 
16th May 2009, 16:44  #8  Link  
Moderator
Join Date: Nov 2001
Location: Netherlands
Posts: 6,312

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

16th May 2009, 17:11  #9  Link 
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. 
16th May 2009, 20:58  #10  Link 
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. 
16th May 2009, 21:29  #11  Link 
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...v01subpix.pdf 
16th May 2009, 23:16  #15  Link 
Registered User
Join Date: Nov 2007
Posts: 246

Smelling CatmullRom36 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... 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. 
17th May 2009, 02:33  #20  Link 
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 catmullrom is the way the tangents are computed at point 2 and 3 . in akima . P1P2 in point 2 in catmullrom. 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. 
Thread Tools  Search this Thread 
Display Modes  

