Find the centers of multiple polygons

I have an nx3 matrix that holds the locations for the vertices of a chain of hexagonal shapes (images 1, 2).

Side-note: The coordinates for the vertices of these hexagons are not in any logical order within the matrix. The coordinates were retrieved from an STL file which may be the reason why.

I'm trying to find the centroid of each of these hexagons and index those six points to its respective centroid. Once I've indexed these, I plan on writing some more code to manipulate the six points farther away from their respective centroids - making the hexagons larger and making my structure thicker [ they are currently too thin ] (image 3).

My problem is the centroid and indexing part. I've tried clustering the hexagons into n/6 clusters using the k-means and subclust functions, but there're always points that do not fall in the center of a hexagon (image 4).

Would anyone know what I should try?

[I've also attached the matrix, in case you would like to play around with the data as well.]

Answer 1

[Image Analyst]

Isn't the centroid of a perfect hexagon just the mean of the x and y coordinates of the vertices?

xCentroid = mean(xVertices);
yCentroid = mean(yVertices);

Just do that for every hexagon and you'll have all the centroids.

[Jeremy]

That part I do get, but the matrix itself isn't in any particular order (which I failed to mention; sorry!!), which makes singling out any one hexagon's coordinates difficult as well as finding the centroid of that hexagon and indexing those hexagon vertices to that centroid, if that makes sense.

I'm trying out other statistical clustering methods to pick out the hexagons and then apply the code you suggested. Thanks!

[Image Analyst]

So it's just a random assortment of vertices where the hexagon groupings only become apparent after all the dots have been plotted?

[Jeremy]

Yes. The data was imported from an STL file, where I presume the random ordering is necessary to construct mesh/surface type objects. Any ideas?

[Image Analyst]

Nope.

Answer 2

[DGM]

If you have an STL file, then you already have connectivity information. You just need to get rid of the longitudinal connections in order to isolate the groups of vertices.

Consider the given tubeplot:

% the triangulated surface
unzip tubeplot.stl.zip % for the forum
T = stlread('tubeplot.stl');
[F V] = t2fv(T); % attached

% preview display
patch('faces',F,'vertices',V,'facecolor','w','edgecolor','k');
view(3); view(-34,28); camlight; 
axis equal; grid on

% isolate short edges
E = edges(T); % edges as indices [idxA idxB]
Ev = cat(3,V(E(:,1),:),V(E(:,2),:)); % edges as coordinates [x y z]
D = sqrt(sum(diff(Ev,1,3).^2,2)); % edge lengths
shortedges = D <= 1.2*prctile(D,20); % just the short edges
E = E(shortedges,:); % get rid of long edges

% plot them
figure
Ev = permute(Ev(shortedges,:,:),[3 1 2]);
plot3(Ev(:,:,1),Ev(:,:,2),Ev(:,:,3))
axis equal; grid on

% find the connected groups of vertex indices
G = graph(E(:,1),E(:,2)); 
figure
plot(G,'nodelabel',1:size(G.Nodes,1),'nodefontsize',6)

% sort them into a cell array of vertex indices
% calculate centroid position for each group
I = G.conncomp;
ngroups = numel(unique(I));
groups = cell(ngroups,1);
centroids = zeros(ngroups,3);
for k = 1:ngroups
    groups{k} = find(I==k);
    centroids(k,:) = mean(V(groups{k}(:),:),1);
end

% plot them again with their centroids
figure
plot3(Ev(:,:,1),Ev(:,:,2),Ev(:,:,3)); hold on
plot3(centroids(:,1),centroids(:,2),centroids(:,3),'k.','markersize',10)
axis equal; grid on

At that point, you have a set of grouped vertex indices and a set of corresponding centroid positions. You can do with that whatever you want. The triangulation should still remain valid unless the radius becomes large WRT the curvature.

Note that in this example STL, the ends of the object are capped using an extra center vertex. In this case, we might need to do some extra work to remove those two points from their corresponding vertex groups, but even if you just blindly offset all the points with respect to the group centroid, those cap centers shouldn't move far, since they should be coincident with the centroid. Not all tubeplots will be capped this way, so I'm not going to solve a problem that likely wasn't representative of OP's model.

The exact parameters used for selecting the short edges may vary depending on the tubeplot's resolutions. As the circumferential and axial edge lengths become similar, it may be difficult to isolate the intended groups. Again, as radius becomes large WRT curvature, the problem becomes more difficult.

I didn't bother making this example period-correct for 2014, so there are some anachronisms.

• Obviously, the triangulation() class didn't exist in 2014, but the TriRep() class did, so very little would need to be changed.
• The built-in STL tools didn't exist until R2018b, but there were tools available on the FEX. Again, little would need to be changed. AFAIK, until 2015, there were no STL decoders on the FEX which returned a unique vertex list. In order to actually detect connectivity, you'd need to prune the vertex list yourself.
• The biggest issue would be the use of the graph() class. There were tools on the FEX for working with undirected graphs, but it would take some work to adapt the code to suit those tools.
• I'm not sure if I've slipped in any implicit array expansion, but that might need to be reverted if I did.

EDIT:

I looked at OP's data to see if each tube was closed with an extra vertex, and they are. I have to wonder if modifying the STL is even necessary. If the original data is available, just generate the STL with the desired diameter to begin with. I used FEX #5562 to generate the example object, and that supports setting the tube radius and the number of circumerential points (e.g. 0.1 and 6 in this case).

Even if the original data isn't available, we can still avoid the need to do any of the vertex offsets. All we need to do is use the calculated centroids, along with other connectivity information from the STL's edges to reconstruct the original data.

[DGM]

For example

Let's say we had some curves in 3D and we generated a tubeplot and saved it as an STL -- but for some reason we lost the original data. I'm going to save the data for sake of comparison, but we're just going to work with the STL.

n = 25; % number of data points per curve
r = 0.1; % tube radius
nc = 6; % number of ring points

% the first object
t = linspace(0,2*pi,n);
curve = [cos(t)-0.1*t; sin(t); 0.2*(t-pi).^2];
[x,y z] = tubeplot(curve,r,nc); %FEX #5562
[F1,V1] = mesh2tri(x,y,z,'f');
data{1,1} = curve.';

% the second object
t = linspace(0,2*pi,n);
curve = [-cos(t)-0.1*t; sin(t)+0.2*t; real(0.2*(t-pi).^1.8)];
[x,y z] = tubeplot(curve,r,nc);
[F2,V2] = mesh2tri(x,y,z,'f');
data{2,1} = curve.';

% concatenate F,V data
F = [F1; F2+size(V1,1)];
V = [V1; V2];

% write the STL and data for the original curves
stlwrite(triangulation(F,V),'tubeplot2.stl')
save('tubedata.mat','data')

% preview display
patch('faces',F,'vertices',V,'facecolor','w','edgecolor','k');
view(3); view(-54,47); camlight; 
axis equal; grid on

So we start with the STL, sort it into connected components (each tube). We then do as we did before, and reduce each tube to its central curve, which should be an estimate of the original data.

% the triangulated surface
T = stlread('tubeplot2.stl');
[F V] = t2fv(T);

% process each object in the model
E = edges(T); % edges as indices [idxA idxB]
G = graph(E(:,1),E(:,2)); 
I = G.conncomp;
ngroups = numel(unique(I));
newdata = cell(ngroups,1);
for k = 1:ngroups
    mask = all(ismember(E,find(I==k)),2);	
    newdata{k} = detubify(E(mask,:),V);
end
view(-54,47) % set the viewing angle

% compare against the original data
load tubedata.mat
error = [immse(data{1},newdata{1}) immse(data{2},newdata{2})]

error = 1×2

1.0e-16 * 0.9586 0.6567

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The results are a decent estimate of the original data. We can then just use tubeplot() to make the whole solid model using whatever parameters we want, without needing to mess with offsets.

% reconstruct the tubeplot using new parameters
r = 0.2; % tube radius
nc = 12; % number of ring points
F = zeros(0,3); V = F;
for k = 1:ngroups
    [x y z] = tubeplot(newdata{k}.',r,nc);
    [F1,V1] = mesh2tri(x,y,z,'f');
    F = [F; F1+size(V,1)]; %#ok<*AGROW>
    V = [V; V1];
end

% write the new STL if desired
stlwrite(triangulation(F,V),'newtubeplot.stl')

% preview display
figure
patch('faces',F,'vertices',V,'facecolor','w','edgecolor','k');
view(3); view(-54,47); camlight; 
axis equal; grid on

function centroids = detubify(E,V)
% CENTROIDS = DETUBIFY(EDGES,VERTICES)
% given edge and vertex lists corresponding to a single
% triangulated tubeplot object, estimate the original data
% by calculating the centroids of the circumferential edges.
% also, plot everything just for sake of demonstration.

    % edges as coordinates [x y z]
    Ev = cat(3,V(E(:,1),:),V(E(:,2),:));
    
    % isolate short edges
    D = sqrt(sum(diff(Ev,1,3).^2,2)); % edge lengths
    shortedges = D <= 1.2*prctile(D,20); % just the short edges
    Es = E(shortedges,:); % get rid of long edges
    Esv = permute(Ev(shortedges,:,:),[3 1 2]);
    
    % i guess the edges need to start at 1 for graph() to work
    offset = min(Es(:)) - 1;
    Es = Es - offset;
    
    % find the connected groups
    G = graph(Es(:,1),Es(:,2));
    I = G.conncomp;
    
    % calculate centroid position for each group
    ngroups = numel(unique(I));
    centroids = zeros(ngroups,3);
    for k = 1:ngroups
        idx = find(I==k) + offset;
        centroids(k,:) = mean(V(idx(:),:),1);
    end
    
    % plot everything
    plot3(Esv(:,:,1),Esv(:,:,2),Esv(:,:,3),'b'); hold on
    plot3(centroids(:,1),centroids(:,2),centroids(:,3),'k')
    axis equal; grid on
end

Again, since we're using connectivity instead of vertex distance alone, we don't have to worry nearly as much about places where tubes are close to each other or close to themselves. Even in a case like the final model, where the objects intersect, the two objects will be processed independently, and we don't need to worry about intersecting rings of vertices being grouped with each other.