In Earth Science research, sometimes we need to construct 3D surfaces from given points, for example: creating the fault surface, locating a subducting slab from earthquake hypocenters, etc.  in a region of interest in X-Y plane.

In this example, we will show how to create a best-fit quadratic surface from given points in 3D using Matlab. The code is written in the following steps:

  1. Input the data 3D points: x, y, z
  2. Calculate the best fitting curve, store the parameters in C matrix.
  3. Define a region of interest xx, yy and using C to get zz.
  4. Plot the xx, yy, zz surface with x, y, z data.
  5. Make comparison (probabilistic density function of misfit between z and zz).

We use getPolygonGrid.m from mathworks.com.

Main program:

fit1fit2

% We start with some some 3d points
data = mvnrnd([0 0 0], [1 -0.5 0.8; -0.5 1.1 0; 0.8 0 1], 50);

x = data(:,1); y = data(:,2); z = data(:,3);


% Make a best-fit quadratic curve from the given points
C = x2fx(data(:,1:2), 'quadratic') \ data(:,3);


% Define an area of interest in X-Y plane by polyon nodes
xv = [-3 -3 0 3 3 0];
yv = [ 0 3 3 0 -3 -3];

% Create a area of interest in x, y plane
inpoints = getPolygonGrid(xv, yv, 1/0.1);
xx = inpoints(:,1);
yy = inpoints(:,2);

% Create the corresponding zz value from the above grid
zz = [ones(numel(xx),1) xx(:) yy(:) xx(:).*yy(:) xx(:).^2 yy(:).^2] * C;
zz = reshape(zz, size(xx));


% plot points and surface
figure('Renderer','opengl')
hold on;

% Plot the points
line(data(:,1), data(:,2), data(:,3), 'LineStyle','none', 'Marker','.', 'MarkerSize',25, 'Color','r')

% Make a triangular surface plot from xx, yy, zz
tri = delaunay(xx,yy); %x,y,z column vectors
trisurf(tri,xx,yy,zz, 'FaceColor','interp', 'EdgeColor','b', 'FaceAlpha',0.2);

% Change view propeties
title('Fitting surface');
grid on; axis tight equal;
view(9,9);
xlabel x; ylabel y; zlabel z;
colormap(cool(64))

% Calculate the misfit between surface and each points:
zfit = [ones(numel(x),1) x y x.*y x.^2 y.^2] * C;
h = -2: 0.01: 2;
misfit = z - zfit;
mu = mean(misfit); sigma = std(misfit);
f = exp(-(h-mu).^2./(2*sigma^2))./(sigma*sqrt(2*pi));
figure
title('Probability function of misfit');
hold on
plot(h,f,'LineWidth',1.5)