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

• I turned my comments into an answer. Jan 14 at 15:49

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.

• I think this is the better answer — 'canonically', if informally, doing line-polytope intersections tends to be more robust through faces than through vertices. Jan 14 at 19:33
• I think I like this approach better. Something I hadn’t considered is that some point inside of a triangle can be closer to some other vertex than any of the triangle’s vertices. This has implications for what “sufficiently small” means and may require checking. Passing a ray through faces should bypass this issue. Jan 15 at 23:13
• @quirkyquark: My implementation intersects the ray $r$ with all surface triangles, and sorts the intersections. That leads to the $\frac{1}{2}(p+q)$ calculation. Jan 15 at 23:51
• Yes, I understood that you looped over them and choose the first intersection after the centroid. I was saying that in my original conception of the problem, something like $\delta$ actually isn’t good enough to always work. Jan 15 at 23:56

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.

• For step 1, you can also choose a random $w$ and find the $v$ which maximizes $v \cdot w$. If the $v$ is unique, and it almost always will be, then you can use that $v$ and $w$. Jan 14 at 19:36