How to plot cirlces around this spots..?

I have this image of blood cells. I want to draw circles around it..

The code is :

im = imread('image.jpg');

im = im(:,:,3);

minR=34; maxR=40; thresh=0.33; delta=12;

origim = im;

if length(size(im))>2

im = rgb2gray(im);

end

maxR2 = 2*maxR;

hough = zeros(size(im,1)+maxR2, size(im,2)+maxR2, maxR-minR+1);

[X Y] = meshgrid(0:maxR2, 0:maxR2);

Rmap = round(sqrt((X-maxR).^2 + (Y-maxR).^2));

Rmap(Rmap<minR | Rmap>maxR) = 0;

edgeim = edge(im, 'canny', [0.1 0.2]);

[Ey Ex] = find(edgeim);

[Cy Cx R] = find(Rmap);

for i = 1:length(Ex);

Index = sub2ind(size(hough), Cy+Ey(i)-1, Cx+Ex(i)-1, R-minR+1);

hough(Index) = hough(Index)+1;

end

twoPi = 0.9*2*pi;

circles = zeros(0,4);

for radius = minR:maxR

slice = hough(:,:,radius-minR+1);

twoPiR = twoPi*radius;

slice(slice<twoPiR*thresh) = 0; 

[y x count] = find(slice)

circles = [circles; [x-maxR, y-maxR, radius*ones(length(x),1), count/twoPiR]];

end

circles = sortrows(circles,-4);

i = 1;

while i<size(circles,1)

j = i+1;

while j<=size(circles,1)

    if sum(abs(circles(i,1:3)-circles(j,1:3))) <= delta

      circles(j,:) = [];

    else

      j = j+1;

    end

end

i = i+1;

end

if nargout==0 % Draw circles

figure, imshow(origim), hold on;

for i = 1:size(circles,1)

    x = circles(i,1)-circles(i,3);

    y = circles(i,2)-circles(i,3);

    w = 2*circles(i,3);

    rectangle('Position', [x y w w], 'EdgeColor', 'red', 'Curvature', [1 1]);

end

hold off;

end

But it is not working.. Can anyone help to solve this..? please help..

[Star Strider]

Image processing isn’t an area of my expertise, but your image is.

What do you mean by ‘spots’? You have one polymorphonuclear leukocyte, and a large number of largely abnormal erythrocytes. The small ‘spots’ within the red cells are Howell-Jolly bodies (remnant DNA normally eliminated in the spleen), and the red cells that do not contain them are spherocytes (smaller than normal red cells lacking the normal central concavity).

Do you want to identify the red cells (and estimate their sizes), the Howell-Jolly bodies, or something else?

[Nimisha]

Yes i want to detect leukocyte and erythrocytes. And need to count them.. i tried using above code.

but it not work.. it seems correct to me.. then why..?

do youknow any hough transform method for it..?

Answer 1

[Dasharath Gulvady]

Nimisha, By executing your code, I could see that the resulting "circles" is a 0-by-1 empty matrix. It looks like the issue is with the values of the variables “thresh”,” minR”, “maxR” which are hardcoded. Refer to my comments below for one of the for-loops, where I think the issue is:

 for radius = minR:maxR

    slice = hough(:,:,radius-minR+1); % why is minR hard coded?
    twoPiR = twoPi*radius ;  

    % values in "slice" are small compared to the value "twoPiR*thresh". Hence the       resulting "slice" will contain all zeros. 
    slice(slice<twoPiR*thresh) = 0; 

     % Since "slice" has all zeros, x and y are empty matrices
    [y x count] = find(slice);

    % Since x and y are empty matrices, "circles" is an empty matrix
    circles = [circles; [x-maxR, y-maxR, radius*ones(length(x),1), count/twoPiR]]; 
    %circles
end

I would suggest taking a look at the hard coded values at the beginning of your code and change them appropriately. Since the code is pretty big, feel free to add some comments in your code which would help in understanding what exactly you are trying to achieve and how you are doing it.

Answer 2

[Image Analyst]

Have you tried imfindcircles? Here is Brett's demo: http://www.mathworks.com/matlabcentral/fileexchange/34365-findcirclesgui (Brett works at The Mathworks)

