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)
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'})
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.
.
