Expand to find inner boundary of object

I have a stack of 3D data that has been split up into slices. Each slice looks something like what is shown here in 'Scatter.'

I need to find the inner and outer boundary of the points. 'boundary()' works really well to find the outer boundary. What I initially did was found all the data points on the outer boundary, subtracted them from the data set. Then tried to use the boundary function again to get the inner boundary. The problem is, the data isn't very clean every time and sometimes outer boundary points aren't on the boundary line, so they get lumped in with the inner circle data. Ultimately, this gives me something like what is shown here in 'Surface' where it jumped the inner boundary out to a point that didn't get subtracted from the data set.

If there was a function like boundary() that instead grows the region out from a centroid that would be ideal and I could find the inner boundary by growing from the centroid. I thought about using a distance minimization method too, but the points are pretty random and I feel like it would still pick up on outer boundary points. Unless anyone has any suggestions about this approach too.

[John BG]

Hi Carolin

this is John BG <mailto:[e-mail] [e-mail]>

can you provide a sample set of points? or a script generating the points?

or should the readers scan the picture attached to your question and attempt answer from there?

appreciating time and attention, awaiting answer

John BG

<mailto:[e-mail] [e-mail]>

[Carolyn]

Attached is the data that goes along with this scatter plot.

[Jan]

And now you want to find the convex hull for the outer shape? Or should it be an alpha shape including the slightly concave parts? What is the wanted solution for the inner shape, which looks noisier? It would be useless to post a solution, when you do not exactly define, what you want.

Answer 1

[Akira Agata]

I've just come up with a good idea. How about clustering data points into 2 groups based on distance from the center point to each data point? Here is my try.

% Read the data file
D = xlsread('SampleData.xls');

% Calculate distance from center for each point
center = mean(D);
dist = vecnorm(D-center,2,2);

% Look at the distribution of the distance
figure
histogram(dist)

% Based on the histogram, set the threshold to separate 
% inner/outer ring with 140.
% If you have Statistics and Machine Learning Toolbox, you can 
% automatically separate by using unsupervised clustering method,
% such as k-means.
idxOuter = dist > 140;
idxInner = ~idxOuter;

% Show the result
figure
scatter(D(idxOuter,1),D(idxOuter,2))
hold on
scatter(D(idxInner,1),D(idxInner,2))
box on
grid on
legend({'Outer ring','Inner ring'})

[Carolyn]

This is great! I didn't even think about this kind of approach. Thank you!

[Image Analyst]

That is a clever approach. But it will work for only certain images. It will not work for images where the inner boundary gets really close to the outer boundary, like at the 7 o'clock position of the binary image you posted in your original post at the top.

Answer 2

[Jan]

Is the "outer boundary" of the points the convex hull? If so, you can determine the points on this hull at first. Afterwards you can transform the positions: Subtract the center of all points, then replace each point by 1/r(i), while r(i) is the distance to the center. This moves the inner points to the outside. Then find the indices of the points on the convex hull again, to find the "inner boundary".

See also: https://www.mathworks.com/help/matlab/ref/alphashape.html, especially the example: fill holes.

[Carolyn]

Yes. The outer boundary is the convex hull though I didn't use the convex hull function to determine the boundary. I have tried using the alphaShape function, but each image I'm trying to process is fairly irregular. For example:

This is the best fit I get after adjusting the radius, but the same radius can't be used for every data set in the stack. The same parameters used on a different image give this:

Could you elaborate on the transformation you're suggesting? What center? The centroid of the boundary? What would be the best way to find this value?

[Jan]

If you provide some input data, I could test some ideas before claiming, that they are working.

With "center" I meant the mean value of the points or the center of mass.

Answer 3

[Image Analyst]

I think this could be a case of someone asking to do something that is what they think they need to do but really don't. So you're starting with a binary image that is a slice of a 3-D volumetric image. And you say the main problem is "I need to find the inner and outer boundary of the points." This is easily done without doing fitting/regression,smoothing, alpha shapes or whatever. You simply call bwboundaries():.

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 = 20;

