MATLAB插值函数浅析：griddata（‘v4’）

  griddata是MATLAB里的一个功能强大的二维及三维插值函数，它提供了以下几种插值式。             

        据表知 ，v4插值法较为独特，根据双调和样条理论进行运算的，而且插值结果极为平滑，所做出来的图相当美观（当然在一些初始点较少的区域内会丧失一定的准确性）。

        我通过解析源码的方式，对v4插值法从编码方面进行分析，下面是MATLAB给出的v4相关的源码及注释。

function vq = gdatav4(x,y,v,xq,yq)
%GDATAV4 MATLAB 4 GRIDDATA interpolation
%   Reference:  David T. Sandwell, Biharmonic spline
%   interpolation of GEOS-3 and SEASAT altimeter
%   data, Geophysical Research Letters, 2, 139-142,
%   1987.  Describes interpolation using value or
%   gradient of value in any dimension.

[x, y, v] = mergepoints2D(x,y,v);

xy = x(:) + 1i*y(:);

% Determine distances between points
d = abs(xy - xy.');

% Determine weights for interpolation
g = (d.^2) .* (log(d)-1);   % Green's function.
% Fixup value of Green's function along diagonal
g(1:size(d,1)+1:end) = 0;
weights = g \ v(:);

[m,n] = size(xq);
vq = zeros(size(xq));
xy = xy.';

% Evaluate at requested points (xq,yq).  Loop to save memory.
for i=1:m
    for j=1:n
        d = abs(xq(i,j) + 1i*yq(i,j) - xy);
        g = (d.^2) .* (log(d)-1);   % Green's function.
        % Value of Green's function at zero
        g(d==0) = 0;
        vq(i,j) = g * weights;        
    end
end

        这段代码写得很好，接下来我们进行逐行分析：

        第一行导入参数，(x , y , v)是已知值及其坐标，（xq , yq）是所需要插值的坐标，vq是返回值。之后的注释是MATLAB开发者所参考的文献。

        第二行代码用到的mergepoints2D函数是开发者内部函数，并没有开放给普通用户使用，我会附在后面，它主要做的工作大体有：

        1、接受三个值，可以理解为x,y以及（x,y)对应的z，并且用reshape函数将x,y,z改写成n*1形式的列向量。

        2、将所接受的值分为实数与虚数，将(x,y)于对应的v合并为新矩阵xyv，格式是（y,x,v)。

        之后的一段代码以一种很简洁的方式求解了每个点之间的欧几里得距离（其实就是距离），代码及解析如下：

xy = x(:) + 1i*y(:);
d = abs(xy - xy.');

        1、虚数化处理，将每一个点（x,y）写成x + y*i格式，命名为xy,注意由于x和y的数量，此时xy是个向量。

        2、matlab里行向量可以减去列向量，其值为，每一行减每一列，所以如果xy是100组数据那么我们现在有了10000组数据，是所有点之间的间距。我想以图论的观点来看，d是一个无向加权邻接矩阵，邻接矩阵是线路图的一种存储方式。

        之后他使用格林函数求解权重：

d = abs(xy - xy.');
g = (d.^2) .* (log(d)-1);   % Green's function.
g(1:size(d,1)+1:end) = 0;
weights = g \ v(:);

        

        1、在matlab里里面对于矩阵X，exp(X)与log(X)都是对每一个元素进行处理

        2、那么对于d里面的任一个值，我们假定它是x，它是点到点的距离，它会经历(x^2)*[log(x) - 1]。

        3、由于在第二步变换中，d矩阵的主对角线会被破坏，所以第三步，便将主对角线回归为0，这里是为了防止在解线性方程组时出现奇异矩阵问题（即避免自相关）。

        4、求解权重就类似于是在求解一个n元一次方程组。

        下一段较为通俗，给出结果矩阵的大小，以及xy虚数矩阵的的非共轭转置。

[m,n] = size(xq);
vq = zeros(size(xq));
xy = xy.';% 虚数序列转置

        这里xy变成了一个横坐标，这是为了之后循环里求解vq做准备。

        之后经过大遍历，便能求解出所有要求的vq，这也体现了这个插值法的弱点，若将（xq,yq）设置得过于离谱，v4依旧会给出插值结果，而其他插值法普遍不会。这里给出我常用的（xq,yq）是(linspace(min(X), max(Y))', linspace(min(X), max(Y)))。

        下面给出mergepoints2D函数源码及部分注释。

function [x, y, v] = mergepoints2D(x,y,v)
% 排序
% Sort x and y so duplicate points can be averaged
% 输入值x,y,z需要是列向量，其实我看后面的代码，即使是矩阵应该也没有关系，reshape函数会将矩阵褛乘向量
%Need x,y and z to be column vectors
sz = numel(x);%求x数组个数
x = reshape(x,sz,1);
y = reshape(y,sz,1);
v = reshape(v,sz,1);%规范化
myepsx = eps(0.5 * (max(x) - min(x)))^(1/3);%这是在计算容差
myepsy = eps(0.5 * (max(y) - min(y)))^(1/3);
%eps是浮点数最小精度范围，为2的负51次方，应用在第60行

% look for x, y points that are indentical (within a tolerance)
% average out the values for these points
if isreal(v)% 如果v是实数
    xyv = matlab.internal.math.mergesimpts([y, x, v], [myepsy, myepsx, Inf], 'average');%这个均值求的是极小范围的
    x = xyv(:,2);
    y = xyv(:,1);
    v = xyv(:,3);%x,y,z建立联系，写入矩阵xyv，写入格式（y，x,v）;
else
    % if z is imaginary split out the real and imaginary parts
    xyv = matlab.internal.math.mergesimpts([y, x, real(v), imag(v)], ...%之后是虚数部分，他将实部和虚部分开求，我没看
        [myepsy, myepsx, Inf, Inf], 'average');
    x = xyv(:,2);
    y = xyv(:,1);
    % re-combine the real and imaginary parts
    v = xyv(:,3) + 1i*xyv(:,4);
end
% give a warning if some of the points were duplicates (and averaged out)
if sz>numel(x)
    warning(message('MATLAB:griddata:DuplicateDataPoints'));
end

        另外给出模仿此函数写的python脚本（好久没有动python了，都不知道python有什么关于插值的第三方库），代码经过python3.9上运行成功。

import numpy as np  
  
def mergepoints2D(x, y, v):  
    points = np.column_stack((y, x, v))  
    return points  
  
def gdatav4(x, y, v, xq, yq):  

    points = mergepoints2D(x, y, v)  
      
    xy = points[:, 1] + 1j * points[:, 0]    
      
    d = np.abs(xy[:, np.newaxis] - xy[np.newaxis, :])   
    g = (d**2) * (np.log(d) - 1)  
    np.fill_diagonal(g, 0)  
    weights = np.linalg.solve(g, v)  
    vq = np.zeros((xq.shape[0], xq.shape[1]))  
      
    for i in range(xq.shape[0]):  
        for j in range(xq.shape[1]):  
            d_q = np.abs(xq[i, j] + 1j * yq[i, j] - xy)  
            g_q = (d_q**2) * (np.log(d_q) - 1)  
            g_q[d_q == 0] = 0   
            vq[i, j] = np.dot(g_q, weights)  
      
    return vq