I'm facing the same problem but with a different trig function and the added benefit that we only consider $-b<x<b$: $f(x) = a \cdot cos(arcsin(\frac{x}{b}))$
There are multiple ways of arriving here, but eventually we end up as you describe looking for an $x$ where the derivative of the distance function is
$$x + f'(x) (f(x)-y1)) -x1 = 0$$
with $(x_1,y_1)$ being the point whose distance we are looking for. This equation is the derivative of the distance function D above. Now to actually substitute and solve this.
First we determine
$$f'(x) = \frac{-a^2x}{b^2 f(x)}$$
Then we substitute
$$x + \frac{-a^2x}{b^2 a \cdot cos(arcsin(\frac{x}{b}))} (a \cdot cos(arcsin(\frac{x}{b}))-y_1) -x_1 = 0$$
and "simplify" to
$$x + \frac{y_1ax-a^2x \cdot cos(arcsin(\frac{x}{b}))}{b^2 \cdot cos(arcsin(\frac{x}{b}))} -x_1 = 0$$
$$x \cdot cos(arcsin(\frac{x}{b}))+ \frac{y_1a}{b^2}x-\frac{a^2 \cdot cos(arcsin(\frac{x}{b}))}{b^2}x -x_1 \cdot cos(arcsin(\frac{x}{b})) = 0$$
$$(x -\frac{a^2}{b^2}x -x_1) \cdot cos(arcsin(\frac{x}{b})) = -\frac{y_1a}{b^2}x$$
At this point we can use the fact that $cos(arcsin(x)) = \sqrt{1-x^2}$ and square both sides:
$$(x -\frac{a^2}{b^2}x -x_1)^2 \cdot (1-(\frac{x}{b})^2) = (-\frac{y_1a}{b^2}x)^2$$
I don't think I need to tell you that working this out results in an absolutely delightful polynomial whose degree turns out to be 4. This means we will obtain four roots when we are finally done, and will then have to validate which of them actually is the true $x$ we were looking for. Mathematically, this is pretty much as far as I could get. Wolfram Alpha (wolframalpha.com) when you give it the generic quartic polynomial ($ax^4+bx^3...$) in fact gives you the four equations for the four roots, but they are horrendously long and in our case, those $a, b, c...$ of course abbreviate terms which themselves consist of sums and differences involving the $a$ and $b$ in our original $f$.
However, a mathematical dead end was no option for me, since my program needs a value for this. In case anyone ever finds themselves in the same situation, I wanted to mention the solutions I found to get to a result. My program is in python, but I hope that whatever you are working in offers implementations of the same approaches:
- We take advantage of our knowledge that the distance function always has exactly one minimum. Use a minimizer algorithm to approximate the minimum of the distance function D. I implemented this with scipy.optimize.minimize_scalar. This worked but running this minimizer for 12000 points does take a while, so I looked for a better option.
optimize.minimize_scalar(
dist_to_vitelline,
args=(x1,y1,a,b),
bounds=(-b,b),
method="Bounded"
)
def dist_to_vitelline(x0,x1,y1,a,b):
return np.sqrt((y1-a * np.cos(np.arcsin((x0 / b))))**2+(x1-x0)**2)
- We instead jump in further along our earlier analysis and instead use an algorithm to find the root of the derivative of the distance function. Scipy offers a specialised function for this, optimize.root_scalar. This was quite a bit faster (maybe 10x or so?) and can maybe be sped up further by specifying
method='newton' and also supplying the second derivative of the distance function as an additional argument fprime=dxdx_dist_to_vitelline.
optimize.root_scalar(
dx_dist_to_vitelline,
args=(x1,x2,a,b),
bracket=(-b,b),
x0=x1
)
def dx_dist_to_vitelline(x0,x1,y1,a,b):
fx0 = a * np.cos(np.arcsin(x0/b))
dx_fx0 = (a**2 / b**2) * (-1 * x0 / fx0)
return x0 + dx_fx0 * (fx0 - y1) - x1