Detecting the number of pixels between connected objects

Question

I have a binary image and I need to find the number of pixels between the connected objects I'm using MATLAB and here's my image

Thx

Answer 1

An interesting problem, but not too difficult to solve. It will just require a lot of code, but each line of code is very simple in its understanding. This is the basic algorithm that I would perform:

1. Find all circles in your image using the Circular Hough Transform - imfindcircles
2. Remove the circles from your image
3. Use regionprops to compute the areas of the remaining square objects and determine the minimum bounding box required to encapsulate the squares where they had circles overlapping on them.
4. These minimum bounding box co-ordinates would be the top left and bottom right co-ordinates of a square if it were to be completely filled. Because there are circles that overlap onto squares, we can find the area of overlap by simply taking the full area of each bounding box and subtract with the detected area. The amount that is left would be areas that were overlapping from the circles.
5. To figure out the area that is being overlapped by the two circles, all we need to do is draw both circles on the image and do logical AND operation. This will find the areas where both circles intersect. Once we have those, we simply count up the pixels that are remaining and that would be the overlapping area for the two circles.

I'm going to assume you have the Image Processing Toolbox, or my method will not work. You also need at least MATLAB R2012a for the Circular Hough Transform to be included with your MATLAB distribution.

---

Step #1

We need to first invert your image so that white corresponds to object pixels and black corresponds to background. Otherwise, regionprops and imfindcircles won't work. I'm going to read in your image directly from StackOverflow:

im = imread('https://i.sstatic.net/WQtgP.png'); %// Already binary!
im = ~im;

Next, let's use the Circular Hough Transform. We would use it like so:

[centres, radii] = imfindcircles(im, [25 60]);

im would be the binary image we just converted and the array [25 60] denotes the minimum and maximum radii of the circles we want to detect. In this case, the minimum radius is 25 while the maximum radius is 60. centres and radii denote the column and row co-ordinates of the centre of each circle and the radius of each circle respectively. The first column of centres is the column co-ordinate and the second column is the row co-ordinate. After running this code, we detect 4 circles (as expected). centres and radii look like this:

centres =

  267.5005   67.5005
  233.5152  200.4808
   83.2850   83.2850
  117.6691  118.0193

radii =

   33.3178
   33.1137
   32.9332
   32.8488

This means that one circle was detected at column 267.50, row 67.50 with a radius of 33.3178 which is read from the first row of each result. You can follow along with the rest of the values in each variable. What's great about the Circular Hough Transform is that it can detect partial circles, which is very nice due to the overlapping of circles and squares. If you want to see the detected results, you can use viscircles like so:

imshow(im);
viscircles(centres, radii, 'DrawBackgroundCircle', false);

I set DrawBackgroundCircle to false because when we draw the detected circles, by default it draws an outline of white over the circle. I don't want you to confuse this with any object pixels, and so I set this flag to false.

This is the figure we get:

Cool!

---

Step #2

The easiest way would be to loop over all of your circles, create a circular grid of co-ordinates and set these pixels in your image to false to remove these from your image. You can do this by using meshgrid to create a 2D grid of row and column locations and index into your image to set them to false. Something like this:

[X,Y] = meshgrid(1:size(im,2), 1:size(im,1));

im_no_circles = im;
for idx = 1 : numel(radii)
    r = radii(idx);
    c = centres(idx,:);

    mask = (X - c(1)).^2 + (Y - c(2)).^2 <= r^2;               
    im_no_circles(mask) = false;
end

This code will take each circle and set the corresponding locations in your image to false. However, what the image looks like now is:

There are some edge artifacts from the circles due to quantization noise. We can remove these pixels by using bwareaopen. Any areas whose areas are below a certain amount, remove those from the image. We can choose something like 50, because those squares are certainly beyond areas of 50 pixels while the spurious pixels aren't. After we do this, let's do a morphological opening with imopen to get rid of any noisy pixels that are connected to each of the squares so that we can truly get the square shapes without any of the circles.

So:

im_no_circles_open = bwareaopen(im_no_circles, 50);
im_open = imopen(im_no_circles_open, strel('square', 5));

