Find minimum enclosing circle

I am using the function minboundcircle from:

https://www.mathworks.com/matlabcentral/fileexchange/34767-a-suite-of-minimal-bounding-objects?s_tid=srchtitle

However, it is not finding the absolute edge which is introducing some error. Is there a way I can mitigate the error and find the absolute edge.

Example:

[Benjamin Thompson]

Since this is code from a third party author on the File Exchange site, you should ask the author directly. You can post your question to the link on File Exchange above.

[John D'Errico]

I cannot know how you used the code, or what you tried to do. So I cannot guess what you did wrong either.

Answer 1

[Walter Roberson]

Use regionprops() to find the Centroid and PixelList . The circle center goes at the centroid. The maximum distance between the Centroid and the PixelList coordinates defines the radius of the circle.

Answer 2

[John D'Errico]

So let me use my own code to solve your problem.

im = imread('IMG_3108.jpg');

[redpixr,redpixc] = find((im(:,:,1) > 140) & (im(:,:,2) < 65) & (im(:,:,3) < 65));

[C,R] = minboundcircle(redpixr,redpixc);
t = linspace(0,2*pi);
X = R*cos(t) + C(1);
Y = R*sin(t) + C(2);
plot(redpixr,redpixc,'b.',X,Y,'r-')
axis equal

Which works fine. Don't forget that the axes there are screwed up as I plotted it, because I used plot, and I plotted pixel coordinates in terms of rows and columns of the original array.

This is probably where your attempt failed. Of course, I'm making only a wild-ass guess here. So you will need to carefully transform the center of that circle back into the correct domain for that image.

[jordan Ashton]

Hello, can you check my recent comment I included my code? Thanks

[jordan Ashton]

Hi John,

I think I may have been unclear. I am using point cloud data. I am uploading a .ply file which contains xyz points. The image above is a picture of my screen.

[jordan Ashton]

Code:

cloud = pcread('data.ply');             
data=cloud.Location;                      
xyz = double(cloud.Location); 
[rpc] = rotatePointCloudAlongZ(xyz, 'x');   
x = rpc(:,1);                             
y = rpc(:,2);
[center,radius] = minboundcircle(x,y,true)
theta = linspace(0,2*pi,100);  
xc = center(1) + radius*cos(theta);
yc = center(2) + radius*sin(theta);
plot(x,y,'ro',xc,yc,'b-')        

% rotatePointCloudAlongZ function

function pc2 = rotatePointCloudAlongZ(pc, direction)
% This function rotates the point cloud of the highest eigenvector  
% to be along z-direction. 
% For the 2nd highest eigenvector, the user can choose to align 
% along x or y direction 

% parameters: 
% inputs: 
% pc: point cloud of size N * 3 size
% direction: 'x' or 'y' (align the 2nd highest eigenvector along 'x' or
% 'y' direction 
% output:
% pc2: rotated point cloud 
%% Bring the point cloud center to the origin

pc = pc - mean(pc);

%% This section align u normal vector along z-direction

% Obtain the eigenvector of the highest eigenvalue
u = pcaEig(pc, 'max');

% Calculate the angles of the normal vector 
[alpha, beta] = unitVectorToAngle(u);

% Align the point cloud along x-axis followed by aligning along z-axis
% YOU CAN REMOVE THE PI IF YOU WANT TO FLIP THE POINT CLOUD ALONG
% Z-DIRECTION 
[~, Ry, Rz] = rotationalMatrix(-alpha, pi-beta);
pc2 = rotatePC(pc, Ry, Rz);

%% This section align v normal vector along x or y direction

switch direction
    case 'x'
        offset = 0;
    case 'y'
        offset = pi/2;
end

% Obtain the eigenvector of the 2nd highest eigenvalue 
v = pcaEig(pc2, 'middle');

% Calculate the angle of the projected v-vector along the xy-plane
% with respect to the x-axis
[alpha, ~] = unitVectorToAngle(v);

% Calculate the rotational matrix for the angle
[~, Ry, Rz] = rotationalMatrix(offset - alpha, 0);

% Rotate the point cloud 
pc2 = rotatePC(pc2, Ry, Rz);

end
function [Rx, Ry, Rz] = rotationalMatrix(alpha, beta)

Rx = [1 0 0; 0 cos(beta) -sin(beta); 0 sin(beta) cos(beta)];
Ry = [cos(beta) 0 sin(beta); 0 1 0; -sin(beta) 0 cos(beta)];
Rz = [cos(alpha) -sin(alpha) 0; sin(alpha) cos(alpha) 0; 0 0 1];

end
function u = pcaEig(pc, magnitude)
%% Obtain the covariance matrix

covariance = cov([pc(:, 1) pc(:, 2) pc(:, 3)]);

%% Calculate the eigenvectors and obtain the normal vector

[V, D] = eig(covariance);
diagonalEigenvalues = diag(D);

%% Output the normal vector 

% Sort the eigenvectors based on size of eigenvalues 
[~, I] = sort(diagonalEigenvalues);
V = V(:, I);

switch magnitude
    case 'max'
        % Choose the eigenvector of the highest eigenvalue
        u = V(:, 3); 
    case 'middle'
        % Choose the eigenvector of the middle eigenvalue
        u = V(:, 2); 
    case 'min'
        % Choose the eigenvector of the lowest eigenvalue
        u = V(:, 1); 
end

end
function [alpha, beta] = unitVectorToAngle(u)
% Rotational angle between the projected u on the xy plane and the x-axis
alpha = atan2(u(2), u(1)); 
% Rotational angle between the u vector and the z-axis
beta = atan2(sqrt(u(1)^2 + u(2)^2), u(3));

end
function pc2 = rotatePC(pc, Ry, Rz)
% Convert the point cloud to 3 * N format
matrix = pc';
% rotation around z axis to align point cloud along x axis 
matrix2 = Rz*matrix;
matrix2 = Ry*matrix2;
% Ouput the point cloud 
pc2 = matrix2';
end

Answer 3

[Image Analyst]

Here's my code. Looks like it works as expected for me and John. What did you do differently? Did you not use viscircles()???

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;
%===============================================================================
% Get the name of the image the user wants to use.
baseFileName = 'IMG_3108.jpg';
folder = pwd;
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
    % The file doesn't exist -- didn't find it there in that folder.
    % Check the entire search path (other folders) for the file by stripping off the folder.
    fullFileNameOnSearchPath = baseFileName; % No path this time.
    if ~exist(fullFileNameOnSearchPath, 'file')
        % Still didn't find it.  Alert user.
        errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
        uiwait(warndlg(errorMessage));
        return;
    end
end
%=======================================================================================
% Read in demo image.
rgbImage = imread(fullFileName);
% Get the dimensions of the image.
[rows, columns, numberOfColorChannels] = size(rgbImage)
% Display image.
subplot(2, 2, 1);
imshow(rgbImage, []);
impixelinfo;
axis on;
caption = sprintf('Original Color Image\n%s', baseFileName);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0.05 1 0.95]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
drawnow;

% Create the mask for the red regions.
[mask, maskedRGBImage1] = createMask(rgbImage);
% Display the color segmentation mask image.
subplot(2, 2, 2);
imshow(mask, []);
title('Color Segmentation Mask Alone', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
axis('on', 'image');
drawnow;
% Now do clean up by hole filling, and getting rid of small blobs.
mask = imfill(mask, 'holes');
mask = bwareafilt(mask, 1); % Take largest blob only.
% Display the color segmentation mask image.
subplot(2, 2, 3);
imshow(mask, []);
title('Final Mask', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
axis('on', 'image');
drawnow;

% Get x and y
boundary = bwboundaries(mask);
x = boundary{1}(:, 2);
y = boundary{1}(:, 1);
subplot(2, 2, 1);

% Get the minimum bounding circle.
hullflag = true;
[center,radius] = minboundcircle(x,y,hullflag)

% Display image.
subplot(2, 2, 4);
imshow(rgbImage, []);
impixelinfo;
axis on;
caption = sprintf('With Boundary and Bounding Circle');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
% Display boundary
hold on;
plot(x, y, 'b-', 'LineWidth', 4);
% Show circle
viscircles(center, radius, 'color', 'r');
uiwait(helpdlg('Done!'));

function [BW,maskedRGBImage] = createMask(RGB)
%createMask  Threshold RGB image using auto-generated code from colorThresholder app.
%  [BW,MASKEDRGBIMAGE] = createMask(RGB) thresholds image RGB using
%  auto-generated code from the colorThresholder app. The colorspace and
%  range for each channel of the colorspace were set within the app. The
%  segmentation mask is returned in BW, and a composite of the mask and
%  original RGB images is returned in maskedRGBImage.
% Auto-generated by colorThresholder app on 11-Apr-2022
%------------------------------------------------------
% Convert RGB image to chosen color space
I = rgb2hsv(RGB);
% Define thresholds for channel 1 based on histogram settings
channel1Min = 0.888;
channel1Max = 0.202;
% Define thresholds for channel 2 based on histogram settings
channel2Min = 0.333;
channel2Max = 1.000;
% Define thresholds for channel 3 based on histogram settings
channel3Min = 0.000;
channel3Max = 1.000;
% Create mask based on chosen histogram thresholds
sliderBW = ( (I(:,:,1) >= channel1Min) | (I(:,:,1) <= channel1Max) ) & ...
    (I(:,:,2) >= channel2Min ) & (I(:,:,2) <= channel2Max) & ...
    (I(:,:,3) >= channel3Min ) & (I(:,:,3) <= channel3Max);
BW = sliderBW;
% Initialize output masked image based on input image.
maskedRGBImage = RGB;
% Set background pixels where BW is false to zero.
maskedRGBImage(repmat(~BW,[1 1 3])) = 0;
end

[jordan Ashton]

Code:

cloud = pcread('data.ply');             
data=cloud.Location;                      
xyz = double(cloud.Location); 
[rpc] = rotatePointCloudAlongZ(xyz, 'x');   
x = rpc(:,1);                             
y = rpc(:,2);

[center,radius] = minboundcircle(x,y,true)
theta = linspace(0,2*pi,100);  
xc = center(1) + radius*cos(theta);
yc = center(2) + radius*sin(theta);

plot(x,y,'ro',xc,yc,'b-')          








function pc2 = rotatePointCloudAlongZ(pc, direction)
    % This function rotates the point cloud of the highest eigenvector  
    % to be along z-direction. 
    % For the 2nd highest eigenvector, the user can choose to align 
    % along x or y direction 
    
    % parameters: 
    % inputs: 
    % pc: point cloud of size N * 3 size
    % direction: 'x' or 'y' (align the 2nd highest eigenvector along 'x' or
    % 'y' direction 
    % output:
    % pc2: rotated point cloud 
    %% Bring the point cloud center to the origin
    
    pc = pc - mean(pc);
    
    %% This section align u normal vector along z-direction
    
    % Obtain the eigenvector of the highest eigenvalue
    u = pcaEig(pc, 'max');
       
    % Calculate the angles of the normal vector 
    [alpha, beta] = unitVectorToAngle(u);
    
    % Align the point cloud along x-axis followed by aligning along z-axis
    % YOU CAN REMOVE THE PI IF YOU WANT TO FLIP THE POINT CLOUD ALONG
    % Z-DIRECTION 
    [~, Ry, Rz] = rotationalMatrix(-alpha, pi-beta);
    pc2 = rotatePC(pc, Ry, Rz);
    
    %% This section align v normal vector along x or y direction
    
    switch direction
        case 'x'
            offset = 0;
        case 'y'
            offset = pi/2;
    end
    
    % Obtain the eigenvector of the 2nd highest eigenvalue 
    v = pcaEig(pc2, 'middle');
    
    % Calculate the angle of the projected v-vector along the xy-plane
    % with respect to the x-axis
    [alpha, ~] = unitVectorToAngle(v);
    
    % Calculate the rotational matrix for the angle
    [~, Ry, Rz] = rotationalMatrix(offset - alpha, 0);
    
    % Rotate the point cloud 
    pc2 = rotatePC(pc2, Ry, Rz);
    
end
function [Rx, Ry, Rz] = rotationalMatrix(alpha, beta)
    
    Rx = [1 0 0; 0 cos(beta) -sin(beta); 0 sin(beta) cos(beta)];
    Ry = [cos(beta) 0 sin(beta); 0 1 0; -sin(beta) 0 cos(beta)];
    Rz = [cos(alpha) -sin(alpha) 0; sin(alpha) cos(alpha) 0; 0 0 1];
    
end
function u = pcaEig(pc, magnitude)
    %% Obtain the covariance matrix
    
    covariance = cov([pc(:, 1) pc(:, 2) pc(:, 3)]);
    
    %% Calculate the eigenvectors and obtain the normal vector
    
    [V, D] = eig(covariance);
    diagonalEigenvalues = diag(D);
    
    %% Output the normal vector 
    
    % Sort the eigenvectors based on size of eigenvalues 
    [~, I] = sort(diagonalEigenvalues);
    V = V(:, I);
    
    switch magnitude
        case 'max'
            % Choose the eigenvector of the highest eigenvalue
            u = V(:, 3); 
        case 'middle'
            % Choose the eigenvector of the middle eigenvalue
            u = V(:, 2); 
        case 'min'
            % Choose the eigenvector of the lowest eigenvalue
            u = V(:, 1); 
    end
    
end
function [alpha, beta] = unitVectorToAngle(u)
    % Rotational angle between the projected u on the xy plane and the x-axis
    alpha = atan2(u(2), u(1)); 
    % Rotational angle between the u vector and the z-axis
    beta = atan2(sqrt(u(1)^2 + u(2)^2), u(3));
    
end
function pc2 = rotatePC(pc, Ry, Rz)
    % Convert the point cloud to 3 * N format
    matrix = pc';
    % rotation around z axis to align point cloud along x axis 
    matrix2 = Rz*matrix;
    matrix2 = Ry*matrix2;
    % Ouput the point cloud 
    pc2 = matrix2';
end

[jordan Ashton]

Hi all, I have attached code but still confused on what I need to do to solve the issue. Data is included in zip file. Thank you

[Image Analyst]

Unrecognized function or variable 'rotatePointCloudAlongZ'.

[jordan Ashton]

Link to 'rotatePointCloudAlongZ':

https://www.mathworks.com/matlabcentral/fileexchange/79472-align-rotate-point-cloud-along-z-direction-based-on-pca

[Image Analyst]

Not planning on dowloading some File Exchange submission to help you. What about my original solution on the red spot? Didn't that contain the spot like it should?

[jordan Ashton]

Thank you for your patience. I have updated code so it should run fine. Yes your example does have te result I am looking for I am just a bit confused becuase I am not working with image data.

[Image Analyst]

Well now we're all confused because you uploaded an image and now you say you're not working with images. Well, anyway, I think we all showed that John's function works.

[jordan Ashton]

I uploaded an image of my screen it is clear from my code that I am using a .ply file with xyz data.