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

  • 2
    $\begingroup$ 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. $\endgroup$
    – Matt D
    Commented Feb 13, 2018 at 21:48
  • 1
    $\begingroup$ 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.) $\endgroup$
    – greybeard
    Commented Feb 28, 2018 at 4:42
  • $\begingroup$ @greybeard I think I can stay with the convex hull thing. I am not sure I get your point. $\endgroup$
    – 0x90
    Commented Feb 28, 2018 at 4:52
  • $\begingroup$ (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$).) $\endgroup$
    – greybeard
    Commented Feb 28, 2018 at 5:26
  • $\begingroup$ 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. $\endgroup$
    – greybeard
    Commented Feb 28, 2018 at 5:28

1 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

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]);


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;

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;

plot3(X,Y,Z, '.', 'MarkerSize',15);
hold on
[k,v] = convhull(X,Y,Z);

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

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);

enter image description here

Demo in 2D:

close all
R=[1, 5];
center=[0, 0];

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);

% 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);

enter image description here

  • 1
    $\begingroup$ With appropriate discretization I think you can convert this to maximum inscribed rectangle - for which there are algorithms linear of around log - I have my own which are ok - dont know the efficiency - the more advanced\efficient are also more complex to understand\implement etc $\endgroup$
    – gpasch
    Commented Nov 24, 2019 at 8:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.