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

Question

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.

Answer 1

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[[1]]]
(* 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.

Answer 2

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

Answer 3

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}]]