Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Surface area of voxels in a 3D binary image

Subject: Surface area of voxels in a 3D binary image

From: Nathan

Date: 31 May, 2011 22:22:02

Message: 1 of 3

I have a 3D matrix which contains clusters of 1's (Each 1 represents a voxel that is designated as a pore in a sample of hydrated cement). I want to find the surface area of the pores in the space. As a simple example, let's say A is such a 3D matrix:

A(:,:,1) =
     0 0 0 0
     0 0 0 0
     0 0 0 0
A(:,:,2) =
     0 0 0 0
     0 1 1 0
     0 0 0 0
A(:,:,3) =
     0 0 0 0
     0 0 0 0
     0 0 0 0

The way that I am visualizing this, those two 1's in the middle of A(:,:,2) are two 1um x 1um x 1um cubes. In this case, I would consider the surface area to be 10 um^2 (or the exterior area of two cubes that share a face).

I have sifted through files on related to the image processing toolbox, and have found nothing that answers my question (I should make the disclaimer that I didn't find one that was obvious to me, at least). I appreciate any help on this!

Subject: Surface area of voxels in a 3D binary image

From: ImageAnalyst

Date: 1 Jun, 2011 02:23:59

Message: 2 of 3

Nathan:
That's a very good question, but I believe the answer is no, there's
no simple built-in function to do that, that I'm aware of. However
you can get it in a few lines by doing a convolution to count up the
exposed faces, then labeling and counting the pixel values. But it's
a little bit tricky/advanced so I thought I'd better do it myself
rather than let you flail away with just a few hints. Here I do it
with your example, although I added a single-voxel blob so that we
would have two regions. Of course the single voxel has an exposed
surface area of 6. Here's the demo for 6-connected blobs. The
extension to 18 or 26 connected blobs should be straightforward in
case you want that.

Instructions: copy the code below, paste into a code editor window and
save as "surfaceArea3D.m" then fix the line breaks (if any were
introduced by the newsreader) and run it.

% IMPORTANT: The newsreader may break long lines into multiple lines.
% Be sure to join any long lines that got split into multiple single
lines.
% These can be found by the red lines on the left side of your
% text editor, which indicate syntax errors, or else just run the
% code and it will stop at the split lines with an error.
%
% Demo function to calculate the exposed surface area
% of a 6-connected 3D binary volume.
function surfaceArea3D()
clc; % Clear the command window.
clear;
workspace;

% Generate sample data set of 2 6-connected blobs.
% It has one blob of 2 voxels and one blob of 1 voxel.
A(:,:,1) =[...
     0 0 0 0
     0 0 0 0
     0 0 0 0]
A(:,:,2) =[...
     0 0 0 0
     0 1 1 0
     0 0 0 0]
A(:,:,3) =[...
     0 0 0 0
     0 0 0 0
     1 0 0 0 ]

% Construct kernel where we can count the
% number of 6-connected neighbor voxels.
conn6Kernel = zeros([3,3,3]);
conn6Kernel(2,2,1) = 1;
conn6Kernel(1,2,2) = 1;
conn6Kernel(2,1,2) = 1;
conn6Kernel(2,3,2) = 1;
conn6Kernel(3,2,2) = 1;
conn6Kernel(2,2,3) = 1;
% Display it in the command window for the user to see.
conn6Kernel

% For each voxel, determine how many 6-connected neighbors it has.
sumOfFaces = convn(A, conn6Kernel, 'same')

% Find number of exposed faces for each voxel.
surfaceArea = 6 * A - sumOfFaces
% Mask out zero voxels that have negative exposed faces.
surfaceArea(surfaceArea<0) = 0

% Now we simply label the volume and sum up the values of each region.
binaryVolume = surfaceArea > 0
cc = bwconncomp(binaryVolume, 6);

% Now it's labeled, so now we measure the PixelValues in each blob.
measurements = regionprops(cc, surfaceArea, 'PixelValues');

% Go through the regions, listing each regions exposed surface area.
numberOfRegions = length(measurements);
for k = 1 : numberOfRegions
% For this region, find out the sum of exposed faces.
thesePixelValues = [measurements(k).PixelValues]
% Sum up the number of exposed faces for each voxel in the blob.
thisRegionsArea = sum(thesePixelValues);
fprintf('The exposed surface area for region %d is %d\n',...
k, thisRegionsArea);
end

% Alert user that we're done.
message = sprintf('Done with demo.\nLook in the command window.');
uiwait(msgbox(message));
%========================================

What you should see:


A =

     0 0 0 0
     0 0 0 0
     0 0 0 0


A(:,:,1) =

     0 0 0 0
     0 0 0 0
     0 0 0 0


A(:,:,2) =

     0 0 0 0
     0 1 1 0
     0 0 0 0


A(:,:,1) =

     0 0 0 0
     0 0 0 0
     0 0 0 0


A(:,:,2) =

     0 0 0 0
     0 1 1 0
     0 0 0 0


A(:,:,3) =

     0 0 0 0
     0 0 0 0
     1 0 0 0


conn6Kernel(:,:,1) =

     0 0 0
     0 1 0
     0 0 0


conn6Kernel(:,:,2) =

     0 1 0
     1 0 1
     0 1 0


conn6Kernel(:,:,3) =

     0 0 0
     0 1 0
     0 0 0


sumOfFaces(:,:,1) =

     0 0 0 0
     0 1 1 0
     0 0 0 0


sumOfFaces(:,:,2) =

     0 1 1 0
     1 1 1 1
     1 1 1 0


sumOfFaces(:,:,3) =

     0 0 0 0
     1 1 1 0
     0 1 0 0


surfaceArea(:,:,1) =

     0 0 0 0
     0 -1 -1 0
     0 0 0 0


surfaceArea(:,:,2) =

     0 -1 -1 0
    -1 5 5 -1
    -1 -1 -1 0


surfaceArea(:,:,3) =

     0 0 0 0
    -1 -1 -1 0
     6 -1 0 0


surfaceArea(:,:,1) =

     0 0 0 0
     0 0 0 0
     0 0 0 0


surfaceArea(:,:,2) =

     0 0 0 0
     0 5 5 0
     0 0 0 0


surfaceArea(:,:,3) =

     0 0 0 0
     0 0 0 0
     6 0 0 0


binaryVolume(:,:,1) =

     0 0 0 0
     0 0 0 0
     0 0 0 0


binaryVolume(:,:,2) =

     0 0 0 0
     0 1 1 0
     0 0 0 0


binaryVolume(:,:,3) =

     0 0 0 0
     0 0 0 0
     1 0 0 0


thesePixelValues =

     5
     5

The exposed surface area for region 1 is 10

thesePixelValues =

     6

The exposed surface area for region 2 is 6

Subject: Surface area of voxels in a 3D binary image

From: Sarah Miller

Date: 4 Dec, 2012 20:02:08

Message: 3 of 3

ImageAnalyst,

I am trying to do a 26 connected version of this code and would like to make sure that I am doing it correctly.

1) For a 26 connected kernel, I would create a 3x3x3 matrix, each entity equal to 1 except the center voxel (2,2,2), which would equal 0, correct?

2) Are there any other adjustments, specifically to (calculations of) sumOfFaces, surfaceArea, cc (=bwconncomp(...))?

Thank you for your help!

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us