# How to find the intersection of a Line with the boundary of a Region?

Given a 2-D Region $\Omega$, i.e.

RegionQ[ω] && RegionDimension[ω] == 2


is True, and a pair of points p1 and p2 such that

RegionMember[ω, p1] && !RegionMember[ω, p2]


is True, I need to find the intersection of the Line[{p1, p2}] with the RegionBoundary[ω]. Ideally the point-intersection nearest to p1 if possible and if not so much expensive.

Ideally, the method should give:

• an exact answer if the region is "exact" and if analytically possible;
• an answer with a Precision and/or Accuracy similar to the input for "numerical" regions;
• an answer p such that RegionMember[RegionBoundary[ω, p] gives True.

None of the above is strictly required; I can accept to obtain an approximate answer even for an exact region; in that case a way to eventually specify the Precision/Accuracy sought would be helpful.

The first basic strategy I tried is based on the following function:

RegionNearest[RegionIntersection[Line[{p1, p2}], RegionBoundary[ω]], p1]


This probably satisfies all the requirements above but is very very slow even for something like 100 pairs of points.

Following another post here I also tried something like:

r = RegionIntersection[Line[{p1, p2}];
Block[{x, y}, {x, y} /. (RegionMember[r, {x, y}] // Reduce // ToRules)]


which is not as robust, but is faster, even if not enough.

The third way I tried is a kind of bisection algorithm:

q = RegionMember[ω];
Module[{n1=N@p1, n2=N@p2, nm},
FixedPoint[(If[q[nm=Mean[{n1,n2}]],n1=nm,n2=nm];nm)&,n1]
]


This approximately works, is fast enough. Unfortunately I don't know how to control the precision of the output and the answer I get often doesn't satisfy the third requirement.

Is there a way to get better results?

I considered also working on a someway discretized version of the region (I think I can control the accuracy of the mesh and so the accuracy of the answer), but I don't know how to start and if it would be useful.

• Can the region be any kind of symbolic region? For special region types, such as polygons / circles / etc. it is easier to devise fast algorithms. – JEM_Mosig Jul 1 '17 at 4:19

## 3 Answers

For many kinds of regions you can use Solve. Here, I work with random polygons as an example.

Generate a random polygon

Ω = Polygon[SortBy[RandomReal[{-1, 1}, {10, 2}], ArcTan[Sequence @@ #] &]]

(* Polygon[{{-0.874595, -0.590699}, {-0.295508, -0.924563}, {0.29796, \
-0.370009}, {0.161826, 0.12431}, {0.871938, 0.913455}, {0.336113,
0.827302}, {-0.202137, 0.952052}, {-0.569595,
0.702745}, {-0.320668, 0.179451}, {-0.571918, 0.150367}}] *)


Generate a random point inside the polygon

p1 = RandomPoint[Ω]
(* {0.0387464, 0.459741} *)


Generate a random point outside the polygon

p2 = RandomPoint[RegionDifference[Rectangle[{-2, -2}, {2, 2}], Ω]]
(* {0.61232, -1.04145} *)


Find the point(s) of intersection

c = Quiet[
Values@Solve[
{x, y} ∈ Line[{p1, p2}] && {x, y} ∈ Line@@Ω,
{x, y}
],
{Solve::ratnz}
]
(* {{0.165393, 0.128274}} *)


Here I convert the Polygon to a Line, because Solve appears to be unable to handle Polygon in Version 11.0.1. I also suppress an unimportant warning message about inexact input.

The generated point obeys your condition

RegionMember[RegionBoundary[Ω], c[]]
(* True *)


Here is a graphical representation:

Graphics[{
Ω,
Dashed, Gray, Arrow[{p1, p2}],
PointSize -> Large,
Red, Point[p1],
Blue, Point[p2],
Orange, Point[c]
}] Using the WorkingPrecision option of Solve and SetPrecision on the input parameters, you can get answers with arbitrary precision.

If multiple points of intersection are found, you may use Nearest[c, p1] to find the one closest to p1.

• Line@@Ω does not close the polygon correctly. I used this: Line[Append[#, First@#] & /@ Ω[]]. – Danvil Jan 12 '18 at 8:05

I've made some other attempts, but I think this is a good solution.

ω = Rectangle[{0, 0}];
RegionQ[ω] && RegionDimension[ω] == 2


True

p1 = {0, 0}; p2 = {1.1, 1.1};
RegionMember[ω, p1] && ! RegionMember[ω, p2]


True

Graphics[{ω, Red, PointSize[0.05], Point[p1], Point[p2]}] q = RegionMember[ω];
Module[{n1 = N@p1, n2 = N@p2, nm},
FixedPoint[(If[q[nm = Mean[{n1, n2}]], n1 = nm, n2 = nm]; nm) &, n1]]

• But how is it different from the solution in the question? – C. E. Jun 30 '17 at 23:08

Just an example:

p = RandomReal[{-1, 1}, {50, 2}];
dm = DelaunayMesh[pts];
a = BoundaryMesh[dm];
pts = MeshCoordinates[a];
mr = MeshRegion[pts, MeshCells[a, 1]];
mc = MeshCells[mr, 1] /. Line[{x_, y_}] :> Line[pts[[{x, y}]]];
il = InfiniteLine[{{0, 0}, {1, 1}}];
int = (RegionIntersection[#, il] & /@ mc) /.
EmptyRegion[_] :> Sequence[]
Show[a, Graphics[{il, Red, PointSize[0.02], int}]] 