<p>It turns out that no, while you can use bicubic Lagrange interpolation for bicubic texture sampling, it isn't the highest quality option, and probably not actually likely to be used.</p>

<p>Cubic hermite splines are a better tool for the job.</p>

<p>Lagrange interpolation will make a curve that passes through the data points, thus preserving C0 continuity, but hermite splines preserve the derivatives at the edges while also passing through the data points, thus preserving C1 continuity and looking much better.</p>

<p>This question has some great info on cubic hermite splines:
<a href="https://dsp.stackexchange.com/questions/18265/bicubic-interpolation">https://dsp.stackexchange.com/questions/18265/bicubic-interpolation</a></p>

<p>Here is the cubic hermite version of the code I posted in the question:</p>

<pre class="lang-cplusplus prettyprint-override"><code>//=======================================================================================
vec3 CubicHermite (vec3 A, vec3 B, vec3 C, vec3 D, float t)
{
    float t2 = t*t;
    float t3 = t*t*t;
    vec3 a = -A/2.0 + (3.0*B)/2.0 - (3.0*C)/2.0 + D/2.0;
    vec3 b = A - (5.0*B)/2.0 + 2.0*C - D / 2.0;
    vec3 c = -A/2.0 + C/2.0;
    vec3 d = B;

    return a*t3 + b*t2 + c*t + d;
}

//=======================================================================================
vec3 BicubicHermiteTextureSample (vec2 P)
{
    vec2 pixel = P * c_textureSize + 0.5;

    vec2 frac = fract(pixel);
    pixel = floor(pixel) / c_textureSize - vec2(c_onePixel/2.0);

    vec3 C00 = texture2D(iChannel0, pixel + vec2(-c_onePixel ,-c_onePixel)).rgb;
    vec3 C10 = texture2D(iChannel0, pixel + vec2( 0.0        ,-c_onePixel)).rgb;
    vec3 C20 = texture2D(iChannel0, pixel + vec2( c_onePixel ,-c_onePixel)).rgb;
    vec3 C30 = texture2D(iChannel0, pixel + vec2( c_twoPixels,-c_onePixel)).rgb;

    vec3 C01 = texture2D(iChannel0, pixel + vec2(-c_onePixel , 0.0)).rgb;
    vec3 C11 = texture2D(iChannel0, pixel + vec2( 0.0        , 0.0)).rgb;
    vec3 C21 = texture2D(iChannel0, pixel + vec2( c_onePixel , 0.0)).rgb;
    vec3 C31 = texture2D(iChannel0, pixel + vec2( c_twoPixels, 0.0)).rgb;    

    vec3 C02 = texture2D(iChannel0, pixel + vec2(-c_onePixel , c_onePixel)).rgb;
    vec3 C12 = texture2D(iChannel0, pixel + vec2( 0.0        , c_onePixel)).rgb;
    vec3 C22 = texture2D(iChannel0, pixel + vec2( c_onePixel , c_onePixel)).rgb;
    vec3 C32 = texture2D(iChannel0, pixel + vec2( c_twoPixels, c_onePixel)).rgb;    

    vec3 C03 = texture2D(iChannel0, pixel + vec2(-c_onePixel , c_twoPixels)).rgb;
    vec3 C13 = texture2D(iChannel0, pixel + vec2( 0.0        , c_twoPixels)).rgb;
    vec3 C23 = texture2D(iChannel0, pixel + vec2( c_onePixel , c_twoPixels)).rgb;
    vec3 C33 = texture2D(iChannel0, pixel + vec2( c_twoPixels, c_twoPixels)).rgb;    

    vec3 CP0X = CubicHermite(C00, C10, C20, C30, frac.x);
    vec3 CP1X = CubicHermite(C01, C11, C21, C31, frac.x);
    vec3 CP2X = CubicHermite(C02, C12, C22, C32, frac.x);
    vec3 CP3X = CubicHermite(C03, C13, C23, C33, frac.x);

    return CubicHermite(CP0X, CP1X, CP2X, CP3X, frac.y);
}
</code></pre>

<p>Here is a picture showing the difference between sampling methods.  From left to right: Nearest Neighbor, Bilinear, Lagrange Bicubic, Hermite Bicubic</p>

