# How to: pixel packs in an image to representative point coordinates (x,y)

Let's have an image (gray-scale or even binary) as shown on the following figure at the left hand side, the goal is to generate a list of points i.e., coordinates in the form of (x,y) for each pack of the dark pixels in the image.
What are the proper image processing tools to do this and where are they available?

1) Here you may find some more details about the problem. (Note the variation in the size of packs)

I may suggest having packs detected to compute the convex-hull boundary for each and then find the representative centroid {see this for details}.

2)
Here is the result produced by application of Distance Transform (suggested by "Libor"). Note my annotations on the figure. The method does not work as it was promising!

3)
Erosion eliminates small packs!

from __future__ import division
from scipy import zeros, ndimage as dsp
from pylab import subplot,plot,matshow,show

img = zeros((30,30))
img[10:14,10:14] = 1
img[16:17,16:17] = 1
img[19:23,19] = 1
img[19,19:23] = 1

subplot(221)
matshow(img,0)

subplot(222)
y = dsp.binary_erosion(img,[[1,1],[1,1]])
matshow(y,0)

subplot(223)
y = dsp.binary_erosion(img,[[0,1,0],[1,1,1],[0,1,0]])
matshow(y,0)

subplot(224)
y = dsp.binary_erosion(img,[[1,1,1],[1,1,1],[1,1,1]])
matshow(y,0)

show()


4)
Well here is a Python (i.e., the language of love :) ) implementation of the labeling idea (also proposed by "Jean-Yves" below):

subplot(221)
l,n = dsp.label(img)
sl = dsp.find_objects(l)
for s in sl:
x = (s[1].start+s[1].stop-1)/2
y = (s[0].start+s[0].stop-1)/2
plot(x,y,'wo')


and the result:

Note that although it is done in Python so quick due to Scipy performance, the background procedure in label function should be an exhausting iteration. This may be considered as a trade-off. So for a while I keep being eager to seek more efficient algorithms. And also note that in the given code above I found the center of geometry so simply while for complex or asymmetric shapes this may cause positioning to be biased. That is it is a work in progress ;).

5)
Here is a complex case (a real image) captured from here on which the labelling proposal applied and you see the results. Note that it took only 0.015 s for whole procedure including labelling and finding the objects. Scipy guys, did very well job, I think. Wow! {right-click on image, click on view image for full resolution}

• Do you know: 1) Number of packs. 2) Max/Min x-size of a pack. 3) Max/Min y-size of a pack. – Spacey Apr 28 '12 at 23:38
• @Mohammad Given is only an image. Thus none of the required information you mentioned is available. – Developer Apr 29 '12 at 1:13
• If you can detect the number of packs, then I would use that as the value of 'k' in a generic k-means algorithm. It will converge to the center of each cluster, given a value of k. However, I am not sure how you will be able to determine the value of 'k'. How are you going to determine that? – Spacey Apr 29 '12 at 2:08
• You need sub-pixel accuracy for distance transform / erosion to detect maxima in the small patches (e.g. 2x2). There are maxima, but cannot be detected since you skip them with your sampling. In case of distance transform, you can upsample your image by factor of two or compute the transform for sub-pixel positions (0.0, 0.5, 1.0, 1.5 ...). In case of erosion, you can implement it using PDE-based morphology (iterative). – Libor Apr 29 '12 at 20:03
• @Developer this question of mine on SO might be of interest: stackoverflow.com/questions/4087919/… – Ivo Flipse Apr 30 '12 at 12:26

Just a naive suggestion: do you know about component labeling?

The technique is about finding chunks of "touching" pixels and assigning them a label, e.g. an integer number. You can then interrogate each chink separately, looking for the pixel that share the same label.

In MATLAB, here is the function that does it trivially: bwlabel

• Thanks for the proposal. I had already approached this way and as you mentioned I updated the post with an implemenation in Python. I emerged some issues regarding use of labeling in the Update 4: – Developer Apr 30 '12 at 10:12
• I upvoted this since it is another approach to solve the problem. Thanks for sharing :) – Developer Apr 30 '12 at 10:16
• Could you clarify the update 4? I don't get what you mean by trade off and exhaustion. – Jean-Yves May 1 '12 at 0:19
• Doing labeling is an iterative job requires lots of computation. However Scipy guys has done excellent job (or some tricks!) to do it almost immediately. I put a real complex case for evaluation and it worked only in 0.015 s. Isn't it awesome? Note that the second part i.e., finding the center of detected packs is still a question. Apparently applying convex-hull for such a mssive set of packs is not recommended. – Developer May 1 '12 at 3:46
• There I disagree: connected-components labeling is NOT expansive. It runs in linear time in most modern implementations, and some even managed to have it run multithreaded. – Jean-Yves May 1 '12 at 18:50

