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

What I did in the end was something like that:

import numpy as np
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt

def plot_random_points_on_circles(N: int, R: list, center: list):
    """
    Generates two sets of N random points on two circles with radii R and centers at center.
    Scales the points and calculates the convex hull of the scaled points.
    Plots the convex hull and the original points.

    :param N: Number of points to generate on each circle
    :type N: int
    :param R: List of radii for the two circles
    :type R: list
    :param center: List of x and y coordinates for the center of the circles
    :type center: list
    """
    x = np.zeros((len(R), N))
    y = np.zeros((len(R), N))
    theta = np.linspace(0, 2 * np.pi, N + 1)
    theta = theta[:-1]

    for i in range(len(R)):
        x[i,:] = R[i] * (np.random.randn() + 1) * np.cos(theta + np.random.rand()) + center[0]
        y[i,:] = R[i] * (np.random.randn() + 1) * np.sin(theta + np.random.rand()) + center[1]

    x = np.concatenate((x[0,:], x[1,:]))
    y = np.concatenate((y[0,:], y[1,:]))
    rSquared = x**2 + y**2
    q = rSquared / max(rSquared)**2
    xx = x / q
    yy = y / q

    k = ConvexHull(np.column_stack((xx, yy)))
    plt.plot(x[k.vertices], y[k.vertices], 'r-', x, y, 'b*')
    plt.show()

N = 20
R = [2, 6]
center = [0, 0]

plot_random_points_on_circles(N, R, center)

Answer 2

If I understand well your question, you want the largest volume inside another volume, without points in common between the two volumes.

The outer volume is built from a subset of the set of points. The obvious solution is to build the inner volume with the rest of points.

A volume from a set of points can be made in several ways. If the volume is not convex, then you need some more info (e.g. minimum angle between faces) because you get starred polytopo or cuasi-convex, or some other shape.

For convex volume I recomend the 3D Delaunay construction, with tetrahedra. The boundary is defined by the faces of "tets" that are not shared with other "tets".

Remove from the full set of points those belonging to the boundary: Each tet in boundary has a fourth point that does not lie on the boundary.

The inner volume is another Delaunay construction. Perhaps you only need the fourth points from the previous boundary-tets.

Answer 3

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.