% Compile voxel images of voids in un-deformed states from x-ray CT images
% Use the x-y-z coordinates of feature centroids identified by DREAM3D, and
% load a 91x91x91 box around that void
%
% Orion L. Kafka
%
% History:
% 7/30/20 Written
%
clear all; close all; clc

% Input locations
meanLineHeight = 200;
specimenName = 'V3-6';
posID = specimenName(end);

mode = 'noncontinuous';

featureFolderPath =...
 ['C:\Pos-',posID,'-scans\',specimenName,'-mats\'];

outputFolderPath = 'C:\Users\olk\Desktop\in-situ-paper\scanDataPlots\';

if strcmp(specimenName,'H2-4')
    expVect = 175:181;
elseif strcmp(specimenName,'V2-4')
    expVect = 191:198;
elseif strcmp(specimenName,'H3-4')
    expVect = 185:189;
elseif strcmp(specimenName,'V3-4')
    expVect = 18:24;
elseif strcmp(specimenName,'H1-4')
    expVect = 137:139;
    mode = 'continuous';
elseif strcmp(specimenName,'V1-4')
    expVect = 70:72;
    mode = 'continuous';
elseif strcmp(specimenName,'H2-2')
    expVect = 106:113;
elseif strcmp(specimenName,'H3-2')
    expVect = 126:130;
elseif strcmp(specimenName,'V2-2')
    expVect = 80:86;
elseif strcmp(specimenName,'V3-2')
    expVect = 156:160;
elseif strcmp(specimenName,'V2-1')
    expVect = 205:211;
elseif strcmp(specimenName,'H2-6')
    expVect = 116:123;
elseif strcmp(specimenName,'H3-6')
    expVect = 150:152;
    mode = 'continuous';
elseif strcmp(specimenName,'V2-6')
    expVect = 51:58;
elseif strcmp(specimenName,'V3-6')
    expVect = 170:172;
    mode = 'continuous';
end

namePrefix = 'Feat';
namePostfix1 = '_0015-0522';
namePostfix2 = '_0523-1044';
namePostfix3 = '_1045-1566';
namePostfix4 = '_1567-2144';

nameExtension = '.csv';

%expVect = 156:160; % V3-2
%expVect = 126:130; % H3-2
%expVect = 80:86; % V2-2
%expVect = 106:113; % H2-2

%namePostfix1 = '_0001-1044';
%namePostfix2 = '_1045-2088';

count = 1;
for expNum = expVect

    featTable1 = importdata(fullfile(featureFolderPath,...
        [namePrefix,padWithZeros(3,expNum),namePostfix1,nameExtension]),',',2);
    featTable2 = importdata(fullfile(featureFolderPath,...
        [namePrefix,padWithZeros(3,expNum),namePostfix2,nameExtension]),',',2);
    featTable3 = importdata(fullfile(featureFolderPath,...
        [namePrefix,padWithZeros(3,expNum),namePostfix3,nameExtension]),',',2);
    featTable4 = importdata(fullfile(featureFolderPath,...
        [namePrefix,padWithZeros(3,expNum),namePostfix4,nameExtension]),',',2);

    numFeatures1 = str2double(featTable1.textdata{1,1});
    numFeatures2 = str2double(featTable2.textdata{1,1});
    numFeatures3 = str2double(featTable3.textdata{1,1});
    numFeatures4 = str2double(featTable4.textdata{1,1});

    numFeatures = numFeatures1 + numFeatures2 + numFeatures3 + numFeatures4;
    
    featData1 = featTable1.data;
    featData2 = featTable2.data;
    featData3 = featTable3.data;
    featData4 = featTable4.data;
    
    % remove features that touch the outter boundary
    featDataInner1 = featData1(featData1(:,end) == 0,:);
    featDataInner2 = featData2(featData2(:,end) == 0,:);
    featDataInner3 = featData3(featData3(:,end) == 0,:);
    featDataInner4 = featData4(featData4(:,end) == 0,:);
    
    % add the z-stack height back in
    featDataInner2(:,13) = featDataInner2(:,13) + 523*0.65;
    featDataInner3(:,13) = featDataInner3(:,13) + 1045*0.65;
    featDataInner4(:,13) = featDataInner4(:,13) + 1567*0.65;
    
    
    % combine top and bottom halves (1-1044 and 1045-2088)
    featDataInner = vertcat(featDataInner1,featDataInner2,featDataInner3,featDataInner4);
    
    % pick up only things bigger than 64 voxels
    featDataInner = featDataInner(featDataInner(:,15) > 125,:);
    
    % sort by volume
%    featDataInner = sortrows(featDataInner,17);
    
    % pick up only the 100 biggest
%    featDataInner = featDataInner(end-99:end,:);
    
    
    centroids1{count} = featDataInner(:,11);
    centroids2{count} = featDataInner(:,12);
    centroids3{count} = featDataInner(:,13);

    AspectRatios0{count} = featDataInner(:,3);
    AspectRatios1{count} = featDataInner(:,4);
    
    EulerAngs0{count} = featDataInner(:,5);
    EulerAngs1{count} = featDataInner(:,6);
    EulerAngs2{count} = featDataInner(:,7);
    
    sizes{count} = featDataInner(:,17);
    
    maxVoidSize(count) = max(featDataInner(:,17));
    
    totalVol(count) = sum(featDataInner(:,17));
    
    count = count + 1;
end

if strcmp(mode,'continuous')
    AspectRatios0{1,2} = vertcat(AspectRatios0{1,2},AspectRatios0{1,3});
    AspectRatios1{1,2} = vertcat(AspectRatios1{1,2},AspectRatios1{1,3});
    totalVol1 = totalVol(end-1) + totalVol(end);
    totalVol = [totalVol(1),totalVol1];
end

figure(1)
plot(totalVol,'O:')
xlabel('Image num')
ylabel('Total void volume')


dcols = ['b','g','m','k','y','c','r','w'];

%{
figure(2)
for jj = 1:max(size(AspectRatios0))
    ptSizes = sizes{1,jj} / maxVoidSize(jj);
    for ii = 1:length(sizes{1,1})
        hold on
        plot(AspectRatios0{1,jj}(ii),AspectRatios1{1,jj}(ii),...
            [dcols(jj),'.'],...,'color',[0.5,0.5,jj/max(size(AspectRatios0))],...
            'markersize',70*ptSizes(ii))
        hold off
    end
end
xlim([0,1])
ylim([0,1])
xlabel('C/A (where C>B>A)')
ylabel('B/A')
%}
%

fg4 = figure(4);
histfit(AspectRatios0{1,1});
hold on
line([mean(AspectRatios0{1,1}),...
      mean(AspectRatios0{1,1})],[0,meanLineHeight],'color',[1,0,0],'linewidth',2,'linestyle','-.')

hs2 = histfit(AspectRatios0{1,max(size(AspectRatios0))});
set(hs2(2),'color',[0,0,0],'linestyle','--')
line([mean(AspectRatios0{1,max(size(AspectRatios0))}),...
      mean(AspectRatios0{1,max(size(AspectRatios0))})],[0,meanLineHeight],...
      'color',[0,0,0],'linewidth',2,'linestyle',':')

hold off
xlabel(featTable1.colheaders(3))
ylabel('Count')
exportgraphics(fg4,fullfile(outputFolderPath,[specimenName,'_AR0.eps']))
legend('Initial state','Initial Normal Density Function', 'Initial Mean', 'Final State', 'Final Normal Density Function', 'Final Mean')


fg5 = figure(5);
histfit(AspectRatios1{1,1});
hold on
line([mean(AspectRatios1{1,1}),...
      mean(AspectRatios1{1,1})],[0,meanLineHeight],'color',[1,0,0],'linewidth',2,'linestyle','-.')

hs2 = histfit(AspectRatios1{1,max(size(AspectRatios1))});
set(hs2(2),'color',[0,0,0],'linestyle','--')
line([mean(AspectRatios1{1,max(size(AspectRatios1))}),...
      mean(AspectRatios1{1,max(size(AspectRatios1))})],[0,meanLineHeight],...
      'color',[0,0,0],'linewidth',2,'linestyle',':')

%histfit(AspectRatios1{1,3});
hold off
xlabel(featTable1.colheaders(4))
ylabel('Count')
exportgraphics(fg5,fullfile(outputFolderPath,[specimenName,'_AR1.eps']))
legend('Initial state','Initial Normal Density Function', 'Initial Mean', 'Final State', 'Final Normal Density Function', 'Final Mean')
%}
%{
figure(6)
bar([sizes{1,1}(end-10:end),sizes{1,2}(end-10:end)]);
%hold on
%bar(sizes{1,2}(end-10:end));
xlabel('Count')
ylabel('Shape Volume')
legend('80','86')
hold off
%

figure(7)
histogram(EulerAngs0{1,1})
hold on
histogram(EulerAngs0{1,2})
histogram(EulerAngs0{1,3})
hold off
xlabel(featTable1.colheaders(5))
ylabel('Count')
legend('106','109','113')

figure(8)
histogram(EulerAngs1{1,1})
hold on
histogram(EulerAngs1{1,2})
histogram(EulerAngs1{1,3})
xlabel(featTable1.colheaders(6))
ylabel('Count')
legend('106','109','113')
hold off

figure(9)
histogram(EulerAngs2{1,1})
hold on
histogram(EulerAngs2{1,2})
histogram(EulerAngs2{1,3})
xlabel(featTable1.colheaders(7))
ylabel('Count')
legend('106','109','113')
hold off
%}
%{    
    figure(1)
    hold on
    histogram(featDataInner(:,3))
    xlabel(featTable1.colheaders(3))
    ylabel('Count')
    legend('106','113')
    %legend('106','107','108','109','110','111','112','113')
    hold off

    figure(2)
    hold on
    histogram(featDataInner(:,4))
    xlabel(featTable1.colheaders(4))
    ylabel('Count')
    hold off
%
    figure(expNum)
    histogram(featDataInner(:,17))
    xlabel(featTable1.colheaders(17))
    ylabel('Count')
%
    figure(4)
    hold on
    histogram(featDataInner(:,6))
    xlabel(featTable.colheaders(6))
    ylabel('Count')
    hold off
    
    figure(5)
    hold on
    histogram(featDataInner(:,7))
    xlabel(featTable.colheaders(7))
    ylabel('Count')
    hold off
%}