Interpolating 1D Gaussian into 2D Gaussian

Question

Let's say I have a 1D Gaussian function. Its length is 600 for that matter.

I want to Interpolate it into 2D Gaussian of the size 600 X 600.

This is the code I wrote (OTFx is the Gaussian Function, 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]);

The problem is I get Overshoot at the end:

How can I prevent it? Is there better Interpolating algorithm for Monotonic Descending Functions?

Answer 1

EDIT: based on your clarification, it's clear what's going on. You are trying to interpolate a function beyond the range of available data -- i.e. you are going from interpolation to extrapolation. Splines are going to result in the overshoot that you are observing. The solution is simply to make sure that your 1D function has values in the interval [min(r), max(r)]. Note that in the original data, max(r) is about 424, while the function you are interpolating is defined on the range [-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)

Answer 2

Leo is right in his diagnosis. I'd like to suggest a simpler (I hope) remedy: to do what you want to do (which is basically rotate a Gaussian around its symmetry axis) and to get a reasonable answer in a 600x600 square you need a Gaussian 600*sqrt(2)=849 pixels long. If you can do that, then all further thttp://stackoverflow.com/questions/2443046/interpolating-1d-gaussian-into-2d-gaussianrickery is not necessary.

Edit: In other words: if you rotate something 600 pixels wide around its center, you get a circle 600 pixels in diameter. You want to cover a 600x600 square. For that you need a circle 849 pixels in diameter, since this is the diagonal of the square.

Answer 3

In the particular case of the gaussian, you may compute the gaussian by using the fact that it is separable:

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

So you still just need to store just OTFx in memory.