But what do you REALLY want to do? I'm sure drawing circles it not it - that's not useful. Do you want the count, the average equivalent circular diameter, the area fraction, or what??? How about doing color segmentation like I show you in my File Exchange. Or maybe you want to do something like this: http://www.mathworks.com/products/demos/image/color_seg_k/ipexhistology.html

[Nimisha]

im = imread('image.jpg');

im = rgb2gray(im);

e = edge(im, 'canny');

radii = 20:1:55;

b = e;

rrange = radii;

% get indices of non-zero features of b

[featR, featC] = find(b);

% set up accumulator array - with a margin to avoid need for bounds

[nr, nc] = size(b);

nradii = length(rrange);

margin = ceil(max(rrange));

nrh = nr + 2*margin; % increase size of accumulator

nch = nc + 2*margin;

h = zeros(nrh*nch*nradii, 1, 'uint32'); % 1-D for now, uint32 a touch

% get templates for circles at all radii - these specify accumulator

% elements to increment relative to a given feature

tempR = []; tempC = []; tempRad = [];

for i = 1:nradii

    r = rrange(i);

    l = round(r/sqrt(2));

if round(sqrt(r.^2 - l.^2)) < l % if crosses diagonal

    l = l-1;

end

% generate coords for 1/8 of the circle, a dot on each row

x0 = 0:l;

y0 = round(sqrt(r.^2 - x0.^2));

% Check for overlap

if y0(end) == l

    l2 = l;

else

    l2 = l+1;

end

% assemble first quadrant

x = [x0 y0(l2:-1:2)];

y = [y0 x0(l2:-1:2)];

% add next quadrant

x0 = [x y];

y0 = [y -x];

% assemble full circle

x = [x0 -x0];

y = [y0 -y0];

end

tR = x;

tC = y;

        tempR = [tempR tR]; %#ok<*AGROW>

    tempC = [tempC tC];

    tempRad = [tempRad repmat(i, 1, length(tR))];

% Convert offsets into linear indices into h - this is similar to sub2ind.

tempInd = uint32( tempR+margin + nrh*(tempC+margin) + nrh*nch*(tempRad-1) );

featInd = uint32( featR' + nrh*(featC-1)' );

% Loop over features

for f = featInd

    % shift template to be centred on this feature

    incI = tempInd + f;

    % and update the accumulator

    h(incI) = h(incI) + 1;

end

% Reshape h, convert to double, and apply options

h = reshape(double(h), nrh, nch, nradii);

    h = h(1+margin:end-margin, 1+margin:end-margin, :);

    margin = 0;

   H = bsxfun(@rdivide, h, reshape(rrange, 1, 1, length(rrange)));

ip = inputParser;

% required

htest = @(h) validateattributes(h, {'double'}, {'real' 'nonnegative' 'nonsparse'});

ip.addRequired('h', htest);

rtest = @(radii) validateattributes(radii, {'double'}, {'real' 'positive' 'vector'});

ip.addRequired('radii', rtest);

% optional

mtest = @(n) validateattributes(n, {'double'}, {'real' 'nonnegative' 'integer' 'scalar'});

ip.addOptional('margin', 0, mtest);

% parameter/value pairs

stest = @(s) validateattributes(s, {'double'}, {'real' 'nonnegative' 'scalar'});

ip.addParamValue('smoothxy', 0, stest);

ip.addParamValue('smoothr', 0, stest);

ip.addParamValue('threshold', [], stest);

nptest = @(n) validateattributes(n, {'double'}, {'real' 'positive' 'integer' 'scalar'});

ip.addParamValue('npeaks', [], nptest);

nhtest = @(n) validateattributes(n, {'double'}, {'odd' 'positive' 'scalar'});

ip.addParamValue('nhoodxy', [], nhtest);

ip.addParamValue('nhoodr', [], nhtest);

ip.parse(h, radii, 'nhoodxy', 15, 'nhoodr', 21, 'npeaks', 10);

params = ip.Results;

% smooth the accumulator - xy

