Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

For a ellipse $E(\theta)$, which owns the parametric equation as follows:

$\begin{cases} x = a \sin\theta + b \cos\theta + c \\ y = d \sin\theta + e \cos\theta + f \\ \end{cases}$

Now, I have a point $P=\{x_p,y_p\}$, to know the relationship between $E(\theta)$ and $P=\{x_p,y_p\}$, I using this method:

$$ \begin{cases} \sin \theta =\frac{(c-x_p)e+b(y_p-f)}{b d-a e}\\ \cos \theta =\frac{(-c+x_p)d+a(f-y_p)}{b d-a e} \end{cases} $$

compute $m=\sin^2 \theta+\cos^2 \theta$

  • if $m=1$, the point $P$ on the ellipse $E(\theta)$
  • if $m>1$, the point $P$ outside the ellipse $E(\theta)$
  • if $m<1$, the point $P$ inside the ellipse $E(\theta)$

However, for the compound curve, for example

  • ellipse segment $E_1(\theta) \qquad \theta \in [0, 2.4798],[5.8629, 2 \pi]$

  • ellipse segment $E_2(\theta) \qquad \theta \in [3.1275, 6.1325]$

  • line segment $L(E_1(2.4798),E_2(3.1275))$


mat1 = {{0., -5., 0}, {-5.2203, 0., 1.7945}}; 
mat2 = {{-0.8583, -4.9384, 0.1765}, {-5.4189, 0.7822, 2.3088}};
θ1 = 2.4798;
θ2 = 3.1275;
Show[
  {ParametricPlot[
     mat1.{Sin[θ], Cos[θ], 1}, {θ, 0, 2.4798}, PlotStyle -> Red], 
   ParametricPlot[
     mat1.{Sin[θ], Cos[θ], 1}, {θ, 5.8629, 2 Pi}, PlotStyle -> Red],
   ParametricPlot[
     mat2.{Sin[θ], Cos[θ], 1}, {θ, 3.1275, 6.1325}],
   Graphics[
     {Line[{mat1.{Sin[θ1], Cos[θ1], 1}, mat2.{Sin[θ2], Cos[θ2], 1}}]}]}, 
  PlotRange -> All]

enter image description here

Now, I have two points $P_1,P_2$, whose coordinate are {{-4.50722, 0.455707}, {5.2748, 2.68502}},respectively.So I would like to know is there a simple method to determine whether $P$ inside the boundary?

enter image description here

share|improve this question
2  
Can't you use RegionMember? – RunnyKine Mar 26 at 6:04
4  
if performance is important, the fastest and most precise approach I suspect will be to use FindRoot to find where/if each segment intersects a line from P to (Infinity,Py). Then just count the roots, odd->inside – george2079 Mar 26 at 11:55
    
@george2079 Could you give me a explanation about (Infinity,Py)? Thanks:) – Shutao TANG Mar 26 at 14:52
up vote 19 down vote accepted

One way you can do this is to first discretize the graphics, then turn it into a region with an interior using DelaunayMesh, and finally using RegionMember:

mat1 = {{0., -5., 0}, {-5.2203, 0., 1.7945}}; 
mat2 = {{-0.8583, -4.9384, 0.1765}, {-5.4189, 0.7822, 2.3088}};
θ1 = 2.4798;
θ2 = 3.1275;
gr=
 Show[
  {ParametricPlot[
     mat1.{Sin[θ], Cos[θ], 1}, {θ, 0, 2.4798}, PlotStyle -> Red], 
   ParametricPlot[
     mat1.{Sin[θ], Cos[θ], 1}, {θ, 5.8629, 2 Pi}, PlotStyle -> Red],
   ParametricPlot[
     mat2.{Sin[θ], Cos[θ], 1}, {θ, 3.1275, 6.1325}],
   Graphics[
     {Line[{mat1.{Sin[θ1], Cos[θ1], 1}, mat2.{Sin[θ2], Cos[θ2], 1}}]}]}, 
  PlotRange -> All];

Now we turn the graphics into a region:

dg = DelaunayMesh@MeshCoordinates@DiscretizeGraphics@gr;

And we can now determine that {0,3} is inside this region:

RegionMember[dg, {0., 3.}]
True

