Expand to find inner boundary of object
46 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
Jan
2018 年 5 月 20 日
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.
採用された回答
Akira Agata
2018 年 5 月 25 日
編集済み: Akira Agata
2018 年 5 月 25 日
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'})
2 件のコメント
Image Analyst
2018 年 5 月 25 日
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.
その他の回答 (3 件)
Jan
2018 年 5 月 18 日
編集済み: Jan
2018 年 5 月 18 日
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 件のコメント
Jan
2018 年 5 月 18 日
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
2018 年 5 月 20 日
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;
3 件のコメント
Image Analyst
2018 年 5 月 20 日
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?
参考
カテゴリ
Help Center および File Exchange で Computer Vision with Simulink についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!