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

example

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
    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 at 21:48
  • 1
    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 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 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 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 at 5:28
up vote 0 down vote accepted

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

enter image description here


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*')

enter image description here


Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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