How do I find a center of mass of imperfect objects in a binary picture?

I have a set of binary images with dashed circular objects (circles within circles, imperfect circles) (as seen in the attached picture). I need to find the coordinates (x,y) of the center of mass of the inner circle most of each object.

Can anyone help me?

Answer 1

[Image Analyst]

Simple.

1. Threshold the image to form a binary image.
2. Do a morphological closing to close gaps
3. Fill holes
4. Throw out small blobs
5. Take convex hull
6. Call regionprops to find centroids.

Try this:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

%===============================================================================
% Read in gray scale demo image.
folder = pwd; % Determine where demo folder is (works with all versions).
baseFileName = 'imperfect circles.jpg';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
	% The file doesn't exist -- didn't find it there in that folder.
	% Check the entire search path (other folders) for the file by stripping off the folder.
	fullFileNameOnSearchPath = baseFileName; % No path this time.
	if ~exist(fullFileNameOnSearchPath, 'file')
		% Still didn't find it.  Alert user.
		errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
		uiwait(warndlg(errorMessage));
		return;
	end
end
rgbImage = imread(fullFileName);
% Display the image.
subplot(2, 3, 1);
imshow(rgbImage, []);
title('Original Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
hp = impixelinfo();

% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(rgbImage)
if numberOfColorChannels > 1
	% It's not really gray scale like we expected - it's color.
	% Use weighted sum of ALL channels to create a gray scale image.
	% 		grayImage = rgb2gray(rgbImage);
	% ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
	% which in a typical snapshot will be the least noisy channel.
	grayImage = rgbImage(:, :, 2); % Take green channel.
else
	grayImage = rgbImage; % It's already gray scale.
end
% Now it's gray scale with range of 0 to 255.

% Display the histogram of the image.
subplot(2, 3, 2);
[counts, binLocations] = imhist(grayImage);
% Suppress bin 1 because it's so tall
counts(1) = 0;
bar(binLocations, counts);
grid on;
title('Histogram of Image', 'FontSize', fontSize, 'Interpreter', 'None');

%------------------------------------------------------------------------------
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
drawnow;

% Binarize the image
% Get the mask where the region is solid.
binaryImage = grayImage > 128;

% Crop off last 2 lines.  For some reason, the next to the last line is all white.  Set them equal to false.
binaryImage(end-1:end, :) = false;	% Blacken last 2 lines.

% Do a morphological closing to connect lines
se = strel('disk', 5, 0);
binaryImage = imclose(binaryImage, se);
% Fill blobs:
binaryImage = imfill(binaryImage, 'holes');
% Display the image.
subplot(2, 3, 3);
imshow(binaryImage, []);
title('Initial Binary Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
hp = impixelinfo();
drawnow;

% Extract only those larger than 1000 pixels in area
binaryImage = bwareaopen(binaryImage, 1000);
% Display the image.
subplot(2, 3, 4);
imshow(binaryImage, []);
title('Closed Binary Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
drawnow;

% Extract only those larger than 1000 pixels in area
binaryImage = bwconvhull(binaryImage, 'objects');
% Display the image.
subplot(2, 3, 5);
imshow(binaryImage, []);
title('Final Binary Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
drawnow;

% Use it to mask the original image.
finalImage = grayImage; % Initialize
finalImage(~binaryImage) = 0; % Erase outside the mask.
% Display the image.
subplot(2, 3, 6);
imshow(finalImage, []);
title('Final, Masked Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');

% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
% Plot the borders of all the clusters on the original grayscale image using the coordinates returned by bwboundaries.
title('Outlines, from bwboundaries()', 'FontSize', fontSize); 
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
hold on;
boundaries = bwboundaries(binaryImage);
numberOfBoundaries = size(boundaries, 1);
for k = 1 : numberOfBoundaries
	thisBoundary = boundaries{k};
	plot(thisBoundary(:,2), thisBoundary(:,1), 'g', 'LineWidth', 2);
end
hold off;

% Find the centroids
props = regionprops(binaryImage, 'Centroid', 'EquivDiameter');
xyCentroids = [props.Centroid];
xCentroids = xyCentroids(1:2:end)
yCentroids = xyCentroids(2:2:end)
% Plot centroids over the image with a large red cross.
hold on;
for k = 1 : length(xCentroids)
	thisX = xCentroids(k);
	thisY = yCentroids(k);
	thisDiameter = props(k).EquivDiameter;
	plot(thisX, thisY, 'r+', 'MarkerSize', thisDiameter, 'LineWidth', 2);
end

[inon hen]

Thank you. The program is excellent but it can not identify the centers of the objects in cases where there is an object connected to another object. As in the attached picture. (I need to find the center of each individual object)

can you help me?

[Image Analyst]

This is not an easy problem to solve. Let's say you have a topographic map of the rock mountains. Now, tell me how many peaks there are. Do you see the problem. Why don't you just use impts() or ginput() and ask the user? You'll probably say "I need it to be automatic", well if you take a year to automate it and it will save your users only a few minutes or hours, then was that time well spent? Do you want to spend your time analyzing your sample for whatever the real reason for your research is, or do you want to spend your time perfecting a 2-D peak analysis algorithm for months just to save a few minutes of analysis time? If I were you, I'd be practical and just locate the centers manually and be done with it.