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

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!

Answer 1

[Will Nitsch]

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)

[Oshikuru]

Thanks! that helps a lot!

[Will Nitsch]

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;