二维插值
二维插值 原理
二维插值原理
二维插值是一种基于已知数据点的二维曲线或曲面估计方法。
它广泛应用于图像处理、地理信息系统、物理模拟等领域。
在二维插值中，我们假设已知的数据点位于一个二维平面上，每个数据点都有一个对应的数值。
我们的目标是通过这些已知数据点，来推断出未知位置上的数值。
常见的二维插值方法包括线性插值、拉格朗日插值和样条插值等。
线性插值是最简单的二维插值方法之一。
它假设在两个相邻数据点之间，数值的变化是线性的。
我们可以通过这两个相邻数据点之间的线段来估计未知位置上的数值。
拉格朗日插值则使用一个多项式来拟合已知数据点。
该多项式会经过所有已知数据点，并通过它们来估计未知位置上的数值。
它的优点是能够完全通过已知数据点来插值，但在高维情况下容易产生过拟合问题。
样条插值是一种基于局部插值的方法。
它通过在每个局部区域上拟合一个低阶多项式来实现插值。
这些局部多项式在相邻区域处满足平滑和连续性条件，从而得到整体平滑的插值结果。
除了上述方法外，还有其他一些二维插值方法，如反距离加权插值和克里金插值等。
总的来说，二维插值通过已知数据点之间的关系来估计未知位置上的数值。
不同的插值方法在计算复杂度、精度和平滑性等方面存在差异，根据具体应用场景的需求，选择合适的插值方法是非常重要的。
interp2d的用法python
interp2d的用法python在Python中，我们经常需要对二维数据进行插值以得到缺失或者非均匀采样的数据点的估计值。
Python的SciPy库提供了interp2d函数，可以通过二维插值方法来估计这些缺失值。
interp2d函数的基本用法如下：```pythonimport numpy as npfrom scipy.interpolate import interp2d# 创建x、y坐标轴x = np.arange(0, 10, 1)y = np.arange(0, 10, 1)# 创建z二维数据z = np.random.rand(10, 10)# 创建插值函数f = interp2d(x, y, z, kind='linear')# 输入新的x、y坐标获取估计值x_new = np.arange(0, 10, 0.5)y_new = np.arange(0, 10, 0.5)z_new = f(x_new, y_new)# 打印估计值print(z_new)```上述代码中，我们首先导入numpy和interp2d函数。
然后创建了一个二维数据z，以及对应的x和y坐标轴。
接下来，我们使用interp2d函数创建了一个插值函数f。
这里使用了线性插值的方法，也可以使用其他的插值方法，如'nearest'最近邻插值、'cubic'三次样条插值等。
然后，我们输入新的x_new和y_new坐标，通过调用插值函数f，得到了估计值z_new。
在实际应用中，interp2d函数的使用非常灵活。
我们可以根据实际情况选择不同的插值方法，也可以根据需要设置其他参数来调整估计值的准确性。
另外，值得注意的是，传递给interp2d函数的x、y坐标轴和z二维数据必须是一维的。
如果传入的是二维数据，需要使用flatten()函数将其转换为一维。
另外，插值函数f返回的估计值是一个二维的数组，其形状与输入的x_new和y_new坐标的网格大小相同。
二维插值算法原理
二维插值算法原理二维插值算法是一种在二维空间中根据已知的数据点来估计未知位置上的数值的算法。
它广泛应用于图像处理、地理信息系统和数值模拟等领域。
其原理是基于数学上的连续性和局部平滑性假设，通过利用已知数据点的信息，对未知位置上的数值进行估计。
二维插值算法的基本思想是根据已知的数据点的数值和位置，构建出一个合适的数学模型。
对于每一个未知位置，通过模型可以预测其数值。
这个模型常常是一个多项式函数或者其它形式的连续函数，以便于能够在整个二维空间中插值。
其中最常见的二维插值算法是双线性插值。
双线性插值法假设每个未知位置上的数值都是由其相邻四个已知点的数值线性插值得到的。
具体而言，假设已知的四个点为A、B、C、D，它们的数值分别为f(A)、f(B)、f(C)、f(D)。
对于未知位置P，可以通过以下公式计算得到其数值f(P)：f(P) = (1 - u)(1 - v) f(A) + u(1 - v) f(B) + (1 - u)v f(C) + uv f(D)其中，u和v是分别表示未知位置在水平和垂直方向上的相对位置的权重。
这种方法实现简单，计算效率高，可以较为准确地插值出未知位置上的数值。
除了双线性插值之外，还有其它一些更复杂的二维插值算法，如三次样条插值、Kriging插值等。
这些算法在不同的应用场景下具有不同的优势。
例如，三次样条插值在处理光滑函数时效果较好，而Kriging插值则适用于处理具有空间相关性的数据。
选择适合的插值算法可以提高插值结果的质量。
在实际应用中，二维插值算法在处理图像、地理数据和模拟结果等方面具有重要意义。
通过插值算法，可以将有限的离散数据转换为连续的函数，从而对未知位置上的数值进行预测和分析。
同时，它也为数据的可视化提供了基础，使得我们能够更直观地理解数据的分布和变化规律。
总之，二维插值算法是一种有指导意义的数学工具，它通过在二维空间中根据有限的已知数据点估计未知位置上的数值。
二维线性插值方法及其在平面温度场计算中的应用
［ ］ 万和 ， 显 峰 ． ２叶 周 Ｂ Ｔ模 式 ， 应 成 为 政 府 垫 资工 程 的合 法 外 不
衣 ［］施 工企 业 管理 ， ０ ６ １１） ２ －２ ． Ｊ． ２０ （ ： ４ ５
［ ］ 歌 ，蒙 鹏 程 ． 于 Ｂ 模 式在 我 国 运 用 的探 讨 ［ ］ 程 建 ３黄 关 Ｔ Ｊ． 工
性， 在投 资建 设 领 域 发 挥 市 场 机 制 的 积极 作 用 。在 城 市 基础 设 施建设 中， 当前 的 政 策 鼓 励 和 允许 采 用 多 种 投 融 资模 式 、 种 多
建 设 方式 进 行城 市 基础 设施 项 目的建 设 和 经 营 。 Ｂ Ｔ建 设 方 式
下 ， 受 到鼓 励 和 支 持 的。 一 是
建 设 程 序 、按 照相 关 制 度 和 法 律 法 规 要 求 通 过 招 投 标 公 平 合 法 地 获 得 政 府 投 资项 目建 设 的 建 设 方式 ，在 当前 的政 策框 架
随 着 投 资 体 制 改 革 的 深 入 ，进 一 步 放 开 社 会 资 金 进 入 的
领域 ， 已经 是 各 级 政 府 的 共 识 ， 以期 充 分调 动 社 会 资 金 的积 极
ｒ１ ２ｘ １ ｒ ３
：
．
Ｐ 为 曲 线 上 ｘ＝ｆ ｘ ） 已 知 点 ， 标 分 别 为 （ １， ） （ ， （ ’两 坐 ｘ’Ｘ 和 Ｘ
（ ２－Ｉ Ｘ２Ｘ２ｇ ×ｇ－ ｘ２ １ （３．１ ｇ ×ｇ Ｘ１ ！（ ３ １ １ ２ ２ Ｘ２ ｘ １ ２ １ － ． ． Ｘ － ） ＂ ｆ ＾ ） － ） － ． ￣
（ ４ Ｘ ｅ ＝０’ Ｘｋ １） －
ｆＸＸ 。价 １，于 Ｘ１３ 等 ２ 一ｌ １ １＝ － 一 － Ｘ １
1二维插值算法与实现
1二维插值算法与实现二维插值算法是一种在二维坐标系上进行插值的技术。
它可以根据已知数据点的值，在未知数据点上推断出一个逼近该值的估计值。
二维插值算法广泛应用于图像处理、地理信息系统、气象学等领域。
最常用的二维插值算法有线性插值和双线性插值。
线性插值算法在二维坐标系上根据已知数据点之间的线性关系进行推断。
双线性插值算法则利用已知数据点周围的四个最近邻数据点之间的线性关系，并根据权重进行加权平均来估计未知数据点的值。
下面将对线性插值和双线性插值算法的实现进行详细介绍。
1.线性插值算法实现：线性插值算法的思想是根据已知数据点之间的线性关系推断未知数据点的值。
假设有两个已知数据点(x1,y1)和(x2,y2)，目标是在这两个点之间的未知坐标(x,y)上估计一个值。
算法的步骤如下：- 计算坐标点x在已知数据点x1和x2之间的相对位置，即插值比例factor = (x - x1) / (x2 - x1)；- 通过线性关系计算该未知坐标上的估计值，即y = y1 + (y2 - y1) * factor。
线性插值算法的实现过程如下所示：```pythondef linear_interpolation(x1, y1, x2, y2, x):factor = (x - x1) / (x2 - x1)y = y1 + (y2 - y1) * factorreturn y```2.双线性插值算法实现：双线性插值算法是在二维坐标系上进行插值的技术，它利用已知数据点周围的四个最近邻数据点之间的线性关系来估计未知数据点的值。
假设已知数据点分别为(x1,y1,v1)、(x1,y2,v2)、(x2,y1,v3)和(x2,y2,v4)，目标是在未知坐标(x,y)上估计一个值。
算法的步骤如下：- 计算坐标点x和y在已知数据点x1、y1、x2和y2所构成的矩形区域内的相对位置，即插值比例factor_x = (x - x1) / (x2 - x1) 和factor_y = (y - y1) / (y2 - y1)；- 分别在x和y方向上进行线性插值得到两个估计值，即v_a = v1+ (v3 - v1) * factor_x 和 v_b = v2 + (v4 - v2) * factor_x；- 在v_a和v_b之间进行线性插值，得到在未知坐标上的估计值 v = v_a + (v_b - v_a) * factor_y。
griddata二维插值
griddata⼆维插值 对某些设备或测量仪器来说，采集的数据点的位置不是规则排列的⽹格结构（可参考），对于这种数据⽤散点图（每个采样点具有不同的值或权重）不能很好的展⽰其内部结构，因此需要对其进⾏插值，⽣成⼀个规则的栅格图像。
可采⽤griddata函数对已知的数据点进⾏插值，数据点（X, Y）不要求规则排列。
下图分别使⽤Nearest、Linear、Cubic三种插值⽅法对数据点进⾏插值。
插值算法都要⾸先⾯对⼀个共同的问题—— 邻近点的选择（利⽤靠近插值点附近的多个邻近点，构造插值函数计算插值点的函数值）。
应尽可能使所选择的邻近点均匀地分布在待估点周围。
因为使⽤适量的已知点对于插值的精度很重要，当已知点过多时，会使插值准确率下降，因为过多的信息量会掩盖有⽤的信息；被选择的邻近点构成的点分布不均匀，若某个⽅位上的数据点过多，或点集的某些数据点位置较集中，或者数据点过远，都会带来较⼤的插值误差。
griddata函数内部会先对已知的数据点进⾏Delaunay三⾓剖分，当完成三⾓剖分完成后，根据每个三⾓⽚顶点的数值，就可以对三⾓形区域内任意点进⾏线性插值（参考LinearNDInterpolator函数）。
对三⾓⽚内的插值点可采⽤质⼼插值，即是对插值点只考虑与该点最邻近周围点的影响，确定出插值点与最邻近的周围点之相互位置关系，求出与周围点的影响权重因⼦，以此建⽴线性插值公式（参考）。
"""Simple N-D interpolation.. versionadded:: 0.9"""## Copyright (C) Pauli Virtanen, 2010.## Distributed under the same BSD license as Scipy.### Note: this file should be run through the Mako template engine before# feeding it to Cython.## Run ``generate_qhull.py`` to regenerate the ``qhull.c`` file#cimport cythonfrom libc.float cimport DBL_EPSILONfrom libc.math cimport fabs, sqrtimport numpy as npimport scipy.spatial.qhull as qhullcimport scipy.spatial.qhull as qhullimport warnings#------------------------------------------------------------------------------# Numpy etc.#------------------------------------------------------------------------------cdef extern from"numpy/ndarrayobject.h":cdef enum:NPY_MAXDIMSctypedef fused double_or_complex:doubledouble complex#------------------------------------------------------------------------------# Interpolator base class#------------------------------------------------------------------------------class NDInterpolatorBase(object):"""Common routines for interpolators... versionadded:: 0.9"""def__init__(self, points, values, fill_value=np.nan, ndim=None,rescale=False, need_contiguous=True, need_values=True):"""Check shape of points and values arrays, and reshape values to(npoints, nvalues). Ensure the `points` and values arrays areC-contiguous, and of correct type."""if isinstance(points, qhull.Delaunay):# Precomputed triangulation was passed inif rescale:raise ValueError("Rescaling is not supported when passing ""a Delaunay triangulation as ``points``.")self.tri = pointspoints = points.pointselse:self.tri = Nonepoints = _ndim_coords_from_arrays(points)values = np.asarray(values)_check_init_shape(points, values, ndim=ndim)if need_contiguous:points = np.ascontiguousarray(points, dtype=np.double)if need_values:self.values_shape = values.shape[1:]if values.ndim == 1:self.values = values[:,None]elif values.ndim == 2:self.values = valueselse:self.values = values.reshape(values.shape[0],np.prod(values.shape[1:]))# Complex or real?self.is_complex = np.issubdtype(self.values.dtype, plexfloating) if self.is_complex:if need_contiguous:self.values = np.ascontiguousarray(self.values,dtype=plex128)self.fill_value = complex(fill_value)else:if need_contiguous:self.values = np.ascontiguousarray(self.values, dtype=np.double) self.fill_value = float(fill_value)if not rescale:self.scale = Noneself.points = pointselse:# scale to unit cube centered at 0self.offset = np.mean(points, axis=0)self.points = points - self.offsetself.scale = self.points.ptp(axis=0)self.scale[~(self.scale > 0)] = 1.0 # avoid division by 0self.points /= self.scaledef _check_call_shape(self, xi):xi = np.asanyarray(xi)if xi.shape[-1] != self.points.shape[1]:raise ValueError("number of dimensions in xi does not match x")return xidef _scale_x(self, xi):if self.scale is None:return xielse:return (xi - self.offset) / self.scaledef__call__(self, *args):"""interpolator(xi)Evaluate interpolator at given points.Parameters----------x1, x2, ... xn: array-like of floatPoints where to interpolate data at.x1, x2, ... xn can be array-like of float with broadcastable shape.or x1 can be array-like of float with shape ``(..., ndim)``"""xi = _ndim_coords_from_arrays(args, ndim=self.points.shape[1])xi = self._check_call_shape(xi)shape = xi.shapexi = xi.reshape(-1, shape[-1])xi = np.ascontiguousarray(xi, dtype=np.double)xi = self._scale_x(xi)if self.is_complex:r = self._evaluate_complex(xi)else:r = self._evaluate_double(xi)return np.asarray(r).reshape(shape[:-1] + self.values_shape)cpdef _ndim_coords_from_arrays(points, ndim=None):"""Convert a tuple of coordinate arrays to a (..., ndim)-shaped array."""cdef ssize_t j, nif isinstance(points, tuple) and len(points) == 1:# handle argument tuplepoints = points[0]if isinstance(points, tuple):p = np.broadcast_arrays(*points)n = len(p)for j in range(1, n):if p[j].shape != p[0].shape:raise ValueError("coordinate arrays do not have the same shape") points = np.empty(p[0].shape + (len(points),), dtype=float) for j, item in enumerate(p):points[...,j] = itemelse:points = np.asanyarray(points)if points.ndim == 1:if ndim is None:points = points.reshape(-1, 1)else:points = points.reshape(-1, ndim)return pointscdef _check_init_shape(points, values, ndim=None):"""Check shape of points and values arrays"""if values.shape[0] != points.shape[0]:raise ValueError("different number of values and points")if points.ndim != 2:raise ValueError("invalid shape for input data points")if points.shape[1] < 2:raise ValueError("input data must be at least 2-D")if ndim is not None and points.shape[1] != ndim:raise ValueError("this mode of interpolation available only for ""%d-D data" % ndim)#------------------------------------------------------------------------------# Linear interpolation in N-D#------------------------------------------------------------------------------class LinearNDInterpolator(NDInterpolatorBase):"""LinearNDInterpolator(points, values, fill_value=np.nan, rescale=False)Piecewise linear interpolant in N dimensions... versionadded:: 0.9Methods-------__call__Parameters----------points : ndarray of floats, shape (npoints, ndims); or DelaunayData point coordinates, or a precomputed Delaunay triangulation.values : ndarray of float or complex, shape (npoints, ...)Data values.fill_value : float, optionalValue used to fill in for requested points outside of theconvex hull of the input points. If not provided, thenthe default is ``nan``.rescale : bool, optionalRescale points to unit cube before performing interpolation.This is useful if some of the input dimensions haveincommensurable units and differ by many orders of magnitude.Notes-----The interpolant is constructed by triangulating the input datawith Qhull [1]_, and on each triangle performing linearbarycentric interpolation.Examples--------We can interpolate values on a 2D plane:>>> from scipy.interpolate import LinearNDInterpolator>>> import matplotlib.pyplot as plt>>> np.random.seed(0)>>> x = np.random.random(10) - 0.5>>> y = np.random.random(10) - 0.5>>> z = np.hypot(x, y)>>> X = np.linspace(min(x), max(x))>>> Y = np.linspace(min(y), max(y))>>> X, Y = np.meshgrid(X, Y) # 2D grid for interpolation>>> interp = LinearNDInterpolator(list(zip(x, y)), z)>>> Z = interp(X, Y)>>> plt.pcolormesh(X, Y, Z, shading='auto')>>> plt.plot(x, y, "ok", label="input point")>>> plt.legend()>>> plt.colorbar()>>> plt.axis("equal")>>> plt.show()See also--------griddata :Interpolate unstructured D-D data.NearestNDInterpolator :Nearest-neighbor interpolation in N dimensions.CloughTocher2DInterpolator :Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D. References----------.. [1] /"""def__init__(self, points, values, fill_value=np.nan, rescale=False):NDInterpolatorBase.__init__(self, points, values, fill_value=fill_value, rescale=rescale)if self.tri is None:self.tri = qhull.Delaunay(self.points)def _evaluate_double(self, xi):return self._do_evaluate(xi, 1.0)def _evaluate_complex(self, xi):return self._do_evaluate(xi, 1.0j)@cython.boundscheck(False)@cython.wraparound(False)def _do_evaluate(self, const double[:,::1] xi, double_or_complex dummy): cdef const double_or_complex[:,::1] values = self.valuescdef double_or_complex[:,::1] outcdef const double[:,::1] points = self.pointscdef const int[:,::1] simplices = self.tri.simplicescdef double c[NPY_MAXDIMS]cdef double_or_complex fill_valuecdef int i, j, k, m, ndim, isimplex, inside, start, nvaluescdef qhull.DelaunayInfo_t infocdef double eps, eps_broadndim = xi.shape[1]start = 0fill_value = self.fill_valueqhull._get_delaunay_info(&info, self.tri, 1, 0, 0)out = np.empty((xi.shape[0], self.values.shape[1]),dtype=self.values.dtype)nvalues = out.shape[1]eps = 100 * DBL_EPSILONeps_broad = sqrt(DBL_EPSILON)with nogil:for i in range(xi.shape[0]):# 1) Find the simplexisimplex = qhull._find_simplex(&info, c,&xi[0,0] + i*ndim,&start, eps, eps_broad)# 2) Linear barycentric interpolationif isimplex == -1:# don't extrapolatefor k in range(nvalues):out[i,k] = fill_valuecontinuefor k in range(nvalues):out[i,k] = 0for j in range(ndim+1):for k in range(nvalues):m = simplices[isimplex,j]out[i,k] = out[i,k] + c[j] * values[m,k]return out#------------------------------------------------------------------------------# Gradient estimation in 2D#------------------------------------------------------------------------------class GradientEstimationWarning(Warning):pass@cython.cdivision(True)cdef int _estimate_gradients_2d_global(qhull.DelaunayInfo_t *d, double *data, int maxiter, double tol,double *y) nogil:"""Estimate gradients of a function at the vertices of a 2d triangulation.Parameters----------info : inputTriangulation in 2Ddata : inputFunction values at the verticesmaxiter : inputMaximum number of Gauss-Seidel iterationstol : inputAbsolute / relative stop tolerancey : output, shape (npoints, 2)Derivatives [F_x, F_y] at the verticesReturns-------num_iterationsNumber of iterations if converged, 0 if maxiter reachedwithout convergenceNotes-----This routine uses a re-implementation of the global approximatecurvature minimization algorithm described in [Nielson83] and [Renka84]. References----------.. [Nielson83] G. Nielson,''A method for interpolating scattered data based upon a minimum norm network''.Math. Comp., 40, 253 (1983)... [Renka84] R. J. Renka and A. K. Cline.''A Triangle-based C1 interpolation method.'',Rocky Mountain J. Math., 14, 223 (1984)."""cdef double Q[2*2]cdef double s[2]cdef double r[2]cdef int ipoint, iiter, k, ipoint2, jpoint2cdef double f1, f2, df2, ex, ey, L, L3, det, err, change# initializefor ipoint in range(2*d.npoints):y[ipoint] = 0## Main point:## Z = sum_T sum_{E in T} int_E |W''|^2 = min!## where W'' is the second derivative of the Clough-Tocher# interpolant to the direction of the edge E in triangle T.## The minimization is done iteratively: for each vertex V,# the sum## Z_V = sum_{E connected to V} int_E |W''|^2## is minimized separately, using existing values at other V.## Since the interpolant can be written as## W(x) = f(x) + w(x)^T y## where y = [ F_x(V); F_y(V) ], it is clear that the solution to# the local problem is is given as a solution of the 2x2 matrix# equation.## Here, we use the Clough-Tocher interpolant, which restricted to# a single edge is## w(x) = (1 - x)**3 * f1# + x*(1 - x)**2 * (df1 + 3*f1)# + x**2*(1 - x) * (df2 + 3*f2)# + x**3 * f2## where f1, f2 are values at the vertices, and df1 and df2 are# derivatives along the edge (away from the vertices).## As a consequence, one finds## L^3 int_{E} |W''|^2 = y^T A y + 2 B y + C## with## A = [4, -2; -2, 4]# B = [6*(f1 - f2), 6*(f2 - f1)]# y = [df1, df2]# L = length of edge E## and C is not needed for minimization. Since df1 = dF1.E, df2 = -dF2.E, # with dF1 = [F_x(V_1), F_y(V_1)], and the edge vector E = V2 - V1,# we have## Z_V = dF1^T Q dF1 + 2 s.dF1 + const.## which is minimized by## dF1 = -Q^{-1} s## where## Q = sum_E [A_11 E E^T]/L_E^3 = 4 sum_E [E E^T]/L_E^3# s = sum_E [ B_1 + A_21 df2] E /L_E^3# = sum_E [ 6*(f1 - f2) + 2*(E.dF2)] E / L_E^3## Gauss-Seidelfor iiter in range(maxiter):err = 0for ipoint in range(d.npoints):for k in range(2*2):Q[k] = 0for k in range(2):s[k] = 0# walk over neighbours of given pointfor jpoint2 in range(d.vertex_neighbors_indptr[ipoint],d.vertex_neighbors_indptr[ipoint+1]):ipoint2 = d.vertex_neighbors_indices[jpoint2]# edgeex = d.points[2*ipoint2 + 0] - d.points[2*ipoint + 0]ey = d.points[2*ipoint2 + 1] - d.points[2*ipoint + 1]L = sqrt(ex**2 + ey**2)L3 = L*L*L# data at verticesf1 = data[ipoint]f2 = data[ipoint2]# scaled gradient projections on the edgedf2 = -ex*y[2*ipoint2 + 0] - ey*y[2*ipoint2 + 1]# edge sumQ[0] += 4*ex*ex / L3Q[1] += 4*ex*ey / L3Q[3] += 4*ey*ey / L3s[0] += (6*(f1 - f2) - 2*df2) * ex / L3s[1] += (6*(f1 - f2) - 2*df2) * ey / L3Q[2] = Q[1]# solvedet = Q[0]*Q[3] - Q[1]*Q[2]r[0] = ( Q[3]*s[0] - Q[1]*s[1])/detr[1] = (-Q[2]*s[0] + Q[0]*s[1])/detchange = max(fabs(y[2*ipoint + 0] + r[0]),fabs(y[2*ipoint + 1] + r[1]))y[2*ipoint + 0] = -r[0]y[2*ipoint + 1] = -r[1]# relative/absolute errorchange /= max(1.0, max(fabs(r[0]), fabs(r[1])))err = max(err, change)if err < tol:return iiter + 1# Didn't converge before maxiterreturn 0@cython.boundscheck(False)@cython.wraparound(False)cpdef estimate_gradients_2d_global(tri, y, int maxiter=400, double tol=1e-6): cdef const double[:,::1] datacdef double[:,:,::1] gradcdef qhull.DelaunayInfo_t infocdef int k, ret, nvaluesy = np.asanyarray(y)if y.shape[0] != tri.npoints:raise ValueError("'y' has a wrong number of items")if np.issubdtype(y.dtype, plexfloating):rg = estimate_gradients_2d_global(tri, y.real, maxiter=maxiter, tol=tol) ig = estimate_gradients_2d_global(tri, y.imag, maxiter=maxiter, tol=tol) r = np.zeros(rg.shape, dtype=complex)r.real = rgr.imag = igreturn ry_shape = y.shapeif y.ndim == 1:y = y[:,None]y = y.reshape(tri.npoints, -1).Ty = np.ascontiguousarray(y, dtype=np.double)yi = np.empty((y.shape[0], y.shape[1], 2))data = ygrad = yiqhull._get_delaunay_info(&info, tri, 0, 0, 1)nvalues = data.shape[0]for k in range(nvalues):with nogil:ret = _estimate_gradients_2d_global(&info,&data[k,0],maxiter,tol,&grad[k,0,0])if ret == 0:warnings.warn("Gradient estimation did not converge, ""the results may be inaccurate",GradientEstimationWarning)return yi.transpose(1, 0, 2).reshape(y_shape + (2,))#------------------------------------------------------------------------------# Cubic interpolation in 2D#------------------------------------------------------------------------------@cython.cdivision(True)cdef double_or_complex _clough_tocher_2d_single(qhull.DelaunayInfo_t *d, int isimplex,double *b,double_or_complex *f,double_or_complex *df) nogil:"""Evaluate Clough-Tocher interpolant on a 2D triangle.Parameters----------d :Delaunay infoisimplex : intTriangle to evaluate onb : shape (3,)Barycentric coordinates of the point on the trianglef : shape (3,)Function values at verticesdf : shape (3, 2)Gradient values at verticesReturns-------w :Value of the interpolant at the given pointReferences----------.. [CT] See, for example,P. Alfeld,''A trivariate Clough-Tocher scheme for tetrahedral data''.Computer Aided Geometric Design, 1, 169 (1984);G. Farin,''Triangular Bernstein-Bezier patches''.Computer Aided Geometric Design, 3, 83 (1986)."""cdef double_or_complex \c3000, c0300, c0030, c0003, \c2100, c2010, c2001, c0210, c0201, c0021, \c1200, c1020, c1002, c0120, c0102, c0012, \c1101, c1011, c0111cdef double_or_complex \f1, f2, f3, df12, df13, df21, df23, df31, df32cdef double g[3]cdef double \e12x, e12y, e23x, e23y, e31x, e31y, \e14x, e14y, e24x, e24y, e34x, e34ycdef double_or_complex wcdef double minvalcdef double b1, b2, b3, b4cdef int k, itricdef double c[3]cdef double y[2]# XXX: optimize + refactor this!e12x = (+ d.points[0 + 2*d.simplices[3*isimplex + 1]]- d.points[0 + 2*d.simplices[3*isimplex + 0]])e12y = (+ d.points[1 + 2*d.simplices[3*isimplex + 1]]- d.points[1 + 2*d.simplices[3*isimplex + 0]])e23x = (+ d.points[0 + 2*d.simplices[3*isimplex + 2]]- d.points[0 + 2*d.simplices[3*isimplex + 1]])e23y = (+ d.points[1 + 2*d.simplices[3*isimplex + 2]]- d.points[1 + 2*d.simplices[3*isimplex + 1]])e31x = (+ d.points[0 + 2*d.simplices[3*isimplex + 0]]- d.points[0 + 2*d.simplices[3*isimplex + 2]])e31y = (+ d.points[1 + 2*d.simplices[3*isimplex + 0]]- d.points[1 + 2*d.simplices[3*isimplex + 2]])e14x = (e12x - e31x)/3e14y = (e12y - e31y)/3e24x = (-e12x + e23x)/3e24y = (-e12y + e23y)/3e34x = (e31x - e23x)/3e34y = (e31y - e23y)/3f1 = f[0]f2 = f[1]f3 = f[2]df12 = +(df[2*0+0]*e12x + df[2*0+1]*e12y)df21 = -(df[2*1+0]*e12x + df[2*1+1]*e12y)df23 = +(df[2*1+0]*e23x + df[2*1+1]*e23y)df32 = -(df[2*2+0]*e23x + df[2*2+1]*e23y)df31 = +(df[2*2+0]*e31x + df[2*2+1]*e31y)df13 = -(df[2*0+0]*e31x + df[2*0+1]*e31y)c3000 = f1c2100 = (df12 + 3*c3000)/3c2010 = (df13 + 3*c3000)/3c0300 = f2c1200 = (df21 + 3*c0300)/3c0210 = (df23 + 3*c0300)/3c0030 = f3c1020 = (df31 + 3*c0030)/3c0120 = (df32 + 3*c0030)/3c2001 = (c2100 + c2010 + c3000)/3c0201 = (c1200 + c0300 + c0210)/3c0021 = (c1020 + c0120 + c0030)/3## Now, we need to impose the condition that the gradient of the spline # to some direction `w` is a linear function along the edge.## As long as two neighbouring triangles agree on the choice of the# direction `w`, this ensures global C1 differentiability.# Otherwise, the choice of the direction is arbitrary (except that# it should not point along the edge, of course).## In [CT]_, it is suggested to pick `w` as the normal of the edge.# This choice is given by the formulas## w_12 = E_24 + g[0] * E_23# w_23 = E_34 + g[1] * E_31# w_31 = E_14 + g[2] * E_12## g[0] = -(e24x*e23x + e24y*e23y) / (e23x**2 + e23y**2)# g[1] = -(e34x*e31x + e34y*e31y) / (e31x**2 + e31y**2)# g[2] = -(e14x*e12x + e14y*e12y) / (e12x**2 + e12y**2)## However, this choice gives an interpolant that is *not*# invariant under affine transforms. This has some bad# consequences: for a very narrow triangle, the spline can# develops huge oscillations. For instance, with the input data## [(0, 0), (0, 1), (eps, eps)], eps = 0.01# F = [0, 0, 1]# dF = [(0,0), (0,0), (0,0)]## one observes that as eps -> 0, the absolute maximum value of the # interpolant approaches infinity.## So below, we aim to pick affine invariant `g[k]`.# We choose## w = V_4' - V_4## where V_4 is the centroid of the current triangle, and V_4' the# centroid of the neighbour. Since this quantity transforms similarly # as the gradient under affine transforms, the resulting interpolant# is affine-invariant. Moreover, two neighbouring triangles clearly# always agree on the choice of `w` (sign is unimportant), and so# this choice also makes the interpolant C1.## The drawback here is a performance penalty, since we need to# peek into neighbouring triangles.#for k in range(3):itri = d.neighbors[3*isimplex + k]if itri == -1:# No neighbour.# Compute derivative to the centroid direction (e_12 + e_13)/2. g[k] = -1./2continue# Centroid of the neighbour, in our local barycentric coordinates y[0] = (+ d.points[0 + 2*d.simplices[3*itri + 0]]+ d.points[0 + 2*d.simplices[3*itri + 1]]+ d.points[0 + 2*d.simplices[3*itri + 2]]) / 3y[1] = (+ d.points[1 + 2*d.simplices[3*itri + 0]]+ d.points[1 + 2*d.simplices[3*itri + 1]]+ d.points[1 + 2*d.simplices[3*itri + 2]]) / 3qhull._barycentric_coordinates(2, d.transform + isimplex*2*3, y, c) # Rewrite V_4'-V_4 = const*[(V_4-V_2) + g_i*(V_3 - V_2)]# Now, observe that the results can be written *in terms of# barycentric coordinates*. Barycentric coordinates stay# invariant under affine transformations, so we can directly# conclude that the choice below is affine-invariant.if k == 0:g[k] = (2*c[2] + c[1] - 1) / (2 - 3*c[2] - 3*c[1])elif k == 1:g[k] = (2*c[0] + c[2] - 1) / (2 - 3*c[0] - 3*c[2])elif k == 2:g[k] = (2*c[1] + c[0] - 1) / (2 - 3*c[1] - 3*c[0])c0111 = (g[0]*(-c0300 + 3*c0210 - 3*c0120 + c0030)+ (-c0300 + 2*c0210 - c0120 + c0021 + c0201))/2c1011 = (g[1]*(-c0030 + 3*c1020 - 3*c2010 + c3000)+ (-c0030 + 2*c1020 - c2010 + c2001 + c0021))/2c1101 = (g[2]*(-c3000 + 3*c2100 - 3*c1200 + c0300)+ (-c3000 + 2*c2100 - c1200 + c2001 + c0201))/2c1002 = (c1101 + c1011 + c2001)/3c0102 = (c1101 + c0111 + c0201)/3c0012 = (c1011 + c0111 + c0021)/3c0003 = (c1002 + c0102 + c0012)/3# extended barycentric coordinatesminval = b[0]for k in range(3):if b[k] < minval:minval = b[k]b1 = b[0] - minvalb2 = b[1] - minvalb3 = b[2] - minvalb4 = 3*minval# evaluate the polynomial -- the stupid and ugly way to do it,# one of the 4 coordinates is in fact zerow = (b1**3*c3000 + 3*b1**2*b2*c2100 + 3*b1**2*b3*c2010 +3*b1**2*b4*c2001 + 3*b1*b2**2*c1200 +6*b1*b2*b4*c1101 + 3*b1*b3**2*c1020 + 6*b1*b3*b4*c1011 +3*b1*b4**2*c1002 + b2**3*c0300 + 3*b2**2*b3*c0210 +3*b2**2*b4*c0201 + 3*b2*b3**2*c0120 + 6*b2*b3*b4*c0111 +3*b2*b4**2*c0102 + b3**3*c0030 + 3*b3**2*b4*c0021 +3*b3*b4**2*c0012 + b4**3*c0003)return wclass CloughTocher2DInterpolator(NDInterpolatorBase):"""CloughTocher2DInterpolator(points, values, tol=1e-6)Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D. .. versionadded:: 0.9Methods-------__call__Parameters----------points : ndarray of floats, shape (npoints, ndims); or DelaunayData point coordinates, or a precomputed Delaunay triangulation. values : ndarray of float or complex, shape (npoints, ...)Data values.fill_value : float, optionalValue used to fill in for requested points outside of theconvex hull of the input points. If not provided, thenthe default is ``nan``.tol : float, optionalAbsolute/relative tolerance for gradient estimation.maxiter : int, optionalMaximum number of iterations in gradient estimation.rescale : bool, optionalRescale points to unit cube before performing interpolation.This is useful if some of the input dimensions haveincommensurable units and differ by many orders of magnitude.Notes-----The interpolant is constructed by triangulating the input datawith Qhull [1]_, and constructing a piecewise cubicinterpolating Bezier polynomial on each triangle, using aClough-Tocher scheme [CT]_. The interpolant is guaranteed to becontinuously differentiable.The gradients of the interpolant are chosen so that the curvatureof the interpolating surface is approximatively minimized. Thegradients necessary for this are estimated using the globalalgorithm described in [Nielson83]_ and [Renka84]_.Examples--------We can interpolate values on a 2D plane:>>> from scipy.interpolate import CloughTocher2DInterpolator>>> import matplotlib.pyplot as plt>>> np.random.seed(0)>>> x = np.random.random(10) - 0.5>>> y = np.random.random(10) - 0.5>>> z = np.hypot(x, y)>>> X = np.linspace(min(x), max(x))>>> Y = np.linspace(min(y), max(y))>>> X, Y = np.meshgrid(X, Y) # 2D grid for interpolation>>> interp = CloughTocher2DInterpolator(list(zip(x, y)), z)>>> Z = interp(X, Y)>>> plt.pcolormesh(X, Y, Z, shading='auto')>>> plt.plot(x, y, "ok", label="input point")>>> plt.legend()>>> plt.colorbar()>>> plt.axis("equal")>>> plt.show()See also--------griddata :Interpolate unstructured D-D data.LinearNDInterpolator :Piecewise linear interpolant in N dimensions.NearestNDInterpolator :Nearest-neighbor interpolation in N dimensions.References----------.. [1] /.. [CT] See, for example,P. Alfeld,''A trivariate Clough-Tocher scheme for tetrahedral data''.Computer Aided Geometric Design, 1, 169 (1984);G. Farin,''Triangular Bernstein-Bezier patches''.Computer Aided Geometric Design, 3, 83 (1986)... [Nielson83] G. Nielson,''A method for interpolating scattered data based upon a minimum norm network''.Math. Comp., 40, 253 (1983)... [Renka84] R. J. Renka and A. K. Cline.''A Triangle-based C1 interpolation method.'',Rocky Mountain J. Math., 14, 223 (1984)."""def__init__(self, points, values, fill_value=np.nan,tol=1e-6, maxiter=400, rescale=False):NDInterpolatorBase.__init__(self, points, values, ndim=2,fill_value=fill_value, rescale=rescale)if self.tri is None:self.tri = qhull.Delaunay(self.points)self.grad = estimate_gradients_2d_global(self.tri, self.values,tol=tol, maxiter=maxiter)def _evaluate_double(self, xi):return self._do_evaluate(xi, 1.0)def _evaluate_complex(self, xi):return self._do_evaluate(xi, 1.0j)@cython.boundscheck(False)@cython.wraparound(False)def _do_evaluate(self, const double[:,::1] xi, double_or_complex dummy):cdef const double_or_complex[:,::1] values = self.valuescdef const double_or_complex[:,:,:] grad = self.gradcdef double_or_complex[:,::1] outcdef const double[:,::1] points = self.pointscdef const int[:,::1] simplices = self.tri.simplicescdef double c[NPY_MAXDIMS]cdef double_or_complex f[NPY_MAXDIMS+1]cdef double_or_complex df[2*NPY_MAXDIMS+2]cdef double_or_complex wcdef double_or_complex fill_valuecdef int i, j, k, m, ndim, isimplex, inside, start, nvaluescdef qhull.DelaunayInfo_t infocdef double eps, eps_broadndim = xi.shape[1]start = 0fill_value = self.fill_valueqhull._get_delaunay_info(&info, self.tri, 1, 1, 0)out = np.zeros((xi.shape[0], self.values.shape[1]),dtype=self.values.dtype)nvalues = out.shape[1]eps = 100 * DBL_EPSILONeps_broad = sqrt(eps)with nogil:for i in range(xi.shape[0]):# 1) Find the simplexisimplex = qhull._find_simplex(&info, c,&xi[i,0],&start, eps, eps_broad)# 2) Clough-Tocher interpolationif isimplex == -1:# outside triangulationfor k in range(nvalues):out[i,k] = fill_valuecontinuefor k in range(nvalues):for j in range(ndim+1):f[j] = values[simplices[isimplex,j],k]df[2*j] = grad[simplices[isimplex,j],k,0]df[2*j+1] = grad[simplices[isimplex,j],k,1]w = _clough_tocher_2d_single(&info, isimplex, c, f, df)out[i,k] = wreturn outView Code 下图中红⾊的是已知采样点，蓝⾊是待插值的栅格点，三⾓形内部栅格点的数值可通过线性插值或其它插值⽅法计算出，三⾓形外部的点可在函数中指定⼀个数值（默认为NaN）。
插值模型
三次样条插值
比分段线性插值更光滑。
y
a
xi-1 xi
bx
在数学上，光滑程度的定量描述是：函数(曲 线)的k阶导数存在且连续，则称该曲线具有k阶光 滑性。
光滑性的阶次越高，则越光滑。是否存在较低
次的分段多项式达到较高阶光滑性的方法？三次 样条插值就是一个很好的例子。
12
三次样条插值
通过全部已知节点,即
再用
计算插值，即
21
第二种（散乱节点）：
y
0
x
22
已知n个节点
其中
互不相同，
构造一个二元函数
通过全部已知节点,即
再用
计算插值，即
返回
23
最邻近插值
y
(x1, y2) (x2,y2)
(x1, y1) (x2, y1)
简记为： f (xi, yj)=f1，f (xi+1, yj)=f2，f (xi+1, yj+1)=f3，f (xi, yj+1)=f4
25
分两片的函数表达式如下： 第一片(下三角形区域)： (x, y)满足
插值函数为：
f (x, y) f1 (f2 f1)(x xi ) (f3 f2 )(y y j )
S(x) {si (x), x [xi1, xi ], i 1,n}
1) si ( x) ai x3 bi x2 ci x di 2) S (xi ) yi (i 0,1,n) 3) S ( x) C 2[ x0 , xn ]
二维插值算法原理
二维插值算法是一种用于在二维空间中估计未知数据点的方法。
它基于已知数据点的值和位置来推断未知数据点的值。
以下是常见的二维插值算法原理之一：双线性插值。
双线性插值是一种基于四个最近邻数据点进行估计的方法。
假设我们有一个二维网格，已知在四个顶点上的数据点的值和位置。
要估计某个位置处的未知数据点的值，双线性插值算法按照以下步骤进行：
1.找到目标位置的最近的四个已知数据点，通常称为左上、右上、左下和右下。
2.计算目标位置相对于这四个已知数据点的相对位置权重。
这可以通过计算目标位置到每个已知数据点的水平和垂直距离，然后根据距离来计算相对权重。
3.根据权重对四个已知数据点的值进行加权平均。
这里的加权平均可以使用线性插值进行计算。
4.得到目标位置的估计值作为插值结果。
双线性插值算法基于以下两个假设：
-在目标位置的附近，插值曲面在水平和垂直方向上是一致的，即呈现线性关系。
-已知数据点之间的变化不会很剧烈，即目标位置与附近已知数据点的值之间存在一定的连续性。
双线性插值算法是一种简单而有效的二维插值方法，适用于平滑、连续变化的数据。
但对于非线性、不规则的数据分布，或者存在边界情况的情况下，可能需要使用其他更复杂的插值算法来获得更准确的估计结果。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
用MATLAB作网格节点数据的插值
z=interp2(x0,y0,z0,x,y,’method’)
被插值点 的函数值
插值 节点
被插值点
插值方法
‘nearest’ 最邻近插值 ‘linear’ 双线性插值 ‘cubic’ 双三次插值 缺省时, 双线性插值
要求x0,y0单调；x，y可取为矩阵，或x取 行向量，y取为列向量，x,y的值分别不能超出 x0,y0的范围。
例：测得平板表面3*5网格点处的温度分别为： 82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 试作出平板表面的温度分布曲面z=f(x,y)的图形。
1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图. 输入以下命令： x=1:5; y=1:3; temps=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86]; mesh(x,y,temps) 2．以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值.
x
O
将四个插值点（矩形的四个顶点）处的函数值依次 简记为： f (xi, yj)=f1，f (xi+1, yj)=f2，f (xi+1, yj+1)=f3，f (xi, yj+1)=f4
分两片的函数表达式如下：
第一片(下三角形区域)： (x, y)满足
插值函数为： f ( x, y) f1 (f 2 f1 )( x x i ) (f 3 f 2 )( y y j )
f=a1+a2/x + + +
f=aebx +
+
-bx f=ae + +
+ +
+ + +
+
+ +
O
x
二维或高维情形的最邻近插值，与被插值点最邻近的 节点的函数值即为所求。 注意：最邻近插值一般不连续。具有连续性的最简单 的插值是分片线性插值。
返回
分片线性插值
y
(xi, yj+1) (xi+1, yj+1) (xi, yj) (xi+1, yj)
通过全部已知节点,即
再用
计算插值，即
第二种（散乱节点）：
y
0
x
已知n个节点
其中 互不相同，
构造一个二元函数
通过全部已知节点,即
再用
计算插值，即
返回
最邻近插值
y
( x1 , y2 )
( x2 , y2 )
( x1 , y1 ) ( x2 , y 1 )
第二片(上三角形区域)：(x, y)满足 y j1 y j y (x x i ) y i x i 1 x i 插值函数为： f ( x, y) f1 (f 4 f1 )( y y j ) (f 3 f 4 )( x x i )
注意：(x, y)当然应该是在插值节点所形成的矩形区 域内。显然，分片线性插值函数是连续的； 返回
1200 X Y 1200 1600 2000 2400 2800 3200 3600 1130 1320 1390 1500 1500 1500 1480 1600 1250 1450 1500 1200 1200 1550 1500 2000 1280 1420 1500 1100 1100 1600 1550 2400 1230 1400 1400 1350 1550 1550 1510 2800 1040 1300 900 1450 1600 1600 1430 3200 900 700 1100 1200 1550 1600 1300 3600 500 900 1060 1150 1380 1600 1200 4000 700 850 950 1010 1070 1550 980
双线性插值
y
( x1 , y2 )
( x2 , y2 )
x
( x1 , y1 ) ( x2 , y1 )
O
双线性插值是一片一片的空间二次曲面构成。 双线性插值函数的形式如下： f ( x, y) (ax b)(cy d) 其中有四个待定系数，利用该函数在矩形的四个顶 点（插值节点）的函数值，得到四个代数方程，正 好确定四个系数。 返回
返回
用MATLAB作散点数据的插值计算
插值函数griddata格式为:
cz =griddata（x，y，z，cx，cy，‘method’） 被插值点 的函数值
插值 节点 被插值点 插值方法
‘nearest’ 最邻近插值 ‘linear’ 双线性插值 ‘cubic’ 双三次插值 'v4'- Matlab提供的插值方法 缺省时, 双线性插值
再输入以下命令:
xi=1:0.2:5; yi=1:0.2:3;
zi=interp2(x,y,temps,xi',yi,'cubic');
mesh(xi,yi,zi) 画出插值后的温度分布曲面图.
例 山区地貌：
在某山区测得一些地点的高程如下表。平面区域为 1200<=x<=4000,1200<=y<=3600) 试作出该山区的地貌图和等高线图。
要求cx取行向量，cy取为列向量。
线性最小二乘拟合 f(x)=a1r1(x)+ …+amrm(x)中 函数{r1(x), …rm(x)}的选取 1. 通过机理分析建立数学模型来确定 f(x)； 2. 将数据 (xi,yi) i=1, …n 作图，通过直观判断确定 f(x)： f=a1+a2x + + + + + f=a1+a2x+a3x2 + + + + + f=a1+a2x+a3x2 + + + + +
二维插值
一、二维插值定义 二、网格节点插值法
最邻近插值 分片线性插值 双线性插值
三、用Matlab解插值问题
网格节点数据的插值
散点数据的插值
返回
二维插值的定义
第一种（网格节点）：
y
O
x
已知 mn个节点
其中 互不相同，不妨设
构造一个二元函数