一位插值、二维插值.

MATLAB中相关命令a aabs 绝对值、模、字符的ascii码值acos 反余弦acosh 反双曲余弦acot 反余切acoth 反双曲余切acsc 反余割acsch 反双曲余割align 启动图形对象几何位置排列工具all 所有元素非零为真angle 相角ans 表达式计算结果的缺省变量名any 所有元素非全零为真area 面域图argnames 函数m文件宗量名asec 反正割asech 反双曲正割asin 反正弦asinh 反双曲正弦assignin 向变量赋值atan 反正切atan2 四象限反正切atanh 反双曲正切autumn 红黄调秋色图阵axes 创建轴对象的低层指令axis 控制轴刻度和风格的高层指令b bbar 二维直方图bar3 三维直方图bar3h 三维水平直方图barh 二维水平直方图base2dec x进制转换为十进制bin2dec 二进制转换为十进制blanks 创建空格串bone 蓝色调黑白色图阵box 框状坐标轴break while 或for 环中断指令brighten 亮度控制c ccapture （3版以前）捕获当前图形cart2pol 直角坐标变为极或柱坐标cart2sph 直角坐标变为球坐标cat 串接成高维数组caxis 色标尺刻度cd 指定当前目录cdedit 启动用户菜单、控件回调函数设计工具cdf2rdf 复数特征值对角阵转为实数块对角阵ceil 向正无穷取整cell 创建元胞数组cell2struct 元胞数组转换为构架数组celldisp 显示元胞数组内容cellplot 元胞数组内部结构图示char 把数值、符号、内联类转换为字符对象chi2cdf 分布累计概率函数chi2inv 分布逆累计概率函数chi2pdf 分布概率密度函数chi2rnd 分布随机数发生器chol cholesky分解clabel 等位线标识cla 清除当前轴class 获知对象类别或创建对象clc 清除指令窗clear 清除内存变量和函数clf 清除图对象clock 时钟colorcube 三浓淡多彩交叉色图矩阵colordef 设置色彩缺省值colormap 色图colspace 列空间的基close 关闭指定窗口colperm 列排序置换向量comet 彗星状轨迹图comet3 三维彗星轨迹图compass 射线图compose 求复合函数cond （逆）条件数condeig 计算特征值、特征向量同时给出条件数condest 范-1条件数估计conj 复数共轭contour 等位线contourf 填色等位线contour3 三维等位线contourslice 四维切片等位线图conv 多项式乘、卷积cool 青紫调冷色图copper 古铜调色图cos 余弦cosh 双曲余弦cot 余切coth 双曲余切cplxpair 复数共轭成对排列csc 余割csch 双曲余割cumsum 元素累计和cumtrapz 累计梯形积分cylinder 创建圆柱d ddblquad 二重数值积分deal 分配宗量deblank 删去串尾部的空格符dec2base 十进制转换为x进制dec2bin 十进制转换为二进制dec2hex 十进制转换为十六进制deconv 多项式除、解卷delaunay delaunay 三角剖分del2 离散laplacian差分demo matlab演示det 行列式diag 矩阵对角元素提取、创建对角阵diary matlab指令窗文本内容记录diff 数值差分、符号微分digits 符号计算中设置符号数值的精度dir 目录列表disp 显示数组display 显示对象内容的重载函数dlinmod 离散系统的线性化模型dmperm 矩阵dulmage-mendelsohn 分解dos 执行dos 指令并返回结果double 把其他类型对象转换为双精度数值drawnow 更新事件队列强迫matlab刷新屏幕dsolve 符号计算解微分方程e eecho m文件被执行指令的显示edit 启动m文件编辑器eig 求特征值和特征向量eigs 求指定的几个特征值end 控制流for等结构体的结尾元素下标eps 浮点相对精度error 显示出错信息并中断执行errortrap 错误发生后程序是否继续执行的控制erf 误差函数erfc 误差补函数erfcx 刻度误差补函数erfinv 逆误差函数errorbar 带误差限的曲线图etreeplot 画消去树eval 串演算指令evalin 跨空间串演算指令exist 检查变量或函数是否已定义exit 退出matlab环境exp 指数函数expand 符号计算中的展开操作expint 指数积分函数expm 常用矩阵指数函数expm1 pade法求矩阵指数expm2 taylor法求矩阵指数expm3 特征值分解法求矩阵指数eye 单位阵ezcontour 画等位线的简捷指令ezcontourf 画填色等位线的简捷指令ezgraph3 画表面图的通用简捷指令ezmesh 画网线图的简捷指令ezmeshc 画带等位线的网线图的简捷指令ezplot 画二维曲线的简捷指令ezplot3 画三维曲线的简捷指令ezpolar 画极坐标图的简捷指令ezsurf 画表面图的简捷指令ezsurfc 画带等位线的表面图的简捷指令f ffactor 符号计算的因式分解feather 羽毛图feedback 反馈连接feval 执行由串指定的函数fft 离散fourier变换fft2 二维离散fourier变换fftn 高维离散fourier变换fftshift 直流分量对中的谱fieldnames 构架域名figure 创建图形窗fill3 三维多边形填色图find 寻找非零元素下标findobj 寻找具有指定属性的对象图柄findstr 寻找短串的起始字符下标findsym 机器确定内存中的符号变量finverse 符号计算中求反函数fix 向零取整flag 红白蓝黑交错色图阵fliplr 矩阵的左右翻转flipud 矩阵的上下翻转flipdim 矩阵沿指定维翻转floor 向负无穷取整flops 浮点运算次数flow matlab提供的演示数据fmin 求单变量非线性函数极小值点（旧版）fminbnd 求单变量非线性函数极小值点fmins 单纯形法求多变量函数极小值点（旧版）fminunc 拟牛顿法求多变量函数极小值点fminsearch 单纯形法求多变量函数极小值点fnder 对样条函数求导fnint 利用样条函数求积分fnval 计算样条函数区间内任意一点的值fnplt 绘制样条函数图形fopen 打开外部文件for 构成for环用format 设置输出格式fourier fourier 变换fplot 返函绘图指令fprintf 设置显示格式fread 从文件读二进制数据fsolve 求多元函数的零点full 把稀疏矩阵转换为非稀疏阵funm 计算一般矩阵函数funtool 函数计算器图形用户界面fzero 求单变量非线性函数的零点g ggamma 函数gammainc 不完全函数gammaln 函数的对数gca 获得当前轴句柄gcbo 获得正执行"回调"的对象句柄gcf 获得当前图对象句柄gco 获得当前对象句柄geomean 几何平均值get 获知对象属性getfield 获知构架数组的域getframe 获取影片的帧画面ginput 从图形窗获取数据global 定义全局变量gplot 依图论法则画图gradient 近似梯度gray 黑白灰度grid 画分格线griddata 规则化数据和曲面拟合gtext 由鼠标放置注释文字guide 启动图形用户界面交互设计工具h hharmmean 调和平均值help 在线帮助helpwin 交互式在线帮助helpdesk 打开超文本形式用户指南hex2dec 十六进制转换为十进制hex2num 十六进制转换为浮点数hidden 透视和消隐开关hilb hilbert矩阵hist 频数计算或频数直方图histc 端点定位频数直方图histfit 带正态拟合的频数直方图hold 当前图上重画的切换开关horner 分解成嵌套形式hot 黑红黄白色图hsv 饱和色图i iif-else-elseif 条件分支结构ifft 离散fourier反变换ifft2 二维离散fourier反变换ifftn 高维离散fourier反变换ifftshift 直流分量对中的谱的反操作ifourier fourier反变换i, j 缺省的"虚单元"变量ilaplace laplace反变换imag 复数虚部image 显示图象imagesc 显示亮度图象imfinfo 获取图形文件信息imread 从文件读取图象imwrite 把图象写成文件ind2sub 单下标转变为多下标inf 无穷大info mathworks公司网点地址inline 构造内联函数对象inmem 列出内存中的函数名input 提示用户输入inputname 输入宗量名int 符号积分int2str 把整数数组转换为串数组interp1 一维插值interp2 二维插值interp3 三维插值interpn n维插值interpft 利用fft插值intro matlab自带的入门引导inv 求矩阵逆invhilb hilbert矩阵的准确逆ipermute 广义反转置isa 检测是否给定类的对象ischar 若是字符串则为真isequal 若两数组相同则为真isempty 若是空阵则为真isfinite 若全部元素都有限则为真isfield 若是构架域则为真isglobal 若是全局变量则为真ishandle 若是图形句柄则为真ishold 若当前图形处于保留状态则为真isieee 若计算机执行ieee规则则为真isinf 若是无穷数据则为真isletter 若是英文字母则为真islogical 若是逻辑数组则为真ismember 检查是否属于指定集isnan 若是非数则为真isnumeric 若是数值数组则为真isobject 若是对象则为真isprime 若是质数则为真isreal 若是实数则为真isspace 若是空格则为真issparse 若是稀疏矩阵则为真isstruct 若是构架则为真isstudent 若是matlab学生版则为真iztrans 符号计算z反变换j j , k kjacobian 符号计算中求jacobian 矩阵jet 蓝头红尾饱和色jordan 符号计算中获得jordan标准型keyboard 键盘获得控制权kron kronecker乘法规则产生的数组l llaplace laplace变换lasterr 显示最新出错信息lastwarn 显示最新警告信息leastsq 解非线性最小二乘问题（旧版）legend 图形图例lighting 照明模式line 创建线对象lines 采用plot 画线色linmod 获连续系统的线性化模型linmod2 获连续系统的线性化精良模型linspace 线性等分向量ln 矩阵自然对数load 从mat文件读取变量log 自然对数log10 常用对数log2 底为2的对数loglog 双对数刻度图形logm 矩阵对数logspace 对数分度向量lookfor 按关键字搜索m文件lower 转换为小写字母lsqnonlin 解非线性最小二乘问题lu lu分解m mmad 平均绝对值偏差magic 魔方阵maple &nb, sp; 运作maple格式指令mat2str 把数值数组转换成输入形态串数组material 材料反射模式max 找向量中最大元素mbuild 产生exe文件编译环境的预设置指令mcc 创建mex或exe文件的编译指令mean 求向量元素的平均值median 求中位数menuedit 启动设计用户菜单的交互式编辑工具mesh 网线图meshz 垂帘网线图meshgrid 产生"格点"矩阵methods 获知对指定类定义的所有方法函数mex 产生mex文件编译环境的预设置指令mfunlis 能被mfun计算的maple经典函数列表mhelp 引出maple的在线帮助min 找向量中最小元素mkdir 创建目录mkpp 逐段多项式数据的明晰化mod 模运算more 指令窗中内容的分页显示movie 放映影片动画moviein 影片帧画面的内存预置mtaylor 符号计算多变量taylor级数展开n nndims 求数组维数nan 非数（预定义）变量nargchk 输入宗量数验证nargin 函数输入宗量数nargout 函数输出宗量数ndgrid 产生高维格点矩阵newplot 准备新的缺省图、轴nextpow2 取最接近的较大2次幂nnz 矩阵的非零元素总数nonzeros 矩阵的非零元素norm 矩阵或向量范数normcdf 正态分布累计概率密度函数normest 估计矩阵2范数norminv 正态分布逆累计概率密度函数normpdf 正态分布概率密度函数normrnd 正态随机数发生器notebook 启动matlab和word的集成环境null 零空间num2str 把非整数数组转换为串numden 获取最小公分母和相应的分子表达式nzmax 指定存放非零元素所需内存o oode1 非stiff 微分方程变步长解算器ode15s stiff 微分方程变步长解算器ode23t 适度stiff 微分方程解算器ode23tb stiff 微分方程解算器ode45 非stiff 微分方程变步长解算器odefile ode 文件模板odeget 获知ode 选项设置参数odephas2 ode 输出函数的二维相平面图odephas3 ode 输出函数的三维相空间图odeplot ode 输出函数的时间轨迹图odeprint 在matlab指令窗显示结果odeset 创建或改写ode选项构架参数值ones 全1数组optimset 创建或改写优化泛函指令的选项参数值orient 设定图形的排放方式orth 值空间正交化p ppack 收集matlab内存碎块扩大内存pagedlg 调出图形排版对话框patch 创建块对象path 设置matlab搜索路径的指令pathtool 搜索路径管理器pause 暂停pcode 创建预解译p码文件pcolor 伪彩图peaks matlab提供的典型三维曲面permute 广义转置pi （预定义变量）圆周率pie 二维饼图pie3 三维饼图pink 粉红色图矩阵pinv 伪逆plot 平面线图plot3 三维线图plotmatrix 矩阵的散点图plotyy 双纵坐标图poissinv 泊松分布逆累计概率分布函数poissrnd 泊松分布随机数发生器pol2cart 极或柱坐标变为直角坐标polar 极坐标图poly 矩阵的特征多项式、根集对应的多项式poly2str 以习惯方式显示多项式poly2sym 双精度多项式系数转变为向量符号多项式polyder 多项式导数polyfit 数据的多项式拟合polyval 计算多项式的值polyvalm 计算矩阵多项式pow2 2的幂ppval 计算分段多项式pretty 以习惯方式显示符号表达式print 打印图形或simulink模型printsys 以习惯方式显示有理分式prism 光谱色图矩阵procread 向maple输送计算程序profile 函数文件性能评估器propedit 图形对象属性编辑器pwd 显示当前工作目录q qquad 低阶法计算数值积分quad8 高阶法计算数值积分(quadl)quit 推出matlab 环境quiver 二维方向箭头图quiver3 三维方向箭头图r rrand 产生均匀分布随机数randn 产生正态分布随机数randperm 随机置换向量range 样本极差rank 矩阵的秩rats 有理输出rcond 矩阵倒条件数估计real 复数的实部reallog 在实数域内计算自然对数realpow 在实数域内计算乘方realsqrt 在实数域内计算平方根realmax 最大正浮点数realmin 最小正浮点数rectangle 画"长方框"rem 求余数repmat 铺放模块数组reshape 改变数组维数、大小residue 部分分式展开return 返回ribbon 把二维曲线画成三维彩带图rmfield 删去构架的域roots 求多项式的根rose 数扇形图rot90 矩阵旋转90度rotate 指定的原点和方向旋转rotate3d 启动三维图形视角的交互设置功能round 向最近整数圆整rref 简化矩阵为梯形形式rsf2csf 实数块对角阵转为复数特征值对角阵rsums riemann和s ssave 把内存变量保存为文件scatter3 三维散点图sec 正割sech 双曲正割semilogx x轴对数刻度坐标图semilogy y轴对数刻度坐标图series 串联连接set 设置图形对象属性setfield 设置构架数组的域setstr 将ascii码转换为字符的旧版指令sign 根据符号取值函数signum 符号计算中的符号取值函数sim 运行simulink模型simget 获取simulink模型设置的仿真参数simple 寻找最短形式的符号解simplify 符号计算中进行简化操作simset 对simulink模型的仿真参数进行设置simulink 启动simulink模块库浏览器sin 正弦sinh 双曲正弦size 矩阵的大小slice 立体切片图solve 求代数方程的符号解spalloc 为非零元素配置内存sparse 创建稀疏矩阵spconvert 把外部数据转换为稀疏矩阵spdiags 稀疏对角阵spfun 求非零元素的函数值sph2cart 球坐标变为直角坐标sphere 产生球面spinmap 色图彩色的周期变化spline 样条插值spones 用1置换非零元素sprandsym 稀疏随机对称阵sprank 结构秩spring 紫黄调春色图sprintf 把格式数据写成串spy 画稀疏结构图sqrt 平方根sqrtm 方根矩阵squeeze 删去大小为1的"孤维"sscanf 按指定格式读串stairs 阶梯图std 标准差step 阶跃响应指令str2double 串转换为双精度值str2mat 创建多行串数组str2num 串转换为数strcat 接成长串strcmp 串比较strjust 串对齐strmatch 搜索指定串strncmp 串中前若干字符比较strrep 串替换strtok 寻找第一间隔符前的内容struct 创建构架数组struct2cell 把构架转换为元胞数组strvcat 创建多行串数组sub2ind 多下标转换为单下标subexpr 通过子表达式重写符号对象subplot 创建子图subs 符号计算中的符号变量置换subspace 两子空间夹角sum 元素和summer 绿黄调夏色图superiorto 设定优先级surf 三维着色表面图surface 创建面对象surfc 带等位线的表面图surfl 带光照的三维表面图surfnorm 空间表面的法线svd 奇异值分解svds 求指定的若干奇异值switch-case-otherwise 多分支结构sym2poly 符号多项式转变为双精度多项式系数向量symmmd 对称最小度排序symrcm 反向cuthill-mckee排序syms 创建多个符号对象t ttan 正切tanh 双曲正切taylortool 进行taylor逼近分析的交互界面text 文字注释tf 创建传递函数对象tic 启动计时器title 图名toc 关闭计时器trapz 梯形法数值积分treelayout 展开树、林treeplot 画树图tril 下三角阵trim 求系统平衡点trimesh 不规则格点网线图trisurf 不规则格点表面图triu 上三角阵try-catch 控制流中的try-catch结构type 显示m 文件u uuicontextmenu 创建现场菜单uicontrol 创建用户控件uimenu 创建用户菜单unmkpp 逐段多项式数据的反明晰化unwrap 自然态相角upper 转换为大写字母v vvar 方差varargin 变长度输入宗量varargout 变长度输出宗量vectorize 使串表达式或内联函数适于数组运算ver 版本信息的获取view 三维图形的视角控制voronoi voronoi多边形vpa 任意精度（符号类）数值w wwarning 显示警告信息what 列出当前目录上的文件whatsnew 显示matlab中readme文件的内容which 确定函数、文件的位置while 控制流中的while环结构white 全白色图矩阵whitebg 指定轴的背景色who 列出内存中的变量名whos 列出内存中变量的详细信息winter 蓝绿调冬色图workspace 启动内存浏览器x x , y y , z zxlabel x轴名xor 或非逻辑yesinput 智能输入指令ylabel y轴名zeros 全零数组zlabel z轴名zoom 图形的变焦放大和缩小ztrans 符号计算z变换。

