smooth 2D array and maintain boundaries
显示 更早的评论
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
2018-11-8
Don't include the edge points in smooth....is that okay?
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
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
2018-11-8
There are simpler way if don't want to preserve approximately boundary data.
I might do something not entirely accurate with the FFT (size, flip, etc...) in the above code, but it doesn't look too bad.
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
2018-11-9
Oh yes, probably there is a bug
inside = ~bdr(:);
should be
inside = ~bmask(:);
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Descriptive Statistics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