This is finally what we get:

Step #3

We now call regionprops like so:

s = regionprops(im_open, 'Area', 'BoundingBox');

s will contain a structure describing each of our square objects. If we examine each attribute, this is what we get:

areas = [s.Area].'
bb = reshape([s.BoundingBox], [], numel(areas)).'

areas =

        4178
        4138
        4489

bb =

  134.5000   50.5000   67.0000   67.0000
  150.5000  200.5000   67.0000   67.0000
  334.5000   59.5000   67.0000   67.0000

Each square area is found in the array areas while bb is an array where each row contains information about each square. Specifically, the first and second elements of each row are the column co-ordinates and row co-ordinates of the top-left corner of each bounding box that is used to fully encapsulate the object. The third and fourth elements are the width and height of each bounding box. As such, this is telling us that we need a box of 67 x 67 pixels for each bounding box to fully encapsulate the object. This is also equal to the total area if the entire bounding box was full.

Step #4

How MATLAB detects the boxes is that it goes from top to bottom, left to right. As such, there are only two boxes that overlapped with any circles, and so those are the first two boxes in the detected result we need to look at. Therefore, for the square and the circle on the top left that were overlapping, the total area would be:

overlap1 = bb(1,3)*bb(1,4) - areas(1)

overlap1  =

311

Remember, we can find the total area of a square by simply multiplying its width and height together. When we subtract this with the actual area taken up by the square without the circle, we get the pixels that were taken by the circle.

Similarly, for the box at the bottom:

overlap2 = bb(2,3)*bb(2,4) - areas(2)

overlap2  =

351

Step #5

Finally, what's left are the actual circles themselves. All we have to do is create a blank image, draw both circles in this image together with a logical AND operator, and find the total area of what was remaining. In the first step, these two overlapping circles correspond to: (x,y) = 83.2850, 83.2850 and (x,y) = 117.6691, 118.0193. These correspond to the last two circles detected in the Circular Hough Transform. Let's get these two circles and create our mask:

As such:

centre1 = centres(3,:);
centre2 = centres(4,:);
radii1 = radii(3);
radii2 = radii(4);
circle1 = (X - centre1(1)).^2 + (Y - centre1(2)).^2 <= radii1^2;
circle2 = (X - centre2(1)).^2 + (Y - centre2(2)).^2 <= radii2^2;
final_two = circle1 & circle2;

If we show this image, this is what we get:

This is visualizing the overlap between the two circles. What's left now is to simply count up the overlap:

overlap3 = sum(final_two(:))

overlap3 = 

515

---

Whew! That was a lot of work!

Edit

You wish to find those intersecting areas like the above with the two circles, but for the rest of the image. There are two more intersecting areas that we need to find. The first one is at the bottom of the image between the circle and square, while the second one is between the square and the circle at the top left.

To get the first one, simply create a circular mask like we did earlier and do a logical AND with the filled-in square that is touching this circle. The circle that's at the bottom is the second detected circle from the output of the Circular Hough Transform. Also, the affected square is the second one detected from regionprops. Therefore:

r = radii(2);
c = centres(2,:);
mask_circle = (X - c(1)).^2 + (Y - c(2)).^2 <= r^2;
mask_square = false(size(im));  
mask_square(floor(bb(2,2)):floor(bb(2,2)) + bb(2,4), floor(bb(2,1)):floor(bb(2,1)) + bb(2,3)) = true;
intersect1 = mask_circle & mask_square;

Here's what the first intersected area looks like:

---

You can apply the same logic above to the other square and circular areas at the top. You just need to select the right square and circle. This would be the fourth detected circle and the first detected square:

r = radii(4);
c = centres(4,:);
mask_circle = (X - c(1)).^2 + (Y - c(2)).^2 <= r^2;
mask_square = false(size(im));  
mask_square(floor(bb(1,2)):floor(bb(1,2)) + bb(1,4), floor(bb(1,1)):floor(bb(1,1)) + bb(1,3)) = true;
intersect2 = mask_circle & mask_square;

This is what we get: