Find shortest path which intersects a point in a polygon

Sorry if I'm explaining this badly, math in English can be a bit troublesome.

I have a polygon, I have a random point inside that polygon. From this point I want a line "drawn" edge-to-edge and to intersect the point, but I want this line to be the shortest possible. See my image below:

The red dot indicates the random point inside polygon. The green dotted line is the shortest path/line (which I'm looking for) The blue vague line is just example of longer lines that does not match the criteria (shortest path of all paths). And, obviously I want the path to intersect the red point.

(My real problem is that I want to find the line AND all the coordinates above that line, but that can be another problem for another day unless someone's feeling really ambitious)

Edit: I want to do this because I want to somewhat (not true physics) simulate the (2D) behavior of cracking a rock and thus want to know what piece of the rock that should split away.

Also, a solution for a convex-polygon would suffice (even though the image shows a non-convex).

• Most of the time it is better to post the real problem. Why don't you o this? – miracle173 Aug 11 '20 at 16:35
• This is an interesting question. I suspect it does not have a nice answer. – Ethan Bolker Aug 11 '20 at 16:35
• I suspect that a simple algorithm can be implemented for convex polygons, but for the general case this problem is quite hard. – K.defaoite Aug 11 '20 at 16:57
• @miracle173 I've updated question to explain why I wanna do this. Please let me know of some other way to do what I wanna do :) – Whyser Aug 11 '20 at 20:19
• @K.defaoite actually I think a convex polygon solution would suffice (the image is misleading, didn't give it much thought about being non-convex when I drew it in paint). Did you have something in mind? – Whyser Aug 11 '20 at 20:23

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$$

• I like code more than I like math expressions (even though I hugely appreciate everyone's effort so far, I haven't been able to try them out yet). Do you think you could either explain the logics parts of the Mathematica code or..a longshot...convert it into Java/C(#)? 😇 – Whyser Aug 12 '20 at 22:55
• @Whyser ASAP, I will try to create a pseudo code to facilitate the conversion to another programming language. – Cesareo Aug 12 '20 at 23:12
• @Whyser Included a geometrical explanation to the presented algorithm. – Cesareo Aug 13 '20 at 8:25
• Thanks alot! This might be a stupid question but..in the Mathematica-code where do you define what lambda (λ) actually is/calculated? If I understand correctly from the added "pseudo code" λ will be 0-1 and that lets say k = 1,2...10, lambda10 woudl be 1, right? – Whyser Aug 13 '20 at 12:17
• The command sol = Solve[Thread[l[p0, lambda, vj] == segs[[k]]], {lambda, Subscript[ lambda, k]}][[1]]; calculates the $\lambda_k$ as well as $\lambda = \lambda_0$ for each $k$ – Cesareo Aug 13 '20 at 12:24

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.

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:

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):

1. Consider doing a Fast Fourier Transformation given an arbitrary curve (if you do not have the equations of the curve).
2. Choose the point of which you need the shortest chord.
3. Make a for-loop and implement it such that it makes a large number of circles with varying radius and a fixed center.
4. 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.
5. If the two tangents are parallel, you've already got the required points so compute the distance.
6. 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:

1. Arbitrary loop (maybe a polygon)
• @Whyser check this for structural implementation. – Anindya Prithvi Aug 13 '20 at 6:20