Интерполяция 1D Гаусса в 2D Гаусса - PullRequest
Новая
       16

Интерполяция 1D Гаусса в 2D Гаусса

2 голосов
/ 14 марта 2010

Допустим, у меня есть 1D функция Гаусса. Его длина 600.

Я хочу интерполировать его в 2D гауссиан размером 600 X 600.

Это код, который я написал (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]);

Проблема в том, что в конце я получаю Overshoot: alt text

Как я могу предотвратить это превышение? Есть ли лучший алгоритм интерполяции для монотонных нисходящих функций?

ПРИМЕЧАНИЕ: Я ищу общее решение для интерполяции одномерной функции в двумерную радиально-симметричную функцию, гауссовский пример только.

Ответы [ 3 ]

5 голосов
/ 14 марта 2010

РЕДАКТИРОВАТЬ: на основе ваших разъяснений, ясно, что происходит. Вы пытаетесь интерполировать функцию за пределами диапазона доступных данных - то есть вы переходите от интерполяции к экстраполяции. Сплайны могут привести к выбросу, который вы наблюдаете. Решение состоит в том, чтобы просто убедиться, что ваша 1D функция имеет значения в интервале [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)

альтернативный текст http://img54.imageshack.us/img54/8255/clipboard01vz.png

2 голосов
/ 15 марта 2010

Лев прав в своем диагнозе.Я хотел бы предложить более простое (я надеюсь) лекарство: чтобы делать то, что вы хотите сделать (то есть вращать гауссиан вокруг своей оси симметрии), и чтобы получить разумный ответ в квадрате 600x600, вам нужен гауссовский 600 * sqrt(2) = 849 пикселей в длину.Если вы можете сделать это, то все дальнейшие t /1795414/interpolyatsiya-1d-gaussa-v-2d-gaussa не нужны.

Редактировать : Другими словами: если вы вращаете что-то шириной 600 пикселей вокруг его центра, вы получите круг600 пикселей в диаметре.Вы хотите покрыть 600x600 квадрат .Для этого вам понадобится круг диаметром 849 пикселей, так как это диагональ квадрата.

1 голос
/ 15 марта 2010

В частном случае гауссиана, вы можете вычислить гауссиан, используя тот факт, что он отделим: 

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

Так что вам все еще нужно хранить только OTFx в памяти.

...