finding intersecting points of two 3D bodies

I have a set of coordinates (stored in vectors Xa, Ya, Za) that describe a 3D body A (including the coordinates "inside"). A second set of coordinates (Xb, Yb, Zb) describes a second 3D body B. I am trying to find the points of A that "intersect" with B, preferably including some margin. I have tried using convhulln and inhull (see code below), but inhull returns false for all elements of A. I also tried delauneyTriangulation and convexHull, but I cannot get that to work either. I would be very grateful for any ideas/comments/corrections!

C = convhulln([Xb Yb Zb]);
W = [Xa, Ya, Za];
in = inhull(W,C);

[Ankit]

function result = checkIntersection(vertices1, faces1, vertices2, faces2)
% vertices1 and vertices2 are N-by-3 arrays of vertex coordinates
% faces1 and faces2 are M-by-3 arrays of indices into vertices1 and vertices2

% Check if any of the edges or faces of the first geometry intersect with the second geometry
for i = 1:size(faces1, 1)
    face1 = vertices1(faces1(i, :), :);
    for j = 1:size(faces2, 1)
        face2 = vertices2(faces2(j, :), :);
        [intersect, ~, ~] = intersectLineTriangle(face1, face2);
        if intersect
            result = true;
            return;
        end
    end
end

% Check if any of the vertices of the first geometry are inside the second geometry
for i = 1:size(vertices1, 1)
    vertex1 = vertices1(i, :);
    [inside, ~] = isPointInMesh(vertex1, vertices2, faces2);
    if inside
        result = true;
        return;
    end
end

% Check if any of the vertices of the second geometry are inside the first geometry
for i = 1:size(vertices2, 1)
    vertex2 = vertices2(i, :);
    [inside, ~] = isPointInMesh(vertex2, vertices1, faces1);
    if inside
        result = true;
        return;
    end
end

result = false;
end

function [intersect, t, u] = intersectLineTriangle(line, triangle)
% line is a 2-by-3 array representing a line segment
% triangle is a 3-by-3 array representing a triangle
% intersect is true if the line intersects the triangle, false otherwise
% t and u are the barycentric coordinates of the intersection point

edge1 = triangle(2,:) - triangle(1,:);
edge2 = triangle(3,:) - triangle(1,:);
normal = cross(edge1, edge2);
if norm(normal) == 0
    % triangle is degenerate
    intersect = false;
    t = NaN;
    u = NaN;
    return;
end
dir = line(2,:) - line(1,:);
w0 = line(1,:) - triangle(1,:);
a = -dot(normal,w0);
b = dot(normal,dir);
if abs(b) < eps
    % line and triangle are parallel
    if a == 0
        % line lies in triangle plane
        intersect = true;
        t = 0;
        u = 0;
        return;
    else
        % line is outside triangle plane
        intersect = false;
        t = NaN;
        u = NaN;
        return;
    end
end
t = a/b;
if t < 0 || t > 1
    % intersection point is outside line segment
    intersect = false;
    u = NaN;
    return;
end
I = line(1,:) + t*dir;
w = I - triangle(1,:);
uu = dot(w,edge1)/dot(edge1,edge1);
vv = dot(w,edge2)/dot(edge2,edge2);
if uu >= 0 && vv >= 0 && uu + vv <= 1
    % intersection point is inside triangle
    intersect = true;

Answer 1

[John D'Errico]

I'm not sure why inhull would have been expected to work, as that is not the purpose of that tool.

Are the objects convex at all? Slightly harder if not, as then you will probably need to learn about alpha shapes.

If your goal is merely to fin where the perimeter boundaries interesct, I'd use a tool like a convexhull (if convex is true) or an alphashape if not. That will help you to create the boundary as a polygon. Then find the intersection points of two polygons. Doug Schwarz's tool, intersections (found on the file exchange) will help you there.

As far as "some margin" is concerned, I have no idea what that means. Fuzzy, vague words can be impossible to answer.

[debbie]

The objects aren't necessarily convex. My goal is to delete or highlight all points of A that intersect with B. What I meant by "margin" is that I want to include those points of A that are close to B, but not "yet" intersecting. Ideally I would like to set a tolerance or something comparable when it comes to testing whether a point intersects with B or not. Sorry for the poor description.

[debbie]

I tried again using alphaShape (see code below) and it works perfectly. Thanks for your help!

shp = alphaShape(Xb, Yb, Zb); % create alphaShape
tf = inShape(shp,Xa, Ya, Za); % test if points of A are inside B

[Ankit]

function result = checkIntersection(vertices1, faces1, vertices2, faces2)

% vertices1 and vertices2 are N-by-3 arrays of vertex coordinates

% faces1 and faces2 are M-by-3 arrays of indices into vertices1 and vertices2

% Check if any of the edges or faces of the first geometry intersect with the second geometry

for i = 1:size(faces1, 1)

face1 = vertices1(faces1(i, :), :);

for j = 1:size(faces2, 1)

face2 = vertices2(faces2(j, :), :);

[intersect, ~, ~] = intersectLineTriangle(face1, face2);

if intersect

result = true;

return;

end

end

end

% Check if any of the vertices of the first geometry are inside the second geometry

for i = 1:size(vertices1, 1)

vertex1 = vertices1(i, :);

[inside, ~] = isPointInMesh(vertex1, vertices2, faces2);

if inside

result = true;

return;

end

end

% Check if any of the vertices of the second geometry are inside the first geometry

for i = 1:size(vertices2, 1)

vertex2 = vertices2(i, :);

[inside, ~] = isPointInMesh(vertex2, vertices1, faces1);

if inside

result = true;

return;

end

end

result = false;

end

function [intersect, t, u] = intersectLineTriangle(line, triangle)

% line is a 2-by-3 array representing a line segment

% triangle is a 3-by-3 array representing a triangle

% intersect is true if the line intersects the triangle, false otherwise

% t and u are the barycentric coordinates of the intersection point

edge1 = triangle(2,:) - triangle(1,:);

edge2 = triangle(3,:) - triangle(1,:);

normal = cross(edge1, edge2);

if norm(normal) == 0

% triangle is degenerate

intersect = false;

t = NaN;

u = NaN;

return;

end

dir = line(2,:) - line(1,:);

w0 = line(1,:) - triangle(1,:);

a = -dot(normal,w0);

b = dot(normal,dir);

if abs(b) < eps

% line and triangle are parallel

if a == 0

% line lies in triangle plane

intersect = true;

t = 0;

u = 0;

return;

else

% line is outside triangle plane

intersect = false;

t = NaN;

u = NaN;

return;

end

end

t = a/b;

if t < 0 || t > 1

% intersection point is outside line segment

intersect = false;

u = NaN;

return;

end

I = line(1,:) + t*dir;

w = I - triangle(1,:);

uu = dot(w,edge1)/dot(edge1,edge1);

vv = dot(w,edge2)/dot(edge2,edge2);

if uu >= 0 && vv >= 0 && uu + vv <= 1

% intersection point is inside triangle

intersect = true;