I have a list of points, which represent positions of peaks in an image. The nature of the image itself shouldn't matter, as I already have the peak positions.

As a simple example, consider the following code, which produces the image below.

points = {{27, 38}, {31, 50}, {33, 56}, {34, 38}, {39, 51}, {39, 
    63}, {40, 31}, {42, 45}, {42, 55}, {47, 27}, {47, 50}, {48, 
    35}, {48, 65}, {49, 43}, {52, 57}, {54, 50}};

img = ReplacePart[ConstantArray[0., {96, 96}], points -> 1.];
img = ImageAdjust@GaussianFilter[ImageRotate@Image@img, 3.];

Show[Image[img, ImageSize -> Large], 
 Graphics[{Red, Circle[#, 4] & /@ points}]]

enter image description here

I'm interested in extracting the region of the image around each point, which will be a ragged array. (I don't want square patches).

I'd like the result to be a list of the $(x,y)$ coordinates of the pixels contained within each region.

For example, I've added some red circles around each point.

  • For the solitary circles, the part I'd like to extract is simply the coordinates of all the pixels within the circle.
  • The tricky part is for the points where the circles overlap, particularly in the middle of the image. In this case I want to extract the pixels within the union of the overlapping circles. For example, these two circles overlap, so I want the coordinates of the pixels in the region below:

enter image description here

MorphologicalComponents is one way to do it, as, of course, instead of drawing the points as Gaussians, I could draw them as circles. This can be extended to retrieve the pixel coordinates within each region.

dImg = ColorNegate@
  Binarize@Rasterize[Graphics[{Disk[#, 3.5] & /@ points}]]
GraphicsRow[{dImg, MorphologicalComponents[dImg] // Colorize}, 
  ImageSize -> Large]

enter image description here

Ideally, I'm looking for a non-image-based approach that uses two inputs.

  1. The list of points (e.g. as given at the top of the question)
  2. The radius of the circles, in pixels

Although if the graphical approach using MorphologicalComponents is the most efficient, then so be it.

  • 1
    $\begingroup$ You should skip the binarize and use a Dilation on the original point image instead. I would have written up an answer, but your reasons for not using an image/integer bases approach are a bit vague. Do you have important reasons why a pixel based approach would not be sufficient? $\endgroup$ – halirutan Dec 17 '14 at 13:44
  • $\begingroup$ Not necessarily an important reason. My thinking was that once I had a numerical list of points, it would be more efficient to work with that rather than converting to an image again. $\endgroup$ – dr.blochwave Dec 17 '14 at 13:48

Here is a graph based solution. However the result (the pixels) will be dependent on resolution, and this solution currently uses the lowest resolution possible.

withinRange[range_][pt1_, pt2_] := Norm[pt1 - pt2] <= 2 range
adjacencyMatrix = Boole@Outer[withinRange[3.5], points, points, 1]  -IdentityMatrix@Length@points;

Some work may be saved by only traversing the upper triangular of the matrix and then reflecting that portion, since the adjacency matrix is symmetric. Get the circles which are connected:

components = ConnectedComponents@AdjacencyGraph@adjacencyMatrix;

Get the pixels:

getPixels[range_][center_] := center + # - {range + 1, range + 1} & /@ Position[DiskMatrix[range], 1]
pixels = Union @@@ Map[Composition[getPixels[3.5], points[[#]] &], components, {2}];

And finally a visualization:

  Join @@ 
   MapIndexed[With[{i = First@#2}, Floor@# -> i & /@ #] &, pixels]


  • 1
    $\begingroup$ Do you feel the lack of built in function to cluster data? It is kind of desired functionality, right? It is strange that user have to create adjacencyMatrix and stuff by himself. Or have I missed something? (+1ed ofc) $\endgroup$ – Kuba Dec 18 '14 at 10:24
  • 1
    $\begingroup$ @Kuba Yes, and I would also like to have a KDTree or a Quadtree to make collision checking more scalable. Other scientific computing packages such as numpy for Python include a KDTree. $\endgroup$ – C. E. Dec 18 '14 at 11:58

I would work purely on the list of points. You have the radius r, so you can easily identify connected regions with

points = {{27, 38}, {31, 50}, {33, 56}, {34, 38}, {39, 51}, {39, 
    63}, {40, 31}, {42, 45}, {42, 55}, {47, 27}, {47, 50}, {48, 
    35}, {48, 65}, {49, 43}, {52, 57}, {54, 50}};
data = ReplacePart[ConstantArray[0., {96, 145}], points -> 1];
radius = 4;

comp = MorphologicalComponents[Image[Dilation[data, DiskMatrix[radius-1]]]];

(see that I used radius-1 because we already have a point!)

Mathematica graphics

Since the labels in the component image always run from 0 (background) to the highest component number, extracting the positions of the each region is as simple as

pos = Position[comp, #] & /@ Range[Max@Flatten[comp]];

Care should be taken when combining graphics due to the different coordinate systems

With[{dim = Reverse@Dimensions[comp]},
   Inset[Colorize[Reverse@comp], {1, 1}, {.5, .5}, dim],
   White, Circle[Reverse@#, 3.5] & /@ points}, 
  PlotRange -> Transpose[{{1, 1}, dim}]]

Mathematica graphics

  • $\begingroup$ I'm not sure about the fact that DiskMatrix forces you to use integer radii. Anyway, not too important $\endgroup$ – Dr. belisarius Dec 17 '14 at 14:23
  • $\begingroup$ @belisarius No, it doesn't force you but no matter what you supply as argument, the output is a bit-matrix. At least this is how I understood it and how it behaves. $\endgroup$ – halirutan Dec 17 '14 at 15:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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