-1
$\begingroup$

This question already has an answer here:

Let's create a rather dense grid of initial conditions

xlim = 4;
ds = 0.05;
data = Flatten[Table[{i, j}, {i, -xlim, xlim, ds}, {j, -xlim, xlim, ds}], 1];
nic = Length[data]
L0 = ListPlot[data, PlotStyle -> {Blue, PointSize[0.001]}];

Then we define a closed two-dimensional curve

C0 = ContourPlot[x^2/4 + y^2/9 == 1, {x, -5, 5}, {y, -5, 5}, 
     PlotPoints -> 100, ContourStyle -> {Black, Thick}];
d0 = C0[[1, 1]];
S0 = ListPlot[d0, PlotStyle -> {Black, Thick}, Joined -> True];

Here I would like to clarify that in my real case scenario we do not know the analytical expression of the closed curve. We only have a list of points such as d0.

P0 = Show[{L0, S0}, AspectRatio -> 1]

enter image description here

Now I would like to create a new list data2 which will contain all the points of the initial data that are inside the closed curve.

I am using version 9.0 in 32bit Win XP.

Any suggestions?

NOTE: This question point-in a polygon test looks very similar. However in my case performance over a large number of points is important.

$\endgroup$

marked as duplicate by J. M. is away plotting Mar 15 '18 at 2:34

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • 1
    $\begingroup$ You could use a point-in-polygon test. $\endgroup$ – J. M. is away Jun 19 '16 at 6:39
  • 1
    $\begingroup$ @J.M. That looks like a duplicate to me; do you disagree? $\endgroup$ – Mr.Wizard Jun 19 '16 at 6:42
  • 1
    $\begingroup$ @Mr.Wizard I saw that. But I think as a general question, easy solutions in V10 would help future visitors. $\endgroup$ – happy fish Jun 19 '16 at 6:49
  • 1
    $\begingroup$ @Vaggelis_Z It would appear that of the existing answers Daniel's is likely most optimized for speed: (9489) -- have you tried that yet? $\endgroup$ – Mr.Wizard Jun 19 '16 at 6:52
  • 1
    $\begingroup$ The question you linked has 11 answers - were all of them too slow? At least one answer includes a timing measurement that (I think) means it can process 1 million points per second. How many points do you have? $\endgroup$ – Niki Estner Jun 19 '16 at 10:12
4
$\begingroup$

First you produce your data

xlim = 4;
ds = 0.25;
data = Flatten[Table[{i, j}, {i, -xlim, xlim, ds}, {j, -xlim, xlim, ds}], 1];
nic = Length[data]

Then you make your contour from arbitrary points and get the polygon out of it

pts = RandomReal[{2, 3}] {Cos[#], Sin[#]} & /@ Range[0, 2 Pi, Pi/20];
contour = ListLinePlot[pts, AspectRatio -> 1,  
                       PlotRange -> {{-xlim, xlim}, {-xlim, xlim}}]
poly = Cases[Normal@contour, Line[x_] :> x, Infinity];

enter image description here

Actually here poly is pts. This method comes handy if you want to get the points from an analytical expression. In worst case if your data is not ordered, use

poly = Sort[pts, ArcTan @@ #1 > ArcTan @@ #2 &];

And then use RegionMember to find the points.

data2 = Select[data, RegionMember[Polygon[poly], #] &];

ListPlot[{data, data2}, 
PlotStyle -> {PointSize[Medium], PointSize[Large]}, AspectRatio -> 1]

enter image description here

For Mathematica9

For MMA9 you have to use Graphics`Mesh`PointWindingNumber which I first found here. In that case,

inPolyQ[poly_, pt_] := Graphics`Mesh`PointWindingNumber[poly, pt] =!= 0
data2 = Select[data, inPolyQ[poly, #] &];

ListPlot[{data, data2}, 
 PlotStyle -> {PointSize[Medium], PointSize[Large]}, AspectRatio -> 1]

enter image description here

For details you can check my answer to Is there a simple strategy to determine whether a point is inside a boundary?

$\endgroup$
  • $\begingroup$ This works only when the curve has a known analytical expression. In my real case, the curve is given only as a set of points. $\endgroup$ – Vaggelis_Z Jun 19 '16 at 9:34
  • $\begingroup$ Just make a close contour anyway. Hold on I am going to give an example with arbitrary point. $\endgroup$ – Sumit Jun 19 '16 at 9:36
  • $\begingroup$ RegionMember is introduced in V10, may need to replace it to run on OP's computer $\endgroup$ – happy fish Jun 19 '16 at 9:50
  • $\begingroup$ @happyfish Sorry, I noticed that. Unfortunately now I don't have M9 in my hand, so I just describe the method for M9 without an example. $\endgroup$ – Sumit Jun 19 '16 at 10:01
  • $\begingroup$ In v9 when I evaluate data2 = Select[data, inPolyQ[poly, #] &]; I get exactly the same points as in data. $\endgroup$ – Vaggelis_Z Jun 19 '16 at 10:14
2
$\begingroup$

In this case, at least, what you seek can be achieved by evaluating the function defining the contour at each of your test points

data1 = Select[data, Function[{xy}, xy[[1]]^2/4 + xy[[2]]^2/9 < 1]];
L1 = ListPlot[data1, PlotStyle -> {Blue, PointSize[0.001]}];
P1 = Show[{L1, S0}, AspectRatio -> Automatic]

enter image description here

$\endgroup$
  • $\begingroup$ No, this is not what I want. This is just a simple example. In my real case the analytical expression of the closed curve is unknown. I have only a set of $(x_0,y_0)$ points which define a closed curve. $\endgroup$ – Vaggelis_Z Jun 19 '16 at 10:39

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