How to find the "interior boundary" / "interior convex hull" for a list of 3D points?

Question

I need some help with writing this algorithm.

For a given set of lines in space, I am trying to find the accessible volume when the origin (reference point) is 0.5,0.5,0.5. Currently, I do the following:

For each line, calculate the distance to the origin (0.5,0.5,0.5). Then, gather all these perpendicular distance points on all the lines into a list.

Now, I would like to calculate the "interior" (neither the boundary nor the convhull), because I want to evaluate the accessible volume for a ball centered at (0.5,0.5,0.5).

---

For example I would like to compute with my algorithm the green (internal line) in this simple example:

---

The configuration:

The closest points from the origin (0.5,0.5,0.5) to the lines

Only the points for whom I want the "internal boundary" be computed. Meaning the shape that bounds all the point either outside of the interior or on its boundary.

Here is the code for which I want something else rather than convhull:

close all
N=30;

S1 = cell(1, N);
for k = 1:N, S1{k} = rand(1, 3); end
S2 = cell(1, N);
for k = 1:N, S2{k} = rand(1, 3); end

M1 = cat(3, S1{:});
M2 = cat(3, S2{:});
M  = permute(cat(1, M1, M2), [1, 3, 2]);
figure
plot3(M(:, :, 1), M(:, :, 2), M(:, :, 3))
hold on
[x,y,z] = sphere;
x=x/100;y=y/100;z=z/100;
plot3(x+0.5,y+0.5,z+0.5)


figure 
hold on

NearestIntersectionPoints = cell(1,N); 
for k = 1:N 
    tmp1 = M(1,k,:); tmp2 = M(2,k,:);
    v1=tmp1(1,:); v2=tmp2(1,:);
    [d, intersection] = point_to_line([0.5,0.5,0.5], v1, v2);
    
    [x,y,z] = sphere;
    x=x/500;y=y/500;z=z/500;
    plot3(x+intersection(1),y+intersection(2),z+intersection(3))
    NearestIntersectionPoints{k} = intersection;

end


MHull = cat(3,NearestIntersectionPoints{:});
X=MHull(:,1,:); Y=MHull(:,2,:); Z=MHull(:,3,:);
X=X(:); Y=Y(:); Z=Z(:);
k = boundary(X,Y,Z);
hold on

plot3(X(k),Y(k),Z(k), 'r-*')


function [d,intersection] = point_to_line(pt, v1, v2)
      a = v1 - v2;
      b = pt - v2;
      d = norm(cross(a,b)) / norm(a);
      theta = asin(norm(cross(a,b))/(norm(a)*norm(b)));
      intersection = v2 + a * cos(theta);
      
end

Answer 1

I would do it like this:

1. tetrahedronize your pointcloud

   so create a mesh consisting of tetrahedrons where no tetrahedron intersect any other or contain any point in it. I do it like this:

   1. structures

   you need list of points,triangles and tetrahedrons. Each triangle need one counter which will tell you if it is used once or twice.

   2. create first tetrahedron

   by 4 nested loops through all points and check if formed tetrahedron does not contain any point inside. If not stop as you found your first tetrahedron. This is O(n^5) but as there are a lot of valid tetrahedrons it will never reach such high runtime... Now just add this tetrahedron to triangle and tetrahedron lists.

   3. find next tetrahedron

   now loop through all triangles that has been used once. for each form tetrahedron by using those 3 points used by it and find 4th point the same way as in #2. Valid tetrahedron must not contain any points in it and also must not intersect any existing tetrahedron in the list.

   To ensure whole volume will be filled without holes you need to prioritize the process by preferring tetrahedrons with more triangles already in list. So first search 4 triangles if no found than 3 etc ...

   For each found valid tetrahedron add it to the lists and look again until no valid tetrahedron can be formed ... The whole process is around O(n^2) so be careful with too many points in pointcloud. Also having normals for triangles stored can speed the tests a lot ...

2. outer boundary

   outer boundary consist of triangles in list which have been used just once

3. interior boundary

   interior gap tetrahedrons should be larger than all the others. So check their size against average size and if bigger they are most likely a gap. So group them together to lists. Each gap have only large tetrahedrons and all of them must share at least one face (triangle). Now just count the triangle usage for each group alone and all the triangles used just once will form your gap/hole/interior boundary/mesh.

If your point density is uniform you can adapt this:

• Finding holes in 2d point sets?

And create a voxel map of point density... voxels with no density are either gap or outer space. This can be used for faster and better selection of interior tetrahedrons.