clear
close all
O = imread('Mock_Sample1_new contrast_100.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
xlabel('Intensity')
ylabel('Counts')
%%
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.168;
lsize=round(5*(22/lpix));
ave=ones(lsize,lsize)/lsize^2;
BWave=conv2(BW,ave,'same');
%%
figure
imagesc(BWave,[0.2 0.8])
axis off
axis image
colormap jet
colorbar
saveas(gcf,'MockLocalVf.png')
%%
%save('MockVf.mat','BW','BW1','maskedImage')
load 'MockVf.mat'
%%
Vfcall=Vf
Vfb(1)=sum(BW1.*BW,'all')/sum(BW1,'all')
%%
StrVfb1=['V_{f} = ',num2str(Vfb(1),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;
BW1a(~BW)=0;
BW0(BW1a)=5;
%%
figure
imagesc(BW0)
axis off
axis image
colormap jet
colorbar('Ticks',[0,5,10],...
         'TickLabels',{'Matrix',StrVfb1,'Fibers'})
title(StrVf)
saveas(gcf,'MockMaskVf.png')
