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.

  • $\begingroup$ I turned my comments into an answer. $\endgroup$
    – Matt F.
    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.

  • $\begingroup$ I think this is the better answer — 'canonically', if informally, doing line-polytope intersections tends to be more robust through faces than through vertices. $\endgroup$ Jan 14 at 19:33
  • $\begingroup$ 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. $\endgroup$ Jan 15 at 23:13
  • $\begingroup$ @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. $\endgroup$ Jan 15 at 23:51
  • 1
    $\begingroup$ 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. $\endgroup$ 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.

enter image description here

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.

  • $\begingroup$ 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$. $\endgroup$
    – Matt F.
    Jan 14 at 19:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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