# Algorithm to extract line-like contour in 3d

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.

• Interesting question! What kind of output are you looking for? How do you want the output to be expressed? The set $\mathcal{C}$ is an infinite set, so it's not possible to output a list of all elements of $\mathcal{C}$. Also, do you have an analytical expression for $F$ and $G$? Are you able to compute gradients for them? – D.W. Jan 15 '18 at 8:32
• @D.W. Thank you! I have no specific requirements for the output. In the end I want to maybe have something like a list/collection of discretized versions of the curves (maybe as graphs?) on this intersection (e.g. a sphere and a cylinder intersection should yield two curves) which I can then further analyze. In my case there are several such disconnected curves. And yes, that is surely possible! I wanted to pose the question in a general way, but my functions are conceptually easy. I will add them as an edit :) Also I forgot to mention: I want my functions with finite domain – Banana Jan 15 '18 at 8:47

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$$

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.

• Exactly what I was looking for, thank you for elaborating! – Banana Jan 16 '18 at 1:33
• I find myself having troubles with the initial points. My intersection consists of many disconnected components. Penalizing is too much of a computational effort (because I trace many curves, the sum becomes too large and has to be reevaluated all the time). I tried a different approach, finding zeros of F, d_z F and the x-component of the intersection curves tangent vector, but this also seems not efficient enough, taking minutes to find initial points even for relatively small domains. Do you know of any further algorithms to find initial points, by any chance? – Banana Jan 24 '18 at 3:24
• @Banana, I'm curious whether you've tried gradient descent starting from a bunch of starting points on a coarse grid? – D.W. Jan 24 '18 at 8:31
• I didn't, because I thought it wouldn't be reliable in finding all components – Banana Jan 24 '18 at 9:48