%   Copyright 2018 Vedrana Andersen Dahl, vand@dtu.dk 
clear
close all
addpath functions

load bonesub
gray_th = 3*10^4;
div = [0.5 30.5,70.5,100.5,130.5,170.5 size(V,3)+0.5]; % layer divisions
midlayer = floor(0.5*(div(1:end-1) + div(2:end))); % middle slice of a layer

%% volume inspection -- use arrows
show_matvol(permute(V,[3,2,1]))
hold on, plot([1;size(V,2)]*ones(1,5),[1;1]*div(2:end-1),'r','LineWidth',2)
text(size(V,2)/2*ones(1,6),midlayer,num2cell(1:6),'Color','r','FontSize',20)
show_matvol(V)

%% still figures
x = 100;
figure, imagesc(squeeze(V(x,:,:))'); axis image, colormap gray
hold on, plot([1;size(V,2)]*ones(1,5),[1;1]*div(2:end-1),'r','LineWidth',2)
text(size(V,2)/2*ones(1,6),midlayer,num2cell(1:6),'Color','r','FontSize',20)
title(['Longitudinal slice, x=',num2str(x)]), xlabel('y'), ylabel('z')

figure
for i=1:numel(midlayer)
    subplot(2,3,i), imagesc(V(:,:,midlayer(i))), axis image, colormap gray
    title(['Region ' num2str(i), ', z=',num2str(midlayer(i))])
    xlabel('y'), ylabel('x')
end

%% binarization
B = medfilt3(medfilt3(V,[3,3,3])>gray_th,[3,3,3]);
show_matvol(B)

%% still figures
figure, imagesc(squeeze(B(x,:,:))'); axis image, colormap gray
hold on, plot([1;size(B,2)]*ones(1,5),[1;1]*div(2:end-1),'r','LineWidth',2)
text(size(V,2)/2*ones(1,6),midlayer,num2cell(1:6),'Color','r','FontSize',20)
title(['Longitudinal slice, x=',num2str(x)]), xlabel('y'), ylabel('z')

figure
for i=1:numel(midlayer)
    subplot(2,3,i), imagesc(B(:,:,midlayer(i))), axis image, colormap gray
    title(['Region ' num2str(i), ', z=',num2str(midlayer(i))])
    xlabel('y'), ylabel('x')
end

%% distance transform
distances = bwdist(~B) - bwdist(B) - B + 0.5; % distance from voxel center to interface
show_matvol(distances)

%% skeletonization of fibres and air 
% criterium: local maxima in any of the 3 planes for larger number of sceleton points)
conn = cat(3,zeros(3,3),[0 1 0; 1 1 1; 0 1 0],zeros(3,3));
skeleton = ((imregionalmax(distances,conn) | imregionalmax(distances,permute(conn,[2 3 1]))...
    | imregionalmax(distances,permute(conn,[3 1 2]))) & B)...
    |  (imregionalmin(distances,conn) | imregionalmin(distances,permute(conn,[2 3 1]))...
    | imregionalmin(distances,permute(conn,[3 1 2]))) & ~B;

% alternative with 3D local maxima, which yields less points
% conn = cat(3,[0 0 0; 0 1 0; 0 0 0],[0 1 0; 1 1 1; 0 1 0],[0 0 0; 0 1 0; 0 0 0]);
% skeleton = (imregionalmax(distances,conn) & B) | (imregionalmax(distances,conn) & ~B);

% cleaning
% skeleton = skeleton & abs(distances)>1;
skeleton([1,end],:,:) = false; skeleton(:,[1,end],:) = false; skeleton(:,:,[1,end]) = false;

%% visualizing single slice with a special colormap
z = 100; % and only a slice
[U,cmap,cbar] = color_skeleton(B(:,:,z),skeleton(:,:,z),distances(:,:,z));
figure, imagesc(U), colormap(cmap), axis image, title(['Slice z=',num2str(z)])
colorbar(cbar{:})

[U,cmap,cbar] = color_skeleton(squeeze(B(x,:,:))',squeeze(skeleton(x,:,:))',squeeze(distances(x,:,:))');
figure, imagesc(U), colormap(cmap), axis image, 
title(['Longitudinal slice, x=',num2str(x)]), xlabel('y'), ylabel('z')
colorbar(cbar{:})

%% visualizing multiple slices using special colormap
[U,cmap,cbar] = color_skeleton(B,skeleton,distances);
figure
for i=1:numel(midlayer)
    subplot(2,3,i), imagesc(U(:,:,midlayer(i)),[0,255]), axis image, colormap(cmap)
    title(['region ' num2str(i), ', z=',num2str(midlayer(i))])
    xlabel('y'), ylabel('x')
end

%% volume slicing using special colormap
show_matvol(U,[0,255]), colormap(cmap); % full volume
colorbar(cbar{:})

%% measuring porosity for a whole volume
volume_all = sum(B(:))/numel(B); % phase volume per block volume
area_all = compute_area(B)/numel(B); % phase surface per block volume

% measures from distance fields
figure
inside_all = fit_distribution(distances(skeleton & B),true);
a = gca; title([{'Thickness for a whole sample'};a.Title.String])
xlim([0 30]), xlabel('Radius (pixels)'), ylabel('Count')

figure
outside_all = fit_distribution(-distances(skeleton & ~B),true);
a = gca; title([{'Spacing for a whole sample'};a.Title.String])
xlim([0 30]), xlabel('Radius (pixels)'), ylabel('Count')

%% performing layered collection of measures

volume_layer = zeros(numel(midlayer),1);
area_layer = zeros(numel(midlayer),1);
inside_layer = zeros(numel(midlayer),5);
outside_layer = zeros(numel(midlayer),5);

fig_in = figure;
fig_out = figure;
for i = 1:numel(midlayer)
    range = ceil(div(i)):floor(div(i+1));
    B_layer = B(:,:,range);
    distances_layer = distances(:,:,range);
    skeleton_layer = skeleton(:,:,range);
    
    volume_layer(i) = sum(B_layer(:))/numel(B_layer);
    area_layer(i) = compute_area(B_layer)/numel(B_layer);
    
    % measures from distance fields
    figure(fig_in), subplot(2,3,i)
    inside_layer(i,:) = fit_distribution(distances_layer(skeleton_layer & B_layer),true);
    legend off, xlim([0 25])
    title(['Thickness for region ' num2str(i)])
    
    figure(fig_out), subplot(2,3,i)
    outside_layer(i,:) = fit_distribution(-distances_layer(skeleton_layer & ~B_layer),true);
    legend off, xlim([0 25])
    title(['Separation for region ' num2str(i)])
end

%% saving results
% fid = fopen('results.tex','w');
% fprintf(fid, 'Measure & Sample & Region 1 & Region 2 & Region 3 & Region 4 & Region 5 & Region 6 \\\\\\hline\n');
% fprintf(fid, 'Bone volume & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[volume_all; volume_layer]);
% fprintf(fid, 'Area & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[area_all; area_layer]);
% fprintf(fid, 'Area per bone volume & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[area_all; area_layer]./[volume_all; volume_layer]);
% fprintf(fid, 'Thickness, mode & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[inside_all(5); inside_layer(:,5)]);
% fprintf(fid, 'Thickness, median & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[inside_all(1); inside_layer(:,1)]);
% fprintf(fid, 'Thickness, mean & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[inside_all(3); inside_layer(:,3)]);
% fprintf(fid, 'Separation, mode & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[outside_all(5); outside_layer(:,5)]);
% fprintf(fid, 'Separation, median & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[outside_all(1); outside_layer(:,1)]);
% fprintf(fid, 'Separation, mean & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[outside_all(3); outside_layer(:,3)]);
% fclose(fid);
% 
% for savefig = [3,4,6,7,9,10,11,13,14,15,16]
%     set(savefig,'Units','Normalized','Position',[0.1 0.1 0.6 0.8])
%     saveas(savefig,['Figure_',num2str(savefig),'.eps'],'epsc');    
% end