将 1D 高斯插值到 2D 高斯答案

【问题标题】：Interpolating 1D Gaussian into 2D Gaussian将 1D 高斯插值到 2D 高斯
【发布时间】：2010-03-14 17:32:13
【问题描述】：

假设我有一个一维高斯函数。它的长度是600。

我想将它内插成 600 X 600 大小的 2D 高斯。

这是我写的代码（OTFx 是高斯函数，OTF - 2d Interpolated Function）：

[x, y] = meshgrid([-300:299], [-300:299]);
r = sqrt((x .^ 2) + (y .^ 2));

OTF = interp1([-300:299], OTFx, r(:), 'spline');
OTF = reshape(OTF, [600, 600]);

问题是我最后得到了 Overshoot：

如何防止这种过冲？单调递减函数有更好的插值算法吗？

注意：我正在寻找一种通用解决方案，用于将 1D 函数内插到 2D 径向对称函数中，高斯只是一个示例。

【问题讨论】：

标签： matlab interpolation

【解决方案1】：

编辑：根据您的说明，很清楚发生了什么。您正在尝试对超出可用数据范围的函数进行插值 - 即您正在从插值到外推。样条曲线将导致您观察到的过冲。解决方案只是确保您的一维函数具有区间 [min(r), max(r)] 中的值。请注意，在原始数据中，max(r) 约为 424，而您要插值的函数定义在 [-300,299] 范围内

% Simulated overshoot, see left figure:
x1d = [-300:299];
[x,y]=meshgrid(x1d,x1d);
r = sqrt(x.^2+y.^2);
gsn1d = exp(-x1d.^2/500);
lowpass = @(x)(x1d > -x & x1d < x);
gsn1dcutoff = ifft(fftshift(lowpass(10).*fftshift(fft(gsn1d))));
plot(gsn1dcutoff)
OTF2d = reshape(interp1(x1d,gsn1dcutoff,r(:),'spline'),[length(x1d),length(x1d)]);
mesh(OTF2d)

% Quick and dirty fix, see right figure:
x1dExtended = linspace(min(x1d*sqrt(2)),max(x1d*sqrt(2)),ceil(length(x1d)*sqrt(2)));
gsn1dE = exp(-x1dExtended.^2/500);
% ^^^ note that this has 600*sqrt(2) points and is defined on the diagonal of your square.   Now we can low-pass filter in the freq. domain to add ripple in space domain:
lowpass = @(x)(x1dExtended > -x & x1dExtended < x);
gsn1dcutoff = -real(ifft(fftshift(lowpass(10).*fftshift(fft(gsn1dE)))));
plot(gsn1dcutoff)
OTF2d = reshape(interp1(x1dExtended,gsn1dcutoff,r(:),'spline'),[length(x1d),length(x1d)]);
mesh(OTF2d)

alt text http://img54.imageshack.us/img54/8255/clipboard01vz.png

【讨论】：

好吧，实际上 OTFx 是从图像估计的，但它与常数和 epsilon 与“高斯”相同。尝试以下方法来复制它。创建一个长度为 21 的高斯。使用带有零填充的 fft 将其转换为长度 600。尝试对其进行插值以具有 2D 高斯，就像将 2D 高斯转换为频域一样。问题是，有没有其他方法可以做到这一点？也许是一个更好的算法，应用一些约束等......这类似于“振铃”效果。谢谢。
查看我的更新答案。我认为其他插值算法不会更好地推断可用数据，但请随时查找 interp1 的帮助并使用不同的算法。要点是：为避免意外，请确保为插值提供足够的数据。
感谢您的回复。也许如果我对这个问题有更多的了解，你可能会为我提供更多帮助。插值函数必须处处为正。通常它应该是单调降序。我从类似于长度为 21 的一维高斯的东西开始，然后将其“FFT”成 600（或者我想要使用零填充的任何长度）。我想要的是频域中最好的 600X600 2D 插值。作为测试用例，我们可以采用长度为 21 的“纯”一维高斯，并尝试将其插值尽可能接近 21X21 二维高斯的 600X600 二维 DFT。谢谢。

【解决方案2】：

Leo 的诊断是正确的。我想建议一个更简单（我希望）的补救措施：做你想做的事（基本上是围绕其对称轴旋转高斯）并在 600x600 正方形中获得合理的答案，你需要一个高斯 600*sqrt (2)=849 像素长。如果你能做到这一点，那么所有进一步的 thttp://stackoverflow.com/questions/2443046/interpolating-1d-gaussian-into-2d-gaussianrickery 就没有必要了。

编辑：换句话说：如果你围绕它的中心旋转 600 像素宽的东西，你会得到一个直径 600 像素的圆。您想要覆盖一个 600x600 正方形。为此，您需要一个直径为 849 像素的圆，因为这是正方形的对角线。

【讨论】：

事实上，这正是我的“快速而肮脏的修复”正在做的事情；它可能需要更明确，因为修复与模拟截止的旧代码混合在一起。感谢您指出这一点。
您基本上建议的是，如果您想要一个具有极对称的 600X600 插值，则插值长度为 849 并且仅裁剪 600X600？谢谢。
@Drazick：可能不会。关键是：如果你围绕它的中心旋转 600 像素宽的东西，你会得到一个直径为 600 像素的圆。你想覆盖一个 600x600 的正方形。为此，您需要一个直径为 849 像素的圆，因为这是正方形的对角线。为了后代，我已将其添加到答案的正文中；-)
我们的想法是尽可能将高斯函数的 2D DFT 填充到 600X600（最初是 21X21）。你会怎么做，将 1D 的 DFT 转换为 849，进行旋转然后裁剪或忍受“过冲”？谢谢。
我怎样才能至少创造一个积极的约束。谢谢。

【解决方案3】：

在高斯的特殊情况下，您可以使用可分离的事实来计算高斯：

OTF2(x,y) = exp( - x^2 - y^2) = exp( - x^2) * exp( - y^2) = OTFx(x) * OTFx(y)

所以您仍然只需要在内存中存储 OTFx。

【讨论】：

我以高斯为例。真正的问题是将一维函数内插到二维径向对称函数中。
好的，我明白了...您应该相应地编辑您的问题。过冲似乎来自镜像 - matlav 的 interp 中是否有一个选项？
哇...有一段时间我在 matlab 中编码 :-)