I am trying to find the radius of an arc in an image.

Hello all,

I have a series of images with circular segments on them. I am trying to determine the radius of the arc and the theoretical center point of the circle. My first attempt has been using 'imfindcircles' with no luck. I suspect the function prefers that the entire circle be filled in and not just a boundary. I am including an example image for reference. The other images look very similar to the first except with smaller arc portions visible. Any help on this would be greatly appreciated.

Answer 1

[Image Analyst]

Try this. I'm using find() to get the x and y of the mask rather than a double for loop like @William Rose. The computed radius is 288.54 pixels.

% Demo by Image Analyst
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 = 16;
markerSize = 40;
%--------------------------------------------------------------------------------------------------------
%    READ IN TEST IMAGE
folder = [];
baseFileName = 'arc_image.png';
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
grayImage = imread(fullFileName);
% 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.
    fprintf('It is not really gray scale like we expected - it is color\n');
    % Extract the blue channel.
    grayImage = grayImage(:, :, 3);
end

It is not really gray scale like we expected - it is color

%--------------------------------------------------------------------------------------------------------
% Display the image.
subplot(2, 2, 1);
imshow(grayImage, []);
impixelinfo;
axis('on', 'image');
title('Original Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
% Update 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);
% Maximize window.
g = gcf;
g.WindowState = 'maximized';
drawnow;

%--------------------------------------------------------------------------------------------------------
% Get mask by thresholding at 116.
lowThreshold = 116;
highThreshold = 254;
% Interactively and visually set a threshold on a gray scale image.
% https://www.mathworks.com/matlabcentral/fileexchange/29372-thresholding-an-image?s_tid=srchtitle
% [lowThreshold, highThreshold] = threshold(lowThreshold, highThreshold, grayImage)
circleMask = grayImage >= lowThreshold & grayImage <= highThreshold;
subplot(2, 2, 2);
imshow(circleMask);
impixelinfo;
axis('on', 'image');
title('Initial Mask Image', 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
% Clean up the initial mask
circleMask = bwareafilt(circleMask, 1);
circleMask = imfill(circleMask, 'holes');
subplot(2, 2, 3);
imshow(circleMask);
impixelinfo;
axis('on', 'image');
title('Final Mask Image', 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;

% Get the x and y coordinates of the circle mask
[y, x] = find(circleMask);
% Fit them to a circle.
[xCenter, yCenter, radius, a] = circlefit(x, y)

xCenter =           218.230585096279

yCenter =           376.506380841752

radius =           288.541071609979

a = 3×1

1.0e+00 * -436.461170192557 -753.012761683504 106125.693080183

% Display the original image in the lower right of the figure.
subplot(2, 2, 4);
imshow(grayImage);
impixelinfo;
axis('on', 'image');
% View the circle over the original image
viscircles([xCenter, yCenter], radius, 'Color', 'red', 'LineWidth', 2);
caption = sprintf('Radius = %.2f pixels', radius);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

drawnow;

%===========================================================================================================================================
function [xCenter, yCenter, radius, a] = circlefit(x, y)
% circlefit(): Fits a circle through a set of points in the x - y plane.
% USAGE :
% [xCenter, yCenter, radius, a] = circlefit(X, Y)
% The output is the center point (xCenter, yCenter) and the radius of the fitted circle.
% "a" is an optional output vector describing the coefficients in the circle's equation:
%     x ^ 2 + y ^ 2 + a(1) * x + a(2) * y + a(3) = 0
% by Bucher Izhak 25 - Oct - 1991
numPoints = numel(x);
xx = x .* x;
yy = y .* y;
xy = x .* y;
A = [sum(x),  sum(y),  numPoints;
    sum(xy), sum(yy), sum(y);
    sum(xx), sum(xy), sum(x)];
B = [-sum(xx + yy) ;
    -sum(xx .* y + yy .* y);
    -sum(xx .* x + xy .* y)];
a = A \ B;
xCenter = -.5 * a(1);
yCenter = -.5 * a(2);
radius  =  sqrt((a(1) ^ 2 + a(2) ^ 2) / 4 - a(3));
end

[William Rose]

@Image Analyst, nice!

[Josh Morsell]

Thanks for the help!

Answer 2

[William Rose]

@Josh Morsell,

I opened your image with Microsoft Paint. I cropped off the whilte border regions and I pasted on two black rectangles, to cover up the "Time" label and to cover up the cluster of points that are not on the arc. I saved the image. The resulting image is attached. In Matlab, I opened the image.

imfile='arc_image1.png';
A=imread(imfile);

Since it looks like gray scale, I made an array which is the red chanel only:

Ar=A(:,:,1);

selected the red channel only. I viewed the gray scale values with histogram, to learn what would be a good threshold value for finding the white and light gray points on the arc,and I zoomed in on the y-axis.

histogram(Ar); ylim([0 1000]);

This shows that the gray and white points have a value of about 140 or higher, so I use 140 as a threshold. I get the image size:

[rows,cols]=size(Ar);

I loop over all the points in the image. If the pixel value exceed the threshold, I save the point's x,y cooridnates. This will give me vectors of the x,y coordinates of all the ponts on the arc.

k=1;
for i=1:rows
    for j=1:cols
        if Ar(i,j)>=140
            xy(k,:)=[i,j]; k=k+1;
        end
    end
end

I realize the code above could be shortened and probably done a lot more efficiently, but it works and it is understandable. Now I have an array containing the arc points, and I want to fit a circle to the points. I wrote a function (attached) that implements the cicle-fit method of Randy Bullock (2006, 2017), described here and here. I also tried using singular value decomposition, but it was not as robust as Bullock's method, when I tried it on other data sets. See W. Gander, G.H. Golub, R. Strebel, "Least squares fitting of circles and ellipses", BIT 34: 558-578, 1994.

[r,xc,yc]=fitcircle1(xy);
fprintf('r=%.2f, xc=%.2f, yc=%.2f.\n',r,xc,yc);

It gives the result shown:

r=288.53, xc=326.51, yc=67.24.

Good luck with your work.