How to extract data of lower surface only in a step profile?

Hi,

I have a 3-D plot of step function and want to extract data of lower surface only. I tried to literally read the data of the four corners of the lower surface and apply the code given below. This code seems to work fine to me, however there are two drawbacks of using this method:

1. this method is not efficient as I am manually reading the values.

2. I have to do the same process after every experiment with different type of steps.

I was wondering if anyone could suggest me any process by which I don't need to manually read the four corners and can give me more efficient result. Please find the attached file of the 3D point data, it has 3 columns corresponding to x, y and z values. Here is the code which I am trying:

function [arr]=step_mask(filename, sheet,range)
  %%Reading the File.
  test = xlsread(filename,sheet,range);
  x=test(:,1);
  y=test(:,2);
  z=test(:,3);
  %%--------------------
  %%Masking the Lower terrace.

    pos = (x >= -12.88) & (x <= 45.5) & (y >= -28.08) & (y <= 12.15) & (z >= 407) & (z <= 460);

    arr = [x(pos), y(pos) ,z(pos)];

      end

Function call:

arr=step_mask('Step_scan01_ex.xlsx','Sheet1', 'A:C');
figure;
plot3(arr(:,1),arr(:,2),arr(:,3),'.'); grid on
xlabel('x(mm)'); ylabel('y(mm)'); zlabel('z(mm)');
title('Masked plot');

Answer 1

[Jim Joy]

Hi Swati,

One way to do this is to get an estimate for the magnitude of the gradient at each point in your dataset, and then filter to find the points with a gradient whose magnitude is below a certain threshold. To do this, you will first need to grid your data. This workflow is illustrated below:

 [xgr,ygr] = meshgrid(min(test(:,1)):1:max(test(:,1)),...
                      min(test(:,2)):1:max(test(:,2)));
 vq = griddata(test(:,1),test(:,2),test(:,3),xgr,ygr);
 [GX,GY] = gradient(vq);
 idx = sqrt(GX.^2 + GY.^2) < 0.6;

You can then extract the points in this region, and use the "boundary" function to find the limits of the 'x' and 'y' coordinates. After this, the "inpolygon" function can be used to determine which of the points in your dataset lie outside of the plateau. This is shown below:

 xPt = xgr(idx);
 xPt = xPt(:);
 yPt = ygr(idx);
 yPt = yPt(:);
 k = boundary(xPt,yPt);
 idx2 = ~inpolygon(test(:,1),test(:,2),test(:,3));
 testFilt = [test(idx2,1) test(idx2,2) test(id2,3)];

I have tested this for the dataset you posted, and it works reasonably well. Noisier data may require some filtering (or a different spacing of the mesh produced by "meshgrid"), but the approach should be reasonably robust.

Best Regards,

Jim

[Jim Joy]

I should add that you can use the boundaries of the plateau to extract the surfaces as well.

Best,

Jim

[Swati Jain]

Thank you for your response and sorry for late reply. I have couple of questions about this code. Please accept my apologies in advance as I have so many questions and I am not very good at MATLAB. Could you please help me to understand this code? First of all I got this error when I ran this code.

Not enough input arguments.

Error in inpolygon (line 64)
if ~isvector(xv) || ~isvector(yv)

Error in mask (line 11)
 idx2 = ~inpolygon(test(:,1),test(:,2),test(:,3));

Other than that I have following questions:

1. I didn't quite understand how "griddata" function is working. I read the documentation and it said that it fits the surface of the of the form v = f(x,y) to the scattered data in the vectors (x,y,v). The griddata function interpolates the surface at the query points specified by (xq,yq) and returns the interpolated values, vq. The surface always passes through the data points defined by x and y. But I don't understand how it is interpolating the data? and the value return in vq is the vlue of z data? It is looking like z value to me. Why it is giving me a matrix of 82x47?

2. How gradient has been calculated from vq w.r.t. x and y? I am assuming vq has height value only so how dz/dx and dz/dy has been calculated?

3. How did you choose the following identity?

idx = sqrt(GX.^2 + GY.^2) < 0.6;

why sqrt(GxX.^2+GY.^2)? and why <0.6?

4. How does the following line work?

xPt = xgr(idx);

5. According to "inpolygon" documentation, inpolygon(xq,yq,xv,yv) returns in indicating if the query points specified by xq and yq are inside or on the edge of the polygon area defined by xv and yv.

As far as I know ~ will give me points lying outside of it. Please correct if my understanding is wrong. It should have 4 inputs and you used 3 and I think this is the place the error comes from.

idx2 = ~inpolygon(test(:,1),test(:,2),test(:,3));

[Swati Jain]

Hi,

I changed the code of line

idx2 = ~inpolygon(test(:,1),test(:,2),test(:,3)); to 
idx2 = ~inpolygon(test(:,1),test(:,2),test(:,3),k);

then I get the result attached in the file. The result I received is not masking the lower surface only. It seems that it is interpolating my data and I don't see any masking to the lower surface.