Note that you can get better performance by creating a RegionMemberFunction and applying it directly to a list of points (it's Listable).

rf = RegionMember[dg]; (* create RegionMemberFunction *)
rf[{{0., 3.}, {0., 2.}, {5., 0.}, {0., 8.}, {1., 1.}}] (* apply the function *)

{True, True, False, False, True}

share|improve this answer
    
Apart from DelaunayMesh[] and RegionMember[](I think DelaunayMesh[] is time-consuming ). I would like to know is there other simple method to do this problem? – Shutao TANG Mar 26 at 6:35
    
@ShutaoTANG, you can also use ConvexHullMesh in place of DelaunayMesh. – RunnyKine Mar 26 at 6:50
3  
@ShutaoTANG Another approach is to use SignedRegionDistance in place of RegionMember and if the result is negative it means the point is inside the region and outside otherwise. – RunnyKine Mar 26 at 6:51
    
OK, thanks a lot:) – Shutao TANG Mar 26 at 6:58

I am going to use Graphics`PolygonUtils`PointWindingNumber from R.M's answer to How to check if a 2D point is in a polygon?. In M10 you can use RegionMember.

The idea is to convert your region into a polygon and then test if the point is inside or not.

mat1 = {{0., -5., 0}, {-5.2203, 0., 1.7945}}; 
mat2 = {{-0.8583, -4.9384, 0.1765}, {-5.4189, 0.7822, 2.3088}};

q1 = 2.4798; q2 = 3.1275;
dq = 0.1; (*use smaller value for better precession*)
seg1 = Table[mat1.{Sin[q], Cos[q], 1}, {q, 0, 2.4798, dq}];
seg2 = Table[mat1.{Sin[q], Cos[q], 1}, {q, 5.8629, 2 Pi, dq}];
seg3 = Table[mat2.{Sin[q], Cos[q], 1}, {q, 3.1275, 6.1325, dq}];
seg4 = {mat1.{Sin[q1], Cos[q1], 1}, mat2.{Sin[q2], Cos[q2], 1}};
path = Join[seg1, seg2, seg3, seg4];

