what is the best way to determine the boundary of a porous object in an 2-D image?

edge detection returns many internal edges which is not desired sample image:

http://imageshack.us/photo/my-images/836/36543017.png

Answer 1

[Image Analyst]

According to your definition of interior, here is code that calculates it. If you want the "reference" area to be the convex hull rather than the entire image, you can use bwconvhull().

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures.
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format longg;
format compact;
fontSize = 20;
% http://imageshack.us/photo/my-images/836/36543017.png/
% Read in a standard MATLAB gray scale demo image.
folder = 'C:\Users\mark\Documents\Temporary';
baseFileName = '36543017.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
  errorMessage = sprintf('Error: %s does not exist in the %s folder.', baseFileName, folder);
  uiwait(warndlg(errorMessage));
  % 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 either.', 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);
if numberOfColorBands > 1
  grayImage = grayImage(:,:,1);
end
% Display the original gray scale image.
subplot(2, 2, 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, 2, 2); 
% bar(pixelCount);
% grid on;
% title('Histogram of original image', 'FontSize', fontSize);
% xlim([0 grayLevels(end)]); % Scale x axis manually.

% Threshold image.
binaryImage = grayImage < 65000;
% binaryImage = imfill(binaryImage, 'holes');
% Display the image.
subplot(2, 2, 2);
imshow(binaryImage, []);
title('Binary Image', 'FontSize', fontSize);

% Close gaps
binaryImage = imclose(binaryImage, true(3));
% Display the image.
subplot(2, 2, 3);
imshow(binaryImage, []);
title('After morphological closing', 'FontSize', fontSize);

% Get the interiors
interiors = imclearborder(~binaryImage);
% Get rid of blobs smaller than 60 pixels.
interiors = bwareaopen(interiors, 60);
% Display the image.
subplot(2, 2, 4);
imshow(interiors, []);
title('Interiors', 'FontSize', fontSize);

% Calculate the porosity
porosity = sum(interiors(:))/numel(interiors);
message = sprintf('The porosity (area fraction) of the interiors is %.4f', porosity);
uiwait(helpdlg(message));

[A]

Yes, i understand this, but plz note before all this I have to find the overal boundary, sth like what I drew manually here. http://imageshack.us/photo/my-images/248/testyj.png/

[Image Analyst]

I did note that. See where I mentioned that you could use bwconvhull()? Here, I'll do it for you. Just add this code to the code above:

% Get the envelope using bwconvhull
envelope = bwconvhull(interiors);
convexHullArea = sum(envelope(:));
boundary = bwboundaries(envelope);
% Display the boundary as an outline in the overlay.
hold on;
plot(boundary{1}(:, 2), boundary{1}(:, 1), 'r-', 'LineWidth', 2);
% Calculate porosity with respect to the envelope.
porosity2 = sum(interiors(:)) / convexHullArea;
message = sprintf('The porosity (area fraction) of the interiors is %.4f', porosity2);
uiwait(helpdlg(message));

Answer 2

[Image Analyst]

Use a better edge detector. You'd get a better answer if you uploaded your image so I could then recommend more appropriate functions than edge(). See this http://www.mathworks.com/matlabcentral/answers/7924-where-can-i-upload-images-and-files-for-use-on-matlab-answers for places where you can upload images. Pick one where I can see it right away and not have to download or signup for anything.

If you have the outer edge intact, and just need to get rid of spurious internal edges, then you can use imfill(), but that will not work if the outer edge has a gap in it. Gaps can be closed with imclose() or edge linking methods.

[Image Analyst]

I see you posted your image now. It looks like something that at one time had nice solid boundaries that were completely closed, but now was imaged in some way such that the boundaries are now broken. What I would do it do first to convert to gray level, then invert and threshold. Then call imclose to close the gaps. Then you can call bwboundaries() or bwperim() though I doubt that is the end of your task. Chances are that once you have the boundaries, you are going to want to do something else with them, like measure their area, count them, find the area fraction, or something. What is your final goal with this? I'm sure it's not just finding boundaries. Unfortunately it looks like your image somehow got some bad jpeg artifacts upon uploading. Nonetheless I did the best I could with a few minutes of programming and this is the result:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures.
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format longg;
format compact;
fontSize = 20;
% http://imageshack.us/photo/my-images/836/36543017.png/
% Read in a standard MATLAB gray scale demo image.
folder = 'C:\Users\A\Documents\Temporary';
baseFileName = '36543017.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);
if numberOfColorBands > 1
  grayImage = grayImage(:,:,1);
end
% Display the original gray scale image.
subplot(2, 2, 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, 2, 2); 
% bar(pixelCount);
% grid on;
% title('Histogram of original image', 'FontSize', fontSize);
% xlim([0 grayLevels(end)]); % Scale x axis manually.

% Threshold image.
binaryImage = grayImage < 65000;
% binaryImage = imfill(binaryImage, 'holes');
% Display the image.
subplot(2, 2, 2);
imshow(binaryImage, []);
title('Binary Image', 'FontSize', fontSize);

% Close gaps
binaryImage = imclose(binaryImage, true(3));
% Display the image.
subplot(2, 2, 3);
imshow(binaryImage, []);
title('After morphological closing', 'FontSize', fontSize);

% Get the interiors
interiors = imclearborder(~binaryImage);
% Get rid of blobs smaller than 60 pixels.
interiors = bwareaopen(interiors, 60);
% Display the image.
subplot(2, 2, 4);
imshow(interiors, []);
title('Interiors', 'FontSize', fontSize);

