Matlab interp2的2D样条外推原理是什么？如何手动实现该功能？

我了解已有很多关于Matlab中样条外推的讨论，我有一个2D场景下该功能表现极佳的示例，希望理解其运行逻辑后，基于Math.Net在C#中实现该能力，以下是我的示例：

matrix = zeros(128,128);
barWidth = 10;

%Calculate Midpoints
midpointX = floor(size(matrix,2)/2) + 1;
midpointY = floor(size(matrix,1)/2) + 1;

matrix(midpointY-barWidth:midpointY+barWidth,:) = 1;

% Windowing (I know it could be shortened...)
distanceMatY = -midpointY+1:midpointY-2;
distanceMatY = abs(repmat(distanceMatY',1,size(matrix,2)));
Factor = 0.5*(cos(pi*distanceMatY/barWidth)+1);
index = distanceMatY > barWidth;
Factor(index) = 0;
matrix = matrix.*Factor;

% Rotate matrix
alpha = 30 * pi/180;
y = -midpointY+1:midpointY-2; y = y'; 
x = -midpointX+1:midpointX-1;     %set new COS for rotation in Midpoint
xRot = x*cos(alpha) - y*sin(alpha) + nMP;  %Determine X and Y matrices for
yRot = x*sin(alpha) + y*cos(alpha) + mMP;  %Roatation with angle alpha

rotMatrix = interp2(matrix,xRot,yRot,'linear');      %Interpolate rotated matrix

我构建了一个128×128的零值矩阵，中间位置有一条横跨全宽的条形。

我对条形施加了汉宁窗（Hann Window）平滑边缘，之后使用双线性插值将矩阵旋转30度。

此时边界外的取值被设为NaN，我可以将所有NaN设为0，但我实际需要的是让条形自动延伸。我也可以选择在旋转前对矩阵做填充，旋转后再裁剪回原输入尺寸。
更简单的方案是直接调用interp2(matrix,xRot,yRot,'spline')，此时任意旋转角度下输出的旋转矩阵都完全符合我的预期。

请问2D样条插值的外推逻辑是如何实现的？有没有方法可以手动实现该功能？

解答 #

2D样条插值外推核心逻辑 #

Matlab interp2的spline模式采用可分离式三次样条实现，外推能力来自三次多项式的特性：

• 插值前会先对原始网格的所有点拟合分段三次多项式，要求相邻分段在拼接点处二阶导数连续
• 对于超出原始网格范围的坐标，直接复用最靠近边界的那一段三次多项式进行求值，不需要额外的截断逻辑，所以边界处的平滑趋势会自然向外延伸，刚好适配你旋转后条形延展的需求
• 双线性插值默认返回NaN的原因是它没有预拟合的多项式，仅依赖相邻四个采样点做加权计算，边界外缺少采样点就无法生成有效值

C# 基于Math.Net的实现方案 #

你可以按以下步骤实现完全对齐Matlab逻辑的2D样条外推能力：

1. 先安装Math.Net.Numerics依赖包，它内置了成熟的一维三次样条实现
2. 2D样条插值采用可分离实现：先按行做一维样条拟合，再按列对行插值的结果做一维样条拟合，两次一维样条拼接就是完整的2D样条能力
3. 外推不需要额外处理，Math.Net的三次样条默认支持边界外的多项式求值，直接传入超出原始网格范围的坐标即可得到外推值

参考实现代码 #

using MathNet.Numerics.Interpolation;
using MathNet.Numerics.LinearAlgebra;

// 初始化原始128x128矩阵，省略汉宁窗赋值逻辑，和你Matlab的逻辑对齐即可
var rawMatrix = Matrix<double>.Build.Dense(128, 128, 0);
// 原始网格坐标
var xGrid = Enumerable.Range(0, 128).Select(v => (double)v).ToArray();
var yGrid = Enumerable.Range(0, 128).Select(v => (double)v).ToArray();

// 预拟合所有行的三次样条
var rowSplines = new IInterpolation[128];
for (int y = 0; y < 128; y++)
{
    rowSplines[y] = CubicSpline.InterpolateNatural(xGrid, rawMatrix.Row(y).ToArray());
}

// 2D样条求值方法，支持x、y超出0~127的外推场景
double GetSplineValue(double x, double y)
{
    // 先对所有行在x坐标处插值，得到一列中间值
    double[] colValues = rowSplines.Select(s => s.Interpolate(x)).ToArray();
    // 对中间值做列方向的样条插值，得到最终结果
    var colSpline = CubicSpline.InterpolateNatural(yGrid, colValues);
    return colSpline.Interpolate(y);
}

// 生成旋转后的xRot、yRot坐标矩阵后，遍历每个点调用GetSplineValue即可得到插值结果

如果需要更高性能，可以提前预计算列方向的样条参数，避免每次求值都重新拟合列样条。

内容的提问来源于stack exchange，提问作者Till