Algorithm to extract line-like contour in 3d

Question

Hello people on the internet,

I'm currently searching for some kind of fast algorithm that allows me to extract curves in three dimensional space that arise as the intersection of two level sets of smooth functions, i.e. \begin{align*}\mathcal{C} = \{\vec{x} \in \mathbb{R}^3 | \; F(\vec{x})= l \; and \; G(\vec{x})=s \} \end{align*}

I was thinking that something similar to marching cubes might work, but not having a computer science background I'm struggling to make out if and how that would look. Does anyone have a suggestion on what I could do?

Edit: I currently use a very inefficient solution where I compute one isosurface with a preimplemented python method on a 3d grid, then check the isosurface for wanted values of the other function (walking along curves in slices of the surface and checking for crossings of the value) and collect all such points. In the end I group them with a clustering algorithm.

Edit2: In my specific case, the functions are defined with a few random parameters for some finite natural $N$, which are $\varphi_n \in [0,2\pi), A_n \in [0,1], \vec{k}_n \in S^2$ (the unit sphere) and are restricted to a three-dimensional cube of some size $a$. \begin{align*}F(\vec{r}) &= \frac{1}{T} \sum_{n=1}^N A_n \cos(2\pi \vec{k}_n \cdot \vec{r} + \varphi_n) \\ G(\vec{r}) &= \frac{\partial_z F}{|| \nabla F ||} \end{align*} for some normalization constant $T$. I want to work with the intersection curves (that I made out to generally be several disconnected ones) for different parameters.

Answer 1

Here is a generic approach:

1. Find a point on the intersection of these two curved planes, i.e., on the intersection of $F(\vec{x})=l$ and $G(\vec{x})=s$. Call that point $\vec{x}_0$.

2. Trace along the intersection curve, starting at $\vec{x}_0$, and continuing as far as you can. Then do the same starting from $\vec{x}_0$ in the opposite direction.

3. Find a new point $\vec{x}_0$ on the intersection of $F(\vec{x})=l$ and $G(\vec{x})=s$ that's not part of what you've traced so far, and go back to step 2.

To instantiate this, we need to know how to implement each of these steps.

To find a point on the intersection of $F(\vec{x})=l$ and $G(\vec{x})=s$, you can use any method for solving a system of nonlinear equations. There are many. If you're not familiar with the subject, one simple approach is to minimize the objective function

$$\Phi(\vec{x})= \|F(\vec{x})-l\|^2 + \|G(\vec{x})-s\|^2$$

using, e.g., gradient descent.

To trace along the intersection curve, you can use the gradient to help you. Suppose that $\vec{x}_0$ is a point on the intersection curve, and expand the Taylor series for $F$ and $G$. We get

$$F(\vec{x}) = F(\vec{x}_0) + \nabla F(\vec{x}_0) \cdot (\vec{x}-\vec{x}_0) + \dots \approx l + \nabla F(\vec{x}_0) \cdot (\vec{x}-\vec{x}_0)$$

and similarly

$$G(\vec{x}) \approx s + \nabla G(\vec{x}_0) \cdot (\vec{x}-\vec{x}_0).$$

So, for $\vec{x}$ to be on the intersection curve, we require $ \nabla F(\vec{x}_0) \cdot (\vec{x}-\vec{x}_0) = 0$ and $\nabla G(\vec{x}_0) \cdot (\vec{x}-\vec{x}_0)=0$ (assuming $x$ is very near $x_0$). Since $\nabla F(\vec{x}_0)$, $\nabla G(\vec{x}_0)$, and $\vec{x}_0$ are known and computable, this gives us two linear equations in three unknowns, so we can solve for $\vec{x}$ using a linear equation solver. In particular, compute the two vectors $\nabla F(\vec{x}_0)$ and $\nabla G(\vec{x}_0)$, then find a unit vector $\vec{d}$ that is orthogonal to both of them. Now you can trace out the intersection curve from $\vec{x}_0$ by preceding in the direction $\pm \vec{d}$; in other words, the line $\vec{x}_1=\vec{x}_0 + c \vec{d}$ is approximately on the intersection curve, when $c$ is close to 0. (If you want, you can generate $\vec{x}_1$ in this way, then "fine-tune" it by using it as the starting point for gradient descent with objective function $\Phi$ to compute $\vec{x}_1$ more precisely and to make sure you stay on the intersection curve.) This gives you a way to step along the intersection curve, one point at a time.

How to implement step 3? This might depend on the specific functions $F,G$. One approach is to modify the objective function $\Phi$ to penalize points that are near the part of the curve you've already found, e.g.,

$$\Phi(\vec{x}) = \|F(\vec{x})-l\|^2 + \|G(\vec{x})-s\|^2 + \sum_i {1 \over \|\vec{x}-\vec{x}_i\|^2},$$

where the $x_i$'s are points you've previously found (or some well-spaced-out subset of them). Another approach might be to do step 1 (gradient descent on $\Phi$, without the penalty term) starting from many different starting points (e.g., spaced on a coarse grid). There may be a better way in your specific situation; I don't know.

I expect this should be more efficient than enumerating all grid points and evaluating $F,G$ at every grid point, especially as the spacing of your grid points gets small. It does require you to be able to compute the gradients of $F,G$, but in your situation that seems to be no problem.