图像插值技术——双线性插值法在图像处理中，如果需要对图像进⾏缩放，⼀般可以采取插值法，最常⽤的就是双线性插值法。

本⽂⾸先从数学⾓度推导了⼀维线性插值和⼆维线性插值的计算过程，并总结了规律。

随后将其应⽤到图像的双线性插值上，利⽤Matlab编程进⾏图像的缩放验证，实验证明，⼆维线性插值能够对图像做出较好的缩放效果。

数学⾓度的线性插值⼀维线性插值假设有⼀个⼀元函数 y=f(x) , 已知曲线上的两点，A 和 B 的坐标分别为 (x0,y0) 、(x1,y1) 。

现在要在A 和 B 之间通过插值计算出⼀个点 P ,若已知 P点的横坐标 x，如何求出 P点的纵坐标 y ?这⾥我们的插值之所以叫做线性插值，就是因为我们假定了 P 点落在 A 点和 B 点的连线上,使得他们的坐标之间满⾜线性关系。

所以，根据初中的知识，可以得到下⾯的等式：y−y0 y1−y0=x−x0 x1−x0这⾥我们令：α=x−x0 x1−x0于是，我们可以得到P点的纵坐标 y 的表达式：y=(1−α)f(x0)+αf(x1)⼆维线性插值⼀维线性插值可以扩展到⼆维的情况。