% Check that user has the specified Toolbox installed and licensed.
hasLicenseForToolbox = license('test', 'image_toolbox');   % license('test','Statistics_toolbox'), license('test','Signal_toolbox')
if ~hasLicenseForToolbox
  % User does not have the toolbox installed, or if it is, there is no available license for it.
  % For example, there is a pool of 10 licenses and all 10 have been checked out by other people already.
  ver % List what toolboxes the user has licenses available for.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end

%===============================================================================
% Read in gray scale demo image.
folder = pwd; % Determine where demo folder is (works with all versions).
baseFileName = 'surface.jpg';
% Get the full filename, with path prepended.
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
rgbImage = imread(fullFileName);
% Display the image.
subplot(2, 3, 1);
imshow(rgbImage, []);
title('Original Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis on;
hp = impixelinfo();

% Get the dimensions of the image.  
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(rgbImage);
if numberOfColorChannels > 1
  % It's not really gray scale like we expected - it's color.
  % Use weighted sum of ALL channels to create a gray scale image.
%   grayImage = rgb2gray(rgbImage); 
  % ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
  % which in a typical snapshot will be the least noisy channel.
  grayImage = rgbImage(:, :, 1); % Take blue channel.
else
  grayImage = rgbImage; % It's already gray scale.
end
% Now it's gray scale with range of 0 to 255.

% Display the image.
subplot(1, 2, 1);
imshow(grayImage, []);
title('Original Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis on;

%------------------------------------------------------------------------------
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% 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;

%------------------------------------------------------------------------------
% Threshold the image.
binaryImage = grayImage > 127;
% Get rid of the white rectangular frame around the image.
binaryImage = imclearborder(binaryImage);
% Display the image.
subplot(1, 2, 2);
imshow(binaryImage, []);
title('Binary Image with Boundaries', 'FontSize', fontSize, 'Interpreter', 'None');
axis on;
hp = impixelinfo();

% Get the outer and inner boundary.
boundaries = bwboundaries(binaryImage);

% Plot outer boundaries on image.
numberOfBoundaries = size(boundaries, 1);
hold on;
for k = 1 : numberOfBoundaries
  thisBoundary = boundaries{k};
  plot(thisBoundary(:,2), thisBoundary(:,1), '-', 'LineWidth', 5);
end
hold off;

Or you fill the shape and call bwboundaries(), then extract the inner black zone and call bwboundaries() again. This will give you the boundary of the outer edge of the central black region, as opposed to the inner edge of the white ring like the above 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 = 20;

% Check that user has the specified Toolbox installed and licensed.
hasLicenseForToolbox = license('test', 'image_toolbox');   % license('test','Statistics_toolbox'), license('test','Signal_toolbox')
if ~hasLicenseForToolbox
  % User does not have the toolbox installed, or if it is, there is no available license for it.
  % For example, there is a pool of 10 licenses and all 10 have been checked out by other people already.
  ver % List what toolboxes the user has licenses available for.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end

%===============================================================================
% Read in gray scale demo image.
folder = pwd; % Determine where demo folder is (works with all versions).
baseFileName = 'surface.jpg';
% Get the full filename, with path prepended.
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
rgbImage = imread(fullFileName);
% Display the image.
subplot(2, 3, 1);
imshow(rgbImage, []);
title('Original Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis on;
hp = impixelinfo();

% Get the dimensions of the image.  
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(rgbImage);
if numberOfColorChannels > 1
  % It's not really gray scale like we expected - it's color.
  % Use weighted sum of ALL channels to create a gray scale image.
%   grayImage = rgb2gray(rgbImage); 
  % ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
  % which in a typical snapshot will be the least noisy channel.
  grayImage = rgbImage(:, :, 1); % Take blue channel.
else
  grayImage = rgbImage; % It's already gray scale.
end
% Now it's gray scale with range of 0 to 255.

% Display the image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis on;

%------------------------------------------------------------------------------
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% 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;

%------------------------------------------------------------------------------
% Threshold the image.
binaryImage = grayImage > 127;
% Get rid of the white rectangular frame around the image.
binaryImage = imclearborder(binaryImage);
% Display the image.
subplot(2, 2, 2);
imshow(binaryImage, []);
title('Binary Image without White Frame', 'FontSize', fontSize, 'Interpreter', 'None');
axis on;
hp = impixelinfo();

% Invert the image so the center zone is white and the outer zone is white.
% Then call imclearborder again to remove the outer zone.
innerZone = imclearborder(~binaryImage);
% Display the image.
subplot(2, 2, 3);
imshow(innerZone, []);
title('Inner Zone', 'FontSize', fontSize, 'Interpreter', 'None');
axis on;
% Dilate by one pixel if you want the pixels along the inside edge of the white rather than the outside edge of the black.

% Get an image of the outer layer of the black inner zone.
% Most probably this IS NOT NEEDED but for what it's worth, it's here.
perimImage = bwperim(innerZone);
% XOR them to get the ring
% perimeterMask = xor(filledImage, erodedImage);
% Display the image.
subplot(2, 2, 4);
imshow(perimImage, []);
title('Perimeter Mask Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis on;
hp = impixelinfo();

% Get the outer boundary.
outerZone = imfill(binaryImage, 'holes');
outerBoundary = bwboundaries(outerZone);
% Get the inner boundary.
innerBoundary = bwboundaries(innerZone);

% Plot outer boundaries on image.
numberOfBoundaries = size(outerBoundary, 1);
hold on;
for k = 1 : numberOfBoundaries
  thisBoundary = outerBoundary{k};
  plot(thisBoundary(:,2), thisBoundary(:,1), 'r-', 'LineWidth', 2);
end
% Plot inner boundaries on image.
numberOfBoundaries = size(innerBoundary, 1);
for k = 1 : numberOfBoundaries
  thisBoundary = innerBoundary{k};
  plot(thisBoundary(:,2), thisBoundary(:,1), 'g-', 'LineWidth', 2);
end
hold off;

If you want the boundaries smoothed for some reason, simply blur them and threshold again before calling bwboundaries():

windowWidth = 55; % Whatever...
kernel = ones(windowWidth) / windowWidth ^ 2;
binaryImage = conv2(double(binaryImage), kernel, 'same') > 0.5;

[Carolyn]

I appreciate this approach, but the binary image is my final output. I start with just a scatter plot of the inner and outer boundaries as shown above. The binary image here is wrong because the boundaries of the inner and outer rings were not found correctly. I would use bwboundaries() if the scatter data were connected in the binary image, but it isn't so I have tried to recreate the boundaries using the boundary() and inpolygon() function. But the boundary() function isn't identifying the inner boundary correctly.

[Image Analyst]

When you say "I have a stack of 3D data that has been split up into slices." how was the stack created? What was the data source that allowed you to create the stack? A 3-D image or something else? And the "stack" itself - you say it's a bunch of points that can be shown with scatter() - so is the stack essentially an N-by-3 stack of (x,y,z) coordinates of those points? How were those points determined from the original source?

[Carolyn]

I started with an STL, then I found the point cloud of edge vertices of the 3D stl. These vertices were given in (xyz) coordinates that I then separated into slices based on the z coordinate.

Answer 4

[Carolyn]

Thanks everyone for the help. I figured out a work around for the problem.

[arnold]

may I ask how?

[Carolyn]

I created two different objects, found the boundary using the boundary function, then used polyshape to find all points within those two boundaries. I then converted that list of points to binary using poly2mask and subtracted the inner object from the outer object to get just the boundary. Hope this helps. Having access to two different objects at the beginning was key and may not work for every situation.