2

I modified a code based on the shown in https://www.mathworks.com/matlabcentral/answers/377838-please-how-can-i-find-the-center-and-the-radius-of-the-inscribed-circle-of-a-set-of-points in order to find the inscribe circle but I do not understand why the image is rotated. Why and how can I solve it?

Code:

url='https://i.pcmag.com/imagery/reviews/00uaCVfzQ4Gsuhmh85WvT3x-4.fit_scale.size_1028x578.v_1569481686.jpg';

Image = rgb2gray(imread(url));   
Image = imcomplement(Image);
fontSize = 10; 
% determine contours
BW = imbinarize(Image);
BW = imfill(BW,'holes');
[B,L] = bwboundaries(BW,'noholes');
k = 1;
b = B{k};
y  = b(:,2);
x  = b(:,1);

 subplot(2, 2, 1);
plot(x, y, 'b.-', 'MarkerSize', 3);
grid on;
title('Original Points', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% Make data into a 1000x1000 image.
xMin = min(x)
xMax = max(x)
yMin = min(y)
yMax = max(y)
scalingFactor = 1000 / min([xMax-xMin, yMax-yMin])
x2 = (x - xMin) * scalingFactor + 1;
y2 = (y - yMin) * scalingFactor + 1;
mask = poly2mask(x2, y2, ceil(max(y2)), ceil(max(x2)));
% Display the image.
p2 = subplot(2, 2, 2);
imshow(mask);
axis(p2, 'on', 'xy');
title('Mask Image', 'FontSize', fontSize);
% Compute the Euclidean Distance Transform
edtImage = bwdist(~mask);
% Display the image.
p3 = subplot(2, 2, 3);
imshow(edtImage, []);
axis(p3, 'on', 'xy');
% Find the max
radius = max(edtImage(:))
% Find the center
[yCenter, xCenter] = find(edtImage == radius)
% Display circles over edt image.
viscircles(p3, [xCenter, yCenter], radius,'Color','g');
% Display polygon over image also.
hold on;
plot(x2, y2, 'r.-', 'MarkerSize', 3, 'LineWidth', 2);
title('Euclidean Distance Transform with Circle on it', 'FontSize', fontSize);
% Display the plot again.
subplot(2, 2, 4);
plot(x, y, 'b.-', 'MarkerSize', 3);
grid on;
% Show the circle on it.
hold on;
% Scale and shift the center back to the original coordinates.
xCenter = (xCenter - 1)/ scalingFactor + xMin
yCenter = (yCenter - 1)/ scalingFactor + yMin
radius = radius / scalingFactor
rectangle('Position',[xCenter-radius, yCenter-radius, 2*radius, 2*radius],'Curvature',[1,1]);
title('Original Points with Inscribed Circle', 'FontSize', fontSize);

Original image: enter image description here

Output image enter image description here

alirazi
  • 159
  • 8

1 Answers1

1

[B,L] = bwboundaries(BW,...) returns in B the row and column values (documentation). That is, the first column of B{k} is y, the second one is x.

After changing this bit of code as follows:

y = b(:,1);
x = b(:,2);

you will notice that the image is upside down! That is because in an image the y-axis increases down (y is the row number in the matrix), whereas in the plot the y-axis increases up (mathematical coordinate system).

The axes where you use imshow in are automatically set to the right coordinate system, but then you do axis(p3, 'on', 'xy');, turning it upside down again. Instead, use

axis(p1, 'on', 'image');

on the axes where you don't use imshow (i.e. the top-left and bottom-right ones).

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
  • Thanks a lot. BTW can I assume that 'MajorAxisLength'/2 is the circumscribed circle radius? – alirazi Apr 15 '20 at 14:49
  • @alirazi: In your code, there is a `radius` variable that contains the radius in pixels. `regionprops` uses `MajorAxisLength` as the larger axis of the ellipse that has the same moments of inertia as the object (it's the best fit ellipse, in some sense of "best"). It does not typically correspond to the diameter of the circumscribed circle. The `MaxFeretDiameter` could correspond to the circumscribed circle, but doesn't necessarily, there are cases where this would fail too. – Cris Luengo Apr 15 '20 at 15:07
  • @alirazi: Circumscribing a circle would involve finding three extreme points on the contour. Three points always define a circle. You'd have to iterate over all possible pairs of points until you find a circle that circumscribes all contour points. – Cris Luengo Apr 15 '20 at 15:08
  • @CrisLuengo thanks for the explanation. A solution can be found in "A suite of minimal bounding objects/minboundcircle" function (https://www.mathworks.com/matlabcentral/fileexchange/34767-a-suite-of-minimal-bounding-objects) – wrek Apr 15 '20 at 16:08
  • @wrek: Thanks for the link. John D'Errico is a trustworthy author, I'm sure those tools are very good. – Cris Luengo Apr 15 '20 at 16:11