## Expand to find inner boundary of object

### Carolyn (view profile)

on 17 May 2018
Latest activity Commented on by Carolyn

### Carolyn (view profile)

on 15 Jan 2019 at 0:52

### Akira Agata (view profile)

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

### John BG (view profile)

on 17 May 2018
Hi Carolin
can you provide a sample set of points? or a script generating the points?
appreciating time and attention, awaiting answer
John BG
Carolyn

### Carolyn (view profile)

on 18 May 2018
Attached is the data that goes along with this scatter plot.
Jan

### Jan (view profile)

on 20 May 2018
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.

R2018a

### Akira Agata (view profile)

on 25 May 2018
Edited by Akira Agata

### Akira Agata (view profile)

on 25 May 2018

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.
% 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

on 25 May 2018
Image Analyst

### Image Analyst (view profile)

on 25 May 2018
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.

on 18 May 2018
Edited by Jan

### Jan (view profile)

on 18 May 2018

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".

Carolyn

### Carolyn (view profile)

on 18 May 2018
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

### Jan (view profile)

on 18 May 2018
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.

### Image Analyst (view profile)

on 20 May 2018

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.
% 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');
% 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
% 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.
% 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');
% 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
% 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
% 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

### Carolyn (view profile)

on 20 May 2018
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

### Image Analyst (view profile)

on 20 May 2018
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

### Carolyn (view profile)

on 24 May 2018
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.

### Carolyn (view profile)

on 21 May 2018

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

arnold

### arnold (view profile)

on 12 Jan 2019 at 4:34