Newer
Older
%% The code is written by Behnaz Pirzamanbein, bepi@dtu.dk last version 2020.08.04
%% summing over 10 pixels along normal direction of the smoothed surface
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
save_folder = 'Results\registration\';
if ~exist(save_folder,'dir')
mkdir(save_folder)
end
%% load the surfaces
samples = {'out_out_CD','out_out_MD','in_in_CD','in_in_MD'};
degree = {'0','45','90','180'};
% if you want to run the measure separately you can use the for loop
% for i = 1:4
% for j = 1:4
sample = samples{i};
dg = degree{j};
inner = load(['Results\layer\inner_converted_layer_',sample,'_',dg]);
outer = load(['Results\layer\outer_converted_layer_',sample,'_',dg]);
flag_layer = 'inner';
load(['Results\point\',flag_layer,'_points_',sample,'_',dg])
%% Contrast enhansed volume
switch answer
case 'No'
Vout = Vol;
case 'Yes'
ce = points.contrast_enhansed_range;
Vout = (Vol-ce(1))/(ce(2)-ce(1))*ce(2);
Vout(Vout < ce(1)) = ce(1);
Vout(Vout > ce(2)) = ce(2);
end
%% smooth the layers
p1 = mean(inner.surface_XYZ,4);
p2 = mean(outer.surface_XYZ,4);
Ba = regularization_matrix_open(size(p1,1),100,5000);
p1_smooth = [];
p2_smooth = [];
for m = 1:size(p1,2)
p1_smooth(:,m,:) = (squeeze(p1(:,m,:))'*Ba')';
p2_smooth(:,m,:) = (squeeze(p2(:,m,:))'*Ba')';
end
Bb = regularization_matrix_open(size(p1,2),1000,5000);
p1_s = [];
p2_s = [];
for n = 1:size(p1,1)
p1_s(n,:,:) = (squeeze(p1_smooth(n,:,:))'*Bb')';
p2_s(n,:,:) = (squeeze(p2_smooth(n,:,:))'*Bb')';
end
clear p1_smooth p2_smooth
%% sum up over 10 step in normal direction inward of the material
no_slice = points.slice_no(end)-points.slice_no(1)+1;
s1_conv = [];
s2_conv = [];
h = waitbar(0,'1','Name','Extracting layer ...');
for s = 1:size(p1,2)
slice_no = points.slice_no(1)+s-1;
% inner surface
range = 5:-1:-10;
normals = compute_normals(p1_s(:,s,1:2),1);
% grid along normals
p_X = p1_s(:,s,1) + normals(:,1)*range;
p_Y = p1_s(:,s,2) + normals(:,2)*range;
interp_method = 'nearest';
s1_conv(:,:,s) = interp2(Vout(:,:,slice_no),p_X,p_Y,interp_method);
% outer surface
range = 10:-1:-10;
normals = compute_normals(p2_s(:,s,1:2),1);
% grid along normals
p_X = p2_s(:,s,1) + normals(:,1)*range;
p_Y = p2_s(:,s,2) + normals(:,2)*range;
interp_method = 'nearest';
s2_conv(:,:,s) = interp2(Vout(:,:,slice_no),p_X,p_Y,interp_method);
waitbar(s/no_slice,h,sprintf('slice: %d',slice_no))
end
close(h)
s1 = squeeze(sum(s1_conv,2));
s2 = squeeze(sum(s2_conv,2));
if flag_save
save([save_folder,'registration_',sample,'_',dg],'s1','s2')
imwrite(rescale(s1),[save_folder,'registration_s1_',sample,'_',dg,'.png'])
imwrite(rescale(s2),[save_folder,'registration_s2_',sample,'_',dg,'.png'])
end
% end
% end
close all