How to find the supremum over all the "good" (interior) polytopes for a given set of 3D points?

Question

Let $S \subset \mathbf{R}^3$ be a set of points in 3D and let $O=(x_0,y_0,z_0)$ be the origin/point of reference.

We consider a convex polytope $P$ good / interior if:

1. $P$ is wholly contained within the interior of the convex hull of $S$: $P \subsetneq \text{convexhull}(S)$.

2. $O$ is contained within the interior of $P$, namely: $O \in P$.

3. none of the points in $S$ are in $P$ namely $S \cap P = \emptyset$.

Out of all good convex polytopes, we want to find either the maximum (supremum) of their volumes, i.e. $\sup \text{Vol}(P)$ where $P$ ranges over good polytopes.

---

Here is a 2D modest illustration (the green polytope might not be really the optimal):

The black dots are the set $S$, the origin $O$ is in the middle. The red lines go through the points in $S$ which define the boundary of $S$ (not a good polytope) and the green lines go through a good polytope.

Is there an efficient algorithm to compute this best convex and non-convex polytopes for a given set $S$ of points and origin $O$?

I would like to formulate this problem as an optimization problem or to come up with some algorithm using triangulation, $\alpha$-shapes or $\beta$-skeleton geometric transformation + convex hull, linear programming, maybe working in the dual space? etc.

One idea which I have is "pumping down" the convex hull of Delaunay Triangulation:

DT $\leftarrow$ delaunayTriangulation($S \cup O$)

CH $\leftarrow$ convexHull(DT)

while Not empty(DT) and $O \notin CH$:

$\quad$ i. previousDT $\leftarrow$ DT

$\quad$ ii. update DT $\leftarrow$ DT $\setminus$ CH

return previousDT

---

EDIT: I think what I am looking for in a one-liner of math symbols is:

$$\text{argmax}_{P} \text\{\text{Vol}(P) \, | \, O \in P \setminus \partial P \, ; \, S\cap P = \emptyset \, ; \, P \subsetneq \text{convexhull}(S) \} $$

Where $\partial P$ is the boundary of a polytope $P$.

Answer 1

What I ended up doing is the following:

Given a set of points in $S \subset \mathbf{R}^3$, and a single point as origin, $\mathbf O \in \mathbf{R}^3$ preform the following:

1. Init empty set $S'$.
2. $r = \max\{ ||\mathbf x - \mathbf O ||^2_2 \, , \,\, \mathbf x \in S\}$
3. Create a mapping $F(x): S\to S'$ where $\mathbf x \in S \to \frac{r^2}{||\mathbf x - \mathbf O\|^2} (\mathbf x- \mathbf O)$
4. $CH_{S'} \leftarrow $ the points in $S'$ which form convex hull of $S'$ (the sphere inverted set $S$)
5. Init the interior points set, $IP$ $\leftarrow$ $F^{-1}(s'\in CH_{S'})$ to find the points in $S$ whose sphere inversion points were in $CH_{S'}$.
6. Use Delaunay triangulation to find the volume of the the convex hull of $IP$

---

Matlab code:

%% Init
N=20;

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

M1 = cat(3, S1{:});
M2 = cat(3, S2{:});
M  = permute(cat(1, M1, M2), [1, 3, 2]);

O=[0.5,0.5,0.5];

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(O, v1, v2);
    NearestIntersectionPoints{k} = intersection;
end

MHull = cat(3,NearestIntersectionPoints{:});
X=MHull(:,1,:); Y=MHull(:,2,:); Z=MHull(:,3,:);
R=[ X(:) Y(:) Z(:)];

X=[X(:)]; Y=[Y(:)]; Z=[Z(:)];
DT = delaunayTriangulation(X,Y,Z);

rSquared = (X-O(1)).^2 + (X-O(2)).^2 + (Z-O(3)).^2;

q = rSquared/max(rSquared);
xx=X./q; yy=Y./q; zz=Z./q;

figure 
plot3(X,Y,Z, '.', 'MarkerSize',15);
hold on
[k,v] = convhull(X,Y,Z);
v
trisurf(k,X,Y,Z,'Facecolor','red','FaceAlpha',0.1)

figure
[k,v] = convhull(xx,yy,zz);

trisurf(k,X,Y,Z,'Facecolor','red','FaceAlpha',0.1)
hold on
plot3(X,Y,Z, '.', 'MarkerSize',15)
uniqueK = unique(k)
DTinterior = delaunayTriangulation(X(uniqueK),Y(uniqueK),Z(uniqueK));
DTexterior = delaunayTriangulation(X,Y,Z);

[kin, vin] = DTinterior.convexHull;
[kext, vext] = DTexterior.convexHull; 

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

---

---

Demo in 2D:

close all
N=20;
R=[1, 5];
center=[0, 0];
x=zeros(length(R),N);
y=zeros(length(R),N);
theta=linspace(0,2*pi,N+1);
theta=theta(1:end-1);

for i=1:length(R)
    x(i,:)=R(i)*(randn+1)*cos(theta+rand) + center(1);
    y(i,:)=R(i)*(randn+1)*sin(theta+rand) + center(2);
end

% Concat all the arrays into one set
x=[x(1,:) x(2,:)]; y=[y(1,:) y(2,:)];
rSquared=x.^2 + y.^2;
q = rSquared/max(rSquared)^2;
xx=x./q; yy=y./q;

k = convhull(xx,yy);
plot(x(k),y(k),'r-',x,y,'b*')

---

---