Talk:Bilinear interpolation

From Wikipedia, the free encyclopedia

Contents

[edit] Rewrite

I think that this article explains the topic in not the best way and I propose a major rewrite.

The idea of explaining the method using the "capital-H" (or capital-I) shape may help less mathematically skilled readers, but overcomplicates it, IMHO. After reading (and before some thinking) I was left with the impression that three separate linear operations were required, but the meaning of the intermediates was not explained.

I would describe the process as a weighted average of the four nearest points, with the weights being determined by the areas of the four rectangles divided by their sum. This sum is the area of the rectangle QQQQ.

Even if you don't like my alternative interpretation, it'd be easier to read if 1/( (x2-x1)(y2-y1) ) were in a sigle term at the start.

thoughts? TomViza 02:24, 19 January 2006 (UTC)

I like your idea of coupling the areas with the weights very much. It is easy to see it as a generalization of the one-dimensional case, where the weights are determined by the distances. I agree that this will be less complicated than the H-shape.
However, explaining bilinear interpolation as successive 1d linear interpolations also has its advantages. It explains the name (you do first an interpolation in the x-direction and then one in the y-direction) and it points to an easy way to do an error analysis.
So, what do you think of the following outline: First give the formula and your area-based explanation and then say that the formula can also be derived via the H-shape. -- Jitse Niesen (talk) 13:27, 19 January 2006 (UTC)
I think it is an excellent formulation for interpretation -- the weights of each corner are proportional to the area of the opposing sub-rectangle. The (x2-x1)(y2-y1) term is the area of the whole rectangle, and, for instance, the (x-x1)(y-y1) is the area of the rectangle opposite Q_22. That would hold to higher dimensions: In trilinear, the weights of each component are proportional to the size of the opposing volume. Back down to the univariate linear, case, the weight of each extreme is proportional to the distance opposite the extreme. Conceptually, it seems equivalent to the Barycentric_coordinates_(mathematics) interpolation often used in Finite_element_analysis.
Also, the 'four nearest points' idea works this way only if the points are on a regular grid. If the points are irregular, one violates the first assumption on the page. Drf5n 19:11, 24 September 2006 (UTC)

[edit] Error in matrix formulation

I think there is an error in the matrix formulation. Product of 2x2 and 1x2?

Yes, you're right. Thanks for the notice. The error is now fixed. -- Jitse Niesen (talk) 05:53, 20 April 2006 (UTC)


[edit] Error in Big equations

On using this page to check some interpolations we noticed that some minus signs should be plus signs. These have been changed. Otherwise this page is excellent. How about something similar in the tri-linear equations page?

Are you sure? The formula used to be:
f(x,y) \approx \frac{f(Q_{11})}{(x_1-x_2)(y_1-y_2)} (x-x_2)(y-y_2) - \frac{f(Q_{21})}{(x_1-x_2)(y_1-y_2)} (x-x_1)(y-y_2)
{} - \frac{f(Q_{12})}{(x_1-x_2)(y_1-y_2)} (x-x_2)(y-y_1) + \frac{f(Q_{22})}{(x_1-x_2)(y_1-y_2)} (x-x_1)(y-y_1).
Substituing (x,y) = (x1,y2), the right-hand side evaluates to
- \frac{f(Q_{12})}{(x_1-x_2)(y_1-y_2)} (x_1-x_2)(y_2-y_1) = f(Q_{12}).
So the minus sign seems correct to me.
It would be great if you could edit the trilinear interpolation page along similar lines. -- Jitse Niesen (talk) 04:09, 17 May 2006 (UTC)


(continued)

We did a test in a spread sheet - with all plus signs the answers were correct - minus signs would have given an incorrect answer. If the f(p) equation is plugged into the the f(R1) and f(R2) equations then minus signs do not seem to occur. We could be wrong though .. :-) I can provide a link to a page with a C code listing showing all plus signs - but I'm not sure if such links are allowed on these pages. I won't change the page. safetycritical 21:38, 17 May 2006 (UTC)

Of course, I could also be wrong. Please post the link, or email the code to me directly (you can find my email address by going to my user page and following the link to my home page). -- Jitse Niesen (talk) 03:43, 18 May 2006 (UTC)

[edit] Better form of equation

Wouldn't it be much more intuitive to replace

f(x,y) \approx f(0,0) \, (x-1)(y-1) - f(1,0) \, x(y-1) - f(0,1) \, (x-1)y + f(1,1) xy.

with

f(x,y) \approx f(0,0) \, (1-x)(1-y) + f(1,0) \, x(1-y) + f(0,1) \, (1-x)y + f(1,1) xy.
You're probably right. I modified the formulas to avoid negative quantities like x − 1. I hope I did not introduce errors in the process. -- Jitse Niesen (talk) 04:50, 22 May 2006 (UTC)
Maybe it's just me, but I find this way easier because in this case the (1 − x) and (1 − y) represent the weights on the pixels, and this way the alternate weights ( (1 − x) and x ) add up to 1.

