將一維高斯插值爲二維高斯

我想把它插補成大小600×600。

的2D高斯這是我寫的代碼（OTFx是高斯函數，OTF - 2D內插功能）：

[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]);

問題是我得到最後的過沖： alt text http://i39.tinypic.com/259et0g.png

我該如何防止它？單調遞減函數有更好的插值算法嗎？

來源

2010-03-14 Royi

編輯：根據你的澄清，很清楚發生了什麼事情。您正嘗試在可用數據的範圍之外插入函數 - 即您正在從插值到外推。樣條曲線將導致您正在觀察的過沖。解決方法是簡單地確保您的一維函數具有區間[min（r），max（r）]中的值。請注意，在原始數據，最大（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

來源

2010-03-14 19:02:11

好的，實際上OTFx是從圖像中估計出來的，然而它與恆定和epsilon到「高斯」是完全相同的。請嘗試以下操作來複制它。創建一個長度爲21的高斯。使用帶有零填充的fft將其轉換爲長度600.嘗試將其插值爲具有二維高斯，因爲您將二維高斯轉換爲頻域。問題是，有沒有其他的方法呢？也許一個更好的算法，應用一些約束等... 這是類似於「振鈴」效果。謝謝。 – Royi

查看我更新的答案。我不認爲其他插值算法在推斷超出可用數據的情況下會更好，但可以隨時查找interp1的幫助並使用不同的算法進行操作。重點是：爲了避免意外，請確保您提供了足夠的插值數據。 –

感謝您的回覆。也許如果我對這個問題有更多的瞭解，你可以幫助我更多。內插函數在任何地方都必須是正數。一般來說它應該是單調下降的。我開始的東西類似於一維高斯的長度21然後「FFT」到600（或任何長度，我想要使用零填充）。我想要的是在頻域中最好的600X600 2D插值。作爲測試案例，我們可以採用長度爲21的「純」1維高斯，並嘗試將其內插到儘可能接近21X21 2D高斯的600X600 2D DFT。謝謝。 – Royi

獅子座是正確的，他的診斷定義。我想提出一個更簡單的（我希望）補救措施：做你想做的事（基本上圍繞它的對稱軸旋轉高斯），並在600x600的正方形中得到合理的答案，你需要一個高斯600 * sqrt （2）= 849個像素長。如果你能做到這一點，那麼所有進一步的thttp：//stackoverflow.com/questions/2443046/interpolating-1d-gaussian-into-2d-gaussianrickery是沒有必要的。

編輯：換句話說：如果在中心周圍旋轉600像素寬的東西，則會得到一個直徑爲600像素的圓。您想要覆蓋600x600 方形。爲此，您需要一個直徑爲849像素的圓，因爲這是正方形的對角線。

來源

2010-03-14 23:22:13 AVB

事實上，這正是我的「快速和骯髒的修復」正在做的事情;它可能需要更加明確，因爲修復與模擬截斷的舊代碼混合在一起。感謝您指出了這一點。 –

你基本上建議的是，如果你想要一個帶有極座標對稱的600X600插值，插值長度爲849，只裁剪600X600？謝謝。 – Royi

@Drazick：可能不是。重點是：如果在中心旋轉600像素寬的東西，則會得到一個直徑爲600像素的圓。你想要覆蓋600x600平方米。爲此，您需要一個直徑爲849像素的圓，因爲這是正方形的對角線。我已經添加了這個答案的身體，爲後人;-) – AVB

高斯的特定情況下，你可以使用的事實，它是可分離計算高斯：

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

所以你還是隻需要只OTFx存儲在內存中。

來源

2010-03-15 09:09:54 meduz

我以高斯爲例。真正的問題是將一維函數插入二維徑向對稱函數。 – Royi

好吧，我明白了......你應該相應地編輯你的問題。超調似乎來自鏡像 - 在matlav的interp中有沒有一個選項？ – meduz

將一維高斯插值爲二維高斯

回答

相關問題