How to fit circles to data points in irregular images

Im trying to write an automated analysis script for microscopic images of white blood cells, meaning that every cell is slightly different and it's radii might vary from 5 to 15 microns..

rgb = imread('1p-40x-400ms-DS-Fi2-3.jpg'); I = rgb2gray(rgb);

hy = fspecial('sobel'); hx = hy; Iy = imfilter(double(I), hy, 'replicate'); Ix = imfilter(double(I), hx, 'replicate'); gradmag = sqrt(Ix.^2 + Iy.^2);

C = zeros(size(gradmag)); [m,n]=size(gradmag); count = 0; Max = max(max(gradmag)); for i = 1:m for j = 1:n if gradmag(i,j) < 20; C(i,j) = 0; % count = count + 1; else C(i,j) = 225; %count = count + 1; end end end

edges = edge(C, 'log');

figure(1) imshow(edges)

by utilizing the code above I am able to generate the image named "original" which is attached. I would like to be able to generate circles and fit them onto the image with as low overlap as possible generating a result something like "wanted".

The end goal is to fill the circles and calculate an estimate of the surface coverage, see "goal".

I have tried utilizing imopen/erode/close without significant improvements, i've also tried watershedding but since the edges are very weak at its points (in the microscopic image) and the "original" image contains objects that isnt "closed".

Answer 1

[Will Nitsch]

So I modified an example on the MathWorks website to roughly accomplish your goal. Ultimately I get a bunch of cells whose measurements are stored in the output of 'regionprops', in the 'cellProps' array. From here you could get centroid locations and you can get the equivalent diameter as well to model a circle as you planned. Then you can use 'viscircles' to see them on your image.

I have some experience in microscope image processing as well, and I have a few suggestions. First, try to remove the reticle from your images. It messes up the dynamic range and makes it a bit harder to process the images. Also, try to correct the image so that the field is relatively flat across it. It can be a lot of work writing automation algorithms for microscopy image processing since a lot of cells are relatively amorphous and because you get a lot of image noise, vignetting/illumination gradients, and sometimes foreign particles. In your case, you have heavy cell occlusion in the center of the image, which isn't impossible to overcome but you need the right algorithm. There are some published online. A book I can recommend is 'Microscope Image Processing' by Wu, Merchant and Castleman. Anyways, here is the algorithm adapted from this example:

https://blogs.mathworks.com/steve/2006/06/02/cell-segmentation/

    rgb = imread('1p40x-400ms-DS-Fi2-3.jpg');

    I = rgb2hsv(rgb);
    % Try using the Value channel of the HSV colorspace
    Iv = I(:,:,3); %better than your grayscale img

    % Kill the reticle bar
    Iv(Iv>0.85) = mean2(Iv); % I wouldn't do this in every case but if you can, get rid of the reticles in your image, they throw off the dynamic range and make things weirder

    se = strel('disk',50);
    % Try getting rid of some of the haze with a top hat filter
    tophatFiltered = imtophat(Iv,se);

    figure
    % Do a little prep work. Gaussian blur, adaptive histeq, and sharpen
    I2 = imsharpen(adapthisteq(imfilter(tophatFiltered,fspecial('gaussian',9,3))));

    % Auto threshold
    bw = im2bw(I2, graythresh(I2));
    imshow(bw)

    bw2 = imfill(bw,'holes');
    bw3 = imopen(bw2, ones(5,5));
    bw4 = bwareaopen(bw3, 40);
    bw4_perim = bwperim(bw4);
    overlay1 = imoverlay(I2, bw4_perim, [.3 1 .3]);
    imshow(overlay1)

    mask_em = imextendedmax(I2, 35/255);
    figure
    imshow(mask_em)

    figure

    mask_em = imclose(mask_em, ones(5,5));
    mask_em = imfill(mask_em, 'holes');
    mask_em = bwareaopen(mask_em, 40);
    overlay2 = imoverlay(I2, bw4_perim | mask_em, [.3 1 .3]);
    imshow(overlay2)

    I_eq_c = imcomplement(I2);

    I_mod = imimposemin(I_eq_c, ~bw4 | mask_em);

    L = watershed(I_mod);
    imshow(label2rgb(L))

   % Use the label matrix for region definition, and Iv for intensity measurements
    cellProps = regionprops(L,Iv,'all');

[Daniel Leppälä]

Thanks alot! with a little bit of work I managed to extract almost exactly what I wanted from that piece of code!

The image supplied in this question is from one of the earlier imaging sessions, in later sessions I have taken identical images with and without the scale for analysis/documentation purposes.