How to fit circles in concave parts of the polygon or binary Image ?

I have written the following code, but now I'm stuck. I need to find radius of curvature, which fits best at the concave parts of the figure. Please help.

if true
  % cod
% Reading the Image
clc;
close all
img = imread('black.jpg');
figure(1);
imshow(img);
title('

Original Image');

% Gray Image
imgGray = rgb2gray(img);
figure(2);
imshow(imgGray);
title('Gray Image');

% Binary Image
imgbw = imbinarize(imgGray,0.36);
figure(3);
imshow(imgbw);
imgbw= imcomplement(imgbw);
BW2 = imfill(imgbw,'holes');
title('Binary Image');

%boundary of original
boundaries1 = bwboundaries(BW2);
thisBoundary1= boundaries1{1}; 
      xo = thisBoundary1(:,2);
      yo = thisBoundary1(:,1);

% smoothening the image
windowSize = 51;
kernel = ones(windowSize) / windowSize ^ 2;
blurryImage = conv2(single(BW2), kernel, 'same');
binaryImage = blurryImage > 0.5; % Rethreshold
figure (4)
imshow(binaryImage)
title('Smoothened Image');

% boundary of smoothened image
boundaries = bwboundaries(binaryImage);
thisBoundary = boundaries{1};
      xb = thisBoundary(:,2);
      yb = thisBoundary(:,1);

%finding centroid of smoothened Image
measurements = regionprops(binaryImage , 'Centroid');
centroidx = measurements(1).Centroid(1);
centroidy = measurements(1).Centroid(2);

% polygon approximation
p_reduced = reducepoly(thisBoundary,0.001);
xp=p_reduced(:,2);
yp=p_reduced(:,1);
  end

I want something like this:

[Narges Kheradmandi]

Hey,

Did you find the answer to your question? I have the same problem.

[Image Analyst]

Yes, it's below. If that's not what you want then start your own question in a new thread (not here) and explain why my code is not good for you.

If you have any more questions, then attach your data and code to read it in with the paperclip icon after you read this:

TUTORIAL: How to ask a question (on Answers) and get a fast answer

[Narges Kheradmandi]

Many thanks for your reply. However, I think your code gives us the fitted circle for each point on the boundary of the abject but we only need the circles which fitted to the curved part of the object(as can be seen in the image that I have uploaded. I need the radius of the green circles in the picture. Do you know how can I extract them?

[Image Analyst]

What is the link to your new question? Have you posted it yet? You can reference this discussion in it, but let's not bug @Deepank Singh with emails about activity on his question anymore. Let's start your own discussion thread.

Answer 1

[Image Analyst]

I'd get the boundary points for the shape, then use the FAQ:

https://matlab.fandom.com/wiki/FAQ#How_can_I_fit_a_circle_to_a_set_of_XY_data.3F

to fit a circle to a sliding window of points. So you get a radius of curvature for every single boundary point for that one window size. Then I'd do it for every window size in a reasonable range. For each point for each window size, I'd compute the residuals (distance from point to fitted circle). So now you know what window size gave the best fit. And you also know the radius at that location for that window size. Conceptually pretty simple. Give it a try, since I don't have any demo code for that all ready to hand to you. Probably within an hour you could get it done.

[Deepank Singh]

Thanks for the reply, I only need to find the radius of curvature only at concave inside points, and I couldn't find the points of concavity, where I need to fit the circle. I can fit the circle using circfit function, but first, I need to find the concave points, that I'm unable to find.

I'm new to MATLAB, Please Don't mind if I sounded silly.

[Image Analyst]

You get the (x,y) and radius from the fit. Use inpolygon() to determine if the center is inside or outside your shape. That will tell you if the curve is convex or concave, respectively, at those points.

[Deepank Singh]

Thanks will update you, if it works. Thanks

[Deepank Singh]

Hey, I tried what you said, but still don't know how to get the best circle,

if true
  % concave portion detection
x=0;
for i=1:length(D(:,1))-4
  %  lineseg1 = [ D(i,1) D(i,2);D(i+1,1) D(i+1,2) ]; 
  %  lineseg2 = [ D(i,1) D(i,2);D(i+2,1) D(i+2,2) ];
  %  [in,out1] = intersect(poly1,lineseg1);
  %  [in,out2] = intersect(poly1,lineseg2);

     % if(isempty(out1))
        % if(isempty(out2)) 
             [xr, yr, r, a]  = circfit(D(i:i+2,1),D(i:i+2,2));
th = 0:pi/50:2*pi;
xunit = r * cos(th) + xr;
yunit = r * sin(th) + yr;
poly0 = polyshape(xunit(1:end-1),yunit(1:end-1));
poly2 = intersect(poly0,poly1);
A1 = area(poly0);
A2 =area(poly2);
overlap = A2/A1;
             if(r<1.0*Inradius && inpolygon(xr,yr,xp,yp) && overlap>0.98) 
                      x=x+1;
                 roc(x,:)= [xr yr r];
             end
             [xr, yr, r, a]  = circfit(D(i:i+3,1),D(i:i+3,2));
             th = 0:pi/50:2*pi;
xunit = r * cos(th) + xr;
yunit = r * sin(th) + yr;
poly0 = polyshape(xunit(1:end-1),yunit(1:end-1));
poly2 = intersect(poly0,poly1);
A1 = area(poly0);
A2 =area(poly2);
overlap = A2/A1;
             if(r<1.0*Inradius && inpolygon(xr,yr,xp,yp) && overlap>0.98)                 
                 x=x+1;
                 roc(x,:)= [xr yr r];
             end
             [xr, yr, r, a]  = circfit(D(i:i+4,1),D(i:i+4,2));
             th = 0:pi/50:2*pi;
xunit = r * cos(th) + xr;
yunit = r * sin(th) + yr;
poly0 = polyshape(xunit(1:end-1),yunit(1:end-1));
poly2 = intersect(poly0,poly1);
A1 = area(poly0);
A2 =area(poly2);
overlap = A2/A1;
             if(r<1.0*Inradius && inpolygon(xr,yr,xp,yp)  && overlap>0.98)                 
                 x=x+1;
                 roc(x,:)= [xr yr r];
             end
       %  end              
    %  end    
end     
  end

This is what I got

[Image Analyst]

The code gives radii all along the border. What is your definition of the "best" circle radii?