13

Take the following parametric equation of an implicit curve as an example:

$$ \left\{\quad \begin{array}{rl} x=& 9 \sin 2 t+5 \sin 3 t \\ y=& 9 \cos 2 t-5 \cos 3 t \\ \end{array} \right. $$

enter image description here

For an arbitrarily given point $P=(x_0,y_0)$, numerically, how to determine whether $P$ is inside the blue region defined by the curve or not?

LCFactorization
  • 2,112

3 Answers3

5

In general, you can determine whether a point lies within a closed curve by computing the winding number of the curve around the point. For a smooth closed curve $\gamma=(\gamma_x, \gamma_y):[0,1]\rightarrow \mathbb{R}^2$ (where $\gamma(0)=\gamma(1)$), the winding number around a point $(x_0, y_0)$ not on the curve is equal to the line integral $$ w_\gamma(x_0, y_0)=\frac{1}{2\pi}\oint_{\gamma}d\theta = \frac{1}{2\pi}\oint_{\gamma}\frac{(\gamma_x(t) - x_0) \gamma_y'(t) - (\gamma_y(t) - y_0) \gamma_x'(t)}{(\gamma_x(t) - x_0)^2+(\gamma_y(t)-y_0)^2}dt. $$ The winding number is an integer, and is zero for a point outside the curve and positive (negative) for a point that is enclosed by the curve in a counterclockwise (clockwise) orientation.

mjqxxxx
  • 43,997
  • 1
    This approach does not work in the case of selfintersecting curve and the curve under consderation is such a curve. Consider the interior point (approximately $(0,5)$) where the curve intersects with the $y$-axis. – user64494 Jun 02 '15 at 17:36
  • 1
    @user64494: It doesn't work for a point on the curve, that's true... Generically, if $(x_0, y_0)$ lies on an interior arc of the curve, then a small random displacement will push it off the curve without changing the signature of the winding number. – mjqxxxx Jun 02 '15 at 18:06
  • @ mjqxxxx : You are not right. Consider the point $(1,0)$ on the unit circumference ${(x,y):x^2+y^2=1}$. – user64494 Jun 02 '15 at 18:11
  • @ mjqxxxx : Also look at that page. – user64494 Jun 02 '15 at 18:37
  • @mjqxxxx I like this answer because it is relatively easy to implement numerically and as you said, considering small perturbations and approximations, points on curves do not really matter. In order to calculate the winding number, do I need to convert the parametric equation into implicit $F(x,y)=0$ form (or polar form) and how if needed? – LCFactorization Jun 02 '15 at 22:57
  • 1
    @ miqxxxx & @LCFactorization : The non-zero-winding number has a serious defect. There may be a lot of points in the interior of a closed self-intersecting curve (not only on the curve) which have zero-winding number: e. g. see an example in my answer. – user64494 Jun 03 '15 at 12:56
3

This is a star domain with respect to the origin by its construction. Because of this reason the nonlinear system $$\begin{cases}x_0 = \lambda(9\sin(2t)+5\sin(3t)),\\y_0= \lambda(9\cos(2t)-5\cos(3t)) \end{cases}$$ has a real solution $\{\lambda,\, t\}$ satisfying the restrictions $\{\lambda \ge 0, \lambda < 1,\, t \ge 0,\, t\le 2\pi\}$ if a point $(x_0,y_0)$ is inside the domain, not belonging to its boundary, and has no such solution in the opposite case. The only working method to solve that system is numeric. I use the DirectSearch package of Maple to this end. The DirectSearch package should be downloaded from Maple Application Center and installed in your Maple.

Let us begin from the point $A(1,-1)$:

 sol := DirectSearch:-SolveEquations([1 = lambda*(9*sin(2*t)+5*sin(3*t)), -1 = lambda*(9*cos(2*t)-5*cos(3*t))], {lambda >= 0, t >= 0, lambda <= 1, t <= 2*Pi}, AllSolutions, solutions = 1);

$$ \left[ \begin {array}{cccc} { 1.32446294720147540\times 10^{-26}}& \left[ \begin {array}{c} -{ 1.09245945623115402\times 10^{-13}} \\{ 3.61932706027801032\times 10^{-14}}\end {array} \right] \\&[\lambda= 0.349055727031787666,t= 1.23742356730272540]&77 \end {array} \right] $$

This output means the numerical solution $\lambda= 0.349055727031787666,t= 1.23742356730272540$ satisfies the first equation with the accuracy $- 1.09245945623115402\times 10^{-13} $ and the second equation up to $ 3.61932706027801032\times 10^{-14} $. The number $ 1.32446294720147540\times 10^{-26}$ equals the sum of the squared residuals. The number $77$ is the number of steps to solve the system. We see that the point $A$ is placed in the interior of the domain as $\lambda= 0.349055727031787666 < 1.$ Now we consider the point $B(5,5)$:

sol := DirectSearch:-SolveEquations([5 = lambda*(9*sin(2*t)+5*sin(3*t)), 5 = lambda*(9*cos(2*t)-5*cos(3*t))], {lambda >= 0, t >= 0, lambda <= 1, t <= 2*Pi}, AllSolutions, solutions = 1);

$$ \left[ \begin {array}{cccc} 0.804364869871176834& \left[ \begin {array}{c} - 0.0365718874608482736\\ 0.896117942526946542\end {array} \right] \\&[\lambda= 0.99999999999280198,t= 0.156602991662910751]&99\end {array} \right] $$

The above output means the system has no real solution as the residuals are quit big. The SolveEquations command can return not only exact solutions of the equation system but also any minimums of function $$F:=(5 - \lambda(9\sin(2t)+5\sin(3t)))^2+(5 - \lambda(9\cos(2t)-5\cos(3t)))^2.$$ If the residuals too large then the solution is not exact solution but it is only solution that minimizes the residuals. We can draw the conclusion that $B$ does not belong to the interior.

Finally, it should be noted that the winding number does not discriminate the interior points in this situation. See that page of William F. Basener, Topology and Its Applications: enter image description here

user64494
  • 6,049
0

An alternative approach is to shoot a ray from the point in question to infinity (for example, along a $y=const.$ line) and see how many times it intersects the curve. For an uneven number of intersections, the point is in the closed region.

For more details, see this answer.

Edit: Oh wait ... I see that this would exclude the inner pentagon of your star. Hum, no idea how to address that. Maybe the answer is still useful for people who look for an algorithm for non self-intersecting curves.

Amos Egel
  • 101