Find shortest path which intersects a point in a polygon
The main problem here is how to handle the non convexity. Given the rock 2D shape as a point sequence
$$ S = \{p_k\}, k = 1,\cdots,n $$
we can construct the segments
$$ s_k = \lambda_k p_k + (1-\lambda_k) p_{k+1},\ \ \ 0 \le \lambda_k\le 1 $$
and $s_n = \lambda_n p_n + (1-\lambda_n) p_1$
Now given a point $p_0$ in the $S$ interior, we define a generic line containing $p_0$ as
$$ L_j = p_0 + \lambda_0 v_j,\ \ \ v_j = (\cos t_j, \sin t_j) $$
and then given a direction $t_j$ we determine all possible intersections between $L_j$ and $\{s_k\}, \ \ k = 1,\cdots n$: thus given a $t_j$ we consider as associated interior distance
$$ d_j = \min{{\lambda_0}_k^+}-\max{{\lambda_0}_k^-} $$
where $\lambda_0^-,\lambda_0^+$ indicates if the intersection result gives a $\lambda \le 0$ or $\lambda \ge 0$ respectively. Finally we register for each $t_j$ the minimum $d_j$ found obtaining this way the result. The sweep made with $t_j$ can be choose to the precision needed.
Follows a MATHEMATICA script to solve with specified precision this problem. Here data is the set of points defining the rock profile, and p0 is the interior point. The algorithm performs a sweep from $0$ to $360$ degree, calculating the shortest distance along all possible intersections.
s[p1_, p2_, lambda_] := lambda p1 + (1 - lambda) p2
l[p0_, lambda_, v_] := p0 + lambda v
v = {Cos[t], Sin[t]};
data = {{0, 2.5}, {2.0, 1.8}, {3, 0.5}, {7.0, 10}, {2, 6.0}, {2.5, 8.0}, {0.5, 7.0}};
p0 = {1, 5};
data = AppendTo[data, data[[1]]];
n = Length[data] - 1;
segs = Table[s[data[[k]], data[[k + 1]], Subscript[lambda, k]], {k, 1, n}];
grp = Graphics[{Red, PointSize[0.02], Point[p0]}];
grd = ListLinePlot[data];
grt = Table[Graphics[Text[k, data[[k]]]], {k, 1, n}];
distmin = Infinity;
jmax = 360;
For[j = 0, j <= jmax, j++, tj = 2 Pi j/jmax;
change = False;
vj = v /. {t -> tj};
minresult = -Infinity;
maxresult = Infinity;
For[k = 1, k <= n, k++,
sol = Solve[Thread[l[p0, lambda, vj] == segs[[k]]], {lambda, Subscript[ lambda, k]}][[1]];
lambda0 = Subscript[lambda, k] /. sol;
If[(0 <= lambda0) && (lambda0 <= 1), result = (lambda /. sol), result = Infinity];
If[result != Infinity,
If[result <= 0, If[result >= minresult, minresult = result; topt = tj; change = True]];
If[result >= 0, If[result <= maxresult, maxresult = result; topt = tj; change = True]]]
];
dist = maxresult - minresult;
If[dist <= distmin, distmin = dist; maxr = maxresult; minr = minresult; tmin = topt]
]
vj = v /. {t -> tmin};
pa = l[p0, minr, vj];
pb = l[p0, maxr, vj];
seg = u pa + (1 - u) pb;
gr2 = ParametricPlot[seg, {u, 0, 1}];
grpa = Graphics[{Red, PointSize[0.02], Point[pa]}];
grpb = Graphics[{Red, PointSize[0.02], Point[pb]}];
Show[grp, grd, grt, grpa, grpb, gr2, Axes -> True]
In the figures, the black dot represents $p_0$ and in dashed red the rupture line.
NOTE
The intersections $L_j\cap s_k$ are performed as
$$ p_0+\lambda_0 v_j = \lambda_k p_k + (1-\lambda_k) p_{k+1} $$
giving
$$ \cases{ \lambda_0 = \frac{x_{k+1}(y_0-y_k)+x_0(y_k-y_{k+1})+x_k(y_{k+1}-y_0)}{(y_{k+1}-y_k)\cos t_j+(x_k-x_{k+1})\sin t_j}\\ \lambda_k = \frac{(y_{k+1}-y_0)\cos t_j+(x_0-x_{k+1})\sin t_j}{(y_{k+1}-y_k)\cos t_j+(x_k-x_{k+1})\sin t_j} } $$
Here to have a feasible intersection it is needed $0\le \lambda_k\le 1$
A nice way for doing it would be to find the perpendicular distances from each side. Following that choose $n \choose 2$ distances and find the set which has the least for both elements. If the chosen sides are parallel and comes out to be favourable, then your answer would be the sum of the distances. Else you may follow what's done below:
Am doing for a simplified case:
You can see from here that $$r_1=P_2 \sec(A-B)\ \text{and}\ r_2=P_1 \sec(B)$$ then minimize $r_1+r_2$ differentiating by the changing angle $B$ (since $A$ is fixed). And yipee, you get your solution.
Note: If the sides (whose distance function is the least) don't appear to converge, just make them virtually converge.
For the graph used and manual testing you might use:
- Polygon Version
If you would like to implement it on a program, you would like to follow this (efficient for large number of sides or even loops):
- Consider doing a Fast Fourier Transformation given an arbitrary curve (if you do not have the equations of the curve).
- Choose the point of which you need the shortest chord.
- Make a for-loop and implement it such that it makes a large number of circles with varying radius and a fixed center.
- Post running the loop, add a condition such that the loop breaks when there are two points where both the loop and the circle have a common tangent.
- If the two tangents are parallel, you've already got the required points so compute the distance.
- If not, make an open triangle with the tangent and implement the method adopted for the polygon, it'll be suffice.
For a sample graph you can use:
- Arbitrary loop (maybe a polygon)
Clearly we cannot expect a closed form solution, but a piecewise formula for the distance from each point in the polygon and an algorithm to manage and find the minimum total distance you require is described in the following steps.
a) Translate the polygon so to bring the red point at the origin of the coordinates
b) Express the sides by the vectorial equation $$ {\bf p}_k = t_k {\bf v}_k + \left( {1 - t_k } \right){\bf v}_{k + 1} \quad \left| \matrix{ \;1 \le k \le n - 1 \hfill \cr \;0 \le t_k \le 1 \hfill \cr} \right. $$
c) Convert the sides equations into polar coordinates
That is
$$
\eqalign{
& \left\{ \matrix{
\rho _k \cos \alpha
= t_k v_k \cos \alpha _k + \left( {1 - t_k } \right)v_{k + 1} \cos \alpha _{k + 1} \hfill \cr
\rho _k \sin \alpha
= t_k v_k \sin \alpha _k + \left( {1 - t_k } \right)v_{k + 1} \sin \alpha _{k + 1} \hfill \cr} \right. \cr
& \quad \quad \Downarrow \cr
& \tan \alpha
= {{t_k \left( {v_k \sin \alpha _k - v_{k + 1} \sin \alpha _{k + 1} } \right) + v_{k + 1} \sin \alpha _{k + 1} }
\over {t_k \left( {v_k \cos \alpha _k - v_{k + 1} \cos \alpha _{k + 1} } \right) + v_{k + 1} \cos \alpha _{k + 1} }} \cr
& \quad \quad \Downarrow \cr
& t_k
= v_{k + 1} {{\sin \left( {\alpha _{k + 1} - \alpha } \right)}
\over {\left( {v_k \cos \alpha _k - v_{k + 1} \cos \alpha _{k + 1} } \right)\sin \alpha
- \left( {v_k \sin \alpha _k - v_{k + 1} \sin \alpha _{k + 1} } \right)\cos \alpha }} \cr
& \quad \quad \Downarrow \cr
& \left\{ \matrix{
t_k (\alpha )
= v_{k + 1} {{\sin \left( {\alpha _{k + 1} - \alpha } \right)}
\over {\left( {v_k \cos \alpha _k - v_{k + 1} \cos \alpha _{k + 1} } \right)\sin \alpha
- \left( {v_k \sin \alpha _k - v_{k + 1} \sin \alpha _{k + 1} } \right)\cos \alpha }} \hfill \cr
\rho _k (\alpha )
= {{\left( {v_k \cos \alpha _k - v_{k + 1} \cos \alpha _{k + 1} } \right)t_k (\alpha ) + v_{k + 1} \cos \alpha _{k + 1} }
\over {\cos \alpha }} \hfill \cr} \right. \cr}
$$
where the meaning of the symbols used is evident.
The expression is a bit complicated but well manageable on computer.
d) Partition of angle intervals
Our scope is to find the minimum of $\rho (\alpha ) +\rho (\alpha + \pi ) $ and the relevant $\alpha$.
The function $\rho (\alpha )$ expressed above is piecewise valid for $\alpha _{k} \le \alpha \le \alpha _{k+1}$.
To cope with our goal we will rearrange the angle intervals as follows.
Starting with the following array
$$
\left( {\matrix{ {\left[ {\alpha _1 ,\alpha _2 } \right)} \cr {\rho _1 (\alpha )} \cr } } \right),
\left( {\matrix{ {\left[ {\alpha _2 ,\alpha _3 } \right)} \cr {\rho _2 (\alpha )} \cr } } \right),
\cdots ,
\left( {\matrix{ {\left[ {\alpha _{n - 1} ,\alpha _n } \right)} \cr {\rho _{n - 1} (\alpha )} \cr } } \right),
\left( {\matrix{ {\left[ {\alpha _n ,\alpha _1 } \right)} \cr {\rho _n (\alpha )} \cr } } \right)
$$
we insert $0 = 2 \pi$ and $\pi$
$$
\left( {\matrix{{\left[ {0,\alpha _1 } \right)} \cr {\rho _n (\alpha )} \cr } } \right),
\left( {\matrix{{\left[ {\alpha _1 ,\alpha _2 } \right)} \cr {\rho _1 (\alpha )} \cr } } \right),
\cdots ,
\left( {\matrix{{\left[ {\alpha _m ,\pi } \right)} \cr {\rho _m (\alpha )} \cr } } \right),
\left( {\matrix{{\left[ {\pi ,\alpha _{m + 1} } \right)} \cr {\rho _m (\alpha )} \cr } } \right),
\cdots ,
\left( {\matrix{{\left[ {\alpha _{n - 1} ,\alpha _n } \right)} \cr {\rho _{n - 1} (\alpha )} \cr } } \right),
\left( {\matrix{{\left[ {\alpha _n ,2\pi } \right)} \cr {\rho _n (\alpha )} \cr } } \right)
$$
At this point we consider the two sections of angle intervals
$$
\left\{ \matrix{
\left[ {0,\alpha _1 } \right),\left[ {\alpha _1 ,\alpha _2 } \right), \cdots ,
\left[ {\alpha _m ,\pi } \right) \hfill \cr
\left[ {\pi ,\alpha _{m + 1} } \right), \cdots ,
\left[ {\alpha _{n - 1} ,\alpha _n } \right),\left[ {\alpha _n ,2\pi } \right) \hfill \cr} \right.
$$
deduct $\pi$ from the values in the second
$$
\left\{ \matrix{
\left[ {0,\alpha _1 } \right),\left[ {\alpha _1 ,\alpha _2 } \right), \cdots ,\left[ {\alpha _m ,\pi } \right) \hfill \cr
\left[ {0,\beta _1 = \alpha _{m + 1} - \pi } \right), \cdots ,
\left[ {\beta _{n - m - 1} ,\beta _{n - m} } \right),\left[ {\beta _{n - m} ,\pi } \right) \hfill \cr} \right.
$$
and then "compenetrate" the $\alpha$ and $\beta$ intervals, i.e. arrange $\alpha _k$ and $\beta _k$ sequentially,
into a congruent set of intervals $ \cdots , \left[ {\gamma _{j},\gamma _{j+1} } \right), \cdots$ to
reach and get the following array
$$
\cdots ,\left( {\matrix{
{\left[ {\gamma _j ,\gamma _{j + 1} } \right)} \cr
{r _{j} (\alpha ) = \rho _u (\alpha ) + \rho _v (\alpha + \pi )} \cr
} } \right), \cdots
$$
Finally we can minimize each $r _{j} (\alpha )$ in its interval and choose the minimum.