假设有⼀个⼆元函数 z=f(x,y) , 已知曲⾯上的四点，A 、B 、C、D的坐标分别为 (x0,y0) 、(x1,y0) 、(x1,y1)、(x0,y1) 。

现在要在A 、B 、C、D之间通过插值计算出⼀个点 P ,若已知 P点的坐标 (x,y)，如何求出 P点的函数值坐标 z ?这⾥我们依旧可以仿照⼀维线性插值，进⾏计算。

假设先计算 y 轴⽅向的插值点 P0 和 P1 ,则根据上⾯的推导过程，且令α=y−y0 y1−y0则， P0 的取值 z0为：z0=(1−α)f(x0,y0)+αf(x0,y1) P1 的取值 z1为：z1=(1−α)f(x1,y0)+αf(x1,y1)再计算 x 轴⽅向的插值点 P，令β=x−x0 x1−x0则 P 的取值 z为：z=(1−β)z0+βz1整理得到下⾯的式⼦：z =(1−β)(1−α)f x 0,y 0+αf x 0,y 1+β(1−α)f x 1,y 0+αf x 1,y 1=(1−β)(1−α)f x 0,y 0+(1−β)αf x 0,y 1+β(1−α)f x 1,y 0+βαf x 1,y 1⼩结由⼀维线性插值过渡到⼆维线性插值，我们发现，⼆者在表达式上有相似的规律：⼀维线性插值：y =f (x )α=x p −x 0x 1−x 0y p =(1−α)f x 0+αf x 1⼆维线性插值：z =f (x ,y )α=x p −x 0x 1−x 0,β=y p −y 0y 1−y 0z p =(1−β)(1−α)f x 0,y 0+(1−β)αf x 0,y 1+β(1−α)f x 1,y 0+βαf x 1,y 1图像中的双线性插值我们可以⽤函数来表⽰⼀幅图像（假设为单通道）。

matlab中二维插值函数用法Matlab是一种强大的数学软件，可以对数据进行各种操作和分析。

在Matlab 中，二维插值函数是一种非常常用的工具，它可以帮助我们在已知数据点的基础上，估计其它位置的数值。

在本文中，我们将一步一步介绍二维插值函数的用法，希望对你有所帮助。

首先，我们需要明确一点，二维插值是用于处理二维平面上的数据，即x和y坐标都是变化的。

插值的目的是在这样的二维数据上，估计其它位置的数值，使得整个平面上的数值分布更加密集。

在Matlab中，二维插值函数主要是interp2函数。

该函数的基本语法如下：Vq = interp2(X, Y, V, Xq, Yq, method)其中，X和Y是已知数据点的x和y坐标，V是对应的数值；Xq和Yq是我们要估计数值的位置，Vq是估计得到的数值；method是插值方法，可以是‘linear’、‘nearest’、‘spline’等。

接下来，我们来演示一个具体的例子，以说明interp2函数的用法。

假设我们有一组二维数据如下：X = [1, 2, 3; 1, 2, 3; 1, 2, 3]Y = [1, 1, 1; 2, 2, 2; 3, 3, 3]V = [0, 1, 0; 1, 4, 1; 0, 1, 0]这里，X和Y分别是数据点的x和y坐标，V是对应的数值。

