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?

[KSSV]

Don't include the edge points in smooth....is that okay?

[Alexander Laut]

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.

Answer 1

[Bruno Luong]

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

[Bruno Luong]

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]

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]

Oh yes, probably there is a bug

inside = ~bdr(:);

should be

inside = ~bmask(:);