[edit] Linearities

Hi there...

The text mentions the fact that "the interpolant is not linear", but have you ever seen a linear interpolant?... Linear in relation to what, after all??

I don't know why this phrase is bothering me... Perhaps because I also want to include a dicussion on the fact that this non-linear interpolant actually generates a linear filter that will be used to process the image (as you all know). Wat do you think?... -- NIC1138 14:18, 25 August 2006 (UTC)

The interpolant is the function f, which is a linear function. So "the interpolant is not linear" means nonlinear with respect to the arguments of f, which are x and y.
If the phrase is bothering you, then perhaps you could try reformulating it? By the way, I don't know that this non-linear interpolant actually generates a linear filter that will be used to process the image. Different editors come from different backgrounds, and I don't do image manipulation. Please do add something about this application. -- Jitse Niesen (talk) 14:54, 25 August 2006 (UTC)
Applying the interpolation is a linear transform of the input variables (Q_**), but the surface between the points is not a plane. Three points specify a plane, and the fourth would require some sort of curvature. If you were interpolating on a triangle by finding the point on the plane defined by the triangle, then you'd have a linear interpolant. Drf5n 01:13, 25 September 2006 (UTC)


[edit] Excel Bilinear Interpolation Function

It would seem useful to include software code for many of these excellent mathematical explanations/algorithms. I've included here an Excel user function I wrote yesterday that should work with nearly any 2D table in Excel. The algorithm is based on the algorithm on the main page. There has been very little testing so take this as a starting point.

