Finding a point inside a surface

Question

I have a triangulation of a surface without boundary in $\mathbb{R}^3$. The triangulation gives a unit normal pointing outwards for each triangle. I need to find some point in the interior of the surface.

I will outline my general strategy. I let $V$ denote the set of vertices of the triangulation and compute a number that I will call the fineness $$ \delta = \min_{x, y \in V \mid x\neq y} \lVert x-y\rVert. $$ I pick some vertex $v$ and consider the point $$ p = v -\frac{\delta}{4}u, $$ where $u$ is some unit vector. It is possible to show that $p$ is closer to $v$ than any other vertex in the triangulation, so it suffices to choose $u$ such that it lies below all of the triangles containing $v$.

If $i$ indexes the triangles containing $v$ as a vertex and $n_i$ is the outer normal of the triangle $i$, choosing $u$ such that $\langle u, n_i\rangle >0$ for all $i$ will put $p$ below the surface.

We can map this onto a linear algebra problem by writing $u$ as a weighted combination of the outer normals: $u = \sum_i w_i n_i$. If all of the inner products of the outer normals were non-negative, we could solve this by taking the $w_i$ to be the components of the Perron–Frobenius eigenvector of the matrix $A_{i,j} = \langle n_i, n_j\rangle$. However, the surface can have regions of high curvature, so $A$ is not always non-negative.

A natural notion of the normal at $v$ is the angle weighted averaged of the normals. This lead me to try taking $w_i = \theta_i$, where $\theta_i$ is angle of the triangle $i$ at $v$. However, this alone does not guarantee $\langle u, n_i \rangle > 0$ for all $i$. This seems like some iterative approach is needed, where the weights corresponding to positive inner products should be decreased and those corresponding to non-positive ones should increase. I think one strategy should be to add to the weights of non-positive terms and then normalize all of the weights so that $\sum_i w_i = 1$. If all inner products are positive, we terminate the iterative procedure. Of course, we would have to make $u$ a unit vector at the end of this.

My questions are: Am I on the right track? Is there a simpler approach? If I am on the right track, what update rule for the weights will guarantee convergence?

Edit: The shapes that I am working with are definitely not convex. The data points are from intestinal organoids, which can roughly be shaped like octopi holding broccoli.

Answer 1

Let me suggest another approach, conceptually simple but maybe not as easy to implement as Matt F.'s algorithm.

Let $P$ be the polyhedron. Let $F$ be any face of $P$ with centroid $p$, and $\vec{n}$ the outward normal vector to $F$. Shoot a ray $r$ from outside of $P$ in direction $-\vec{n}$ through $p$. Ignore all intersections until the ray reaches $p$. The ray $r$ must now penetrate to the interior of $P$ on the immediate other side of $F$. Track $r$ until it hits another face of $P$. (It might hit a vertex or an edge of the face.) Let this first hit-point be $q$.

Then $(p+q)/2$ is necessarily strictly interior to $P$.

This method requires careful ray-triangle intersection (assuming all faces are triangles). I implemented this as part of code to answer point-in-polyhedron queries in Section 7.5 of Computational Geometry in C.

Answer 2

A short answer is: $1$) Pick a plane which intersects the surface in exactly one vertex $v$. $2$) From $v$ go slightly inward to some $p$, going so slightly that you can’t run into the planes connecting other vertices.

In more detail:

$1$) Consider all vertices with maximal $z$ coordinate, i.e. where $v \cdot \hat{k}$ is maximal.

• Choose one such $v$ on the convex hull of such points.
• Choose a horizontal vector $h$ such that $v\cdot h$ is maximal among all the vertices with maximal $z$ coordinate.
• Choose $w$ of the form $\hat{k}+\delta h$ so that $v \cdot w$ is maximal among all the vertices of the surface.

The diagram shows all this in one less dimension. The desired plane is the one through $v$ perpendicular to $w$.

$2$) Let $m$ be the mean of the inward unit normals over all the faces that include $v$. Choose $\epsilon$ such that $$(v \cdot w) + \epsilon(m \cdot w) > \max_{v'\neq v}(v' \cdot w)$$ Then $p=v+\epsilon m$ is a point inside the surface, since $m$ points inward, and none of the planes connecting other points go that far in the $w$ direction.