现在我们想要在位置(1.5, 1.5)处进行插值，估计其数值。

我们可以使用interp2函数来完成这个任务：Vq = interp2(X, Y, V, 1.5, 1.5, 'linear')这里，我们使用了‘linear’插值方法，得到的估计值Vq为2。

这就是在位置(1.5, 1.5)处的估计数值。

除了‘linear’插值方法，我们还可以尝试其它方法，比如‘nearest’和‘spline’。

这些方法在实际应用中有着不同的特点，需要根据具体情况选择合适的方法。

例如，‘nearest’方法会选择最近的数据点的数值作为估计值，‘spline’方法则会使用样条插值进行估计。

MATLAB 中的曲线拟合和插值在大量的使用领域中，人们经常面临用一个分析函数描述数据(通常是测量值)的任务。

对这个问题有两种方法。

在插值法里，数据假定是正确的，要求以某种方法描述数据点之间所发生的情况。

这种方法在下一节讨论。

这里讨论的方法是曲线拟合或回归。

人们设法找出某条光滑曲线，它最佳地拟合数据，但不必要经过任何数据点。

图11.1说明了这两种方法。

标有'o'的是数据点；连接数据点的实线描绘了线性内插，虚线是数据的最佳拟合。

11.1 曲线拟合曲线拟合涉及回答两个基本问题：最佳拟合意味着什么？应该用什么样的曲线？可用许多不同的方法定义最佳拟合，并存在无穷数目的曲线。

