%% EPS 207, 2011 Spring, Code by Juhyuk Moon
% Instructor: Burkhard Militzer
% GSI: Dino Bellugi

%%
% Image processing tool box & voidconst.m file are required. 
% The total size of original tomography images is huge (~4 GPa).
% If you want to re-run this code with the datafiles (250 Mb), please send
% an email. I will send you via email. 
% With the datafiles, please start this code from line "42"
%%
clear all; close all
% get directory listing_images
% 'images/*tif' '1KSH'
% '2KSHA/*tif'
% '9KSC/*tif'
listing=dir('10KSA/*tif');
% get number of files matching the query
nFiles = numel(listing);
% print all the file names using a loop
for n=1:nFiles
fprintf([listing(n).name '\n']);
end
%Is = image sample
Is = imread(['10KSA' '/' listing(1).name]);
[nRows nCols nLyrs] = size(Is);
nPix = nRows * nCols;

% target resize 
resize = 494;
cRow = round(nRows/2);
sRow = cRow-resize/2;
cCol = round(nCols/2);
sCol = cCol-resize/2;

% resize original image to {resize by resize}
for k = 1:nFiles;
  I = imread(['10KSA' '/' listing(k).name]);
  I2{k} = I(sRow:sRow+resize-1,sCol:sCol+resize-1);
end

%%%% resized dataset should be imported before this step

% for 3D image of void
[rRows rCols] = size(I2{1});
VM3=false(rRows,rCols,nFiles);

for k = 1:nFiles;
    VMs=false(rRows,rCols);
    I=I2{k};
    % call voidconst function
    VM = voidconst(I);
    VMs(VM(:,:)==1)=1;
    VM3(:,:,k)=VMs;
end

% 3D image re-construction
v=VM3;
p = patch(isosurface(v,0.5));
set(p,'FaceColor','red','EdgeColor','none');
daspect([1 1 1])
view(3); axis tight
camlight 
lighting gouraud

% void distribution calculation
CC = bwconncomp(VM3,18)
voiddata = regionprops(CC, 'basic');
void_vols = [voiddata.Area];
void_vols = void_vols*(0.0087^3);

% void distribution plot void volume - void number
figure, plot(void_vols)
title('void volume distribution');
xlabel('void number');
ylabel('void volume (mm^3)');

% void distribution plot cumulative void count - void volume
n=round(max(void_vols));
separators = [0:0.00001:0.4];
myHist = histc(void_vols, separators);
myHist_n = myHist ./ sum(myHist);
myDistr = cumsum(myHist_n);
figure, semilogx(separators,myDistr,'-ro');
title('cumulative void count as a function of void volume');
axis([10^(-4.5) n 0 1.05])
xlabel('log(void volume,(um^3))');
ylabel('cumulative void count');

% void distribution plot void count - void volume
figure, semilogx(separators,myHist,'-ro');
title('number of voids as a function of void volume');
axis([10^(-4.5) n 0 160])
xlabel('log(void volume,(um^3))');
ylabel('number of voids');



