Finding the radius at regular intervals of theta for an irregular shaped object

For my research project, I need to find the radius at regular intervals of theta for an irregularly shaped object (for this particular case, a kidney bean shape). I can read in the image and convert it to black and white and then use the edge function. But I cannot figure out a way to find the radii like I need too. Please and thank you for your help.

Answer 1

[Image Analyst]

You don't need the edge function. Simply threshold your image, call bwboundaries, call regionprops to get the centroid, then run around the boundary to find the distance from the centroid to each boundary pixel, and the angle. Then if you need evenly spaced angles, and don't have any more than 1 boundary pixel at each angle, then you can use interp1 to get the angles at the values you want. Understand? If not, here's a demo where I did most everything for you, except the final interpolation part (not sure why that is needed though - it's probably not needed).

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures if you have the Image Processing Toolbox.
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

% Read in a standard MATLAB gray scale demo image.
folder = fullfile(matlabroot, '\toolbox\images\imdemos');
baseFileName = 'moon.tif';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
  % File doesn't exist -- didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, '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
grayImage = imread(fullFileName);
% Get the dimensions of the image.  
% numberOfColorBands should be = 1.
[rows columns numberOfColorBands] = size(grayImage);
% Display the original gray scale image.
subplot(2, 3, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Give a name to the title bar.
set(gcf,'name','Demo by ImageAnalyst','numbertitle','off') 

% Let's compute and display the histogram.
[pixelCount grayLevels] = imhist(grayImage);
subplot(2, 3, 2); 
bar(pixelCount);
grid on;
title('Histogram of original image', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

binaryImage = grayImage > 15;
% Clean it up a bit.
binaryImage = bwareaopen(binaryImage, 1000); % Remove small objects
binaryImage = imfill(binaryImage, 'holes'); % Fill holes.
subplot(2, 3, 3);
imshow(binaryImage, []);
title('Binary Image', 'FontSize', fontSize);

% Measure/find the centroid
measurements = regionprops(binaryImage, 'Centroid');
centroid = measurements.Centroid; % [row, column], NOT [x,y]

% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
% Plot the borders of all the objects on the original grayscale image using the coordinates returned by bwboundaries.
subplot(2, 3, 1);
title('Original image with Outlines, from bwboundaries()', 'FontSize', fontSize);
hold on;
boundaries = bwboundaries(binaryImage); % in row, column order, not x,y
numberOfBoundaries = size(boundaries, 1)
for k = 1 : numberOfBoundaries
  thisBoundary = boundaries{k};
  x = thisBoundary(:,2);
  y = thisBoundary(:,1);
  plot(x, y, 'g', 'LineWidth', 2);
end
hold off;

% Calculate the angles in degrees
deltaY = thisBoundary(:,1) - centroid(1);
deltaX = thisBoundary(:,2) - centroid(2);
angles = atand(deltaY ./ deltaX);

% Calculate the distances.
distances = sqrt((thisBoundary(:,1) - centroid(1)).^2 + (thisBoundary(:,2) - centroid(2)).^2);

% Plot distance vs. angle.
subplot(2, 3, 4:6);
plot(angles, distances, 'b-');
grid on;
title('Distance as a function of angle', 'FontSize', fontSize);
xlabel('Angle', 'FontSize', fontSize);
ylabel('Distance', 'FontSize', fontSize);

[Jacob Lundeen]

I am not a very good coder, but I think you are on the right track, but the code you provided me is not getting what I need. The code is getting the outside edge of the image I have, but that isn't what I need.

http://postimg.org/image/jdb46usup/

This is what it is giving me, but I need the kidney bean shape in the middle.

[Image Analyst]

Well up load your original image. You didn't do that yet so I had to take a guess and just used the standard demo image of the moon. It looks like you will need to change your thresholding statement so that your bean is white. But that will also make that boundary white also so you'll need to call imclearborder() to get rid of that white frame around the image.

[Jacob Lundeen]

This is the original image:

http://postimg.org/image/hqk2a4ysp/

I've tried changing the thresholding statement to make the bean white, but considering in the original image the bean is already black...

[Image Analyst]

Did you try it this way?

% Threshold the image.
binaryImage = grayImage < 50;
% Get rid of outside frame.
binaryImage = imclearborder(binaryImage);
% Get rid of small particles.
binaryImage = bwareaopen(binaryImage, 1000);

[Jacob Lundeen]

I = imread('cutout3.png');

level = graythresh(I);

BW = im2bw(I,level);

This is the code I'm using to turn my original image into the grayImage. If I use your code from here, I just get a black square.

[Jacob Lundeen]

size_BW = size(BW);

for ii = 1:size_BW(1)
    for jj = 1:size_BW(2)
        if BW(ii,jj) == 1
            BW(ii,jj) = 0;
        else
            BW(ii,jj) = 1;
        end
    end
end

BI = imclearborder(BW);
BI = bwareaopen(BI,1000);

I just ran this and I got a white kidney bean with a full black field!

[Image Analyst]

Your code is inefficient. Use my code. Is your image a uint8 gray scale image? Probably not - maybe it's color and you need to call rgb2gray() first. What does

whos I

say? By the way, don't use "I" as a variable name - it looks too much like 1 and l (lower case L).

[Jacob Lundeen]

I is a uint8 class image. im2bw turns it into a logical class. When I used your code all I got was a black square.

[Image Analyst]

Let me see your code. My code works fine with a binary (logical) image. That's what this line does:

binaryImage = grayImage > 15;

Obviously you have to adjust the threshold from 15 to 50 but this is what I showed you later when I got your original image:

binaryImage = grayImage < 50;

So everything must work. Let me see your code so I can see what you did wrong.

[Image Analyst]

Jacob, I downloaded your image, made the threshold change I told you to make, and it worked just fine. Not sure what you did, but here's what I did and it works great:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures if you have the Image Processing Toolbox.
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

% Read in a standard MATLAB gray scale demo image.
folder = 'C:\Users\Jacob\Documents\Temporary';
baseFileName = 'cutout3.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
  % File doesn't exist -- didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, '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
grayImage = imread(fullFileName);
% Get the dimensions of the image.  
% numberOfColorBands should be = 1.
[rows, columns, numberOfColorBands] = size(grayImage);
% Display the original gray scale image.
subplot(2, 3, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Give a name to the title bar.
set(gcf,'name','Demo by ImageAnalyst','numbertitle','off') 

% Let's compute and display the histogram.
[pixelCount, grayLevels] = imhist(grayImage);
subplot(2, 3, 2); 
bar(pixelCount);
grid on;
title('Histogram of original image', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

% Threshold the image.
binaryImage = grayImage < 50;
% Get rid of outside frame.
binaryImage = imclearborder(binaryImage);
% Get rid of small particles.
binaryImage = bwareaopen(binaryImage, 1000); % Remove small objects
binaryImage = imfill(binaryImage, 'holes'); % Fill holes.
subplot(2, 3, 3);
imshow(binaryImage, []);
title('Binary Image', 'FontSize', fontSize);

% Measure/find the centroid
measurements = regionprops(binaryImage, 'Centroid');
centroid = measurements.Centroid; % [row, column], NOT [x,y]

% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
% Plot the borders of all the objects on the original grayscale image using the coordinates returned by bwboundaries.
subplot(2, 3, 1);
title('Original image with Outlines, from bwboundaries()', 'FontSize', fontSize);
hold on;
boundaries = bwboundaries(binaryImage); % in row, column order, not x,y
numberOfBoundaries = size(boundaries, 1)
for k = 1 : numberOfBoundaries
  thisBoundary = boundaries{k};
  x = thisBoundary(:,2);
  y = thisBoundary(:,1);
  plot(x, y, 'g', 'LineWidth', 2);
end
hold off;

% Calculate the angles in degrees
deltaY = thisBoundary(:,1) - centroid(1);
deltaX = thisBoundary(:,2) - centroid(2);
angles = atand(deltaY ./ deltaX);

% Calculate the distances.
distances = sqrt((thisBoundary(:,1) - centroid(1)).^2 + (thisBoundary(:,2) - centroid(2)).^2);

% Plot distance vs. angle.
subplot(2, 3, 4:6);
plot(angles, distances, 'b-');
grid on;
title('Distance as a function of angle', 'FontSize', fontSize);
xlabel('Angle', 'FontSize', fontSize);
ylabel('Distance', 'FontSize', fontSize);

[Jacob Lundeen]

This is the code I'm using to convert to grayscale:

level = graythresh(grayImage);
grayImage = im2bw(grayImage,level);

When I do that using your code, the "binaryImage" comes out to a giant black square. If I don't put that code in, the program breaks at the imhist line because the image file isn't the correct file type.

[Image Analyst]

Don't do that. If you have an RGB image just use

if numberOfColorBands > 1
  grayImage = grayImage(:,:,2); % Take the green channel.
end

Strange that I didn't have to do that with the image you uploaded.

[Jacob Lundeen]

Ok, that worked. Thank you so much for your help so far. Two questions: when looking at the original image, you see how the kidney bean in inside the metal "cutout", the center needs to be from the edges of that cutout (this is for research with radiation treatment). From how the code reads, when it is turned into the binary image, the Centroid is the center of the kidney bean, correct? Also, the angles array goes from -89 to +89, and I need it to go 0 to 359, is there a quick way to convert it or should I just setup a loop?

[Image Analyst]

Correct. I'm sure you can figure out the angle problem.

[Jacob Lundeen]

Ok, what would be the easiest way to make the code use the center point of the entire cutout and not the center of the kidney bean?

[Image Analyst]

threshold to find the cutout, then call regionprops. Same concept, you're just getting a different foreground object. Do you understand that regionprops makes measurements on the "white" regions in the binary image? So do whatever you have to do to get the region(s) you want to be white/1/true. Then call regionprops().

[Jacob Lundeen]

I'm figuring it out, slowly. We don't get taught a majority of these image manipulations and such, so I'm learning on the go.