6

I have a project where I need to know the exact minimal distance between a point $(e, f)$ and a sine wave $y = a + b\cdot\sin(cx+d)$

Is there any way of calculating this? If not, is there a way to approximate this?

Thanks in advance!

Edit

Following Ross' and Ben's answers, The distance between a point $(e,f)$ and a sine wave $y = a + b\cdot\sin(cx+d)$ can be calculated by defining it as a distance between two points: $$D=\sqrt{(x-e)^2+(a+b\cdot\sin(cx+d)-f)^2}$$ The $x$ where the distance is minimal can then be calculated by $E=D^2$ and then solving $E'=0$ $$E'= 2(bc\cdot\cos(cx+d)(a+b\cdot\sin(cx+d)-f)-e+x)=0$$ $$x=e-bc\cdot\cos(cx+d)(a+b\cdot\sin(cx+d)-f)$$ This last formula can only be answered exactly in certain cases.

4 Answers4

2

I have had success approximating the distance to the sine curve by using a Chebyshev polynomial of the first kind. These are polynomials that satisfy the following equation:

$$T_n(\cos(x)) = \cos(nx) $$

Somewhat surprisingly, functions that satisfy this equality are polynomials, and the degree is equal to the value $n$. The intuition behind the equality is that if you draw a sine wave onto the surface of a half-cylinder and project the drawing to a plane, you will have the plot for a Chebyshev polynomial. The following diagram from this source (Creative Commons Attribution 4.0 International) illustrates this transformation.

We can use the function $T_2(x) = 2x^2 - 1$ to help us approximate the distance to the sine function. Because the Chebyshev polynomials are more easily related to the cosine curve, we'll get the distance to the cosine curve instead. This is just a shifted sine, so we just need a change of coordinates.

Given a point $(s,t)$ we would like to find a point $(x,y)$ on the cosine curve with frequency $f$. We imagine that our point $(s,t)$ is sitting on the surface of the half-cylinder, together with one cycle of our cosine curve. We can map the point down to the plane, then get the closest point on the aforementioned Chebyshev polynomial $T_2$. Finding the closest point on a parabola is a solved problem in the closed form (requires the solving of a cubic equation), so we can avoid iterative methods.

Once we have the closest point on the parabola, we can map it back to the surface of the half-cylinder. This will be our approximately closest point on the surface of the cosine curve. Care must be taken when mapping back to the half-cylinder, since the closest point on the parabola may be outside the diameter of the cylinder. This can be fixed by clamping its x coordinate to the bounds. Because of this clamping, we also must recalculate the resulting y position to be on the cosine curve. The following pseudocode computes the approximate closest point:

s = sin(s*f)/f;
(x, y) = closest point on parabola 2(f*x)^2 - 1 to (s, t)
x = asin(clamp(x*f,-1.,1.))/f
y = -cos(2f*x)

Unfortunately, this estimate really only works well when near the peak of the cosine curve, and breaks down at the two opposite peaks of the cycle and begins overestimating. To avoid this problem, we can take two estimates for the two peaks on either side of the given point $(s,t)$, then take the closest point of the two.

One thing to note is because we are using $T_2$, the resulting distance will be for the curve $\cos(2fx)$. We get a doubling of the frequency. This can be resolved by simply using a change of variables.

The resulting approximation is very good, but we can use a single Newton's method step to finalize the approximation to make it even more accurate. I am doubtful that doing more than one optimization step will improve things.

I have written an article on my website describing in further detail this method. I have also implemented it on the website shadertoy.

Blackle Mori
  • 121
  • 4
2

Let your point be $(e,f)$ so we don't reuse $x,y$. An approximate approach is as follows: First, find the limits of the half wave of interest. Let's say we are above the curve. You want the local minimum nearest $e$ and the local maximum on the other side of $e$. The perpendicular at a point has slope that is the inverse reciprocal of the derivative. The perpendicular from the maximum will be vertical and on one side of $(e,f)$, the perpendicular from the minimum will be on the other side of $(e,f)$ Call up your favorite one-dimensional root finder to find the $x$ value where the perpendicular goes through $(e,f)$. We have bracketed the root, so it should be easy to find. Now find the distance from $(x,y)$ to $(e,f)$

Ross Millikan
  • 374,822
  • @Macavity: thanks. I have changed to $e,f$ – Ross Millikan Oct 04 '13 at 18:04
  • Thanks for the quick answer! So if I am right, you said to take one half of a sine wave, calculate its derivative and then check at which x the value of the derivative is equal to the peripendicular of the slope of an imaginary line between point (e, f) and the x and the value y at that x? – Wouter Raateland Oct 04 '13 at 18:34
  • That is correct. As the sign of the error changes, a root finder should get there very quickly. They are in many numeric packages and described in any numerical analysis text. – Ross Millikan Oct 04 '13 at 19:43
0

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:

  1. 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)
  1. 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
Armatus
  • 101
0

Following Ross' answer, I will use $(e,f)$ to denote your point.

First, some simplifications:

  • The minimal distance between $(e,f)$ and $y = a + b\sin(cx + d)$ is the same as the minimal distance between $(e,f - a)$ and $y = b\sin(cx + d)$. So if I can solve the problem with $a = 0$, I can easily solve it for other cases as well.
  • Similarly, the minimal distance between $(e,f)$ and $y = b\sin(cx + d)$ is the same as the minimal distance between $(e + d/c,f)$ and $y = b\sin cx$, as long as $c \not= 0$. If $c = 0$ then the problem is not very difficult.
  • Again, assuming $c \not= 0$, the minimal distance between $(e,f)$ and $y = b\sin cx$ is the same as $1/c$ times the minimal distance between $(ec,fc)$ and $y = bc \sin x$.

So it's enough to solve the problem in the form "find the minimal distance between $(e,f)$ and $y = a\sin x$".

The points on the sine curve are of the form $(x,a\sin x)$. The distance to such a point is \[D = \sqrt{(x - e)^2 + (a\sin x - f)^2}.\] When this reaches a minimum, its derivative will be zero. As Jack M points out in the comments, you can throw away the square root, if you like: $D$ will reach its minimum exactly when $D^2$ does.

I've no patience for algebra so I leave the rest to you (unless you really get stuck), but I suspect you'll not end up with a neat expression, because there'll be too many trig functions in there. There are various ways to proceed from there, but I don't currently have time to go into them.

Ben Millwood
  • 14,211
  • Thanks, this will get me there! I will update my question when I finished this piece of algebra. – Wouter Raateland Oct 04 '13 at 18:57
  • 2
    Do you need to keep the square root on $D$? If you minimize $D^2$ you'll have minimized $D$. – Jack M Oct 04 '13 at 19:00
  • Having both $x$ and $\sin x$ will (almost certainly) prevent an algebraic solution-you are into numerics. The choice between minimization and root finding is not easy. Root finding may be more accurate in $x$ because it avoids an $\epsilon^2$ problem, but the distance should be equally accurate. – Ross Millikan Oct 04 '13 at 19:46