I have a self-intersecting polygon defined by $n$ points on the plane:

A <- t(matrix(c(
 0,   0,
 3,   0,
 3,  -1,
 2,  -1,
 2,   2,
 3,   2,
 3,   1,
 0,   1,
 0,   2,
 1,   2,
 1,  -1,
 0,  -1,
 0,   0), nrow=2));

The number of points is less $20$, coordinates are integer. You can see four points of self-intersections in $(1,1)$, $(2,1)$, $(2,0)$, $(1,0)$ on the figure:

plot(A, col='red', type= 'l', xlim=c(min(A[,1]),max(A[,1])),
      ylim=c(min(A[,2]),max(A[,2])), xlab='x', ylab='y'); 
points(A, col='black', pch = 22); 

enter image description here

Question I am looking for an algorithm to define the number and coordinates of self-intersections points.

How to find the number and coordinates of self-intersections points?

  • $\begingroup$ So are the points of your polygon $(0,0)$, $(3,0)$, $(3,-1)$, $(2,-1)$, $(2,2)$, $(3,2)$, ($3,1)$, $(0,1)$, $(0,2)$, $(1,2)$, $(1,-1)$, $(0,-1$, and $(0,0)$ again? $\endgroup$ – Hagen von Eitzen Apr 23 '19 at 15:56
  • 1
    $\begingroup$ Sweep line algorithm, Bentley-Ottmann algorithm $\endgroup$ – Jean-Claude Arbaut Apr 23 '19 at 16:00
  • $\begingroup$ @hagenvoneitzen, a moving point on the polygon must start to move backward to (0,0). $\endgroup$ – Nick Apr 23 '19 at 16:03
  • $\begingroup$ If you have on hand a segment-segment intersection algorithm, then there is an easy $O(n^2)$ algorithm: Intersect each segment with every other. If want speed and willing to code a more complex algorithm, then as Jean-Claude Arbaut and Sharat Chandrasekhar suggest, go with Bentley-Ottmann. $\endgroup$ – Joseph O'Rourke Apr 23 '19 at 22:34
  • 1
    $\begingroup$ If you have no more than 100 vertices, then the $O(n^2)$ algorithm that I’ve outlined below works well and is also applicable in 3- (or for that matter, $n-$) dimensions. It should not involve more than about 30 lines of code. $\endgroup$ – Sharat V Chandrasekhar Apr 24 '19 at 3:27

An approach that will work for a general self-intersecting polygon in 3-dimensional space is as follows:

For each adjacent vertex pair ${i,j}$, loop over every other set of adjacent vertex pairs ${m,n} : m\ne i\ne j,\hspace {0.5 cm} n\ne i\ne j$ and invoke the following procedure:

The position vector of any point along the line joining vertices $i$ and $j$ is

${\bf \vec x} = {\bf \vec x_i} + \xi L_{ij} \bf \vec q \tag{EQ. 1}$

and likewise that of any point along the line joining vertices $m$ and $n$ is

${\bf \vec x} = {\bf \vec x_m} + \eta L_{mn}\bf \vec r \tag{EQ. 2}$

where $$ L_{ij}=\|{\bf \vec x_j -\bf \vec x_i}\|\\ L_{mn}=\|\bf \vec x_n -\bf \vec x_m\|$$ and $${\bf \vec q}=\frac{1}{L_{ij}}({\bf \vec x_j} -{\bf \vec x_i})\\ {\bf \vec r}=\frac{1}{L_{mn}}({\bf \vec x_n -\bf \vec x_m})$$

At the intersection, we have

${\bf \vec p} + {\bf \vec q}L_{ij}\xi -{\bf \vec r} L_{mn}\eta =0 \tag{EQ. 3}$

where ${\bf \vec p}={\bf \vec x_i} -{\bf \vec x_m}$. Now you can find a unit vector $\bf\vec s$ orthogonal to $\bf\vec r$ as a linear combination of the unit vectors $\bf\vec q$ and $\bf\vec r$ as follows:

$\bf\vec s = \alpha \bf\vec q + \beta \bf\vec r\tag{EQ. 4}$

From the condition $\bf\vec s\cdot\bf\vec r=0$ we have

$\alpha \bf\vec q\cdot\bf\vec r +\beta=0\tag{EQ. 5}$

and from the stipulation that $\bf\vec s$ be a unit vector, we have

$\alpha^2 \bf\vec q\cdot\bf\vec q +\beta^2\bf\vec r\cdot\bf\vec r+2\alpha\beta\bf\vec q\cdot\bf\vec r=1\tag{EQ. 6}$

From EQS. 5 and 6, $\beta$ can be eliminated giving an equation for $\alpha$ whereupon only the positive root is relevant. Following this, $\beta$ is given by EQ. 5, and $\bf\vec s$ from EQ. 4

Now, taking the inner product of EQ. 3 with $\bf\vec s$ gives

${\bf \vec s \cdot\bf \vec p} + L_{ij}(\bf \vec s \cdot\bf \vec q)\xi =0\tag{EQ. 7}$

from which you can evaluate $\xi$. If $\xi<0$ or $\xi>1$, then you abort the procedure (the intersection is outside the lines joining the vertices $i$ and $j$), and proceed to the next pair.

Otherwise, calculate $\bf\vec x$ from EQ. 1 and then evaluate $\eta$ from

$\eta= \frac{1}{L_{mn}}(\bf\vec x-\bf\vec x_m)\cdot\bf\vec r \tag{EQ. 8}$

If $0\le\eta\le1$, then $\bf\vec x$ is an intersection.

This is a brute-force approach of $O(n^2)$. A more efficient alternative may be the Bentley-Ottman Algorithm

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.