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)
