smooth 2D array and maintain boundaries

12 次查看（过去 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 个评论
KSSV
KSSV 2018-11-8
Don't include the edge points in smooth....is that okay?
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 个评论
显示 1更早的评论
Alexander Laut
Alexander Laut 2018-11-9
编辑：Alexander Laut 2018-11-9
It looks pretty good and I'll have to spend a bit of time to figure out what's going on, but I notice that if the boundaries are all zero, then there is an issue with reshape.
To me the suspect looks to be the sparse call but I'll have to tinker around to understand the full picture
In the meantime, I just set my border to eps instead of zero and the algorithm seems to work! Thanks!
Bruno Luong
Bruno Luong 2018-11-9
Oh yes, probably there is a bug
inside = ~bdr(:);
should be
inside = ~bmask(:);

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 2-D and 3-D Plots 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by