How do I get two circles from a thick, noisy arc?

I have a noisy picture and i need to get two circles like this picture

[Rik]

Have a read here and here. It will greatly improve your chances of getting an answer.

What have you tried already?

[bassem ltaief]

i have binirayzed and filtering with canny and i start with poly curve fitting

[Rik]

Since you want to fit a circle, why did you try to fit a polynomial? There is no polynomial that will fit a circle. If you want to fit a circle, you essentially want to find the center coordinates and the radius that will minimize the distance between the points and the resulting circle (which you can do by substracting the distance between each point and the centroid from the radius).

If you attach the binary image I might be able to give you a bit more specific guidance.

[bassem ltaief]

this is the binary image

i don't have the centre of circle and the raduis.

from that picture i need to complite the circle and after that i can calcule the centre and the raduis of the two circules

[Rik]

You want to find those three values. You want to fit data to a function. How would you do that?

[bassem ltaief]

That band should be limited by two concentric circles part.

i need to fit a circular band to it when i have a complet circle i can calculate the diameter of each circle and the centre of this two circles have the same coordinates.

[Rik]

You could first fit a circle to all data, and then calculate the distance to the center for every point. A histogram should be able to show the two radii you are interested in.

[Bjorn Gustavsson]

@bassem: Why haven't you tried the solution outlined below? There's more information in your original image than the binarized version, why not use that information to something useful instead of throwing it away? The image appears to be from some kind of interferometric system with circular fringes, the intensity of those fringes have a known intensity-variation that depends on the finess of the interferometer, but for simplicity (or if you cannot be bothered to look it up) a sinc-squared would to as a first approximation. Fitting the centroid and width of that wouldn't be overly complicated - as outlined. Once that's done you have a couple of file-exchange tools for circle-fitting: circle-fitting at FEX

[bassem ltaief]

no is not like you think.

i try it but i have error and i work in fix it but in same time i search to find any other solution with @rik

[Bjorn Gustavsson]

OK "it is not like you think". Could you, pretty please, upload the original image, then I'll get you a solutino.

[bassem ltaief]

Answer 1

[Bjorn Gustavsson]

This is how I'd go about it:

(Perhaps 0, apply some low-pass filtering to reduce the speckle and make image smoother.)

1, find the row-indices for the intensity maxima:

[I_max,idx_max] = max(Img,[],2);

2, fit a sinc^2 to each row:

fitfcn = @(x,pars) pars(1)*sinc((x-pars(2))/pars(3)).^2;
errfcn = @(pars,I,x) sum((I-fitfcn(x,pars)).^2);
resfcn = @(pars,I,x) (I-fitfcn(x,pars));
for i1 = 1:numel(I_max)
  parsFMS(i1,:) fminsearch(@(pars) errfcn(pars,Img(i1,:),1:size(Img,2)),...
                           [I_max(i1),idx_max(i1),10]);
  % Here you could just as well use lsqnonlin and resfcn instead, which might be better...
end

If that works you then have to decide if you want the first zeros from the main-lobe of the sinc^2 or something like the half-widths, but that is for you to chose.

3, (likely) the points to the right and left "first zeroes" will not necessarily form neat arcs, but you can use some of the circle-fitting functions on the file exchange to fit two circles to those sets of points.

HTH

[Bjorn Gustavsson]

Here's a full step-by-step version:

fitfcn = @(x,pars) pars(1)*sinc((x-pars(2))/pars(3)).^2 + pars(4);
resfcn = @(pars,I,x)(I-fitfcn(x,pars));
Im = imread('image.png');
Im = double(Im);
fIm = imfilter(Im,conv2(ones(9)/9^2,ones(15)/15^2),'symmetric');
% If imfilter not available a similar enough filter can be done with conv2
[I_max,idx_max] = max(fIm,[],2);
lsqopts = optimoptions('lsqnonlin');
lsqopts.Display = 'off';
for i1 = 1:numel(I_max)
  parsLSQ(i1,:) = lsqnonlin(@(pars) resfcn(pars,fIm(i1,idx_max(i1)+[-200:200]),idx_max(i1)+[-200:200]),...
                           [I_max(i1),idx_max(i1),10,5e3],...
                           [0 idx_max(i1)-100 0 -1e4],...
                           [20*max(Im(:)) idx_max(i1)+100 150 2e4],...
                           lsqopts);
end
%% illustrate one fit, row 800
figure
i1 = 800;
plot(idx_max(i1)+[-200:200],fIm(i1,idx_max(i1)+[-200:200]),'b','linewidth',2)
hold on
plot(idx_max(i1)+[-200:200],Im(i1,idx_max(i1)+[-200:200]),'b.')
plot(idx_max(i1)+[-200:200],fitfcn(idx_max(i1)+[-200:200],parsLSQ(i1,:)),'r-','linewidth',2)
plot(parsLSQ(i1,2),parsLSQ(i1,1)+parsLSQ(i1,4),'gs','linewidth',2) 
plot(parsLSQ(i1,2)-parsLSQ(i1,3),parsLSQ(i1,4),'gs','linewidth',2) 
plot(parsLSQ(i1,2)-parsLSQ(i1,3),parsLSQ(i1,4),'gs','linewidth',2)
%% Points on image, centre, left and right first minima
figure
imagesc(Im)
hold on
plot(parsLSQ(:,2),1:size(Im,1),'c.')
plot(parsLSQ(:,2)-parsLSQ(:,3),1:size(Im,1),'m.')
plot(parsLSQ(:,2)+parsLSQ(:,3),1:size(Im,1),'r.')
%% Circle-fits
[xcM ycM RM] = circfit(parsLSQ(:,2)-parsLSQ(:,3),1:numel(idx_max));
[xcP ycP RP] = circfit(parsLSQ(:,2)+parsLSQ(:,3),1:numel(idx_max));
[xcC ycC RC] = circfit(parsLSQ(:,2),1:numel(idx_max));             
%% Plot the circle-centres
plot(xcC,ycC,'c.','markersize',18)                    
axis auto
axis equal
plot(xcP,ycP,'r.','markersize',18)
plot(xcM,ycM,'m.','markersize',18)

There's more to be done. Firstly: check that the fit actually converges, check the help and documentation for lsqnonlin to do that. When the output flag for convergence indicates that the fit is not optimal, just restar the optimization from the current best point, this can be put into a while-loop. Another way to handle this is to simply remove those results from the circle-fitting.

Use a more appropriate line-profile than a sinc^2.

Use improved filtering of the image.

But these are for you to tailor.

HTH

[bassem ltaief]

thank you.

But i don't have the Optimization Toolbox

[Bjorn Gustavsson]

You're welcome.

[Bjorn Gustavsson]

That's a bummer. But you can use the very convenient fminsearchbnd (or the newer minimize) function. fminsearchbnd expects the sum-of-squared residuals instead of the residuals (as lsqnonlin does). You should be able to simply replace the resfcn with the errfcn:

errfcn = @(pars,I,x) sum((I-fitfcn(x,pars)).^2);

and lsqnonlin with fminsearchbnd (it takes the same type of par0,parmin,parMAX input as lsqnonlin).

HTH

[Rik]

Comment posted as answer by bassem ltaief:

I have this code but it does not work correctly I need help

for k=1:length(B1)
      boundary = B1{k};
      plot(boundary(:,2), boundary(:,1), 'r','LineWidth',2);
      distances = sqrt((boundary(:,2)-xc).^2 + (boundary(:,1)-yc).^2);
end
% Now find the max of the min distances.
minDistances = min(distances);
maxDistances = max(distances);

[bassem ltaief]

@Rik can you help me ti fix this

[Bjorn Gustavsson]

To get help you need to make it clear in what way it does not work the way you want.

[bassem ltaief]

I have this code which calculates the distance between a coordinate point (xc, yc) with all the border then returns the min and the max distance.

but is don't work correctly

[Bjorn Gustavsson]

Ok, you still haven't explained in what way it doesn't do what you want. One thing for you to think about is what happens with the variable distances for each value of k? Do you keep all the distances from previous iterations, or does it become ovewritten? Do you want the max and min distance for the boundary for each value of k? If so, could you do that by moving the two last lines into the loop and assign the k-th elements of minDistances and maxDistances to the max ans min values of distances?

[bassem ltaief]

i have this parte of code that show the boundaries of the picture ( j ) and after that i need to know what is minDistances and maxDistances from that point to the bordure.

