2
$\begingroup$

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 center...it 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}]; 
n++];]

ListPlot[list, PlotStyle -> {PointSize[.005]}, AspectRatio -> 1]
Improve this question
$\endgroup$
4
  • 3
    $\begingroup$ Maybe you could show the code for generating the point pattern? $\endgroup$
    – bill s
    Jan 23, 2014 at 4:21
  • 1
    $\begingroup$ Did one of the answers in your earlier question address that problem? If so, consider selecting it as the answer. $\endgroup$
    – bobthechemist
    Jan 23, 2014 at 12:33
  • $\begingroup$ 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. $\endgroup$
    – Confused-cius
    Jan 23, 2014 at 13:09
  • 1
    $\begingroup$ 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. $\endgroup$
    – Szabolcs
    Jan 23, 2014 at 20:36

1 Answer 1

Reset to default
8
$\begingroup$

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

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Confused-cius
    Jan 23, 2014 at 21:12
  • $\begingroup$ @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. $\endgroup$
    – Szabolcs
    Jan 23, 2014 at 21:22
  • $\begingroup$ @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.) $\endgroup$
    – Confused-cius
    Jan 24, 2014 at 7:09
  • 1
    $\begingroup$ @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. $\endgroup$
    – user484
    Jan 24, 2014 at 7:36
  • $\begingroup$ @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...) $\endgroup$
    – Confused-cius
    Jan 24, 2014 at 14:20

Your Answer

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

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