所以，从这里开始，我们走向何方？正如它证实的那样，当最佳拟合被解释为在数据点的最小误差平方和，且所用的曲线限定为多项式时，那么曲线拟合是相当简捷的。

数学上，称为多项式的最小二乘曲线拟合。

如果这种描述使你混淆，再研究图11.1。

虚线和标志的数据点之间的垂直距离是在该点的误差。

对各数据点距离求平方，并把平方距离全加起来，就是误差平方和。

这条虚线是使误差平方和尽可能小的曲线，即是最佳拟合。

最小二乘这个术语仅仅是使误差平方和最小00.20.40.60.81-2024681012xy =f (x )Second O rder C urv e Fitting图11.1 2阶曲线拟合在MATLAB 中，函数polyfit 求解最小二乘曲线拟合问题。

为了阐述这个函数的用法，让我们以上面图11.1中的数据开始。

» x=[0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1]; » y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];为了用polyfit ，我们必须给函数赋予上面的数据和我们希望最佳拟合数据的多项式的阶次或度。

如果我们选择n=1作为阶次，得到最简单的线性近似。

matlab中二维插值函数Matlab中的二维插值函数是一种用于在给定有限数据点的情况下估算未知数据点的方法。

它可以用于图像处理、数据分析、数值计算等多个领域。