[B1,L1,N1] = bwboundaries(j,'noholes');
hold on
L1=(L1>0).*L1;
for k=1:length(B1)
   boundary = B1{k};
     plot(boundary(:,2), boundary(:,1), 'r','LineWidth',2);
end

[Bjorn Gustavsson]

You can do something like this:

for k=1:length(B1)
      boundary = B1{k};
      plot(boundary(:,2), boundary(:,1), 'r','LineWidth',2);
      distances = sqrt((boundary(:,2)-xc).^2 + (boundary(:,1)-yc).^2);
      minDistances(k) = min(distances);
      maxDistances(k) = max(distances);
end

Did you manage to get my code-snippet to work? That is a far better way of going about this than any binarized idea.

Answer 2

[Image Analyst]

Here is how I did it:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures if you have the Image Processing Toolbox.
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;

%===============================================================================
% Read in gray scale image.
folder = pwd;
baseFileName = 'arc.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
  % Didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.  numberOfColorBands should be = 1.
[rows, columns, numberOfColorBands] = size(grayImage);
% If it's RGB instead of grayscale, convert it to gray scale.
if numberOfColorBands > 1
	grayImage = rgb2gray(grayImage);
end
% Display the original image.
subplot(2, 2, 1);
imshow(grayImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('Original Image : %s', baseFileName);
title(caption, 'FontSize', fontSize);
impixelinfo;
% Enlarge figure to full screen.
g = gcf;
g.WindowState = 'maximized';
g.NumberTitle = 'off';
g.Name = 'Demo by Image Analyst'
drawnow;

%===============================================================================
% Segment the image.
% Display the histogram
subplot(2, 2, 2);
imhist(grayImage);
grid on;
xlabel('Gray Level', 'FontSize', fontSize);
ylabel('Pixel Count', 'FontSize', fontSize);

% Specify a threshold.  The image is 16 bits in the range 0-65535 so the threshold will be high.
threshold = 11000;
% Place line on histogram at the mean.
xline(threshold, 'Color', 'r', 'LineWidth', 2);

% Create a binary image
mask = grayImage > threshold;
mask = imfill(mask, 'holes');
subplot(2, 2, 3);
imshow(mask);
title('Initial Mask', 'FontSize', fontSize);
drawnow;

% Now call imdilate to get them to merge together.
se = true(7, 1); % Structuring element is a vertical bar.
mask = imclose(mask, se);
% Extract the largest blob only.
mask = bwareafilt(mask, 1);
subplot(2, 2, 4);
imshow(mask);
title('Final Mask with Red Circles Overlaid', 'FontSize', fontSize);
drawnow;

%===============================================================================
% Find left and right edges.
% Scan down the mask getting the left column and right column of the blob.
leftCol = zeros(1, columns);
rightCol = zeros(1, columns);
for row = 1 : rows
	thisRow = mask(row, :);
	t = find(thisRow, 1, 'first');
	if ~isempty(t)
		leftCol(row) = t;
		rightCol(row) = find(thisRow, 1, 'last');
	end
end

%===============================================================================
% Fit the data to a circle.
x = 1 : rows;
% Get rid of lines with no white on them
noData = leftCol == 0;
x(noData) = [];
leftCol(noData) = [];
rightCol(noData) = [];
% Fit the left data to a circle.
[ycLeft, xcLeft, leftRadius, a] = circfit(x, leftCol);
viscircles([xcLeft, ycLeft], leftRadius);
% Fit the right data to a circle.
[ycRight, xcRight, rightRadius, a] = circfit(x, rightCol);
viscircles([xcRight, ycRight], rightRadius);
msgbox('Done!  Thank you Image Analyst!');

%==========================================================================================
function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT  Fits a circle in x,y plane
%
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R.  A is an optional
% output describing the circle's equation:
%
%   x^2+y^2+a(1)*x+a(2)*y+a(3)=0
         
% by Bucher izhak 25/oct/1991

n=length(x);  xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;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;
xc = -.5*a(1);
yc = -.5*a(2);
R  =  sqrt((a(1)^2+a(2)^2)/4-a(3));
end

If you're absolutely sure the radius on the right and left side should be the same (because you don't think, theoretically, it should taper down at the top and bottom), then I would just assign them both the greater of the right and left radius so that it takes on the "flatter" circle.

[Bjorn Gustavsson]

