How to measure surface area layer by layer using Matlab

i am using following code to find area of layers. I am trying but masking somewhere some mistake. Any help from community will be highly appreciated.

Attached file are the xyz-coordinates. I measure the surface area of one layer using bellow code. This coding is combining all layers. I want to find area of each layer individually. Regards to all for their cooperation and guidance.

vertices = load(' attached file ')
vertices = round(sortrows(vertices,2));
Vx = round(vertices(:,1));
Vy = round(vertices(:,2));
Vz = round(vertices(:,3));
G_surface_area =0
layer = []
surface_area = 0

for i = min(Vy):10:max(Vy) % This for loop, i am using to measure by 10 spacing. 
   minVy =  i

        for j = i:i + 9   % This measure is for layer to find its area. 
           
            layer_new = vertices(any(Vy== j,2),:);      
            layer = [layer; {layer_new}];
        end 
             
               layer = cell2mat(layer)                   
                Lx = layer(:,1);
                Ly = layer(:,2);
                Lz = layer(:,3);
                surface_area = polyarea(Lx,Lz)
       
           G_surface_area = G_surface_area + surface_area
            
end

Answer 1

[KSSV]

data = importdata("Cone_20000_points.txt") ; 
data = round(data) ; 
x = data(:,1) ; y = data(:,2) ; z = data(:,3) ; 
zi = unique(z) ;   
N = length(zi) ; 
A = zeros(N,1) ; 
tol = 10^-5 ; 
for i = 1:N
  idx = abs(z-zi(i))<=tol ; 
  % layer coordinates 
  xl = x(idx) ; yl = y(idx) ; zl = z(idx) ; 
  Layer = [x(idx) y(idx) z(idx)] ;   % these are the Layer coordinates 
  % get boundary of layer 
  id = boundary(xl,yl) ; 
  Boundary = [xl(id) yl(id) zl(id)] ;   % Boundary coordinates 
  A(i) = polyarea(xl(id),yl(id)) ; 
end

[M.S. Khan]

Hi Dear KSSV, thanks for your guidacne. The code is showing following error.

Array indices must be positive
integers or logical values.

Error in untitled13 (line 12)
  xi = x(idx) ; yi = y(idx) ; zi =  z(idx) ; 
  
 

You have used in the code i.e. tol = 10^-5 ; what does it mean, could you plz explain me your coding. i will be very thankful. Regards

[KSSV]

There was a typo error..now I have modified code...check now.

[M.S. Khan]

Hi KSSV, Hope and pray, you are fine. can i get it layer by layer? why have you used 10^-5? plz guide.

[KSSV]

It is layer by layer...note that the layers are along the XY plane and each layer have constant z. To capture the layer along z-axes, it needs to get the z values which all remain same. So using tolerance and differrence the indices are obtained. You can also z=zi(i) instead.

[M.S. Khan]

Do you mean, it will be a xy-layer with thickness of z-values, right. i am trying to validates the results. Thanks KSSV for all your cooperation. Regards

[M.S. Khan]

Could you plz revisit it. Using debugging, i checked and its repeating the same areas i.e. 18192. Regards

[KSSV]

You can plot the coordinates of Layer..using plot3 to check it. Note that there is no thickness....the layer is a surface..it will be a circle if your cone is perfect.

[M.S. Khan]

yes i checked it using scatter3(), it is showing a surface of layer but the same layer is repeated again and agains with same area i.e. 18192

[KSSV]

Edited the code....yes there is a bug in the code. I have used zi as unique elements and also a layer z coordinates. This creates a problem. Now check it.

[M.S. Khan]

Hi KSSV, thank for you all your professional support. I got the total area of the cone using codel like: Area = Area + A(i) = 83856 which is very high number. There is some problem in the code.

Do you think, its correct? What will be unit e.g. mm^2 or cm^2? What is the default unit?

[KSSV]

You should be knowing the unit....area will be in the units which x,y,z are in. If you want surface area there is a formula for the cone right? The present layer by area and syrface area will be different. They won't be equal.

[M.S. Khan]

They will not be exactly equal because they will be approximated, right.

They results have huge difference.

surface_area = pi *78 *78 + pi * 78 * 188.4 = 65279.78

surface_area area using our code = 314400

-----------------------------------------------------------------------------------------------------------------------------

Formula for cone's volume = 1/3 *pi * r^2 * h

where r means radius and h means height of the cone.

Volume of the cone = 1/3 * pi * 78 * 78 * 188.4 = 1,200,324.64

