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. Doom9's Forum on the theory of the Spline16/36/64 resizers
 Register FAQ Calendar Search Today's Posts Mark Forums Read  15th May 2009, 23:54 #1  |  Link Wilbert Moderator Join Date: Nov 2001 Location: Netherlands Posts: 6,341 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*x^3 + a*x^2 + a*x + a 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(-k/2+1) = y0 y(-k/2+2) = y1 y(-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.   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.   16th May 2009, 02:56   #3  |  Link
MfA
Registered User

Join Date: Mar 2002
Posts: 1,075
Quote:
 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.   16th May 2009, 07:31   #4  |  Link
Registered Developer

Join Date: Sep 2006
Posts: 9,137
@Wilbert, good job! Although I barely understand most of what you wrote... Quote:
 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!   16th May 2009, 16:05   #5  |  Link
Wilbert
Moderator Join Date: Nov 2001
Location: Netherlands
Posts: 6,341
@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.   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.   16th May 2009, 16:35 #7  |  Link Wilbert Moderator Join Date: Nov 2001 Location: Netherlands Posts: 6,341 @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.   16th May 2009, 16:44   #8  |  Link
Wilbert
Moderator Join Date: Nov 2001
Location: Netherlands
Posts: 6,341
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 panosrc.zip (157.0 KB, 271 views)   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.   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.   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   16th May 2009, 22:27 #12  |  Link lych_necross ZZZzzzz...   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?   16th May 2009, 22:47   #13  |  Link
Egh
Registered User

Join Date: Jun 2005
Posts: 630
Quote:
 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.   16th May 2009, 22:53 #14  |  Link Wilbert Moderator Join Date: Nov 2001 Location: Netherlands Posts: 6,341 @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).   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.   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.   17th May 2009, 00:37 #17  |  Link IanB Avisynth Developer   Join Date: Jan 2003 Location: Melbourne, Australia Posts: 3,173 @Wilbert, Pretty cool stuff Now we finally have our definition.   17th May 2009, 01:32   #18  |  Link
Egh
Registered User

Join Date: Jun 2005
Posts: 630
Quote:
 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?    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.   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.   Thread Tools Search this Thread Show Printable Version Email this Page Search this Thread: Advanced Search Display Modes Linear Mode Switch to Hybrid Mode Switch to Threaded Mode 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 Rules
 Forum Jump User Control Panel Private Messages Subscriptions Who's Online Search Forums Forums Home Announcements and Chat     General Discussion     News     Forum / Site Suggestions & Help General     Decrypting     Newbies     DVD2AVI / DGIndex     Audio encoding     Subtitles     Linux, Mac OS X, & Co Capturing and Editing Video     Avisynth Usage     Avisynth Development     VapourSynth     Capturing Video     DV     HDTV / DVB / TiVo     NLE - Non Linear Editing     VirtualDub, VDubMod & AviDemux     New and alternative a/v containers Video Encoding     (Auto) Gordian Knot     MPEG-4 ASP     MPEG-4 Encoder GUIs     MPEG-4 AVC / H.264     High Efficiency Video Coding (HEVC)     New and alternative video codecs     MPEG-2 Encoding     VP9 and AV1 (HD) DVD, Blu-ray & (S)VCD     One click suites for DVD backup and DVD creation     DVD & BD Rebuilder     (HD) DVD & Blu-ray authoring     Advanced authoring     IFO/VOB Editors     DVD burning Hardware & Software     Software players     Hardware players     PC Hard & Software Programming and Hacking     Development     Translations

All times are GMT +1. The time now is 08:56.

 Doom9.org - Archive - Top