Finding a point inside a surface 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.
 A: 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.
A: 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.