<p><a href="https://i.sstatic.net/h0Ot0.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/h0Ot0.png" alt="enter image description here"></a></p>

CC BY-SA 3.0
2015-08-05T14:51:27.773Z
Although all cubic splines are, in a sense, equivalent, it's probably conceptually easier to use Catmull-Rom splines. e.g.  https://www.cs.cmu.edu/~462/projects/assn2/assn2/catmullRom.pdf
Simon F
CC BY-SA 3.0
2015-08-05T15:57:32.830Z
Do you think the tau parameter helps or hinders in this case?  I might be wrong but i feel like catmull rom is too "user defined" (and must be tuned), whereas the hermite spline attempts to just use information from the data that's there.  It seems like cubic hermite is closer to a "ground truth", which i guess would be something like an ideal sinc filter.  What do you think though?
Alan Wolfe
CC BY-SA 3.0
2015-08-06T09:24:46.650Z
I don't see how Catmull-Rom is "user defined". Once you have a sequence of 4 contiguous points, P[i-1], P[i], P[i+1], P[i+2] (4x4 for 2D case) the curve segment is defined between P[i] and P[i+1] and is C1 continuous with the neighbouring segments. 

A sinc filter is fine for audio but not video. See Mitchell & Netravali: https://www.cs.utexas.edu/~fussell/courses/cs384g-fall2013/lectures/mitchell/Mitchell.pdf   
IIRC Catmull-Rom is a special case of the family of filters they propose, but I think that filter is an approximating curve so, unlike C-R, might not go through the original points.
Simon F
CC BY-SA 3.0
2015-08-06T14:04:06.173Z
That's how he hermite spline works except that the catmull rom spline has an additional parameter tau (tension) that is user defined.  Also, sinc does apply to video, DSP is DSP :P
Alan Wolfe
CC BY-SA 3.0
2015-08-06T14:29:47.557Z
I must admit, I've never seen a tension parameter associated with Catmull Rom splines before, but then I've really only learned of them via Foley & van Dam (et al) or, say, Watt & Watt which, AFAICR, make no mention of such.

Actually, having said that, given there are four constraints - i.e. the curve has to pass through 2 points and have two defined tangents** at those points and it's a cubic - I'm at a bit of a loss as to how there's any more degrees of freedom to support a tension parameter ....