(*Make a continuous path from all the points*)
spath = FindCurvePath[path]//First; 
boundary = path[[#]] & /@ spath;

inPolyQ[poly_, pt_] := Graphics`PolygonUtils`PointWindingNumber[poly, pt] =!= 0


Manipulate[Graphics[Line[boundary], 
 PlotLabel -> Text[Style[StringForm["Point`` is ``", p, 
  If[inPolyQ[boundary, p], "Inside ", "Outside"]], Bold, Italic]]],
 {{p, {0, 0}}, Locator}]

enter image description here

share|improve this answer
    
Thanks @MarcoB for the edit. You left a � after the final bracket. So I am taking the advantage ;) – Sumit Mar 26 at 12:43
    
sorry about that! That's what I get for doing edits on my silly tablet! :-) I've got to look into this PointWindingNumber thing. It sounds very cool! Thanks for digging it out. – MarcoB Mar 26 at 12:54

Perhaps a bit of an overkill. Imagine the points on the boundary and the points you input as electric charges:

mat1 = {{0., -5., 0}, {-5.2203, 0., 1.7945}}; 
mat2 = {{-0.8583, -4.9384, 0.1765}, {-5.4189, 0.7822, 2.3088}};
θ1 = 2.4798;
θ2 = 3.1275;
g=Show[
  {ParametricPlot[
     mat1.{Sin[θ], Cos[θ], 1}, {θ, 0, 2.4798}, PlotStyle -> Red], 
   ParametricPlot[
     mat1.{Sin[θ], Cos[θ], 1}, {θ, 5.8629, 2 Pi}, PlotStyle -> Red],
   ParametricPlot[
     mat2.{Sin[θ], Cos[θ], 1}, {θ, 3.1275, 6.1325}, PlotStyle -> Red],
   Graphics[
     {Thick, Red, Line[{mat1.{Sin[θ1], Cos[θ1], 1}, mat2.{Sin[θ2], Cos[θ2], 1}}]}]}, 
  PlotRange -> All];

step[form_] := 
 Function[# - 
    10^-4 Sum[
      Normalize[ext[[i]] - #]/Norm[ext[[i]] - #]^2, {i, 
       Length@ext}]] /@ form
relax[object_, t_] := First@{
   ext = RandomPoint[DiscretizeGraphics@g, 10^3];
   fin = NestList[step, object, t]; 
   Print@ListAnimate[
     Function@Show[g, #, Axes -> False] /@ 
      Function@ListPlot[#, PlotStyle -> Gray] /@ fin];
   Table[{Part[First@fin, i], 
      If[EuclideanDistance[Part[Last@fin, i], Mean@ext] - 
         EuclideanDistance[Part[First@fin, i], Mean@ext] < 0, 
       Text["Inside"], Text["Outside"]]}, {i, Length@First@fin}
     ] // MatrixForm
   }

Usage

Note that increasing t yields more accurate results.

relax[{{0, 0}, {2, 0}, {5, 6}}, 100]

enter image description here

share|improve this answer

the direct approach:

mat1 = {{0., -5., 0}, {-5.2203, 0., 1.7945}};
mat2 = {{-0.8583, -4.9384, 0.1765}, {-5.4189, 0.7822, 2.3088}};
\[Theta]1 = 2.4798;
\[Theta]2 = 3.1275;
f[1, t_] := mat1.{Sin[t], Cos[t], 1};
f[2, t_] := mat1.{Sin[t], Cos[t], 1};
f[3, t_] := mat2.{Sin[t], Cos[t], 1};
f[4, t_] := 
            t mat1.{Sin[\[Theta]1], Cos[\[Theta]1], 1} +  
      (1 - t) mat2.{Sin[\[Theta]2], Cos[\[Theta]2], 1};
f[1, "range"] = Sequence @@ {0, \[Theta]1};
f[2, "range"] = Sequence @@ {5.8629, 2 Pi};
f[3, "range"] = Sequence @@ {\[Theta]2, 6.1325};
f[4, "range"] = Sequence @@ {0, 1};

p = RandomReal[{-5, 5}, {2}]
Show[{Graphics@{PointSize[.025], Point[p], Arrow[{p, p + {10, 0}}]},
  Table[ParametricPlot[f[k, t], Evaluate@{t, f[k, "range"]}, 
    PlotStyle -> Red], {k, 4}]}, PlotRange -> All]

enter image description here

count the crossings, if odd we are inside.

OddQ@Total@Table[ Length[t /. Quiet@NSolve[f[k, t][[2]] == #[[2]] && 
          Less @@ Riffle[{f[k, "range"]}, t] && 
          f[k, t][[1]] > #[[1]] , t]] , {k, 4}] &@p

True

process a bunch of points:

p = RandomReal[{-8, 8}, {20, 2}];
pin = Select[p, 
   OddQ@Total@Table[ Length[t /. Quiet@NSolve[f[k, t][[2]] == #[[2]] && 
            Less @@ Riffle[{f[k, "range"]}, t] && 
            f[k, t][[1]] > #[[1]] , t]] , {k, 4}] &];
Show[{Graphics@{PointSize[.025], Point[p], Red, PointSize[.02], 
    Point[pin]},
  Table[ParametricPlot[f[k, t], Evaluate@{t, f[k, "range"]}, 
    PlotStyle -> Red], {k, 4}]}, PlotRange -> All]

enter image description here

If you have functions that NSolve can't handle you can use other methods such as FindRoot, but you need a scheme to find all the roots, and then you might loose the performance advantage over other methods.

Also to make this robust you need to deal with cases where your projected line exactly hits a junction between piecewise curves. (do a Union on the roots for example )

share|improve this answer
    gra = Show[{ParametricPlot[
         mat1.{Sin[\[Theta]], Cos[\[Theta]], 1}, {\[Theta], 0, 2.4798}, 
         PlotStyle -> Red], 
        ParametricPlot[
         mat1.{Sin[\[Theta]], Cos[\[Theta]], 1}, {\[Theta], 5.8629, 2 Pi},
          PlotStyle -> Red], 
        ParametricPlot[
         mat2.{Sin[\[Theta]], Cos[\[Theta]], 1}, {\[Theta], 3.1275, 
          6.1325}], 
        Graphics[{Line[{mat1.{Sin[\[Theta]1], Cos[\[Theta]1], 1}, 
            mat2.{Sin[\[Theta]2], Cos[\[Theta]2], 1}}]}]}, 
       PlotRange -> All];
    RegionMember[
 ConvexHullMesh@MeshCoordinates@DiscretizeGraphics@gra, {0, 3}]

True

share|improve this answer
    
Obviously, {0,3} is inside the boundary. Why your answer is False?:) – Shutao TANG Mar 26 at 6:14
    
@ShutaoTANG Well just misunderstand your English.I thought you jedge whether or not in the boundry but inside. – yode Mar 26 at 6:48

Not sure how you would do that in Mathematica, but the simplest way to check if a point is inside a closed shape, is to draw a line through the point.
Then count intersections with the border of the shape from the point going outwards.
If the number is odd, it's inside.

share|improve this answer
1  
This is more of a comment than answer. If you would like to expand on it and write an actual answer, feel free! Otherwise, please post this is a comment. – march Mar 26 at 17:34
    
Also, quoting from @george2079's comment from 5 hours ago: "Then just count the roots, odd->inside." You might want to turn this into code if you can. – thedude Mar 26 at 17:36

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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