Bilinear interpolation
Bilinear interpolation is a method of two-dimensional interpolation that estimates the value of a function at an arbitrary point within a rectangular grid by computing a weighted average based on the known values at the four surrounding grid points. It achieves this by performing linear interpolation successively along each axis—first in one direction (e.g., x) to generate intermediate values, then in the other direction (e.g., y)—resulting in a smooth approximation that is linear in both variables when the other is fixed.[1][2] The technique assumes the function is defined on a regular Cartesian grid and produces a continuous surface that passes exactly through the grid points, though it may introduce some smoothing or blurring in regions of high curvature due to its piecewise bilinear nature.[2] To compute the interpolated value at a point (x, y) where x_0 \leq x \leq x_1 and y_0 \leq y \leq y_1, the method first interpolates along the x-direction at the bottom and top grid lines to obtain values at (x, y_0) and (x, y_1), using weights proportional to the relative distances (x - x_0)/(x_1 - x_0) and (x_1 - x)/(x_1 - x_0); it then interpolates these results along the y-direction with similar weights (y - y_0)/(y_1 - y_0) and (y_1 - y)/(y_1 - y_0).[1] This process yields an explicit formula equivalent to f(x,y) = (1 - u)(1 - v) f(x_0, y_0) + u(1 - v) f(x_1, y_0) + (1 - u) v f(x_0, y_1) + u v f(x_1, y_1), where u = (x - x_0)/(x_1 - x_0) and v = (y - y_0)/(y_1 - y_0), ensuring the estimate lies within the convex hull of the corner values.[3] Bilinear interpolation is widely applied in fields requiring efficient 2D data resampling, such as image processing for scaling or rotating raster images, where it balances computational simplicity with reasonable quality by avoiding artifacts like aliasing in nearest-neighbor methods while being faster than higher-order alternatives like bicubic interpolation.[4] In computer graphics, it serves as a core component for texture mapping, blending colors from a 2D texture onto 3D surfaces by interpolating texel values during rasterization, often implemented in hardware for real-time rendering.[5] Additional uses include geographic information systems (GIS) for interpolating elevation or environmental data on grids and numerical simulations where quick approximations of scattered data are needed without excessive computational overhead.[1]Fundamentals
Definition and Motivation
Linear interpolation in one dimension serves as a foundational concept, estimating the value of a function between two known points through a weighted average that reflects their relative distances.[6] Bilinear interpolation extends this approach to two dimensions, providing an estimate for the value of a function f at an arbitrary point (x, y) inside a rectangular grid by incorporating the known function values at the four nearest grid points: (x_1, y_1), (x_1, y_2), (x_2, y_1), and (x_2, y_2).[7] This technique constructs a smooth surface patch over the rectangle formed by these points, ensuring the approximation is linear along both coordinate axes.[8] The motivation for bilinear interpolation stems from the limitations of univariate methods when dealing with two-dimensional data, such as scattered measurements on a plane in scientific simulations or terrain models. In fields like computer graphics and numerical analysis, it enables efficient approximation of continuous surfaces from discrete grid data, facilitating tasks like rendering or data visualization without requiring complex higher-order computations.[9] It is fundamentally based on repeated linear interpolation along each dimension.[8]Comparison to Univariate Interpolation
Bilinear interpolation extends the principles of univariate linear interpolation from one dimension to two dimensions, adapting the method to handle data arranged on a rectangular grid. In univariate linear interpolation, a value is estimated at a point along a straight line segment connecting exactly two known data points, resulting in a precise linear approximation that is exact at the endpoints and forms a straight line between them. By contrast, bilinear interpolation employs four known values at the corners of a rectangle surrounding the target point, performing successive one-dimensional linear interpolations first along one axis (e.g., x) to generate intermediate values, and then along the other axis (e.g., y) to obtain the final estimate. This separable structure allows bilinear interpolation to approximate a function over a 2D domain while leveraging the simplicity of univariate operations, but it requires the data to be structured in a grid-aligned rectangle for accurate application.[1][2] One key advantage of bilinear interpolation over simpler methods like nearest-neighbor selection lies in its ability to produce smoother transitions across the interpolated region. Nearest-neighbor interpolation merely assigns the value of the closest grid point, often resulting in abrupt, blocky artifacts that preserve no continuity between adjacent areas. Bilinear interpolation, building on univariate linear techniques, generates gradual value changes by weighting contributions from all four surrounding points, yielding a more continuous and visually coherent output suitable for applications involving spatial data.[10][11] Despite these benefits, bilinear interpolation has limitations when compared to univariate linear interpolation, particularly in handling non-rectangular or irregularly shaped data domains. Univariate linear interpolation remains exact along any straight line connecting its two points, faithfully reproducing linear variations without distortion. In two dimensions, however, bilinear interpolation can introduce shearing or warping effects, as the separable process assumes axis-aligned independence that may not hold for skewed or curved underlying distributions, leading to approximations that deviate from the true surface. Geometrically, this manifests as the interpolated values lying on a hyperbolic paraboloid—a saddle-shaped surface patch ruled by straight lines in both directions—rather than the flat plane or straight line of univariate interpolation, which can exaggerate distortions in anisotropic data.[12][13][14]Computation
Repeated Linear Interpolation
Bilinear interpolation can be understood and implemented as a two-stage process of repeated linear interpolation, first along one coordinate axis and then along the other, to estimate the value of a function f(x, y) at an arbitrary point (x, y) within a rectangular cell bounded by known grid points (x_1, y_1), (x_2, y_1), (x_1, y_2), and (x_2, y_2). This method leverages univariate linear interpolation, which computes a weighted average between two points based on their relative distance.[2][3] In the first stage, linear interpolation is applied horizontally along the rows at fixed y = y_1 and y = y_2 to obtain intermediate values at the desired x-coordinate. The intermediate value at (x, y_1) is given by f(x, y_1) = (1 - t) f(x_1, y_1) + t f(x_2, y_1), where t = \frac{x - x_1}{x_2 - x_1} is the normalized distance along the x-direction, assuming x_1 < x < x_2. Similarly, the intermediate value at (x, y_2) is f(x, y_2) = (1 - t) f(x_1, y_2) + t f(x_2, y_2). This step creates two linear approximations spanning the x-interval at the bottom and top edges of the cell.[1][15] In the second stage, linear interpolation is performed vertically along the column at the fixed x using the two intermediate values to yield the final estimate: f(x, y) = (1 - u) f(x, y_1) + u f(x, y_2), where u = \frac{y - y_1}{y_2 - y_1} is the normalized distance along the y-direction, assuming y_1 < y < y_2. The order of horizontal and vertical stages can be reversed without affecting the result, as the operations commute in this linear setting.[16][17] Intuitively, this repeated process constructs a bilinear surface by first forming horizontal linear "slices" across the cell and then blending those slices linearly in the vertical direction, effectively creating a ruled surface that varies linearly in both parameters.[2][3] For implementation, the algorithm can be expressed in pseudocode as follows:This straightforward procedure requires only four known values and a few arithmetic operations, making it computationally efficient for grid-based resampling.[15]function bilinear_interpolate(x, y, x1, y1, x2, y2, f11, f21, f12, f22): if x1 == x2 or y1 == y2: return error // degenerate case t = (x - x1) / (x2 - x1) u = (y - y1) / (y2 - y1) // Clamp t and u to [0,1] if outside cell t = max(0, min(1, t)) u = max(0, min(1, u)) // Horizontal interpolation f_y1 = (1 - t) * f11 + t * f21 // f11 = f(x1,y1), f21 = f(x2,y1) f_y2 = (1 - t) * f12 + t * f22 // f12 = f(x1,y2), f22 = f(x2,y2) // Vertical interpolation return (1 - u) * f_y1 + u * f_y2function bilinear_interpolate(x, y, x1, y1, x2, y2, f11, f21, f12, f22): if x1 == x2 or y1 == y2: return error // degenerate case t = (x - x1) / (x2 - x1) u = (y - y1) / (y2 - y1) // Clamp t and u to [0,1] if outside cell t = max(0, min(1, t)) u = max(0, min(1, u)) // Horizontal interpolation f_y1 = (1 - t) * f11 + t * f21 // f11 = f(x1,y1), f21 = f(x2,y1) f_y2 = (1 - t) * f12 + t * f22 // f12 = f(x1,y2), f22 = f(x2,y2) // Vertical interpolation return (1 - u) * f_y1 + u * f_y2
Polynomial and Matrix Forms
Bilinear interpolation can be formulated as a polynomial of the form f(x, y) = a + b(x - x_1) + c(y - y_1) + d(x - x_1)(y - y_1), where (x_1, y_1), (x_2, y_1), (x_1, y_2), and (x_2, y_2) are the coordinates of the four surrounding grid points, and the coefficients a, b, c, and d are determined by substituting the known function values f_{11} = f(x_1, y_1), f_{21} = f(x_2, y_1), f_{12} = f(x_1, y_2), and f_{22} = f(x_2, y_2) into the equation and solving the resulting linear system.[18] This polynomial form arises naturally from the requirement that f(x, y) is linear in x for any fixed y and linear in y for any fixed x. Starting with linear interpolation in the x-direction at the two y-levels yields expressions that are affine in x with coefficients depending linearly on y; substituting and expanding produces the cross term d(x - x_1)(y - y_1), confirming the bilinear structure.[8] The coefficients are solved explicitly as follows: a = f_{11}, b = (f_{21} - f_{11}) / (x_2 - x_1), c = (f_{12} - f_{11}) / (y_2 - y_1), and d = \frac{f_{22} - f_{21} - f_{12} + f_{11}}{(x_2 - x_1)(y_2 - y_1)}.[19] In matrix notation, the polynomial can be expressed using the basis vector for the global coordinates as f(x, y) = [1, x, y, xy] M \begin{bmatrix} f_{11} \\ f_{12} \\ f_{21} \\ f_{22} \end{bmatrix}, where M is the 4×4 transformation matrix depending on x_1, x_2, y_1, y_2 obtained by inverting the system mapping the basis to the grid point evaluations. This matrix form is equivalent to the repeated linear interpolation approach, as expanding the weighted average (1-t)(1-u)f_{11} + t(1-u)f_{21} + (1-t)uf_{12} + tuf_{22} (with t = (x - x_1)/(x_2 - x_1) and u = (y - y_1)/(y_2 - y_1)) and collecting terms yields the same polynomial expression after algebraic simplification.[8]Unit Square Parameterization
In bilinear interpolation over a general rectangular domain defined by corners at coordinates (x_1, y_1), (x_2, y_1), (x_2, y_2), and (x_1, y_2), the process begins by normalizing the evaluation point (x, y) to parameters s and t within the unit square [0,1] \times [0,1]. This transformation is given bys = \frac{x - x_1}{x_2 - x_1}, \quad t = \frac{y - y_1}{y_2 - y_1}.
Such normalization maps the arbitrary rectangle to a canonical unit square, facilitating standardized computations independent of the specific grid spacing.[20] With the point normalized to (s, t), the interpolated value f(s, t) is computed as a weighted average of the function values at the four corners of the unit square, denoted f_{11} = f(0,0), f_{12} = f(1,0), f_{21} = f(0,1), and f_{22} = f(1,1):
f(s,t) = (1-s)(1-t) f_{11} + s(1-t) f_{12} + (1-s) t f_{21} + s t f_{22}.
This expression arises from applying linear interpolation first along one parameter (e.g., s) and then the other (e.g., t), yielding a bilinear surface.[21] The formula represents a convex combination of the corner values, where the weights (1-s)(1-t), s(1-t), (1-s)t, and st sum to 1 and are nonnegative for s, t \in [0,1], ensuring the result lies within the convex hull of the corner values. This interpretation underscores the method's role as a direct weighted mean, preserving bounds and providing a smooth transition across the domain.[21] The unit square parameterization offers key advantages, including a reduction in the number of parameters by decoupling the interpolation weights from the physical coordinates, which simplifies implementation and avoids repetitive solving of systems for each evaluation. It is particularly prevalent in computer graphics, where normalized coordinates s and t (often called texture coordinates) are standard in APIs like OpenGL for efficient texture mapping and sampling via bilinear interpolation.[20][22]
Properties
Mathematical Characteristics
Bilinear interpolation exhibits linearity in each variable separately, meaning that fixing one coordinate results in a linear function with respect to the other. Overall, the method produces an affine function across the rectangular domain, exactly preserving linear functions evaluated along the grid lines. This property ensures that if the input values at grid points follow an affine relationship, such as f(m, n) = a m + b n + c, the interpolated output matches the exact affine surface u(x, y) = a x + b y + c. The interpolation is C^0 continuous across adjacent rectangular patches, guaranteeing that the function values are seamless and continuous everywhere within the domain. However, it lacks differentiability at the boundaries between patches, where gradients exhibit discontinuities due to the piecewise nature of the construction. This results in smooth value transitions but abrupt changes in slope at grid edges, limiting higher-order smoothness.[23] As a convex combination of the four corner values, bilinear interpolation weights the inputs with non-negative coefficients that sum to one, specifically using normalized distances from the point to the corners. Consequently, the output value always lies within the range bounded by the minimum and maximum of the corner values, preserving convexity and avoiding overshoots beyond the local data extrema.[24] For any four given values at the corners of a rectangle, there exists a unique bilinear function that interpolates these points exactly. This uniqueness stems from the bilinear form being the lowest-degree polynomial (degree at most one in each variable) capable of satisfying the interpolation conditions without additional constraints.[25]Error and Approximation Bounds
Bilinear interpolation provides a second-order approximation for smooth functions, with the error scaling as O(h^2), where h denotes the uniform grid spacing, due to the truncation of higher-order terms in the Taylor expansion beyond the linear components in each direction.[26] This order of accuracy arises because the method exactly reproduces linear functions in x and y, but deviates for quadratic and higher terms, limiting its precision for functions with significant curvature.[27] The maximum approximation error can be bounded using Taylor expansions around the grid points, with the bound depending on the second partial derivatives of f over the local cell. This highlights the dependence on second partial derivatives, emphasizing that error grows with the function's local curvature or non-linearity. In worst-case scenarios, such as regions of high curvature or anisotropic data, bilinear interpolation can introduce noticeable distortions; for instance, circular features may appear as pillow-shaped artifacts due to the method's separable linear averaging, which unevenly smooths edges and corners.[28] Compared to bicubic interpolation, bilinear offers faster computation but reduced accuracy for smooth data, as bicubic incorporates higher-order terms to achieve O(h^4) error, resulting in lower root-mean-square errors (e.g., up to 12.6% improvement in reconstructed image quality).[29] This trade-off makes bilinear suitable for real-time applications where speed outweighs precision demands.[30]Applications
Image Resizing and Processing
Bilinear interpolation plays a central role in image resizing by enabling the scaling of digital images either upward (upsampling) or downward (downsampling) while preserving visual continuity. In upsampling, it mitigates the jagged, blocky effects associated with pixel duplication in nearest-neighbor methods by computing new pixel values as weighted averages of the four surrounding input pixels, resulting in smoother gradients and edges. This approach is particularly valuable for applications requiring natural-looking enlargements without introducing harsh discontinuities. For downsampling, bilinear interpolation facilitates anti-aliasing through the averaging of neighboring pixel intensities, which attenuates high-frequency details that might otherwise produce aliasing artifacts such as moiré patterns or flickering edges.[31] By blending contributions from adjacent pixels, it ensures that reduced-resolution images retain a coherent representation of the original content, though at the potential cost of subtle detail loss. In practical implementations, such as the OpenCV library'scv2.resize function, bilinear interpolation serves as the default method (via the INTER_LINEAR flag) and excels at smoothing transitions during zoom operations, making it suitable for real-time video processing or interactive scaling. Similarly, Adobe Photoshop employs bilinear interpolation as one of its resampling options in the Image Size dialog, where it provides efficient edge smoothing for general-purpose resizing tasks.[32]
A key advantage of bilinear interpolation lies in its computational efficiency, demanding only three multiplications and six additions/subtractions per output pixel, which supports its widespread adoption in raster graphics pipelines for low-latency operations.[14] However, it can introduce minor blurring as an artifact, softening fine details during upsampling; this trade-off is generally beneficial compared to the pronounced aliasing seen in nearest-neighbor interpolation.[33]
To illustrate, the following pseudocode demonstrates bilinear interpolation for a simple 2x upsampling resize, where input coordinates are normalized to the unit square around each set of four source pixels:
This normalization to the unit square streamlines boundary handling and computation.for each output pixel position (i, j) in the enlarged image: src_x = (i * 0.5) # Scale factor for 2x resize src_y = (j * 0.5) x1 = floor(src_x) y1 = floor(src_y) x2 = min(x1 + 1, input_width - 1) y2 = min(y1 + 1, input_height - 1) dx = src_x - x1 dy = src_y - y1 # Weighted average output[i, j] = (1 - dx) * (1 - dy) * input[x1, y1] + dx * (1 - dy) * input[x2, y1] + (1 - dx) * dy * input[x1, y2] + dx * dy * input[x2, y2]for each output pixel position (i, j) in the enlarged image: src_x = (i * 0.5) # Scale factor for 2x resize src_y = (j * 0.5) x1 = floor(src_x) y1 = floor(src_y) x2 = min(x1 + 1, input_width - 1) y2 = min(y1 + 1, input_height - 1) dx = src_x - x1 dy = src_y - y1 # Weighted average output[i, j] = (1 - dx) * (1 - dy) * input[x1, y1] + dx * (1 - dy) * input[x2, y1] + (1 - dx) * dy * input[x1, y2] + dx * dy * input[x2, y2]