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

Question

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); 
grid() 

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?

Answer 1

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