本文将介绍Matlab中常用的二维插值函数及其应用。

Matlab中的二维插值函数包括interp2、griddata和scatteredInterpolant等。

这些函数可以根据给定的数据点和相应的值，生成一个二维插值函数对象。

通过该对象，可以对未知数据点进行插值计算。

其中，interp2函数是最常用的二维插值函数之一。

它可以根据给定的数据点和相应的值，通过线性插值、样条插值或三次样条插值等方法，计算未知数据点的值。

interp2函数的使用方法如下：```Vq = interp2(X,Y,V,Xq,Yq)```其中，X和Y是二维网格的坐标，V是对应网格点的值，Xq和Yq 是待插值的点的坐标，Vq是插值得到的值。

通过interp2函数，可以将离散的数据点插值为连续的曲面，实现对未知数据点的估算。

除了interp2函数，griddata函数也是常用的二维插值函数之一。

它可以根据给定的散点数据和相应的值，生成一个二维插值函数对象。

griddata函数的使用方法如下：```F = griddata(x,y,v,xq,yq)```其中，x和y是散点数据的坐标，v是对应散点数据的值，xq和yq 是待插值的点的坐标，F是插值得到的函数对象。

通过griddata函数，可以实现对离散数据点的插值计算。

scatteredInterpolant函数也是常用的二维插值函数之一。

它可以根据给定的散点数据和相应的值，生成一个二维插值函数对象。

scatteredInterpolant函数的使用方法如下：```F = scatteredInterpolant(x,y,v)```其中，x和y是散点数据的坐标，v是对应散点数据的值，F是插值得到的函数对象。

通过scatteredInterpolant函数，可以实现对散点数据的插值计算。

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函数对已知的数据点进⾏插值，数据点（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）。

二维样条插值函数和二维牛顿插值二维样条插值函数和二维牛顿插值是在二维平面上进行插值的两种常见方法。

它们在科学计算、图像处理、地理信息系统等领域都有广泛的应用。

本文将对二维样条插值函数和二维牛顿插值进行详细介绍，并对它们的使用方法和特点进行比较。

首先，我们先来了解二维样条插值函数。

样条插值是一种光滑插值方法，它通过将插值函数拆分成多个小段，每个小段被称为样条。

二维样条插值函数是在二维平面上进行插值时使用的方法。