Public Function Interp2D(ByVal x As Double, ByVal y As Double, ByVal rgTable As Range) As Double

   '**** Interp2D BILINEAR INTERPOLATION FUNCTION ****
   '(From http://en.wikipedia.org/wiki/Bilinear_interpolation)
   'In mathematics, bilinear interpolation is an extension of linear interpolation
   'for interpolating functions of two variables. The key idea is to perform linear
   'interpolation first in one direction, and then in the other direction.
   '
   'Suppose that we want to find the value of the unknown function f at
   'the point P = (x, y). It is assumed that we know the value of f at the four
   'points Q11 = (x1, y1), Q12 = (x1, y2), Q21 = (x2, y1), and Q22 = (x2, y2).
   '
   'We first do linear interpolation in the x-direction. This yields
   '
   ' f(R1) = (x2 - x)/(x2 - x1) * f(Q11) + (x - x1)/(x2 - x1)*f(Q21)
   '   where R1 = (x, y1)
   'and
   ' f(R2) = (x2 - x)/(x2 - x1) * f(Q12) + (x - x1)/(x2 - x1)*f(Q22)
   '   where R2 = (x, y2)
   '
   'We then proceed by interpolating in the y-direction.
   '
   ' f(P) = (y2 - y)/(y2-y1) * f(R1) + (y - y1)/(y2 - y1) * f(R2)
   '   where P = (x, y)
   '
   'Constraints:
   '   x is a Double and xHeaderRowMinValue <= x <= xHeaderRowMaxValue
   '   y is a Double and yHeaderColumnMinValue <= y <= yHeaderColumnMaxValue
   '   rgTable is an Excel range, specified in the normal Excel way
   '       The range should be a rectangle including the xHeader Row,
   '        the yHeaderRow, and all of the data table.
   '   xHeaderRow values increase to the right
   '   yHeaderColumn values increase downward
   
   Const minPositiveDouble As Double = 4.94065645841247E-324
   Const maxPositiveDouble As Double = 1.7976931348623E+308
   Const Epsilon As Double = 0.0001  'Tolerable error allowed
   
   Dim rgXHdr, rgYHdr, rgXYData As Range
   Dim xHdr(), yHdr(), xyData() As Double  'Array sizes still unknown
   Dim xHdrMin, xHdrMax, yHdrMin, yHdrMax As Double
   
   Dim x1, x2, y1, y2, ixX1, ixX2, ixY1, ixY2, fQ11, fQ12, fQ21, fQ22, fR1, fR2, fP As Double
   Dim row, col As Range
   Dim i, j As Integer
   Dim t As Variant
       
   'Get xHeaderRow, and set up array to hold values
   Set rgXHdr = rgTable.Offset(0, 1).Resize(1, rgTable.columns.Count - 1)
   ReDim xHdr(rgXHdr.columns.Count - 1)    'Size the array
   
   'Find min/max limits for xHeaderRow, and read values into array
   '  Assumes minimum xHeader value is 0, ie no negative numbers
   xHdrMin = maxPositiveDouble
   xHdrMax = 0
   i = 0
   For Each col In rgXHdr.columns
       If xHdrMax < col.Cells(1, 1).Value Then
           xHdrMax = col.Cells(1, 1).Value
       End If
       
       If xHdrMin > col.Cells(1, 1).Value Then
           xHdrMin = col.Cells(1, 1).Value
       End If
       
       xHdr(i) = col.Cells(1, 1).Value
       i = i + 1
   Next
   'Get yHeaderColumn, and set up array to hold values
   Set rgYHdr = rgTable.Offset(1, 0).Resize(rgTable.rows.Count - 1, 1)
   ReDim yHdr(rgYHdr.rows.Count - 1)    'Size the array
   
   'Find min/max limits for yHeaderColumn, and read values into array
   '  Assumes minimum yHeader value is 0, ie no negative numbers
   yHdrMin = maxPositiveDouble
   yHdrMax = 0
   i = 0
   For Each row In rgYHdr.rows
       If yHdrMax < row.Cells(1, 1).Value Then
           yHdrMax = row.Cells(1, 1).Value
       End If
       
       If yHdrMin > row.Cells(1, 1).Value Then
           yHdrMin = row.Cells(1, 1).Value
       End If
       
       yHdr(i) = row.Cells(1, 1).Value
       i = i + 1
   Next
   
   'Verify that x & y are in proper range,
   '  after allowing for Epsilon differences in values
   If x < xHdrMin And x >= xHdrMin - Epsilon Then
           x = xHdrMin
   End If
   If x > xHdrMax And x <= xHdrMax + Epsilon Then
           x = xHdrMax
   End If
   If y < yHdrMin And y >= yHdrMin - Epsilon Then
           y = yHdrMin
   End If
   If y > yHdrMax And y <= yHdrMax + Epsilon Then
           y = yHdrMax
   End If
   
   If x < xHdrMin Or x > xHdrMax Or y < yHdrMin Or y > yHdrMax Then
       MsgBox ("Interp2D: x or y value out of range: x= " & x & " y= " & y)
       End
   End If
       
   'Get xyData, and set up array to hold values
   Set rgXYData = rgTable.Offset(1, 1).Resize(rgTable.rows.Count - 1, rgTable.columns.Count - 1)
   ReDim xyData(rgXYData.columns.Count - 1, rgXYData.rows.Count - 1)    'Size the 2D array
   
   'Read Excel table values into xyData array
   i = 0
   For Each col In rgXYData.columns
       j = 0
       For Each row In col.rows
           xyData(i, j) = row.Cells(1, 1).Value
           j = j + 1
       Next
       i = i + 1
   Next
   
   'Find x1, x2, and corresponding indices ixX1, ixX2 in xHdr
   i = 0
   x1 = xHdrMin
   x2 = xHdrMin
   ixX1 = 0
   ixX2 = 0
   For Each t In xHdr
       If x = t Then
           x1 = t
           x2 = t
           ixX1 = i
           ixX2 = i
           Exit For
       ElseIf x > t Then
           x1 = t
           x2 = xHdr(i + 1)
           ixX1 = i
           ixX2 = i + 1
       End If
       i = i + 1
   Next
   
   'Find y1, y2, and corresponding indices ixY1, ixY2 in yHdr
   i = 0
   y1 = yHdrMin
   y2 = yHdrMin
   ixY1 = 0
   ixY2 = 0
   For Each t In yHdr
       If y = t Then
           y1 = t
           y2 = t
           ixY1 = i
           ixY2 = i
           Exit For
       ElseIf y > t Then
           y1 = t
           y2 = yHdr(i + 1)
           ixY1 = i
           ixY2 = i + 1
       End If
       i = i + 1
   Next
   
   'Get fQ11, fQ12, fQ21, fQ22
   fQ11 = xyData(ixX1, ixY1)
   fQ12 = xyData(ixX1, ixY2)
   fQ21 = xyData(ixX2, ixY1)
   fQ22 = xyData(ixX2, ixY2)
   
   'Get f(R1) and f(R2)
   If x1 = x2 Then
       fR1 = fQ11
       fR2 = fQ12
   Else
       fR1 = ((x2 - x) / (x2 - x1)) * fQ11 + ((x - x1) / (x2 - x1)) * fQ21
       fR2 = ((x2 - x) / (x2 - x1)) * fQ12 + ((x - x1) / (x2 - x1)) * fQ22
   End If
   
   'Get f(P)
   If y1 = y2 Then
       fP = fR1
   Else
       fP = ((y2 - y) / (y2 - y1)) * fR1 + ((y - y1) / (y2 - y1)) * fR2
   End If
       
   Interp2D = fP
   

End Function