smooth 2D array and maintain boundaries

16 ビュー (過去 30 日間)
Alexander Laut
Alexander Laut 2018 年 11 月 7 日
編集済み: Bruno Luong 2018 年 11 月 9 日
I have a 2D array which has some noisy data. I need to smoothen it but preserve the values towards the edges.
I have played around with smoothdata and it works but it seems that it only preserves values along the edges of one dimension. For my purposes, it's important that the edges of the array are unchanged after the smoothing process though are included in that they would affect the smoothing of adjacent points inwards of the array.
Is there a clever way in using smoothdata in both directions to get my result or another alternative?
  2 件のコメント
Alexander Laut
Alexander Laut 2018 年 11 月 8 日
Not quite, if I do that I may get sharp discontinuities towards the boundaries. I'd like for the smoothing algorithm to include the edges as means in weighting the adjacent points in a smooth transition.
In a radical example, say I was describing Z = X^2+Y^2, a paraboloid, but wanted it's edges to intersect with the Z = 0 plane. I could just set the edges to zero, but then I'd have a sharp step from the edge to the inside of the surface.
I'd like the resulting surface to look sort of like a pillow in that you'd still see the bowl of the paraboloid towards the center, but rather than diverge towards the edges, it'd be progressively dragged back to the Z = 0 plane, for which I've defined the edge.

サインインしてコメントする。

採用された回答

Bruno Luong
Bruno Luong 2018 年 11 月 8 日
編集済み: Bruno Luong 2018 年 11 月 9 日
One way (bug edited version):
% Generate fake data
A = peaks(100);
A = A(26:75,21:81); % crop so we get varying boundary data
[m,n] = size(A);
% Add noise
A(2:m-1,2:n-1) = A(2:m-1,2:n-1) + 0.5*randn(m-2,n-2);
% A = A+ 0.5*randn(m,n);
% Solve Laplace's equation
[m,n] = size(A);
x = 2:m-1;
y = 2:n-1;
ifun = @(x,y) sub2ind([m,n],x(:),y(:));
[X,Y] = ndgrid(x,y);
iXY = ifun(X,Y);
I = repmat(iXY ,[5 1]);
J = [iXY;
ifun(X-1,Y);
ifun(X+1,Y);
ifun(X,Y-1);
ifun(X,Y+1) ];
V = [4,-1,-1,-1,-1].*ones(size(iXY));
M = sparse(I,J,V(:),m*n,m*n);
bmask = true(size(A));
bmask(2:m-1,2:n-1) = false;
bdr = bmask.*A;
y= -M*bdr(:);
inside = ~bmask(:);
M = M(inside,inside);
Ainside = M\y(inside);
Ainside = reshape(Ainside,m-2,n-2);
A0 = A;
A0(2:m-1,2:n-1) = Ainside;
% B is equal to 0 on bounddary
B = A - A0;
% Filter inner part using FFT
Bi = B(2:m-1,2:n-1);
fftBi = fft(fft(Bi,[],1),[],2);
px = 10; % filter parameter in x: smaller value, stronger filtering
xf = Gaussfilter(m,px);
py = 10; % filter parameter in y
yf = Gaussfilter(n,py);
fftBi = xf(:) .* fftBi .*yf(:)';
Bif = real(ifft(ifft(fftBi,[],2),[],1));
B(2:m-1,2:n-1) = Bif;
% Put together
Af = B+A0;
% Check
close all
ax1=subplot(1,2,1);
surf(A)
title('Original data');
ax2=subplot(1,2,2);
surf(Af)
title('Filtered data');
linkprop([ax1, ax2],{'CameraUpVector', 'CameraPosition', 'CameraTarget', ...
'XLim', 'YLim', 'ZLim'});
function f = Gaussfilter(m,p)
mid = floor((m-1)/2);
f = exp(-((0:mid)/p).^2);
f = [f, flip(f(2:mid-mod(m,2)))];
end
  3 件のコメント
Bruno Luong
Bruno Luong 2018 年 11 月 9 日
Oh yes, probably there is a bug
inside = ~bdr(:);
should be
inside = ~bmask(:);

サインインしてコメントする。

その他の回答 (0 件)

タグ

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by