How do I obtain the area of each contour on a contour3 plot?

I have written a code that overlaps a 3D scan with some external data from the object being analysed, resulting in a plot similar to the one showed below

I am trying to obtain the area of each contour, but so far I have not had any luck.

To obtain the contourplots, I used the following:

clear all , close all, clc, 

%processing code that selected a specific area of the 3D scan where the
%sample was located
[x1, y1, z_data] = process_3d_data(sample_file);

figure; 
t = tiledlayout(1,1);
ax1 = axes(t);
contour3(ax1, x1, y1, z_data, linspace(0.05, 0.8, 15), 'ShowText', 'off'); 

The sample_file originally comes in an .asc format, but to be able to upload it I have changed it into a .mat format,

Answer 1

[Star Strider]

I wrote some utility routines to get information from contour plots a while ago, and tweaked them for this project. The ap[proach is straightforward if not always easy to describe. The trapz function will calculate the area of a closed contour if its arguments are consistent. Here, that requires a bit of effort, since the components of the contours are not arranged in a way that trapz can integrate them correctly as they exist. First, the individual contours are borken up into non-consecutive pieces, making this impossible using the data as presented. To correct for that, the code calculates the centroids of the contours, and then uses that information to calculate the distances from the centroids (each contour has a different centroid), and the angles from the centroid to the particular element of the contour. The code then sorts both of these by angle, then uses that information to calculate the (x,y) coordinates of the contours in a way that trapz can work with them. The only problem is the outer (lowest) contour. Since it is fractured, the data for it may not be reliable.

Incomplete contours use the axis boundaries (actually, the existing (x,y) at those boundaries) to complete them, so if you do not want the incomplete contour areas to appear in the ‘Results’ table, the easiest way to do that is to delete them from the ‘Levels’ vector. That applies to the lowest level, and those that run off the right edge. My code does not automatically check for those, and other than searching for gaps in the ‘y’-vectors closest to the median or mean ‘y’ value with the highest ‘x’ value, it is not obvious to me how to detect and remove them.

I added ‘LevelVolume’ in the event you want that. It simply multiplies ‘LevelArea’ by the associated ‘Levels’ value for each of them.

The code is relatively straightforward —

% type('process_3d_data.m')                                                       % View Function File


clear all , close all, clc, 

load('data.mat')                                                                % Retrieve 'data' From .mat file
writematrix(data,'data.txt')                                                    % Write .txt File
sample_file = 'data.txt';                                                       % Supply Required Variable
% whos

%processing code that selected a specific area of the 3D scan where the
%sample was located
[x1, y1, z_data] = process_3d_data(sample_file);

figure; 
% t = tiledlayout(1,1);
% ax1 = axes(t);
% contour3(ax1, x1, y1, z_data, linspace(0.05, 0.8, 15), 'ShowText', 'off');
contour3(x1, y1, z_data, linspace(0.05, 0.8, 15), 'ShowText', 'off');
view(0,90)                                                                      % Top View
axis('equal')

% return

figure
[M,H] = contour3(x1, y1, z_data, linspace(0.05, 0.8, 15), 'ShowText', 'off');  % Get Outputs

Levels = H.LevelList;
BufferLevels = buffer(Levels,5)                                             % Show All Levels (Informational)

BufferLevels = 5×3

0.0500 0.3179 0.5857 0.1036 0.3714 0.6393 0.1571 0.4250 0.6929 0.2107 0.4786 0.7464 0.2643 0.5321 0.8000

for k = 1:numel(Levels)
    idx = find(M(1,:) == Levels(k));                                        % Return The 'idx' Value For This Level
    ValidV = rem(M(2,idx),1) == 0;                                          % Checks To Be Certain That A Specific Level Value In 'M(1,:)' Is A Level VAlue And Not A Point On The Associated Contour
    StartIdx{k,:} = idx(ValidV);                                            % Start Index For That Contour (There May Be Several Unconnected Vectors)
    VLen{k,:} = M(2,StartIdx{k});                                           % Associated Vector Length
end

% StartIdx
% VLen

figure
hold on
for k1 = 1:numel(Levels)
    % fprintf('----- k1 = %2d,  Level = %.3f --------------------------------------------------\n', k1, Levels(k1))
    xvc = [];
    yvc = [];
    for k2 = 1:numel(StartIdx{k1})
        idxv = StartIdx{k1}(k2)+1 : StartIdx{k1}(k2)+VLen{k1}(k2);          % Index For Contour 'k1'
        xv = M(1,idxv);                                                     % X-Vector
        xvc = [xvc xv];                                                     % Concatenated X-Vectors
        yv = M(2,idxv);                                                     % Y-Vector
        yvc = [yvc yv];                                                     % Concatenated X-Vectors
        % plot(xv, yv)                                                        % Draw Contour Or Contour Segment (Check)
    end
    if numel(xvc) > 1
        Cx(k1) = trapz(xvc, xvc.*yvc) / trapz(xvc, yvc);                    % Calculate Centroid X-Coordinate
        Cy(k1) = trapz(yvc, xvc.*yvc) / trapz(yvc, xvc);                    % Calculate Centroid Y-Coordinate
        r = hypot(yvc-Cy(k1),xvc-Cx(k1));                                   % Radius From Centroid
        a = atan2(yvc-Cy(k1),xvc-Cx(k1));                                   % Angle With Respect To Centroid
        [as,idx] = sort(a,'descend');                                       % Sort Angles
        rs = r(idx);                                                        % Match Radius Values
        xs = rs .* cos(as);                                                 % Reconstrructed X-Coordinates For This Contour
        ys = rs .* sin(as);                                                 % Reconstrructed Y-Coordinates For This Contour
        LevelArea(k1,:) = trapz(xs, ys);                                    % Area For That Level
        Levelv(k1,:) = Levels(k1);                                          % This Level (Not All Levels Can Be Calculated)
        plot(xs+Cx(k1), ys+Cy(k1))                                          % Plot Reconstructed Contour
    end
end
% plot(Cx, Cy, '+k')
hold off
axis('equal')
title('Contours Using Reconstructed Data')

LevelVolume = Levelv .* LevelArea;
Result = table(Levelv, LevelArea, LevelVolume, 'VariableNames',{'Level','Level Area','Level Volume'})

Result = 13×3 table
     Level     Level Area    Level Volume
    _______    __________    ____________

       0.05     0.42034        0.021017  
    0.10357      95.771          9.9191  
    0.15714      90.544          14.228  
    0.21071      85.198          17.952  
    0.26429      79.727          21.071  
    0.31786      73.773          23.449  
    0.37143      67.055          24.906  
      0.425      59.901          25.458  
    0.47857      51.263          24.533  
    0.53214      42.112           22.41  
    0.58571      32.216          18.869  
    0.63929      21.247          13.583  
    0.69286      8.5663          5.9352  

EDIT — (12 Mar 2024 at 08:15)

It later occurred to me that there were some problems with my original code. These changes correct those errors. I also added more comments to document how the code works.

.

[Goncalo Costa]

You're always a great help. Thanks for that.

[Star Strider]

As always, my pleasure!

Thank you for posting an extremely interesting problem!