You can also run a distance transform on the image, then detect local maxima (searching for pixels with highest/lowest value of all pixels in 3x3 pixel patch - can be larger depending on expected minimum distance between the original blobs).

Note that for detecting features of size 1-3 pixels, you need to double your sampling frequency (either upscale the source image or perform distance transform/erosion with sub-pixel accuracy).

UPDATE:

Both distance transform and erosion approaches assumes that the features you are detecting are convex. Something with U shape, for example, may fire in your detector several times.

A more elaborate method for such segmentation is based on level sets and active contours. It starts with a large closed curve which is iteratively adapted to your features. This method has been used on my university for counting cells and detecting chromosomes in microscope images.

• You can see the result of my implementation of your idea in the section Updates:2. It seems doesn't work as expected :( – Developer Apr 29 '12 at 12:36
• OK, I have extended my answer. – Libor Apr 29 '12 at 20:13
• I upvoted this since it is another approach to solve the problem. Thanks for sharing :) – Developer Apr 30 '12 at 10:16

One option would be to apply repeated morphological erosion to the image until it is fully eroded. At that point, each of the blobs shown above would be reduced down to a single pixel; you could take the locations of those pixels to be the list of points that you're looking for.

• Some points: 1) how to stop iteration? 2) it may be disappear those packs are smaller. 3) it look very computation intensive as iteration required for the whole data set. Note the image could be large. – Developer Apr 28 '12 at 7:18
• I quickly implemented your idea in its current form but it didn't work. It eliminates small packs. All my previous comments apply. – Developer Apr 28 '12 at 7:20
• I'm not sure what you mean by "small packs." Are any shown in the example image? – Jason R Apr 28 '12 at 13:25
• Well, as shown the packs are in different sizes and shapes and as I applied iterative erosion (Python:Scipy:Spatial:Erosion) the smaller (or narrower) ones were disappeared. Note that in actual case the pack in the input image could be in the range of one pixel to very large. – Developer Apr 28 '12 at 15:07
• @Developer Do you know the number of 'packs' ahead of time? For example, in here, you have 6. – Spacey Apr 28 '12 at 16:24

I am afraid that whatever way you pick to do this, it is not going to be straightforward because to assign the targets to clusters you are going to have to go through the image (at least once).

I suppose that getting the points is the easier problem of the two (you are probably already applying some form of thresholding for example).

To recover the clusters that points are grouped into and do it fast, you could create a quad-tree structure that holds "chains of connected pixels" that are somewhere around the area of the quad-tree cell.

In this way, you could iterate through the image and once you come across a pixel that is a target, "push" it in your quad-tree structure.

This "push" operation would start an iterative process which would return the cell (in other words the specific area of the image) where the particular pixel lies in. You could then iterate through all the pixel chains that are assigned to that cell and try to "push" (again) the pixel to the pixel chain. A pixel chain accepts a new pixel if it is at least 1-pixel near to any of its already assigned pixels. If no pixel chain "accepts" the new pixel, then create a new pixel chain into this cell and assign the new pixel to it.

The quad-tree here is a way to limit your search for the nearest pixel-chain and it would be required if your image is big and the targets numerous so that the pixel-chain push operations are done fast. If you know that you are not going to be dealing with a large number of targets then you could even skip the quad-tree and maintain a simple "list of pixel chains". Every time you come across a "target" you could iterate through the lists and try "pushing" the pixel on them. If no list "admits" the pixel then create a new list and assign the pixel in.

Whichever way you choose to do it, at the end of this process you are going to have a set of connected "pixel-chains" (your clusters) which you can then pass to another part of your program that will be dealing with estimating their location. This can be convex-hull, model fitting (for example ellipsoid or other) or simply the mean / median of the x,y coordinates.

I hope this helps.

• I upvoted this since it is another approach to solve the problem. Thanks for sharing :) – Developer Apr 30 '12 at 10:15