Search code examples
matlabimage-processinggeometrypixelbinary-image

Count black pixels within a circle

I have one vector of radii and second vector of hundreds of [X,Y] coordinates. For each possible coordinate-radius pair I have count all black pixels within a circle (whose center is placed in the coordinate) on an input binary image.

What is the fastest way to do it? My only idea is to iterate through every pixel of an image, check the circle equation and then the pixel color but it doesn't seem much optimized for several hundred such operations.

Solution

  • Matlab is great for working with images thanks to the matrix syntax. It does also work with indices so most time you can avoid "iterating through pixels" (although sometimes you'll still have to).

    Instead of checking all the pixels within each circle, and having to detect how many pixels were counted twice, another approach is to create a mask, the same size of you image. Blank this mask for each of your circles (so overlapping pixels are only 'blanked' once), then apply the mask on your original picture and count the remaining illuminated pixels.

    For an example, I have to take some sample data, the image:

    load trees
BW = im2bw(X,map,0.4);
imshow(BW)

    treesBW

    And 20 random point/circle coordinates (you can change the number of points and the min/max radii easily):

    %// still building sample data
s = size(BW) ;
npt = 20 ; Rmin=5 ; Rmax=20 ; %// problem constants

x = randi([1 s(2)]   ,npt,1); %// random X coordinates
y = randi([1 s(1)]   ,npt,1); %//        Y
r = randi([Rmin Rmax],npt,1); %// radius size between 5 to 20 pixels.

    Then we build your custom mask :

    %// create empty mask with enough overlap for the circles on the border of the image
mask = false( s+2*Rmax ) ; 

%// prepare grid for a sub-mask of pixels, as wide as the maximum circle
xmask = -Rmax:Rmax ; 
[xg,yg] = ndgrid(xmask,xmask) ;
rg = sqrt( (xg.^2+yg.^2) ) ;    %// radius of each pixel in the subgrid

for ip=1:npt
    mrow = xmask+Rmax+y(ip) ;  %// calc coordinates of subgrid on original mask
    mcol = xmask+Rmax+x(ip) ;  %// calc coordinates of subgrid on original mask

    cmask = rg <= r(ip) ;      %// calculate submask for this radius
    mask(mrow,mcol) = mask(mrow,mcol) | cmask ; %// apply the sub-mask at the x,y coordinates
end
%// crop back the mask to image original size (=remove border overlap)
mask = mask(Rmax+1:end-Rmax,Rmax+1:end-Rmax) ;
imshow(mask)

    bubbles

    Then you just apply the mask and count :

    %% // Now apply on original image
BWm = ~BW & mask ; %// note the ~ (not) operator because you want the "black" pixels
nb_black_pixels = sum(sum(BWm)) ;
imshow(BWm)

nb_black_pixels =
        5283

    treeMasked