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

3 views (last 30 days)
bassem ltaief on 18 Nov 2020
I have a noisy picture and i need to get two circles like this picture

Show 8 older comments
bassem ltaief on 18 Nov 2020
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 on 18 Nov 2020
OK "it is not like you think". Could you, pretty please, upload the original image, then I'll get you a solutino.
bassem ltaief on 19 Nov 2020

Bjorn Gustavsson on 18 Nov 2020
Edited: Bjorn Gustavsson on 18 Nov 2020
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

Show 8 older comments
Bjorn Gustavsson on 26 Nov 2020
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 on 26 Nov 2020
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 on 26 Nov 2020
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.

Image Analyst on 3 Dec 2020
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
% 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;
subplot(2, 2, 3);
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.
% Extract the largest blob only.
subplot(2, 2, 4);
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);
% Fit the right data to a circle.
[ycRight, xcRight, rightRadius, a] = circfit(x, rightCol);
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.

#### 1 Comment

Bjorn Gustavsson on 3 Dec 2020
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 = 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:

R2019a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!