how to fit ellipse to an image with cylinderical shape

Hi

I have to fit an ellipse to an image with a cylinder(see the attached image). after using the edge detection, i could recognize the ellipse shape but i couldn't find a way how to fit an ellipse on to it, or in other words i couldn't extract that part from rest of the image.

any suggestions?

Thanks!

[Adam Danz]

If you share the code that inputs original_image and produces edge_image, I might be able to play around with it later.

Answer 1

[Image Analyst]

First of all, don't do an edge detection. I don't know why people always think edge detection is the first step in any image processing algorithm. In your case it's definitely NOT. You should do a simple threshold, then scan the image to find the right edge, then get the semi axes from the values. Assuming it's perfectly aligned, try something like (untested)

mask = imbinarize(grayImage);
[rows, columns, numColors] = size(grayImage);
rightColumns = zeros(rows, 1);
for row = 1 : rows
    thisRow = mask(row, :);
    t = find(thisRow, 1, 'last');
    if ~isempty(t)
        rightColumns(row) = t;
    end
end
% Then get rid of any outliers that you may get if you're not perfectly aligned (for you to do).
% Find top and bottom (first and last elements of right rightColumns).
topRow = find(rightColumns, 1, 'first')
bottomRow = find(rightColumns, 1, 'last')
% Get ellipse semi-axes, a and b
b = (bottomRow - topRow) / 2
nonZeroIndexes = rightColumns > 0;
xCenter = min(rightColumns(nonZeroIndexes)
yCenter = mean([topRow, bottomRow])
a = max(rightColumns(nonZeroIndexes)) - min(rightColumns(nonZeroIndexes))
% Now you have the equation of an ellipse: ((x-xCenter)/a)^2 + ((y-yCenter)/b)^2 = 1

[Malik Malik]

Thank you! I have tried to implement this but seems like "t" is always empty and never get any value. am i making any mistake?

clear all;


originalImage = imread('original_image.png');

hFig = figure;
imshow(originalImage, []);
axis on;
title('Original Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);

	grayImage = originalImage;

mask = imbinarize(grayImage);
[rows, columns, numColors] = size(grayImage);
rightColumns = zeros(rows, 1);
for row = 1 : rows
    thisRow = mask(row, :);
    t = find(thisRow, 1, 'last');
    if ~isempty(t)
        rightColumns(row) = t;
    end
end
% Then get rid of any outliers that you may get if you're not perfectly aligned (for you to do).
%Fill holes

mask = imfill(mask, 'holes');
% Get rid of anything touching the edge of the image

mask = imclearborder(mask);

% Get rid of anything less than 7000 pixels
mask = bwareaopen(mask, 7000);

% Find top and bottom (first and last elements of right rightColumns).
topRow = find(rightColumns, 1, 'first');
bottomRow = find(rightColumns, 1, 'last');
% Get ellipse semi-axes, a and b
b = (bottomRow - topRow) / 2;
nonZeroIndexes = rightColumns > 0;
xCenter = min(rightColumns(nonZeroIndexes));
yCenter = mean([topRow, bottomRow]);
a = max(rightColumns(nonZeroIndexes)) - min(rightColumns(nonZeroIndexes));
% Now you have the equation of an ellipse: ((x-xCenter)/a)^2 + ((y-yCenter)/b)^2 = 1

x = xCenter + xCenter*cos(t);
 y = yCenter + yCenter*sin(t);




 plot(x,y);
 axis equal

[Image Analyst]

I think the binarization was not good. Try it manually. Use my interactive thresholding app if you want: It's in my File Exchange

% This demo characterizes the shape of the end of a rod/dowel.
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 = 15;

%--------------------------------------------------------------------------------------------------------
%    READ IN IMAGE
folder = pwd;
baseFileName = 'original_image.png';
% 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
% Read in the image from disk.  If storedColorMap is not empty, it's an indexed image with a stored colormap.
[grayImage, storedColorMap] = imread(fullFileName);
if ~isempty(storedColorMap)
	grayImage = ind2rgb(grayImage, storedColorMap);
end
% 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(grayImage);
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(grayImage);
	% 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 = grayImage(:, :, 2); % Take green channel.
end
% Display the image.
hFig = figure;
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
hFig.WindowState = 'maximized'; % May not work in earlier versions of MATLAB.
drawnow;

subplot(2, 2, 2);
imhist(grayImage);
grid on;
title('Histogram of Original Grayscale Image', 'FontSize', fontSize, 'Interpreter', 'None');

%--------------------------------------------------------------------------------------------------------
%    SEGMENTATION
% Binarize the image
% binaryImage = imbinarize(grayImage);
binaryImage = grayImage > 40;
% Take the largest blob only.
binaryImage = bwareafilt(binaryImage, 1);
% Fill holes, if any.
binaryImage = imfill(binaryImage, 'holes');
% Display the binary image.
subplot(2, 2, 3);
imshow(binaryImage, []);
title('Binary Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');

[rows, columns, numColors] = size(grayImage);
rightColumns = zeros(rows, 1);
for row = 1 : rows
    thisRow = binaryImage(row, :);
    t = find(thisRow, 1, 'last');
    if ~isempty(t)
        rightColumns(row) = t;
    end
end

% Then get rid of any outliers that you may get if you're not perfectly aligned (for you to do).
% Find top and bottom (first and last elements of right rightColumns).
topRow = find(rightColumns, 1, 'first')
bottomRow = find(rightColumns, 1, 'last')
% Get ellipse semi-axes, a and b
b = (bottomRow - topRow) / 2
nonZeroIndexes = rightColumns > 0;
xCenter = min(rightColumns(nonZeroIndexes))
yCenter = mean([topRow, bottomRow])
a = max(rightColumns(nonZeroIndexes)) - min(rightColumns(nonZeroIndexes))
% Now you have the equation of an ellipse: ((x-xCenter)/a)^2 + ((y-yCenter)/b)^2 = 1

[Malik Malik]

Hi

Thanks again, i have tried to improve the binarized image but "t" still came out empty.

regards

[Image Analyst]

Of course t will be empty for rows where there are no white pixels. Don't worry about that. My code contains this code which takes care of that:

    if ~isempty(t)
        rightColumns(row) = t;
    end    

Did you actually run my code, or just look at it and try to copy over some bits and pieces of it here and there into your existing code?

[Malik Malik]

Yes i run your code. and got the value for a and b but to calculate x and y for elipse i need "t" , i think.

x = xCenter + xCenter*cos(t);
y = yCenter + yCenter*sin(t);

[Image Analyst]

That's a different t. That t is the angle. And they're the wrong equations. I think they should be

t = linspace(0, 360, 500);
x = xCenter + a * cosd(t);
y = yCenter + b * sind(t);

Try that. I've got a meeting with the Mathworks over the next two hours. If it doesn't work, let me know and I'll look at it after my Mathworks board meeting.

[Malik Malik]

It does generate the ellipse but its bigger and not fit over the ellipitical part of the image. See the attached results.

Thanks for the help. any suggestions?

Answer 2

[Hiro Yoshino]

I am not an expert - please see the following page and hopefully there are some similar examples for you.

https://jp.mathworks.com/help/images/image-segmentation.html

If you are fancy to app, the try this out:

https://jp.mathworks.com/help/images/segment-image-using-auto-cluster.html

[Malik Malik]

I don't want to select the points by myself. i want to do is that after edge detection ellipitical shape detected seperated from rest of the image as ellipse.