**Unless you mean the tangents can be scaled?
Simon F
CC BY-SA 3.0
2015-08-06T14:34:42.063Z
The tension parameter is shown in that first paper you linked.  And yeah I just mean, the fact that there's a tunable parameter means it must be tuned :p. There is a value of tau you can use such that you end up with a regular cubic hermite spline, and I think it might be "1" but not 100% sure on that.
Alan Wolfe
CC BY-SA 3.0
2015-08-07T08:25:01.567Z
Re that first paper: That serves me right for pasting in the first link I found that seemed suitable :-)    Again, I think I've only ever seen the tangent at a segment end point set to the (assuming I've done the matrix maths correctly) 1/2 the difference of the two surrounding neighbours.
Simon F
CC BY-SA 3.0
2015-08-07T13:36:17.957Z
Ah OK! Anyways yeah, catmull rom / hermite is really freaking awesome I totally agree. Really neat way to interpolate between data points and have C1 continuity at the edges. It really simplifies things (:
Alan Wolfe
CC BY-SA 3.0
2015-08-04T23:00:18.923Z
2015-08-10T10:40:04.633Z
2017-04-13T12:47:29.743Z
Alan Wolfe
<p>I just implemented some interpolated texture sampling by sampling the 4x4 nearest pixels then doing Lagrange interpolation across the x axis to get four values to use Lagrange interpolation on across the y axis.</p>

<p>Is this the same as bicubic interpolation or is it different?  Or are there different kinds of bicubic interpolation, and this is just one of them perhaps?</p>

<p>Webgl Shadertoy implementation here and relevant GLSL (WebGL) code below: <a href="https://www.shadertoy.com/view/MllSzX" rel="nofollow">https://www.shadertoy.com/view/MllSzX</a></p>

<p>Thanks!</p>

<pre class="lang-cplusplus prettyprint-override"><code>float c_textureSize = 64.0;

float c_onePixel = 1.0 / c_textureSize;
float c_twoPixels = 2.0 / c_textureSize;

float c_x0 = -1.0;
float c_x1 =  0.0;
float c_x2 =  1.0;
float c_x3 =  2.0;

//=======================================================================================
vec3 CubicLagrange (vec3 A, vec3 B, vec3 C, vec3 D, float t)
{
    return
        A * 
        (
            (t - c_x1) / (c_x0 - c_x1) * 
            (t - c_x2) / (c_x0 - c_x2) *
            (t - c_x3) / (c_x0 - c_x3)
        ) +
        B * 
        (
            (t - c_x0) / (c_x1 - c_x0) * 
            (t - c_x2) / (c_x1 - c_x2) *
            (t - c_x3) / (c_x1 - c_x3)
        ) +
        C * 
        (
            (t - c_x0) / (c_x2 - c_x0) * 
            (t - c_x1) / (c_x2 - c_x1) *
            (t - c_x3) / (c_x2 - c_x3)
        ) +       
        D * 
        (
            (t - c_x0) / (c_x3 - c_x0) * 
            (t - c_x1) / (c_x3 - c_x1) *
            (t - c_x2) / (c_x3 - c_x2)
        );
}

//=======================================================================================
vec3 BicubicTextureSample (vec2 P)
{
    vec2 pixel = P * c_textureSize + 0.5;

    vec2 frac = fract(pixel);
    pixel = floor(pixel) / c_textureSize - vec2(c_onePixel/2.0);

    vec3 C00 = texture2D(iChannel0, pixel + vec2(-c_onePixel ,-c_onePixel)).rgb;
    vec3 C10 = texture2D(iChannel0, pixel + vec2( 0.0        ,-c_onePixel)).rgb;
    vec3 C20 = texture2D(iChannel0, pixel + vec2( c_onePixel ,-c_onePixel)).rgb;
    vec3 C30 = texture2D(iChannel0, pixel + vec2( c_twoPixels,-c_onePixel)).rgb;

    vec3 C01 = texture2D(iChannel0, pixel + vec2(-c_onePixel , 0.0)).rgb;
    vec3 C11 = texture2D(iChannel0, pixel + vec2( 0.0        , 0.0)).rgb;
    vec3 C21 = texture2D(iChannel0, pixel + vec2( c_onePixel , 0.0)).rgb;
    vec3 C31 = texture2D(iChannel0, pixel + vec2( c_twoPixels, 0.0)).rgb;    

    vec3 C02 = texture2D(iChannel0, pixel + vec2(-c_onePixel , c_onePixel)).rgb;
    vec3 C12 = texture2D(iChannel0, pixel + vec2( 0.0        , c_onePixel)).rgb;
    vec3 C22 = texture2D(iChannel0, pixel + vec2( c_onePixel , c_onePixel)).rgb;
    vec3 C32 = texture2D(iChannel0, pixel + vec2( c_twoPixels, c_onePixel)).rgb;    

    vec3 C03 = texture2D(iChannel0, pixel + vec2(-c_onePixel , c_twoPixels)).rgb;
    vec3 C13 = texture2D(iChannel0, pixel + vec2( 0.0        , c_twoPixels)).rgb;
    vec3 C23 = texture2D(iChannel0, pixel + vec2( c_onePixel , c_twoPixels)).rgb;
    vec3 C33 = texture2D(iChannel0, pixel + vec2( c_twoPixels, c_twoPixels)).rgb;    

    vec3 CP0X = CubicLagrange(C00, C10, C20, C30, frac.x);
    vec3 CP1X = CubicLagrange(C01, C11, C21, C31, frac.x);
    vec3 CP2X = CubicLagrange(C02, C12, C22, C32, frac.x);
    vec3 CP3X = CubicLagrange(C03, C13, C23, C33, frac.x);

    return CubicLagrange(CP0X, CP1X, CP2X, CP3X, frac.y);
}
</code></pre>

CC BY-SA 3.0
2015-08-04T18:17:21.623Z
You could post the relevant shader code here in case of bitrot, no?
joojaa
CC BY-SA 3.0
2015-08-04T18:21:29.027Z
we should have some prettier code markup for shader code, i'll post on meta if someone hasn't beaten me to it!
Alan Wolfe
CC BY-SA 3.0
2015-08-05T22:24:19.397Z
Is that a specific shader language not available in the list of languages covered by our syntax highlighting?
trichoplax is on Codidact now
CC BY-SA 3.0
2015-08-05T22:37:41.687Z
I'm not sure.  It's just GLSL (from webgl to be exact!).  I just did 4 spaces before each line of code, not sure if there's a better way to mark it up...
Alan Wolfe
<texture><interpolation>
Is a cubic Lagrange interpolation tensor product the same as bicubic interpolation?
Alan Wolfe
# Is a cubic Lagrange interpolation tensor product the same as bicubic interpolation?

(Asked by: username_0 on 2015-08-04)

I just implemented some interpolated texture sampling by sampling the 4x4 nearest pixels then doing Lagrange interpolation across the x axis to get four values to use Lagrange interpolation on across the y axis.

Is this the same as bicubic interpolation or is it different? Or are there different kinds of bicubic interpolation, and this is just one of them perhaps?

Webgl Shadertoy implementation here and relevant GLSL (WebGL) code below: <https://www.shadertoy.com/view/MllSzX>

Thanks!

```cplusplus
float c_textureSize = 64.0;

float c_onePixel = 1.0 / c_textureSize;
float c_twoPixels = 2.0 / c_textureSize;

float c_x0 = -1.0;
float c_x1 =  0.0;
float c_x2 =  1.0;
float c_x3 =  2.0;

//=======================================================================================
vec3 CubicLagrange (vec3 A, vec3 B, vec3 C, vec3 D, float t)
{
    return
        A * 
        (
            (t - c_x1) / (c_x0 - c_x1) * 
            (t - c_x2) / (c_x0 - c_x2) *
            (t - c_x3) / (c_x0 - c_x3)
        ) +
        B * 
        (
            (t - c_x0) / (c_x1 - c_x0) * 
            (t - c_x2) / (c_x1 - c_x2) *
            (t - c_x3) / (c_x1 - c_x3)
        ) +
        C * 
        (
            (t - c_x0) / (c_x2 - c_x0) * 
            (t - c_x1) / (c_x2 - c_x1) *
            (t - c_x3) / (c_x2 - c_x3)
        ) +       
        D * 
        (
            (t - c_x0) / (c_x3 - c_x0) * 
            (t - c_x1) / (c_x3 - c_x1) *
            (t - c_x2) / (c_x3 - c_x2)
        );
}

//=======================================================================================
vec3 BicubicTextureSample (vec2 P)
{
    vec2 pixel = P * c_textureSize + 0.5;

    vec2 frac = fract(pixel);
    pixel = floor(pixel) / c_textureSize - vec2(c_onePixel/2.0);

    vec3 C00 = texture2D(iChannel0, pixel + vec2(-c_onePixel ,-c_onePixel)).rgb;
    vec3 C10 = texture2D(iChannel0, pixel + vec2( 0.0        ,-c_onePixel)).rgb;
    vec3 C20 = texture2D(iChannel0, pixel + vec2( c_onePixel ,-c_onePixel)).rgb;
    vec3 C30 = texture2D(iChannel0, pixel + vec2( c_twoPixels,-c_onePixel)).rgb;

    vec3 C01 = texture2D(iChannel0, pixel + vec2(-c_onePixel , 0.0)).rgb;
    vec3 C11 = texture2D(iChannel0, pixel + vec2( 0.0        , 0.0)).rgb;
    vec3 C21 = texture2D(iChannel0, pixel + vec2( c_onePixel , 0.0)).rgb;
    vec3 C31 = texture2D(iChannel0, pixel + vec2( c_twoPixels, 0.0)).rgb;    

    vec3 C02 = texture2D(iChannel0, pixel + vec2(-c_onePixel , c_onePixel)).rgb;
    vec3 C12 = texture2D(iChannel0, pixel + vec2( 0.0        , c_onePixel)).rgb;
    vec3 C22 = texture2D(iChannel0, pixel + vec2( c_onePixel , c_onePixel)).rgb;
    vec3 C32 = texture2D(iChannel0, pixel + vec2( c_twoPixels, c_onePixel)).rgb;    

    vec3 C03 = texture2D(iChannel0, pixel + vec2(-c_onePixel , c_twoPixels)).rgb;
    vec3 C13 = texture2D(iChannel0, pixel + vec2( 0.0        , c_twoPixels)).rgb;
    vec3 C23 = texture2D(iChannel0, pixel + vec2( c_onePixel , c_twoPixels)).rgb;
    vec3 C33 = texture2D(iChannel0, pixel + vec2( c_twoPixels, c_twoPixels)).rgb;    

    vec3 CP0X = CubicLagrange(C00, C10, C20, C30, frac.x);
    vec3 CP1X = CubicLagrange(C01, C11, C21, C31, frac.x);
    vec3 CP2X = CubicLagrange(C02, C12, C22, C32, frac.x);
    vec3 CP3X = CubicLagrange(C03, C13, C23, C33, frac.x);

    return CubicLagrange(CP0X, CP1X, CP2X, CP3X, frac.y);
}

```

**Tags:** texture, interpolation

> **Comments:**
> **username_1**: You could post the relevant shader code here in case of bitrot, no?  
> **username_0**: we should have some prettier code markup for shader code, i'll post on meta if someone hasn't beaten me to it!  
> **username_2**: Is that a specific shader language not available in the list of languages covered by our syntax highlighting?  
> **username_0**: I'm not sure.  It's just GLSL (from webgl to be exact!).  I just did 4 spaces before each line of code, not sure if there's a better way to mark it up...  

---

## 1 answers

### Answer by username_0 on 2015-08-04. Score: 8 | Accepted answer

It turns out that no, while you can use bicubic Lagrange interpolation for bicubic texture sampling, it isn't the highest quality option, and probably not actually likely to be used.

Cubic hermite splines are a better tool for the job.

Lagrange interpolation will make a curve that passes through the data points, thus preserving C0 continuity, but hermite splines preserve the derivatives at the edges while also passing through the data points, thus preserving C1 continuity and looking much better.

This question has some great info on cubic hermite splines:
<https://dsp.stackexchange.com/questions/18265/bicubic-interpolation>

Here is the cubic hermite version of the code I posted in the question:

```cplusplus
//=======================================================================================
vec3 CubicHermite (vec3 A, vec3 B, vec3 C, vec3 D, float t)
{
    float t2 = t*t;
    float t3 = t*t*t;
    vec3 a = -A/2.0 + (3.0*B)/2.0 - (3.0*C)/2.0 + D/2.0;
    vec3 b = A - (5.0*B)/2.0 + 2.0*C - D / 2.0;
    vec3 c = -A/2.0 + C/2.0;
    vec3 d = B;

    return a*t3 + b*t2 + c*t + d;
}

//=======================================================================================
vec3 BicubicHermiteTextureSample (vec2 P)
{
    vec2 pixel = P * c_textureSize + 0.5;

    vec2 frac = fract(pixel);
    pixel = floor(pixel) / c_textureSize - vec2(c_onePixel/2.0);

    vec3 C00 = texture2D(iChannel0, pixel + vec2(-c_onePixel ,-c_onePixel)).rgb;
    vec3 C10 = texture2D(iChannel0, pixel + vec2( 0.0        ,-c_onePixel)).rgb;
    vec3 C20 = texture2D(iChannel0, pixel + vec2( c_onePixel ,-c_onePixel)).rgb;
    vec3 C30 = texture2D(iChannel0, pixel + vec2( c_twoPixels,-c_onePixel)).rgb;

    vec3 C01 = texture2D(iChannel0, pixel + vec2(-c_onePixel , 0.0)).rgb;
    vec3 C11 = texture2D(iChannel0, pixel + vec2( 0.0        , 0.0)).rgb;
    vec3 C21 = texture2D(iChannel0, pixel + vec2( c_onePixel , 0.0)).rgb;
    vec3 C31 = texture2D(iChannel0, pixel + vec2( c_twoPixels, 0.0)).rgb;    

    vec3 C02 = texture2D(iChannel0, pixel + vec2(-c_onePixel , c_onePixel)).rgb;
    vec3 C12 = texture2D(iChannel0, pixel + vec2( 0.0        , c_onePixel)).rgb;
    vec3 C22 = texture2D(iChannel0, pixel + vec2( c_onePixel , c_onePixel)).rgb;
    vec3 C32 = texture2D(iChannel0, pixel + vec2( c_twoPixels, c_onePixel)).rgb;    

    vec3 C03 = texture2D(iChannel0, pixel + vec2(-c_onePixel , c_twoPixels)).rgb;
    vec3 C13 = texture2D(iChannel0, pixel + vec2( 0.0        , c_twoPixels)).rgb;
    vec3 C23 = texture2D(iChannel0, pixel + vec2( c_onePixel , c_twoPixels)).rgb;
    vec3 C33 = texture2D(iChannel0, pixel + vec2( c_twoPixels, c_twoPixels)).rgb;    

    vec3 CP0X = CubicHermite(C00, C10, C20, C30, frac.x);
    vec3 CP1X = CubicHermite(C01, C11, C21, C31, frac.x);
    vec3 CP2X = CubicHermite(C02, C12, C22, C32, frac.x);
    vec3 CP3X = CubicHermite(C03, C13, C23, C33, frac.x);

    return CubicHermite(CP0X, CP1X, CP2X, CP3X, frac.y);
}

```

Here is a picture showing the difference between sampling methods. From left to right: Nearest Neighbor, Bilinear, Lagrange Bicubic, Hermite Bicubic

[![enter image description here](https://i.sstatic.net/h0Ot0.png)](https://i.sstatic.net/h0Ot0.png)

> **Comments:**
> **username_3**: Although all cubic splines are, in a sense, equivalent, it's probably conceptually easier to use Catmull-Rom splines. e.g.  https://www.cs.cmu.edu/~462/projects/assn2/assn2/catmullRom.pdf  
> **username_0**: Do you think the tau parameter helps or hinders in this case?  I might be wrong but i feel like catmull rom is too "user defined" (and must be tuned), whereas the hermite spline attempts to just use information from the data that's there.  It seems like cubic hermite is closer to a "ground truth", which i guess would be something like an ideal sinc filter.  What do you think though?  
> **username_3**: I don't see how Catmull-Rom is "user defined". Once you have a sequence of 4 contiguous points, P[i-1], P[i], P[i+1], P[i+2] (4x4 for 2D case) the curve segment is defined between P[i] and P[i+1] and is C1 continuous with the neighbouring segments.   A sinc filter is fine for audio but not video. See Mitchell & Netravali: https://www.cs.utexas.edu/~fussell/courses/cs384g-fall2013/lectures/mitchell/Mitchell.pdf    IIRC Catmull-Rom is a special case of the family of filters they propose, but I think that filter is an approximating curve so, unlike C-R, might not go through the original points.  
> **username_0**: That's how he hermite spline works except that the catmull rom spline has an additional parameter tau (tension) that is user defined.  Also, sinc does apply to video, DSP is DSP :P  
> **username_3**: I must admit, I've never seen a tension parameter associated with Catmull Rom splines before, but then I've really only learned of them via Foley & van Dam (et al) or, say, Watt & Watt which, AFAICR, make no mention of such.  Actually, having said that, given there are four constraints - i.e. the curve has to pass through 2 points and have two defined tangents** at those points and it's a cubic - I'm at a bit of a loss as to how there's any more degrees of freedom to support a tension parameter ....  **Unless you mean the tangents can be scaled?  
> **username_0**: The tension parameter is shown in that first paper you linked.  And yeah I just mean, the fact that there's a tunable parameter means it must be tuned :p. There is a value of tau you can use such that you end up with a regular cubic hermite spline, and I think it might be "1" but not 100% sure on that.  
> **username_3**: Re that first paper: That serves me right for pasting in the first link I found that seemed suitable :-)    Again, I think I've only ever seen the tangent at a segment end point set to the (assuming I've done the matrix maths correctly) 1/2 the difference of the two surrounding neighbours.  
> **username_0**: Ah OK! Anyways yeah, catmull rom / hermite is really freaking awesome I totally agree. Really neat way to interpolate between data points and have C1 continuity at the edges. It really simplifies things (:  

---