To the best of my understanding the original image is an image of some kind of interference-fringes. In that case they should be concentric circles. Diffraction-patterns from linear gratings typically have shapes of the form:

where and (k is the wavenumber, b and a are the width of the slits and their spacing in the grating respectively) and θ is the angle to the optical axis. Since the fringes are circular the integrals become trickier, and Bessel-functions typically appears, this is what typically appears when one has a Fabry-Perrot (sp?) etalon. Since the circles have reasonably large radius compared to the widths, I thought it safe to simply fit for a sinc^2 - this way one does not throw away the information that is in the image intensity, and the OP can find the most suitable function to model the intensity-variation if necessary.

%% Fitting and Error-function
fitfcn = @(x,pars) pars(1)*sinc((x-pars(2))/pars(3)).^2 + pars(4);
errfcn = @(pars,I,x)sum((I-fitfcn(x,pars)).^2);
%% Image-reading and low-pass filtering
Im = imread('image.png');
Im = double(Im);
fIm = imfilter(Im,conv2(ones(9)/9^2,ones(15)/15^2),'symmetric');
% If imfilter not available a similar enough filter can be done with conv2
[I_max,idx_max] = max(fIm,[],2);
%% Fitting row-by-row:
for i1 = 1:numel(I_max)
  parsLSQ(i1,:) = minimize(@(pars) errfcn(pars,fIm(i1,idx_max(i1)+[-200:200]),idx_max(i1)+[-200:200]),...
                           [I_max(i1),idx_max(i1),10,5e3],...
                           [],[],[],[],...
                           [0 idx_max(i1)-100 0 -1e4],...
                           [20*max(Im(:)) idx_max(i1)+100 150 2e4]);
end
%% illustrate one fit, row 800
figure
i1 = 800;
subplot(2,2,1)
plot(idx_max(i1)+[-200:200],fIm(i1,idx_max(i1)+[-200:200]),'b','linewidth',2)
hold on
plot(idx_max(i1)+[-200:200],Im(i1,idx_max(i1)+[-200:200]),'b.')
plot(idx_max(i1)+[-200:200],fitfcn(idx_max(i1)+[-200:200],parsLSQ(i1,:)),'r-','linewidth',2)
plot(parsLSQ(i1,2),parsLSQ(i1,1)+parsLSQ(i1,4),'gs','linewidth',2) 
plot(parsLSQ(i1,2)-parsLSQ(i1,3),parsLSQ(i1,4),'gs','linewidth',2) 
plot(parsLSQ(i1,2)-parsLSQ(i1,3),parsLSQ(i1,4),'gs','linewidth',2)
%% Points on image, centre, left and right first minima
subplot(2,2,2)
imagesc(Im)
hold on
plot(parsLSQ(:,2),1:size(Im,1),'c.')
plot(parsLSQ(:,2)-parsLSQ(:,3),1:size(Im,1),'m.')
plot(parsLSQ(:,2)+parsLSQ(:,3),1:size(Im,1),'r.')
subplot(2,1,2)
imagesc(Im)
hold on
%% Circle-fits
[xcM ycM RM] = circfit(parsLSQ(:,2)-parsLSQ(:,3),1:numel(idx_max));
[xcP ycP RP] = circfit(parsLSQ(:,2)+parsLSQ(:,3),1:numel(idx_max));
[xcC ycC RC] = circfit(parsLSQ(:,2),1:numel(idx_max));             
%% Plot the circle-centres and circle-sectors
plot(xcC,ycC,'c.','markersize',18)                    
axis auto
axis equal
plot(xcP,ycP,'r.','markersize',18)
plot(xcM,ycM,'m.','markersize',18)
phi_some = linspace(-10,20.5,301)*pi/180;
plot(xcP+RP*cos(phi_some),ycP+RP*sin(phi_some),'r');
plot(xcM+RM*cos(phi_some),ycM+RM*sin(phi_some),'m');
plot(xcC+RC*cos(phi_some),ycC+RC*sin(phi_some),'c');
subplot(2,2,2)
plot(xcP+RP*cos(phi_some),ycP+RP*sin(phi_some),'k--');
plot(xcM+RM*cos(phi_some),ycM+RM*sin(phi_some),'k--');
plot(xcC+RC*cos(phi_some),ycC+RC*sin(phi_some),'k--'); 

Which gives this result: