MATLAB Answers

How to fit circles to data points in irregular images

28 views (last 30 days)
Daniel Leppälä
Daniel Leppälä on 27 Apr 2017
Edited: Will Nitsch on 1 Sep 2017
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".

Answers (1)

Will Nitsch
Will Nitsch on 1 May 2017
Edited: Will Nitsch on 1 Sep 2017
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');
  1 Comment
Daniel Leppälä
Daniel Leppälä on 2 May 2017
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!