Download Interpolation of Continuous Signals: Voronoi Diagrams, Delaunay Triangulation, and Methods and more Lecture notes Computer Graphics in PDF only on Docsity!
Interpolation and Filtering
Data is often discretized in time and/or space. We have only a finite number of sample points, i.e., the continuous signal is only known at few points (data points). But in general data is also needed between this points. By interpolation we obtain a representation that matches the function at the data points. An evaluation at any other point is possible. We can reconstruct the signal at points that are not sampled. For this reconstruction some assumptions are necessary. Often we arrogate smooth functions.
1 Voronoi Diagrams and Delaunay
Triangulation
Given are irregularly distributed points without connectivity information. The problem is to obtain connectivity to find a "good" triangulation. For a set of points there are many possible triangulations. A measure for the quality of a triangulation is the aspect ratio of the so-defined triangles. Long and thin triangles are to be avoided. Scattered data triangulation for 2D: A triangulation of a set of data points S =s 0 , s 1 , , sm ∈ℝ 2 consists of
- Vertices (0D) = S
- Edges (1D) connecting two vertices
- Faces (2D) connecting three vertices A triangulation must satisfy the following criteria:
- ∪ faces=conv S , i.e. the union of all faces including the boundary is the convex hull of all vertices.
- The intersection of two triangles is either empty, or a common vertex, or a common edge, or a common face (tetrahedra). The following triangulations are not valid. Non valid Triangulations The image shows three different triangulations which are not valid. The triangulation a) is not valid because it “contains a hole”. In b) two faces overlap. The triangulation c) contains two vertices which form a “T” (T-vertices). To get a triangulation from scattered data the Delaunay Triangulation, which is tightly connected to the Voronoi diagram, can be used.
1.1 The Voronoi diagram
Given is a set of points X^ ={x 1 ,^ ,^ xn }^ from ℝ d and a distance function dist(x,y).
The Voronoi diagram Vor(X) contains for each point xi a cell V^ ^ xi ^ with V xi ={x∣dist x , xi dist x , x (^) j ∀ j≠i } For each sample every point within a Voronoi region is closer to it than to every other sample. A voronoi diagram The picture shows an example for a Voronoi diagram. For a point xi the corresponding Voronoi cell is given by the intersection of all half spaces h^ xi ,^ x^ j ^ ,^ j≠i^ : V xi =∩ (^) j≠i h xi , x (^) j h xi , x (^) j (^) is seperated by the perpendicular bisector between xi and x (^) j. h xi , x (^) j contains xi. Voronoi cells are convex. Construction of voronoi cells The image illustrates the construction of a half space and the corresponding Voronoi cell.
1.1.1 The Delaunay triangulation
The Delaunay graph Del(X) is the geometric dual of the Voronoi diagram Vor(X). The points in X are nodes. Two nodes xi and x^ j are connected if the Voronoi cells V^ ^ xi and V^ ^ X^ j ^ share the same edge. Delaunay cells are convex. The Delaunay triangulation is the triangulation of the Delaunay graph. A Delaunay graph The picture shows an example for a Delaunay graph. The Delaunay triangulation in 2D: Three points xi ,^ x^ j ,^ xk in X belong to a face from Del(X) if no further point lies inside the circle around xi ,^ x^ j ,^ xk. Two points xi ,^ x^ j form an edge if there is a circle around xi ,^ x^ j that does not contain a third point from X. For each triangle the circumcircle does not contain any other sample. The smallest angle and the ratio of radius of incircle radius of circumcircle is maximized. The triangulation is unique
walk clockwise along convex polygon of triangulation until the tangent points, inserting edges between pi and polygon points update convex hull endfor
- Bowyer-Watson algorithm The Bowyer-Watson algorithm builds the Delaunay triangulation from scattered points in one pass. [Watson-1981-CDD] [Bowyer-1981-CDT] The idea of this algorithm is the incremental insertion of points into the triangulation:
- Start with initial triangulation which covers the domain (e.g. two triangles of bounding box)
- Incremental insertion of points into the triangulation
- All triangles whose circumcircles contain the inserted point are removed
- The resulting cavity is triangulated by linking the inserted point to all vertices of the cavity boundary
- The cavity is star-shaped: Edges from the location of the newly inserted point The Bowyer Watson algorithm The picture illustrates one step of the Bowyer Watson algorithm. The algorithm: determine the super triangle that encompasses all vertices add super triangle vertices to the end of the vertex list add the super triangle to the triangle list for (each point in the vertex list) calculate the triangle circumcircle center and radius insert new point if (new point lies in a circumcircle) add the three triangle edges to the edge buffer remove the triangle from the triangle list delete multiple specified edges from the edge buffer, which leaves the edges of the enclosing polygon add all triangles formed of the point and the enclosing polygon remove all triangles from the triangulation that use the super triangle vertices and remove their vertices from the vertex list The following Applet visualizes the Voronoi diagram and the Delaunay triangulation: DelaunayApplet/DelaunayApplet.html Other techniques, e.g. Radial sweep, Intersecting halfspaces, Divide and conquer (merge- based or split-based), exist but are not discussed here.
1.2 Univariate Interpolation
Univariate interpolation means the interpolation for one variable.
1.2.1 Taylor Interpolation
For the taylor interpolation we use the basis functions mi=^ x i with i ∈ℕ 0 (monom basis) P m ={1, x , x 2 , , x m } is an^ m+1 -dimensional vector space of all polynomials with maximum degree m.
The task is to find coefficients ci^ with^ f^ =∑
i ci⋅x i
. This is the general approach, for xi we can use any other basis function, but if we use the monom basis the interpolation problem can be solved with the Vandermond matrix. The samples are represented by: f^ ^ x^ j =^ f^ j ∀^ j=^1 n The interpolation problem is given by: V ⋅ c = f with the Vandermond matrix V (^) ij = xi j− 1 1 x 1 x 1 2 x 1 n− 1 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ 1 xn xn 2 xn n− 1 ⋅ c 1 c 2 ⋮ cn = f (^1) f (^2) ⋮ f (^) n Properties of the Taylor interpolation:
- Unique solution
- Numerical problems / inaccuracies
- Complete system has to be solved again if a single value is changed
- Not intuitive Example: Given are 3 samples: f(1)=2, f(2)=5, f(4)= This leads to the following linear system of equations. 1 1 1 1 2 4 1 4 16 ∣ 2 5 3 ⇒ 1 1 1 0 1 3 0 0 6 ∣^ 2 3 − 8 ⇒ c 1 =− 11 3 , c 2 = 7 , c 3 =− 4 3 The interpolated function: f^ ^ x^ =−^ 4 3 x 2 7 x− 11 3
1.2.2 Generic interpolation problem
Given are n sampled points X ={xi }⊆⊆ℝ d with function values f^ i. The n -dimensional function space n d has the basis {∮i= 1 n } We search coefficients ci^ with^ f^ =∑ i ci⋅∮i x The samples are represented by: f^ ^ x^ j =^ f^ j ∀^ j=^1 n
We have to solve the linear system of equations M⋅ c =^ f^ with^ M^ ji=∮i ^ x^ j ^ ,^ ci =ci ,^ f^ j =^ f^ j
Note: The number of points n determines dimension of vector space (= degree of polynomials)
B-Spline basis functions. The picture shows B-Spline basis functions.
1.2.5 Piecewise linear interpolation
This is the most simple approach (except for nearest-neighbor sampling). The main advantage is that it is fast to compute. The piecewise linear interpolation is often used in visualization applications. Given are data points ^ x 0 ,^ y 0 ^ ,^ ,^ xn ,^ yn Piecewise linear interpolation The image illustrates the computation of the value f(x) for a point x between two points (^) xi , xi 1. For any point x with xi^ x^ xi 1 evaluate f^ ^ x^ =^1 −u^ yiuyi 1 where u= x− xi xi 1 − xi ∈ [0,1].
1.3 Differentiation on grids
First Approach Idea: Replace differential by “finite differences”. Note that approximating the derivate by f ' x = df dx f x causes subtractive cancellation and large rounding errors for small h. f ' x ≈ f xh− f x h Second Approach Approximate / interpolate (locally) by differentiable function and differentiate this function.
1.3.1 Finite differences on uniform grids with grid size h
(1D case)
Finite differences on uniform grids with grid size h The picture shows data values on a uniform 1D-grid with size h
- Forward differences: f ' xi = f xi 1 − f xi h
- Backward differences: f ' xi = f xi − f xi− 1 h
- Central differences: f ' xi = f xi 1 − f xi− 1 2 h Error estimation:
- Forward / backward differences are first order.
- Central differences are second order.
1.3.2 Finite differences on non-uniform rectilinear grids
Forward and backward differences as for uniform grids with xi 1 − xi= h xi− xi− 1 = h Finite differences on non-uniform grids The illustration shows the differences (^) h and (^) h. The central differences are given by the Taylor expansion around the point xi. f xi 1 = f xi hf ' xi h 2 2 f ' ' xi f xi− 1 = f xi − hf ' xi h 2 2 f ' ' xi
ij x , y= (^) j x ⋅l y
The tensor product is: f x , y= ∑
i=1, j= 1 n , m ij x , ycij This means that we gain one bivariate basis function with two variables out of two univariate basis functions with one variable. The task is to solve a linear system of equations for the unknown coefficients cij. An extension to k dimension is done in the same way. Example: Given are 4 samples: f(0,0)=2, f(0,1)=0.5, f(1,0)=3, f(1,1)= We choose the monom basis for (^) j and l : 1 =1, 2 = x , l=1, 2 =x ⇒ 11 = 1 , 12 = x , 21 = y , 22 = xy This leads to the following system of equations: 2 =c 11 1 2 =c 11 c 21 3 =c 11 c 12 1 =c 11 c 12 c 21 c 22 ⇒ c 11 = 2 , c 12 = 1 , c 21 =− 3 2 , c 22 =− 1 2 The interpolated function is given by: f x , y=− 1 2 xy x− 3 2 y 2
1.4.1 Bilinear interpolation on a rectangle
- Tensor product for two linear interpolations
- 2D local interpolation in a cell
- Known solution of the linear system of equations for the coefficients cij
- Four data points ^ xi ,^ y^ j ^ ,^ ,^ xi 1 ,^ y^ j 1 ^ with scalar values f^ kl =^ f^ ^ xk ,^ yl
- Bilinear interpolation of points (x,y) with xi ^ x^ xi 1 and yi^ y^ yi 1 Bilinear interpolation on a rectangle The picture shows 4 grid points, the corresponding data values and an interpolated data value between these points. f x , y= 1 −[ 1 − f (^) i , j f (^) i1, j ][ 1 − f (^) i , j 1 f (^) i1, j 1 ] = 1 − f (^) j f (^) j 1 with:
f (^) j= 1 − f (^) i , j f (^) i1, j f (^) j 1 = 1 − f (^) i.j 1 f (^) i1, j 1 and local coordinates: = x− xi xi 1 − xi , = y− yi yi 1 − yi , , ∈[0,1] Bilinear interpolation on a rectangle The image shows 4 data values which form a rectangle and a interpolated data value with its local coordinates α and β. The equation above for f(x,y) can be rewritten as follows: f x , y= 1 − 1 − f (^) i , j 1 − f (^) i1, j 1 − f (^) i , j 1 f (^) i1, j 1 The point to be interpolated divides the rectangle into 4 areas. The corresponding data value can be computed by weighting the 4 given data values with area of the corresponding opposite area. Bilinear interpolation on a rectangle The picture shows how the point to be interpolated divides the rectangle into 4 areas. Note: Bilinear interpolation is not linear. This is illustrated in the animation BilinearInterpolation.divx. You see a plane which is generated by interpolating between 4 points. In the animation you see how the shape of the plane changes when the points are moved. Maya-Animation: BilinearInterpolation.divx Trilinear interpolation on a 3D uniform grid:
- Straightforward extension of bilinear interpolation
- Three local coordinates , ,
- Known solution of the linear system of equations for the coefficients cij
- Trilinear interpolation is not linear! Extension to higher order of continuity:
- Piecewise cubic interpolation in 1D
- Piecewise bicubic interpolation in 2D
- Piecewise tricubic interpolation in 3D
- Based on Hermite polynomials
Barycentric coordinates from area/volume considerations for d=2 and d= The image shows examples for barycentric coordinates from area/volume considerations for d=2 and d=3. Barycentric interpolation in a triangle: Geometrically, barycentric coordinates are given by the ratios of the area of the whole triangle and the subtriangles defined by x and any two points of x 1 ,^ x 2 ,^ x 3. Vol x 1 , x 2 , x 3 =det x 1 x 2 x 3 y 1 y 2 y 3 1 1 1 =± 2 Area x 1 , x 2 , x 3 1 = Vol x , x 2 , x 3 Vol x 1 , x 2 , x 3 So we end up with x =∑ i
i xi and ∑
i i= (^1). Barycentric interpolation in a triangle The image illustrates the barycentric interpolation in a triangle.
1.4.2.3 Interpolation in a generic quadrilateral
The main application for this approach are curvilinear grids. The problem is to find a parameterization for arbitrary quadrilaterals. Interpolation in a generic quadrilateral The picture illustrates the problem of finding a parameterization for arbitrary quadrilaterals. The mapping from rectangular domain to quadratic domains is known: Bilinear interpolation on a rectangle. x 12 = 1 ⋅x 1 1 − 1 ⋅x 2 x 34 = 1 ⋅x 4 1 − 1 ⋅x 3 x= 2 ⋅x 12 1 − 2 ⋅x 34
1 ∈[0,1] , 2 ∈[0,1] Computing the inverse of is more complicated. There are two possibilties:
- Analytically solve quadratic system for 1 , 2
- Numerical solution by Newton iteration The final value is given by: f^ = 2 ⋅ 1 ⋅f^1 ^1 − 1 ⋅f^2 ^1 − 2 ⋅ 1 f^4 ^1 − 1 ^ f^3 The Jacobi matrix J is given by: J ij = ∂ i ∂ (^) j J ij describes direction and speed of position changes of when (^) j are varied Newton iteration: start with seed points as start configuration, e.g. (^) i= 1 2 while ∥x− 1 , 2 , 3 ∥ compute (^) J 1 , 2 , 3 transform (^) x in coordinate system (^) J : x= J 1 , 2 , 3 − 1 ⋅ x− 1 , 2 , 3 update i=i x , i Other primitive cell types are possible. Examples: Other primitive cells The illustration shows a prism where the interpolation is done twice barycentric and once linear and a pyramid where the interpolation is done first bilinear on base face and then linear.
1.4.2.4 Inverse distance weighting
The Shepard interpolation [D. Shepard, A two-dimensional interpolating function for irregularly spaced data. Proc. ACM. nat. Conf., 517--524, 1968] was originally developed for scattered data. Interpolated values: f^ ^ x^ =∑ i i x f (^) i The sample points are the vertices of the cell. Basis functions: i ^ x^ =^
∥x−^ xi∥
− p
∑∥x−^ x^ j∥
− p Define values at sample points: f^ ^ xi :=^ f^ i=lim^ x xi^ f^ ^ x
- Projection into subspaces
- Often a mapping to a sub set of the original values is chosen Selection :
- Selection of data according to logical conditions (predicates)
- Example:
- Height field 2 d^^1 v^ with data (x,y,h)
- (^) Dσ ={ x , y , h∣ x^2 y^2 5 km∧h 1 km} Slicing:
- Example: 2D cutting surface (slice) through a 3D volume Slicing The picture shows a cube with a cutting surface.
1.7 Fourier Transform
The fourier transformation is often used for image processing, especially for filtering. Here, the image is converted from the spatial domain to the frequency domain by using the fourier transformation. In the frequency domain the image is multiplied by a filter (e.g. Gauss-filter, box filter, etc.). Afterwards the image will be transformated back to the spatial domain. In the spatial domain a signal ht is given by the amplitude value as a function of time. The analogous representation H v in the frequency domain is a function of the frequency v^. Via the fourier transform these two representations can be converted into each other.
- forward transform: H v=∫ −∞ ∞ ht e − 2 ivt dt
- inverse transform: ht =∫ −∞ ∞ H ve 2 ivt dv
The convolution is defined by: g∗ht =∫
−∞ ∞ g ⋅h t −d The convolution theorem : g∗ht ⇔G v⋅H v ,i.e. a convolution in the time domain corresponds to a multiplication in the frequency domain. Examples:
Examples for functions in the time domain and the corresponding functions in the frequency domain The above images show different functions with their representations in the frequency domain. In applications mostly discrete Fourier transform, which are based in a discrete signal, are used.
1.8 Sampled Signals
Assume that a signal ht is band limited with frequencies smaller than B. The so called Nyquist frequency is defined as vNyq=^2 B The signal can be discretized with a constant step size ^ t=^ 1 vNyq = 1 2 B Sampled signal: h^ j=h^ j⋅^ t^ If only a finite interval j=0..n-1 is used, periodicity is assumed. Sampling Theorem (Shannon 1949): If H f = 0 for all (^) ∣v∣B= vNyq 2 , then ht is uniquely given by the samples hi :
A low pass filter The image shows a low pass filter with a limit frequency in time and frequency domain A high pass filter emphasizes features, e.g. edges. High pass filter The picture shows a high pass filter in time and frequency domain Reconstruction issues : The measurements mt^ ^ of the original signal s^ t^ ^ are based on a point-spread function pt−ti ^ , not on the ideal delta function t−ti . A convolution in coordinate space corresponds to a multiplication in the frequency space: m t =∫ −∞ ∞ p t− s d ⇔ M v= P v S v An Applet for Fourier Analysis can be found at: http://www.gris.uni- tuebingen.de/projects/grdev/applets/fourier/html/index.html
