MATLAB: Please , how can I find the center and the radius of the inscribed circle of a set of points

Automated Driving ToolboxMATLAB

Please , how can I find the center and the radius of the inscribed circle of a set of points

Best Answer

daeze, see if this is what you want. I solved it numerically by converting the data to a digital image, then using the Euclidean distance transform and finding the center of it. The max of the EDT is the center of the circle and its value there is the radius. If you need more than 3 digits of precision then you can increase the size of the image from 1000 to 10000 or larger.
% Initialization / clean-up code.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 25;
x=[2.6381318;-18.66985584;-39.97784349;-69.8090262;-99.64020891;-92.131167993;-84.62315095;-57.4301004;-30.23704914;-3.65279788;22.93145337;34.09278024;45.2541071;23.94611945]
y=[75.28822303;63.72102973;52.15383644;43.02184173;33.88984702;19.48158871;5.07333039;5.07333039;5.07333039;5.07333039;5.07333039;25.77251893;46.4717064;60.87996471]
subplot(2, 2, 1);
plot(x, y, 'b.-', 'MarkerSize', 30);
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);
% Display polygon over image also.
hold on;
plot(x2, y2, 'r.-', 'MarkerSize', 30, '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', 30);
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);
Related Question