% Label the image
labeledImage = bwlabel(interiors);
measurements = regionprops(labeledImage, 'Area', 'Solidity');
allAreas = [measurements.Area]
allSolidities = [measurements.Solidity];
% Get rid of weird-shaped objects
% Get a list of the blobs that meet our criteria and we need to keep.
allowableSolidityIndexes = (allSolidities > 0.8);
keeperIndexes = find(allowableSolidityIndexes);
% Extract only those blobs that meet our criteria, and
% eliminate those blobs that don't meet our criteria.
% Note how we use ismember() to do this.
keeperBlobsImage = ismember(labeledImage, keeperIndexes);
% Re-label with only the keeper blobs kept.
[labeledImage, numberOfBlobs] = bwlabel(keeperBlobsImage, 8);     % Label each blob so we can make measurements of it
binaryImage = labeledImage > 0;
figure;
subplot(2, 2, 1);
imshow(binaryImage, []);
title('Filtered by Solidity', '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') 

% Find boundaries by bwperim
% NOTE: THIS WILL LOOK BROKEN ON THE SCREEN DUE TO SUBSAMPLING/SHRINKING THE IMAGE
boundaries1 = bwperim(binaryImage);
subplot(2, 2, 2);
imshow(boundaries1, []);
title('Boundaries via bwperim()', 'FontSize', fontSize);

% Find the boundaries via bwboundaries.
% First display the original image.
subplot(2, 2, 3);
imshow(grayImage, []);
title('Original Image (Again)', 'FontSize', fontSize);
subplot(2, 2, 4);
imshow(grayImage, []);
title('Boundaries via bwboundaries()', 'FontSize', fontSize);
% 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 coins on the original grayscale image using the coordinates returned by bwboundaries.
hold on;
boundaries = bwboundaries(binaryImage);
numberOfBoundaries = size(boundaries, 1);
for k = 1 : numberOfBoundaries
  thisBoundary = boundaries{k};
  plot(thisBoundary(:,2), thisBoundary(:,1), 'r-', 'LineWidth', 3);
end
hold off;

[A]

Thanks for answer, but it didnt work, I guess need to do some coding instead of using functions directly. Is there standard methods for this, Something like finding the closest nonzero neighbor pixels, etc...BTW, in imclose what is the concept of second variable _imclose(*, ???)_ The purpose is finding the porosity of area inside the boundary

[Image Analyst]

What does "doesn't work" mean? I gave you some sample code that found the insides of the shapes. It does work at doing that, and, given no clear guidance, that's what I decided to do. It was not meant to be a complete turnkey solution for all your images. I don't know what the void space is and what the solid space is. All you've given me is a binary image - no gray scale image for context. What is the boundary? Is it each bean-shaped thing? So you want the porosity inside each little bean shape? Or is the boundary the boundary of the entire group of shapes, and the insides of the individual shapes are porous (void space) and everything else (individual boundaries plus the interstitial space) is solid? Please clarify.

If a boundary is each small individual boundary, and you want the porosity inside each boundary, then I found the insides for you. But inside each is just uniform white. How can I find porosity in that case? Again, I don't know what's void space and what's solid.

If you want the envelope of the outside you can use the union convex hull, performed by bwconvexhull() of the Image Processing Toolbox, or you can use a restricted convex hull (I have code for that) or alpha shapes, both of which will allow the outer shape to bend into bays to allows the boundary to more closely follow the shapes at the outer part of the group. Or you could run imclose with a big kernel, call imfill() and then call bwboundaries. It just depends on what you consider the "outside boundary" of your region to be. The definition is not so clear cut if your region is the group of shapes.

You can look up imclose(). You'll see the second argument is the shape of the structuring element. You can tailor it to certain shapes, and I tried that, taking vertical and horizontal slivers, but a 3 by 3 box seemed to work best.

Answer 3

[Walter Roberson]

regionprops() and extract the boundary pixels list?

[Image Analyst]

Or bwboundaries() or bwperim() assuming you have the "right" image to start with.

Answer 4

[Walter Roberson]

Perhaps what you want to do is use something like John D'Errico's Alpha Shapes routines to wrap a boundary around the collection of points you have. After you have that boundary, regionprops can tell you the porosity.

[A]

Do you know a good literature to start reading about alpha-shape; is it sth like connecting close points with interpolation,,,,?

[Walter Roberson]

http://en.wikipedia.org/wiki/Alpha_shape

[Image Analyst]

http://biogeometry.duke.edu/software/alphashapes/index.html

http://cgm.cs.mcgill.ca/~godfried/teaching/projects97/belair/alpha.html

http://www.mathworks.com/matlabcentral/fileexchange/28851-alpha-shapes

http://www.mathworks.com/matlabcentral/fileexchange/6760-ashape-a-pedestrian-alpha-shape-extractor

[A]

Thanks guy*s*, the comment & documents was very useful

[Image Analyst]

But aren't you going to define porosity or share the gray scale image from where your binary image came from?

[A]

porosity was the goal. the problem was I had to separate inside the porous object from outside. For this I need a proximate of boundary

[Image Analyst]

Which is the porous part? The thick boundaries, or the interiors of them?

[A]

both inside and outside of thick boundary,it is not visible in current pic cuz of single threshholding.

[Image Analyst]

So, the thick boundary is solid, and both inside of that and outside of that is porous. Correct? But what do you want to know? The area fraction of porous pixels? That would be simply

porosity = sum(binaryImage(:))/numel(binaryImage).

[A]

the fraction of pores inside the outer boundary(not thick lines). So first I should define the boundary then I can use your command line. right?

[Image Analyst]

See code in my new answer.