Step2_import_bw_slices.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB CODES ACCOMPANYING QUAN ET AL. (2021) PAPER
% CODES CALCULATE POROSITY ON PROCESSED X-RAY CT IMAGES
% 
% STEP2: IMPORT BLACK/WHITE IMAGES 
% REFER TO README.MD FOR INSTRUCTION
%
% CITE AND CREDIT:
% SUN ET AL. (2021). POWDER TECHNOLOGY, 388:496-504.
% HTTPS://DOI.ORG/10.1016/J.POWTEC.2021.05.006
% 
% TESTED ON MATLAB VERSION 2018(a) OR NEWER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; clear; close all;
%dbstop  if error

%% Enter the ouput 3D binary .mat file name
% output_filename = 'Loose'; % for example 1
output_filename = ['Steel']; % for example 2

%%
[fileName, pathName, ~] =...
    uigetfile(...
    {'*.jpg;*.tif;*.tiff;*.png;*.gif;*.bmp;','All Image Files'},...
    'Select stack images to be processed',...
    'MultiSelect','on'...
    );

%{
Determine the number of files to be processed; >= 0 and <= 2^16-1
    %}
    nFiles = uint16( size(fileName, 2) );
    
    %{
Convert file name from cell to string.
    %}
    fileName = char(fileName);
    
    %{
Get the pixel resolution of each image
    %}
    img = imread( [pathName fileName(1,:)] );
    imgSize = size(img);
    %{
Initiate a 3D matrix for the stack images
    %}
    if ( isa(img, 'double') )
        BW3D = zeros( imgSize(1), imgSize(2), nFiles );
    elseif ( isa(img, 'single') )
        BW3D = zeros( imgSize(1), imgSize(2), nFiles, 'single' );
    elseif ( isa(img, 'uint16') )
        BW3D = zeros( imgSize(1), imgSize(2), nFiles, 'uint16' );
    elseif ( isa(img, 'uint8') )
        BW3D = zeros( imgSize(1), imgSize(2), nFiles, 'uint8' );
    else
        BW3D = false( imgSize(1), imgSize(2), nFiles );
    end
    
    %{
Serial-reading of images to assign values in BW3D matrix
    %}
    for iFile=1:nFiles
        BW3D(:,:,iFile) = imread( [pathName fileName(iFile,:)] );
    end
    %
BW = BW3D;
    
    %{ 
     Binarize the image to logical format  (if the image is not logical type) and crop the margin
    %}
if ~isa(BW,'logical')
    BW = imbinarize(BW3D); 
end
    stats = regionprops3(BW,'Image','BoundingBox');
    BW=stats.Image{1,1};
    %%
%{
    save the 3D binary matrix  (Select the file name!)
   
%}
    save ([output_filename,'_packing.mat'],'BW');  % 

%{
    Visualize your packing by z-slices 
   Thanks to Maysam Shahedi for the visualization function.
   (https://www.mathworks.com/matlabcentral/fileexchange/41334-imshow3d)
%}
imshow3Dfull(BW)