How can I fit both left and right curve with a circle and get the equation of it?

1 view (last 30 days)
How can I find a circle fitting for both left and right curve and get the equation of the fitting circle? I need to know the tangent of the point that circle and horizontal line intersects. Thanks!

Accepted Answer

Will Nitsch
Will Nitsch on 31 Aug 2017
The Iin in the script below is the image you gave. See my comments for information on what I did. I basically just solved for lines that were equidistant to the top, bottom and middle points on each curve. Pretty fun problem to solve!Not 100% correct (see the image below, slightly off in the y direction). The script below solves for the center and radius for each circle relative to the location (1,1) in the original image, so the left center point location is negative.
Iin = imbinarize(rgb2gray(im2double(imread('Screen Shot 201708-28 at 3.26.41 pm.png'))),.01);
% create a solid region of the image
I2 = imdilate(Iin,strel('disk',10));
I2 = imerode(I2,strel('disk',10));
% pad so we have some area in the image to put the circle centers
I2 = padarray(I2,[0,600]);
figure
imshow(I2)
%%get just left and right edges of the region
Ileft = zeros(size(I2));
Iright = zeros(size(I2));
for m = 1:size(I2,1)
idx = find(I2(m,:)==1);
if(length(idx)>0)
Ileft(m,idx(1)) = 1;
Iright(m,idx(end)) = 1;
end
end
% clip off the curves at the top or bottom of each side
Ileft(1:100,:) = 0;
Iright(1:100,:) = 0;
Ileft(530:end,:) = 0;
Iright(530:end,:) = 0;
%%pad left and right images for hough transform so center of circle should be in the image
%Ileft = padarray([zeros(size(Ileft)),Ileft],[500,500]);
%Iright = padarray([Iright,zeros(size(Iright))],[500,500]);
%%get top, bottom and mid points for each line
[yL, xL] = find(Ileft == 1);
[yR, xR] = find(Iright == 1);
topLidx = find(yL == min(yL));
if(length(topLidx)>1) % in case there are more than 1 index
topLidx = topLidx(1);
end
topRidx = find(yR == min(yR));
if(length(topRidx)>1) % in case there are more than 1 index
topRidx = topRidx(1);
end
botLidx = find(yL == max(yL));
if(length(botLidx)>1) % in case there are more than 1 index
botLidx = botLidx(1);
end
botRidx = find(yR == max(yR));
if(length(botRidx)>1) % in case there are more than 1 index
botRidx = botRidx(1);
end
midlengthL = floor((yL(botLidx)-yL(topLidx))/2+yL(topLidx));
midlengthR = floor((yR(botRidx)-yR(topRidx))/2+yR(topRidx));
midLidx = find(yL == midlengthL); % rightmost point
midRidx = find(yR == midlengthR); % leftmost point
topLxy = [xL(topLidx),yL(topLidx)]; % top left point
topRxy = [xR(topRidx),yR(topRidx)]; % top right point
botLxy = [xL(botLidx),yL(botLidx)]; % bottom left point
botRxy = [xR(botRidx),yR(botRidx)]; % bottom right point
midLxy = [xL(midLidx),yL(midLidx)]; % mid left point
midRxy = [xR(midRidx),yR(midRidx)]; % mid right point
%%bisect line between top and bottom with line from middle (center point will be on this line)
distL = sqrt((botLxy(1)-topLxy(1))^2+(botLxy(2)-topLxy(2))^2); % distance between top and bottom on left
unitVecBTL = [(botLxy(1)-topLxy(1)) (botLxy(2)-topLxy(2))]/distL; % unit vector between top and bottom on left
bisectLxy = topLxy+unitVecBTL*distL/2; % midpoint of line between top and bottom points on left
unitVecBisectL = (bisectLxy - midLxy)/sqrt((bisectLxy(1)-midLxy(1))^2+(bisectLxy(2)-midLxy(2))^2); % unit vector from middle of left curve to midpoint top to bottom line (bisector)
distR = sqrt((botRxy(1)-topRxy(1))^2+(botRxy(2)-topRxy(2))^2); % distance between top and bottom on right
unitVecBTR = [(botRxy(1)-topRxy(1)) (botRxy(2)-topRxy(2))]/distR; % unit vector between top and bottom on right
bisectRxy = topRxy+unitVecBTR*distR/2; % midpoint of line between top and bottom points on right
unitVecBisectR = (bisectRxy - midRxy)/sqrt((bisectRxy(1)-midRxy(1))^2+(bisectRxy(2)-midRxy(2))^2); % unit vector from middle of right curve to midpoint top to bottom line (bisector)
%%solve left circle
radius = 1;
initFlag = 1;
errorMinFlag = 0;
% increment radius along bisecting line from midpoint until this point is equidistant
% to the top, middle and bottom points on curve.
while(~errorMinFlag)
midBiL = midLxy + unitVecBisectL*radius;
distTopBisectL = sqrt((midBiL(1)-topLxy(1))^2+(midBiL(2)-topLxy(2))^2);
distBotBisectL = sqrt((midBiL(1)-botLxy(1))^2+(midBiL(2)-botLxy(2))^2);
error = abs(distTopBisectL + distBotBisectL - 2*radius);
if(initFlag == 1)
initFlag = 0;
else
if(error > oldError)
radius = radius - 1;
errorMinFlag = 1;
end
end
radius = radius + 1;
oldError = error;
end
centerLxy = midLxy + unitVecBisectL*radius;
radiusL = radius;
%%solve right circle
clear oldError
clear error
radius = 1;
initFlag = 1;
errorMinFlag = 0;
% increment radius along bisecting line from midpoint until this point is equidistant
% to the top, middle and bottom points on curve.
while(~errorMinFlag)
midBiR = midRxy + unitVecBisectR*radius;
distTopBisectR = sqrt((midBiR(1)-topRxy(1))^2+(midBiR(2)-topRxy(2))^2);
distBotBisectR = sqrt((midBiR(1)-botRxy(1))^2+(midBiR(2)-botRxy(2))^2);
error = abs(distTopBisectR + distBotBisectR - 2*radius);
if(initFlag == 1)
initFlag = 0;
else
if(error > oldError)
radius = radius - 1;
errorMinFlag = 1;
end
end
radius = radius + 1;
oldError = error;
end
centerRxy = midRxy + unitVecBisectR*radius;
radiusR = radius;
%%show circles
figure
imshow(Iin)
% correct center points x values for padding
centerLxy(1) = centerLxy(1)-600;
centerRxy(1) = centerRxy(1)-600;
% (now center points are relative to pixel 1,1 in input image)
viscircles([centerLxy(1),centerLxy(2)],radiusL)
viscircles([centerRxy(1),centerRxy(2)],radiusR)
  2 Comments
Will Nitsch
Will Nitsch on 1 Sep 2017
As a side note, I tried using imfindcircles using just the "Ileft" and "Iright" images (dilated versions too) but the circular hough transform wasn't really picking up the circles.
Also, if you want to detect the tangent at the circle edge points that intersect the horizontal lines, you can just figure out the angle of the circle from the center point of the circle and the top and bottom pixels, for the angle to the top and bottom horizontal lines respectively. Seems like this should do the trick on the angle calculations too:
circle_tangent_left = atand(-(topLxy(2)-centerLxy(2))/((topLxy(1)-600)-(centerLxy(1)+600)))+90;
circle_tangent_right = 180-atand((topRxy(2)-centerRxy(2))/((topRxy(1)-600)-(centerRxy(1)+600)))+90;

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!