Getting average from surf

I have the white light interferometry data of a laser ablated surface in .txt format and I have used the following code to generate the surface profile in matlab. The task is to take the average in the Z direction (height) for the whole surface so that I can make an estimation of the ablation depth (please refer to the image). I would really appreciate any advice on how to achieve this task. Thanks in advance!

clear;

clc;

clear all;

close all;

%Loading the matrix

X=load('5J-1.txt');

M=X;

Xm = 887; % measuring field in µm

Ym = 662; % measuring field in µm

% Interpolating the matrix for missing values

valid= ~isnan(M);
M(valid);   
%Produce Logical-Matrix

good= find(valid);
%Column vector with positions

M(~valid) = interp1(good, M(valid), find(~valid),'nearest');

M2=M;

f1=figure(1)

h1=surf(M2);

set(h1,'LineStyle','none')

shading interp;view([0,0]);

% Filter for refinement 

FX = [1 1 2 1 1;1 1 4 1 1;1 1 4 1 1; 1 1 4 1 1; 1 1 2 1 1];
FX = FX./sum(FX(:));

M3 = imfilter(M2,FX);

f2=figure(2)

h2=surf(M3);

set(h2,'LineStyle','none')

% Histogram of depth values

f3=figure(3)

histogram(M3);

Answer 1

[Suhan]

Since you have X, Y, Z data,you do not need to depend on surf function to find out what the average height (z-coordinate) of the surface is. However, since the surface is non-linear, you cannot simply use MEAN.

There are two methods to accomplish this, depending on the whether or not there are singularities in the surface fit.

If there are discontinuities, the best way to estimate the average height of the surface is to sample a large number of points on the surface and divide by the total number of points sampled. This is similar to a Monte Carlo estimation. The following code shows how do to this.

num_samples = 100;
x_samples = linspace(min(X), max(X), num_samples);
y_samples = linspace(min(Y), max(Y), num_samples);
[xx, yy]     = meshgrid(x_samples, y_samples);
z_samples = fitted_surface(xx(:), yy(:));
average_z = mean(z_samples);

In the above code "fitted_surface" is the surface returned by Curve Fitting Toolbox and "X" and "Y" are the X and Y data.

If the surface fit does not have singularities, the INTEGRAL2 function can be used to numerically evaluate the volume under the surface. So you can have a measure of the average height as the volume of the surface divided by the area projected onto the XY plane. The following code demonstrates how to do this

my_fun = @(x, y) fitted_surface(x, y);
volume = integral2(my_fun, min(X), max(X), min(Y), max(Y));
average_height = volume  / ((max(X) - min(X)) * (max(Y) - min(Y)));

In the above code "fitted_surface" is the surface fit returned by Curve Fitting Toolbox. The anonymous function "my_fun" is necessary to pass the surface fit into the INTEGRAL2 function. The division divides the volume by the area the surface projects onto the XY plane.