clear
close all
O = imread('Case1-Panorama - correlation matching.tif');
%O=load('CropImage.mat')

%% data fix 
%%rect2=[33391 1093 2175 5919]
%rect2=[5000 1 200 7558]
%I = O(rect2(2):rect2(2)+rect2(4),rect2(1):rect2(1)+rect2(3)); % removing invalid areas
I = O(:,:); % removing invalid areas
I = uint8(I/256); % reducing memory use

%% investigatin bias in the image
%figure,
%subplot(2,2,1), imagesc(I), title('big image')
%subplot(2,2,2), plot(mean(I,2),1:size(I,1)), axis ij tight, title('row-wise mean')
%subplot(2,2,3), plot(1:size(I,2),mean(I,1)), axis ij tight, title('column-wise mean')

%% cropping out
p = round(size(I)/2); % center of crop-out
r = 1000; % width of crop-out
C = I((1-r/2:r/2)+p(1),(1-r/2:r/2)+p(2)); % crop-out

%% investigating different thresholding methods
h = hist(C(:),0:255);
th = otsuthresh(h)*255;

%figure
%subplot(3,3,1)
%imagesc(C), axis image, colormap gray, title('input image')

%subplot(3,3,2)
%BW_adaptive = imbinarize(C, 'adaptive');
%imagesc(BW_adaptive), axis image, title('adaptive')

%T = adaptthresh(C);
%subplot(3,3,3), imagesc(T*255), colorbar, axis image, title('adaptive th')

%subplot(3,3,4)
%bar(0:255,h), hold on, plot(th*[1,1],[0,max(h)]), axis tight square
%title('hist and otsu')

%subplot(3,3,5)
%imagesc(C>th), axis image, title('otsu')

%subplot(3,3,6)
%BW_global = imbinarize(C, 'global');
%imagesc(BW_global), axis image, title('global')

%T = adaptthresh(C,0.7,'NeighborhoodSize',[301,301]);
%subplot(3,3,7), imagesc(T*255), colorbar, axis image, title('adaptive th fixed')
%subplot(3,3,8), imagesc(C>T*255), axis image, title('adaptive fixed')

%% Output to paper
figure,
imagesc(C), axis image, colormap gray, title('input image')
figure,
bar(0:255,h), hold on, 
plot(th*[1,1],[0,max(h)]), axis tight square
figure,
imagesc(C>th), axis image, title('otsu')

%% computing volume fraction
%Vf = sum(BW_global(:))/numel(BW_global(:));

%% whole image
%%T = adaptthresh(I,0.7,'NeighborhoodSize',[301,301]);
%%%%BW = I>T*255;
BW = I>th;
figure, imagesc(BW), axis image, title('Otsu')
Vf = sum(BW(:))/numel(BW(:));

%% Vf variation gennem tykkelsen
%Vf2=sum(BW,2)/size(BW,2);
%Vf1=sum(BW,1)/size(BW,1);
%%
StrVf=['V_{f}^{mean}=',num2str(Vf,3)];
%figure,
%subplot(2,2,1), imagesc(I), title('big image')
%subplot(2,2,3), plot(1:size(BW,2),Vf1),  title(['column-wise ',StrVf])
%subplot(2,2,2), plot((Vf2),(1:size(BW,1))),  title(['row-wise ',StrVf])
%set(gca, 'YDir','reverse')
%%
%figure,
%plot(1:size(BW,2),Vf1)
%hold on
%plot(1:size(BW,2),movmean(Vf1,200),'-r')
%title(['column-wise ',StrVf])
%%
%figure,
%plot((Vf2),(1:size(BW,1)))
%hold on
%plot(movmean(Vf2,200),(1:size(BW,1)),'-r')
%title(['row-wise ',StrVf])
%set(gca, 'YDir','reverse')
%%
figure
imagesc(BW)
axis image
colorbar
%%
lpix=0.52725;
lsize=round(5*(16/lpix));
ave=ones(lsize,lsize)/lsize^2;
BWave=conv2(BW,ave,'same');
%%
figure
imagesc(BWave,[0.5 1])
axis off
axis image
colormap jet
colorbar
%%
%save('Case1Vf.mat','BW','BW1','BW2','BW3','BW4','BW5','BW6','BW7','BW8','maskedImage1','maskedImage2','maskedImage3','maskedImage4','maskedImage5','maskedImage6','maskedImage7','maskedImage8')
load 'Case1Vf.mat'
%%
Vfcall=Vf
Vfb(1)=sum(BW1.*BW,'all')/sum(BW1,'all')
Vfb(2)=sum(BW2.*BW,'all')/sum(BW2,'all')
Vfb(3)=sum(BW3.*BW,'all')/sum(BW3,'all')
Vfb(4)=sum(BW4.*BW,'all')/sum(BW4,'all')
Vfb(5)=sum(BW5.*BW,'all')/sum(BW5,'all')
Vfb(6)=sum(BW6.*BW,'all')/sum(BW6,'all')
Vfb(7)=sum(BW7.*BW,'all')/sum(BW7,'all')
Vfb(8)=sum(BW8.*BW,'all')/sum(BW8,'all')
%%
StrVfb1=['V_{f} = ',num2str(Vfb(1),3)];
StrVfb2=['V_{f} = ',num2str(Vfb(2),3)];
StrVfb3=['V_{f} = ',num2str(Vfb(3),3)];
StrVfb4=['V_{f} = ',num2str(Vfb(4),3)];
StrVfb5=['V_{f} = ',num2str(Vfb(5),3)];
StrVfb6=['V_{f} = ',num2str(Vfb(6),3)];
StrVfb7=['V_{f} = ',num2str(Vfb(7),3)];
StrVfb8=['V_{f} = ',num2str(Vfb(8),3)];
%%
Vfbmean=mean(Vfb)
Vfbstd=std(Vfb)
StrVf=['V_{f} = ',num2str(Vf,3),'; V_{f}^b = ',num2str(Vfbmean,3),'\pm',num2str(Vfbstd,3)];
%%
BW0=BW*10;
BW1a=BW1;
BW2a=BW2;
BW3a=BW3;
BW4a=BW4;
BW5a=BW5;
BW6a=BW6;
BW7a=BW7;
BW8a=BW8;
BW1a(~BW)=0;
BW2a(~BW)=0;
BW3a(~BW)=0;
BW4a(~BW)=0;
BW5a(~BW)=0;
BW6a(~BW)=0;
BW7a(~BW)=0;
BW8a(~BW)=0;
BW0(BW1a)=2;
BW0(BW2a)=3;
BW0(BW3a)=4;
BW0(BW4a)=5;
BW0(BW5a)=6;
BW0(BW6a)=7;
BW0(BW7a)=8;
BW0(BW8a)=9;
%%
figure
imagesc(BW0)
title(StrVf)
axis off
axis image
colormap jet
colorbar('Ticks',           [0, 2,      3,      4,      5,      6,      7,      8,      9,      10],...
         'TickLabels',{'Matrix',StrVfb1,StrVfb2,StrVfb3,StrVfb4,StrVfb5,StrVfb6,StrVfb7,StrVfb8,'Fibers'})
saveas(gcf,'Case1Vf.png')