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

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$.

• You could look into alpha shapes, you do have to adjust the 'alpha' parameter, but I imagine if you have a measure of how good a given polytope is you could find it with a few rounds of bisection. The good thing about alpha shapes is that you just need the delaunay triangulation, so it is reasonably efficient. – Matt D Feb 13 '18 at 21:48
• The conditions seem to keep changing. In revision 26, you define a convex good polyhedron to go on and ask for best convex and non-convex. Removing the convex hull seems a reasonable first step. A pity that convert to polar and add, in increasing distance from the origin, every point that doesn't make the polyhedron non-convex (one of the points formerly added interior) doesn't work: you might want to construct an example showing that. (I'm down to 3 (hull) + 1 (shared) + 2 (suboptimal) + 3 points.) – greybeard Feb 28 '18 at 4:42
• @greybeard I think I can stay with the convex hull thing. I am not sure I get your point. – 0x90 Feb 28 '18 at 4:52
• (I have difficulties following your use of set notation. You don't (explicitly) state $P \subset S$. For $P \subsetneq \text{convexhull}(S)$, it suffices to have $P = \text{convexhull}(S)$ (or to add one point (non-convex, if from $S$, too) or to take one point away) - I think you want $P \cap \text{convexhull}(S) = \emptyset$ (and, consequently, condition 3 to read none of the points in $S$ is contained within the interior of $P$).) – greybeard Feb 28 '18 at 5:26
• not sure I get your point the point of my 1st comment is that I can't see from your sketch that this problem isn't trivial. – greybeard Feb 28 '18 at 5:28

What I ended up doing is the following:

Given a set of points in $\mathbf{R}^3$ $S$, and a single point as origin, $\mathbf O$ do 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*')