Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
% 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