-1
$\begingroup$

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.

Improve this question
$\endgroup$
17
  • 1
    $\begingroup$ You could use a point-in-polygon test. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Jun 19, 2016 at 6:39
  • 1
    $\begingroup$ @J.M. That looks like a duplicate to me; do you disagree? $\endgroup$
    – Mr.Wizard
    Commented Jun 19, 2016 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$
    – vapor
    Commented Jun 19, 2016 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
    Commented Jun 19, 2016 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
    Commented Jun 19, 2016 at 10:12

2 Answers 2

Reset to default
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?

Improve this answer
$\endgroup$
9
  • $\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
    Commented Jun 19, 2016 at 9:34
  • $\begingroup$ Just make a close contour anyway. Hold on I am going to give an example with arbitrary point. $\endgroup$
    – Sumit
    Commented Jun 19, 2016 at 9:36
  • $\begingroup$ RegionMember is introduced in V10, may need to replace it to run on OP's computer $\endgroup$
    – vapor
    Commented Jun 19, 2016 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
    Commented Jun 19, 2016 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
    Commented Jun 19, 2016 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

Improve this answer
$\endgroup$
1
  • $\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
    Commented Jun 19, 2016 at 10:39

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