它的主要思想是将数据点之间的曲线进行拟合，得到一个平滑的曲线。

二维样条插值函数的插值曲线不仅能通过给定的数据点，还能将曲线的一阶和二阶导数保持连续。

这保证了插值曲线的光滑性和拟合性能。

在二维样条插值函数中，常用的有三次样条插值和二次样条插值。

三次样条插值利用三次多项式求解，保证了曲线的光滑性和连续性。

而二次样条插值只利用二次多项式，比三次样条插值的计算量要小，但相应地也降低了插值曲线的平滑性。

使用哪种样条插值函数要根据具体的应用场景和需求来选择。

接下来，我们来介绍二维牛顿插值。

二维牛顿插值是利用多项式进行插值的一种方法。

它通过构造一个多项式函数来拟合给定的数据点。

多项式的次数取决于数据点的个数。

二维牛顿插值的关键在于寻找插值节点和系数。

在二维牛顿插值中，插值节点的选取非常重要。

合理选择插值节点可以提高插值函数的拟合性能。

常见的选择方法有均匀节点选择和非均匀节点选择。

均匀节点选择是将二维平面均匀地划分成若干个小区域，然后在每个小区域内选择一个节点作为插值节点。

非均匀节点选择则根据数据点的密集程度来灵活选择节点，从而提高插值函数的拟合精度。

二维牛顿插值可以通过Newton的插值公式推导得到插值函数。

利用插值节点和数据点的差商来表示多项式的系数，从而得到插值函数。

同时，牛顿插值函数的计算速度相对较快，不需要拟合样条的复杂过程。

综上所述，二维样条插值函数和二维牛顿插值在二维平面上进行插值都有各自的特点和应用场景。

matlab中二维插值函数用法摘要：一、Matlab中二维插值函数简介1.插值函数的作用2.Matlab中的二维插值函数二、二维插值函数的用法1.函数语法2.函数参数3.示例三、二维插值函数的应用1.图像处理2.数据插补3.数据预测正文：Matlab中二维插值函数用法Matlab是一款功能强大的科学计算软件，在工程和科学研究中应用广泛。

其中，二维插值函数是Matlab中的一个重要工具，用于实现二维数据的插值和拟合。

本文将介绍Matlab中二维插值函数的用法及其在图像处理、数据插补和数据预测等方面的应用。

一、Matlab中二维插值函数简介插值函数是一种数学方法，用于在已知数据点之间计算新数据点的值。

Matlab提供了多种二维插值函数，包括最邻近插值（nearest）、线性插值（linear）和立方插值（cubic）等。

这些插值函数可以用于图像处理、数据分析和科学计算等领域。

二、二维插值函数的用法1.函数语法在Matlab中，二维插值函数的语法如下：`Z = interp2(X, Y, Z, Xnew, Ynew)`其中，`X`、`Y`和`Z`分别为输入数据的x、y和z值，`Xnew`和`Ynew`为需要插值的新数据点的x、y坐标。

函数返回插值后的新数据点`Znew`。

2.函数参数- `X`：输入数据点的x坐标（1xN向量）- `Y`：输入数据点的y坐标（1xN向量）- `Z`：输入数据点的z坐标（1xN向量）- `Xnew`：新数据点的x坐标（1xM向量）- `Ynew`：新数据点的y坐标（1xM向量）3.示例假设我们有一个数据点集，如下所示：```matlabX = [1, 2, 3, 4, 5];Y = [1, 4, 9, 16, 25];Z = [2, 5, 8, 11, 14];```现在，我们想要在新数据点集`Xnew`和`Ynew`上进行插值：```matlabXnew = [1.5, 2.5, 3.5, 4.5];Ynew = [1.5, 2.5, 3.5, 4.5];```可以使用以下代码进行插值：```matlabZnew = interp2(X, Y, Z, Xnew, Ynew);```插值后的新数据点集`Znew`如下：```matlabZnew =2.50005.00008.000011.000014.0000```三、二维插值函数的应用1.图像处理在图像处理中，二维插值函数可以用于插值和缩放图像。