if params.smoothxy > 0

    [m, hsize] = gaussmask1d(params.smoothxy);

    % smooth each dimension separately, with reflection

    h = cat(1, h(hsize:-1:1,:,:), h, h(end:-1:end-hsize+1,:,:));

    h = convn(h, reshape(m, length(m), 1, 1), 'valid');

    h = cat(2, h(:,hsize:-1:1,:), h, h(:,end:-1:end-hsize+1,:));

    h = convn(h, reshape(m, 1, length(m), 1), 'valid');

end

% smooth the accumulator - r if params.smoothr > 0

    sigma = params.smoothhr;

    hsize = ceil(2.5*sigma);  % reasonable truncation

x = (-hsize:hsize) / (sqrt(2) * sigma);

m = exp(-x.^2);

m = m / sum(m); % normalise

    h = cat(3, h(:,:,hsize:-1:1), h, h(:,:,end:-1:end-hsize+1));

    h = convn(h, reshape(m, 1, 1, length(m)), 'valid');

end

% set threshold

if isempty(params.threshold)

    params.threshold = 0.5 * max(h(:));

end

if isempty(params.nhoodxy) && isempty(params.nhoodxy)

    % First approach to peak finding: local maxima

    % find the maxima

    maxarr = imregionalmax(h);

    maxarr = maxarr & h >= params.threshold;

    % get array indices
    peakind = find(maxarr);

    [y, x, rind] = ind2sub(size(h), peakind);

    peaks = [x'; y'; radii(rind)];

    % get strongest peaks

    if ~isempty(params.npeaks) && params.npeaks < size(peaks,2)

        [~, ind] = sort(h(peakind), 'descend');

        ind = ind(1:params.npeaks);

        peaks = peaks(:, ind);

    end

else

    % Second approach: iterative global max with suppression

    if isempty(params.nhoodxy)

        params.nhoodxy = 1;

    elseif isempty(params.nhoodr)

        params.nhoodr = 1;

    end

    nhood2 = ([params.nhoodxy params.nhoodxy params.nhoodr]-1) / 2;

    if isempty(params.npeaks)

        maxpks = 0;

        peaks = zeros(3, round(numel(h)/100));  % preallocate

    else

        maxpks = params.npeaks;  

        peaks = zeros(3, maxpks);  % preallocate

    end

    np = 0;

    while true

        [vr, r] = max(h);

[vc, c] = max(vr);

[v, k] = max(vc);

c = c(1, 1, k);

r = r(1, c, k);

        % stop if peak height below threshold

        if v < params.threshold || v == 0

            break;

        end

        np = np + 1;

        peaks(:, np) = [c; r; radii(k)];

        % stop if done enough peaks

        if np == maxpks

            break;

        end

        % suppress this peak

        r0 = max([1 1 1], [r c k]-nhood2);

        r1 = min(size(h), [r c k]+nhood2);

        h(r0(1):r1(1), r0(2):r1(2), r0(3):r1(3)) = 0;

    end 

    peaks(:, np+1:end) = [];   % trim

end

% adjust for margin

if params.margin > 0

    peaks([1 2], :) = peaks([1 2], :) - params.margin;

end

% We draw the circles found on the image, using both the positions and the

% radii stored in the peaks array.

imshow(im);

hold on;

for peak = peaks

    r = peak(3);

    l = round(r/sqrt(2));

if round(sqrt(r.^2 - l.^2)) < l % if crosses diagonal

    l = l-1;

end

% generate coords for 1/8 of the circle, a dot on each row

x0 = 0:l;

y0 = round(sqrt(r.^2 - x0.^2));

% Check for overlap

if y0(end) == l

    l2 = l;

else

    l2 = l+1;

end

% assemble first quadrant

x = [x0 y0(l2:-1:2)];

y = [y0 x0(l2:-1:2)];

% add next quadrant

x0 = [x y];

y0 = [y -x];

% assemble full circle

x = [x0 -x0];

y = [y0 -y0];

    plot(x+peak(1), y+peak(2), 'g-');

end

hold off

% For above code i got following resultant image.. but it detects wrong circles..? please give solution...

[Image Analyst]

Just attach the script with the paper clip icon, or else read this: http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup