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.



Finding a point inside a surface
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](http://cs.smith.edu/%7Ejorourke/books/compgeom.html).


