Find discontinuous areas in geometry

Question

Chenglin Li 2023-6-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1980574-find-discontinuous-areas-in-geometry

评论： Chenglin Li 2023-6-15

There is an stl file of a concave surface, known the vertices and joint indexes of each face, as well as the normal vectors of each face, to obtain the discontinuous region, that is, the concave face and the vertex information contained in the face, as shown in the figure below.

My idea is to find the concave part by judging that the Angle between the normal vectors of the adjacent faces is less than or equal to 90 degrees, and output the concave faces, but my program outputs all the faces, I don't know why

[VF] = stlread("ASCao.stl");
% normal is for every outward facing normal vector
[normals] = lin_compute_mesh_normal(VF.Points,VF.ConnectivityList);
% %Plot the mesh with the normals:
lin_plot_face_normal(VF.Points,VF.ConnectivityList,normals);
faces = VF.ConnectivityList;
[verticesOnFace] = findVerticesOnFace(2, VF.Points, faces);
% Find the concave face
concaveFaces = [];
for i = 1 : size(VF.ConnectivityList, 1)
    % Gets the face ids of the current and adjacent faces
    faceId = i;
    adjFaceIds = vertexAttachments(faces, faceId);
    
    % Determine whether the current face is recessed
    isConcave = false;
    for j = 1 : size(adjFaceIds, 2)
        if dot(normals(faceId,:), normals(adjFaceIds(j),:)) <= 0
            % The Angle between the normal vectors of adjacent faces is less than or equal to 90 degrees 
            % (that is, both sides are less inclined, and there may be a concave part).
            isConcave = true;
            break;
        end
    end
    
    % If the current face is concave, it is added to the result array
    if isConcave
        concaveFaces = [concaveFaces; faceId];
    end
end
% Output the ID of the concave face
disp(concaveFaces);
function [attachments] = vertexAttachments(faces, faceId)
% Computes the face ids of all faces adjacent to the specified face
numFaces = size(faces, 1);
adjFaceIds = [];
faceVertices = faces(faceId,:);
for i = 1 : numFaces
    sameVertexCount = 0;
    for j = 1 : 3
        if ismember(faces(i,j), faceVertices)
            sameVertexCount = sameVertexCount + 1;
        end
        
        if sameVertexCount == 2
            adjFaceIds(end+1) = i;
            break;
        end
    end
end
% Returns the face ID of all faces adjacent to the specified face
attachments = adjFaceIds(adjFaceIds ~= faceId);
end
function [verticesOnFace] = findVerticesOnFace(faceId, vertices, faces)
% Find all vertex coordinates on the specified face ID
faceVertices = faces(faceId,:);
verticesOnFace = [];
for i = 1 : size(vertices, 1)
    if ismember(vertices(i,:), vertices(faceVertices,:), 'rows')
        verticesOnFace(end+1,:) = vertices(i,:);
    end
end
end