Finding self-kissing points on a plane curve?

Question

I have a curve in the complex plane given by $$ f(t) = \sum_k r_k\exp(2\pi\mathrm{i}(t+\varphi_k)p_k). $$ Some of the parameters are specially chosen: $r_k>0$, $\sum_k r_k=1$, $p_k\in\mathbb{Z}$, $\mathrm{gcd}(p\ldots)=1$, $\varphi_k\in\mathbb{R}$. These restrictions are probably not very important for this question.

When I plot the curve, it will sometimes have some self-kissing points: points $t_1$ and $t_2$ where $f(t_1)=f(t_2)$ and $f'(t_1)\parallel f'(t_2)$. A self-kissing point doesn't have to be a self-crossing: at the self-intersection with parallel tangent vectors, the normal vectors can be pointing the opposite ways.

Can anybody suggest a good method for finding such points? I can't seem to think of a real-valued function that changes sign at the self-kissing point. I would like to be able to find these points once every frame in an animation, so efficiency can be an issue. If you can suggest a reference for this, that helps too.

Edit Direct root finding doesn't work because such a point would be a multiple root and the curve might come close to touching itself without doing so to numerical accuracy. Optimizing a function like $|f(t_1)-f(t_2)|^2$ doesn't work: it has many spurious critical points. Splitting the interval $[0,1]$ into $N$ parts and treating the curve as approximately quadratic inside each interval also doesn't work because there are two variables in $f(t_1)-f(t_2)$, and so requires $O(N^2)$ work, which is too much. (Although when looking for self-intersection points instead, this can work because of efficient line segment intersection algorithms.)

Answer 1

There exists a clever algorithm for doing this, based on ideas due to Cayley, using the Bézout matrix. I learned this from [1], which describes an algorithm for a generic version of this problem.

If we want to find all the critical points of the function $\|f(s)-f(t)\|^2$ (viewing $f$ as a function $\mathbb{R}\to\mathbb{R}^2$), a sufficient condition is that $$ P_s(t) := (f(s)-f(t))^\top f'(s) = 0, \qquad Q_s(t) := (f(s)-f(t))^\top f'(t) = 0, $$ which is a pair of simultaneous trigonometric equations in $s$ and $t$.

Form the Bézout matrix as the matrix with entries $B_{ij}(s)$ defined by the expansion $$ \frac{P_s(x)Q_s(y) - P_s(y)Q_s(x)}{x-y} = \sum_{i,j} B_{ij}(s) \phi_i(x) \phi_j(y). $$ Here $\phi_i$ are a suitable basis for the function $f,P,Q$. Typically they are either monomials, or Chebyshev polynomials, but presumably can also be elementary trigonometric functions.

Those values of $s$ for which the exists a solution $t$ to the pair of equations $P_s(t) = Q_s(t) = 0$ cause the matrix $B_{ij}(s)$ to be singular. In other words, if we expand the matrix $B_{ij}(s)$ in the basis $\phi_i(s)$ as $$ B_{ij}(s) = \sum_k B_{ij}^{(k)}\phi_k(s), $$ for a solution $s$ to exist it is necessary that $$ \det B_{ij}(s) = 0, $$ which is a polynomial eigenvalue problem (the polynomial's coefficients in the basis $\phi$ are the matrices $B_{ij}^{(k)}$). This is a well-known problem, and has standard, although possibly difficult algorithms, such as the polyeig function in Matlab. In particular it can turned into a linear eigenvalue problem for $s$ by constructing the block-companion matrix (if $\phi$ are monomials) or the block-colleague matrix (if $\phi$ are Chebyshev), and solving that linearized version of the PEP.

It may also be helpful to partition the $(s,t)$-square into some number of subsquares and use the Bézout matrix approach on individual subsquares, as this can substantially lower the degree of the polynomial eigenvalue problem and its linearized version.

References

1. Computing the common zeros of two bivariate functions via Bézout resultants by Yuji Nakatsukasa, Vanni Noferini, Alex Townsend