MATLAB Answers

0

Expand to find inner boundary of object

Asked by Carolyn on 17 May 2018 at 22:47
Latest activity Answered by Carolyn on 21 May 2018 at 0:23

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.

  3 Comments

Hi Carolin

this is John BG jgb2012@sky.com

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

jgb2012@sky.com

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

Jan
on 20 May 2018 at 13:12

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.

Sign in to comment.

3 Answers

Answer by Jan
on 18 May 2018 at 10:26
Edited by Jan
on 18 May 2018 at 10:28

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.

  2 Comments

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
on 18 May 2018 at 19:43

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.

Sign in to comment.


Answer by Image Analyst
on 20 May 2018 at 14:22

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;

  2 Comments

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.

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?

Sign in to comment.


Answer by Carolyn on 21 May 2018 at 0:23

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

  0 Comments

Sign in to comment.


Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today