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

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.

  • $\begingroup$ 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. $\endgroup$
    – JEM_Mosig
    Commented Jul 1, 2017 at 4:19

3 Answers 3


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[
    {x, y} ∈ Line[{p1, p2}] && {x, y} ∈ Line@@Ω,
    {x, y}
(* {{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:

  Dashed, Gray, Arrow[{p1, p2}],
  PointSize -> Large,
  Red, Point[p1],
  Blue, Point[p2],
  Orange, Point[c]

graphical representation

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.

  • 2
    $\begingroup$ Line@@Ω does not close the polygon correctly. I used this: Line[Append[#, First@#] & /@ Ω[[1]]]. $\endgroup$
    – Danvil
    Commented Jan 12, 2018 at 8:05

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

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


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


Graphics[{ω, Red, PointSize[0.05], Point[p1], Point[p2]}]

enter image description here

q = RegionMember[ω];
Module[{n1 = N@p1, n2 = N@p2, nm}, 
 FixedPoint[(If[q[nm = Mean[{n1, n2}]], n1 = nm, n2 = nm]; nm) &, n1]]
  • $\begingroup$ But how is it different from the solution in the question? $\endgroup$
    – C. E.
    Commented Jun 30, 2017 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}]]

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.