Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

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

Following an earlier post here, I have generated a point pattern on the unit square $[0,1]\otimes[0,1]$, which is described by a piecewise constant intensity function constructed so that (in expectation) $2/3$ of the points will fall within a circle of radius $\sqrt{0.1}$ centered on $(0.5,0.5)$. This piecewise intensity function manifests as a visible change in the density of points, forming an "edge" or "boundary" along this circle.

Now, using only the point pattern and no knowledge of the intensity function underlying the process, I am interested in how one could detect this edge, get a formula describing it (assuming we don't know it's a circle with given radius and has to be general) and superimpose the inferred boundary on a ListPlot illustrating the point pattern. My initial ideas are dividing the x and y axes into bins, or getting a CDF in the x and in the y directions and taking partial derivatives in both directions.

Here's the code for the list plot:

n = 0;
list = {};
r := RandomReal[];
coin := RandomChoice[{True, False}];

While[n < 10000, {x, y} = {r, r};
 If[((x - .5)^2 + (y - .5)^2 <= .1 || coin), AppendTo[list, {x, y}]; 

ListPlot[list, PlotStyle -> {PointSize[.005]}, AspectRatio -> 1]
share|improve this question
Maybe you could show the code for generating the point pattern? – bill s Jan 23 '14 at 4:21
Did one of the answers in your earlier question address that problem? If so, consider selecting it as the answer. – bobthechemist Jan 23 '14 at 12:33
Ryan, if you're interested in detecting or identifying features in a general point pattern, you might want to look at kernel principal component analysis as a possible technique. Depending on the your choice of kernel and the type of structure you are interested in detecting, you can often find linear boundaries in the space of principal components that correspond to your non-linear boundaries in the actual pattern. If I have time tonight, I might write an answer demonstrating the technique. – Confused-cius Jan 23 '14 at 13:09
You have already asked several questions, so you must know how this site works. Please do make a little bit of effort to format your posts for readability. – Szabolcs Jan 23 '14 at 20:36

Let's make a set of points:

pts = Select[RandomReal[{-1, 1}, {1000, 2}], Norm[#] < 1 &];
ListPlot[pts, AspectRatio -> Automatic]

enter image description here

Construct the density and plot it, including a contour line at the manually selected value of 0.1, all in one go:

g = SmoothDensityHistogram[pts, Mesh -> {{.1}}, PlotRange -> 1.5 {{-1, 1}, {-1, 1}}]

enter image description here

Use Manipulate to find a good value for the contour, if you like. Use SmoothKernelDistribution manually to set the bandwidth ("smoothening").

Extract and plot that curve:

curve = Cases[Normal[g], Line[pts_] :> pts, Infinity];
ListLinePlot[curve, AspectRatio -> Automatic]

enter image description here

share|improve this answer
Szabolics, this doesn't quite answer the poster's question. He has points on the whole unit square, but with higher density on a circular region within the square. Depending on the relative densities of the points in the two regions, this has the potential to confound your approach somewhat -- Poisson noise can lead to spurious clustering -- and you'll also potentially have to contend with boundary effects from the edge of the square. – Confused-cius Jan 23 '14 at 21:12
@Confused-cius You're right, I didn't remember correctly what he did before. There are points outside of the circle too. But it doesn't invalidate the approach. Take e.g. pts = Join[Select[RandomReal[{-1.5, 1.5}, {3000, 2}], Norm[#] < 1 &], RandomReal[{-1.5, 1.5}, {500, 2}]];, and it works the same way. You are right that if the density difference is smaller then the curve will break up. – Szabolcs Jan 23 '14 at 21:22
@Szabolics, I suspect that your method will depend on both the ratio of the intensities of the "features" versus the "background", and on the total number of points. (Try looking for the feature with only 100 points.) – Confused-cius Jan 24 '14 at 7:09
@Confused-cius: Sure, and try looking for the feature with only 10 points and even a human expert will have trouble... You need a sufficient number of points before the problem is meaningful. – Rahul Jan 24 '14 at 7:36
@RahulNarain: I agree. I just think he's made his job "too easy" by picking a high contrast between foreground and background and choosing $O(10^{3.5})$ points before going hunting. I'm just curious as to how small $N$ can be before the problem is no longer "meaningful" (either with this technique, or with any other). 300 points? Probably, but perhaps not with kernel density estimation (it's too data hungry). 100 points? Maybe. 30 points? Probably not. I'll experiment a bit and post my findings in an actual answer sometime this weekend. (Pesky "work" keeps getting in the way...) – Confused-cius Jan 24 '14 at 14:20

Your Answer


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.