Volume of the cone measured using our code = surface area * height = 314400*188.4 = 59232960

[KSSV]

No what the code calculates is not surface area.......the code is calculating the area of every circle in the layer of cone. How it can be a surface area?

[M.S. Khan]

Really!

if we combine all areas of all circles, will it not be surface area of the cone.

Please guide me how can we proceed then. Also, i was just googling the conversion process of relationship between pixels and mm units. It also needs some extra calculations.

[KSSV]

Yes...the area of each circle added up will not give you surface area. What you have to do is, get all the boundary points, which from a cone and calculate the area of each cell and add them up. You need to proceed like this:

clc; clear all ;
data = importdata("Cone_20000_points.txt") ; 
data = round(data) ; 
x = data(:,1) ; y = data(:,2) ; z = data(:,3) ; 
zi = unique(z) ;   
N = length(zi) ; 
A = zeros(N,1) ; 
tol = 10^-5 ; 
Bounadry = cell(N,1) ; 
for i = 1:N
  idx = abs(z-zi(i))<=tol ; 
  % layer coordinates 
  xl = x(idx) ; yl = y(idx) ; zl = z(idx) ; 
  Layer = [x(idx) y(idx) z(idx)] ;   % these are the Layer coordinates 
  % get boundary of layer 
  id = boundary(xl,yl) ; 
  Boundary{i} = [xl(id) yl(id) zl(id)] ;   % Boundary coordinates 
  A(i) = polyarea(xl(id),yl(id)) ; 
end
Boundary = cell2mat(Boundary(:)) ;   % this will give only a surface cone/ hollow cone 
bx = Boundary(:,1) ; by = Boundary(:,2) ; bz = Boundary(:,3) ; 

Apply delaunayTriangulation to [bx by bz]. It will mesh the cone with triangles. Get the area of each triangle and sum them. This will give you surface area.

[KSSV]

Other option is from the points we can get the radius and height of the cone and use the formula for surface area of cone.

[M.S. Khan]

Dear KSSV, the code is giving some erros as following. i am trying to remove errors.

Error using cat
Dimensions of arrays being concatenated are not consistent.

Error in cell2mat (line 75)
            m{n} = cat(2,c{n,:});

Error in untitled4 (line 20)
Boundary = cell2mat(Boundary) ;   

The used of braces in Boundary{i} = [xl(id) yl(id) zl(id)] is making problem.

[KSSV]

Boundary = cell2mat(Boundary(:));

Answer 2

[M.S. Khan]

Hi Dear KSSV, along with your code, i followed following procedure ( hints from Mr. Uday) to generate the traingulated mesh of the cone. Then i estimated the area of triangulated cone. could you please guide, am i correct??

bx = Boundary(:,1) ; by = Boundary(:,2) ; bz = Boundary(:,3) ; 
T = delaunay(bx, by)
TO = triangulation(T, bx(:),by(:),bz(:)); %creates a 3-D triangulation representation with the point coordinates specified as column vectors X, Y, and Z
trimesh(TO) %plots the triangular mesh
%It will mesh the cone with triangles. 
%Get the area of each triangle and sum them. 
% This will give you surface area. 
%% How to find area
area = 0; %initialize the area
%loop over all the desired small triangles and accumulate the areas
for i = 1:size(TO.ConnectivityList,1)   %the number of rows gives the number of triangles produced
    a =  TO.Points(TO.ConnectivityList(i,:),:); %this gives the 3 vertices of the ith triangle
    if any(a(:,3) > max(bz))
        continue;  %if the triangle has a vertex whose z - coordinate is > 0.4 ignore it
    end
    %extract the three vertices
    p1 = a(1,:);  
    p2 = a(2,:);
    p3 = a(3,:);
    format compact
    area = area + 0.5 * norm(cross(p2-p1,p3-p1));  %calculate the area of the small triangle and add with previous result
end
area

[KSSV]

Yes this should work...This will give you the surface area.

trimesh(TO)

The above mesh was all cone? Actually I don't have access to MATLAB. I work on Octave and many functions are missing in it. Else I would have tried for complete working code.

[M.S. Khan]

Yes, the mesh of triangles for cone. You are my best friend. I wish to have your contact email to ask for your guidance in future if possible. Regards for all your kid guidance.

[KSSV]

That's good... hope the present answer and your formula obtained value are close. You can mail me through my author page or send me linkedin request, we can keep in touch.

[M.S. Khan]

This area is in pixels, right. Do you have any idea how to convert them into mm^2.