Search code examples
matlabinterpolationsplinepoint-clouds

How can I interpolate and smooth 3D point clouds in MATLAB?

I am hoping to get some help with a point cloud data processing problem I am faced with. Basically, I have lots of point cloud data that is patchy and noisey. My aim is therefore to approximate data where it is missing in the “patchy areas” and apply some form of light smoothing to filter the noise.

My first attempt to solve this was the interpolation methods in MATLAB. This was executed as follows and provided good results, in that the interpolated Z points across the working XY grid looks like the shape I am expecting.

% Load Point Cloud:
Point_Cloud  = importdata(‘Point_Cloud_1.txt')
x            = Point_Cloud(1,:)';
y            = Point_Cloud(2,:)';
z            = Point_Cloud(3,:)';

% Interpolate inspection points:
Density = 300;
[X,Y]   = meshgrid(linspace(min(x), max(x), Density), linspace(min(y), max(y), Density));
F       = scatteredInterpolant(x, y, z, 'natural','linear');
Z       = F(X,Y);
Int_PC  = [reshape(X,Density^2,1) reshape(Y,Density^2,1) reshape(Z,Density^2,1)];
Int_PC(any(isnan(Int_PC{i}),2),:) = [];  

% Plot results:
scatter3(x, y, z, 20, 'r', 'fill'); % Original data
hold on;
scatter3(Int_PC(:,1), Int_PC(:,2), Int_PC(:,3), 20, 'r', 'fill'); % Interpolated data

The problem with this is that noisey data is used to compute the interpolant, F, so the above code only solves the patchy data problem.

I then considered spline fitting using the Curve Fitting Toolbox. The thin-plate smoothing spline seems to make sense as it accepts scattered (non-gridded) data and it does not interpolate all points, thereby smoothing noise. The code for this is below. When executed, the results are disappointing because the fit between the original points and the computed surface is very poor (beyond what is needed to smooth noise).

Spline = tpaps([x;y],z);
fnplt(Spline)

Thin plate spline example

Could anybody suggest any solutions?

Thanks.

Solution

  • One proposition, using the Savitzky-Golay Filter:

    So

    1. Use this filter
    2. Use your code to complete the missing part.

    EXAMPLE

    %We build the noisy 3D point cloud

[X,Y] = meshgrid(0:0.1:4.9,0:0.1:4.9);
Z = sin(X)+cos(Y)+0.5*rand(50,50);

% The smoothing (with sgolayfilt(Z,order,length of the moving windows))
t1 = sgolayfilt(Z.',2,25); %smoothing over the x-axis
t2 = sgolayfilt(Z,2,25); %smoothing over the y-axis
t  = (t1.'+t2)/2; %smoothing in booth directions

surf(X,Y,t)

    Before smoothing

    before

    After smoothing

    after

    EDIT

    The same operation but with scattered data:

    %random X and Y

X = rand(100,1)*5;
Y = rand(100,1)*5;

Z = sin(X)+cos(Y);

%Calculate the interpolant
F =  scatteredInterpolant(X,Y,Z);

%Fix the grid size
gs = 0.1; %grid size

tx = min(X(:)):gs:max(X(:));
tz = min(Y(:)):gs:max(Y(:));

%Scattered X,Y to gridded x,y
[x,y] = meshgrid(tx,ty);

%Interpolation over z-axis
z = F(x,y);

t1 = sgolayfilt(z.',2,5);
t2 = sgolayfilt(z,2,5);

t  = (t1.'+t2)/2;

surf(X,Y,t)
colormap winter