-
Notifications
You must be signed in to change notification settings - Fork 1
/
DEMO_02_registration.m
119 lines (102 loc) · 4.76 KB
/
DEMO_02_registration.m
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
% This demo script will register distance maps to a reference shape.
% Registration is done in two steps: in step 1, all distance maps
% are registered to a chosen reference shape and in step 2, all distance
% maps are registered to the mean shape after step 1.
%
% This step requires elastix to be installed and recognised as an external
% command.
% NOTE: DEMO_01_alignment_and_distmap.m should have been run before
% running this script.
clear
dataFolder = fullfile('example_data');
%% Registration step 1: Register distance maps to reference scan
% The result of this step is the first estimate of the mean shape and the
% mapping from the chosen reference model to all other models.
filename.refDistmap = fullfile(dataFolder,'distmap','abs','01_abs.nii.gz');
filename.refModel = fullfile(dataFolder,'surface','aligned','01.stl');
% Set parameter files for Elastix: similarity transform followed by a
% b-spline transform with a grid spacing of 10x10x10 mm
parfile = {'parfiles/parameters_BSpline10.txt'};
results01 = register_distance_maps(...
fullfile(dataFolder,'distmap','abs'),...
filename.refDistmap,...
filename.refModel,...
parfile,...
'transformFolder',fullfile(dataFolder,'registration','step01','transform'),...
'resultFolder',fullfile(dataFolder,'registration','step01'));
%% Step 2: registration to the mean shape of step 1
% Make the mean shape from step 1 the reference shape for step 2 of
% registration. Make distance maps of the mean shape of step 1.
if exist(fullfile(dataFolder,'registration','step02'),'dir')~=7
mkdir(fullfile(dataFolder,'registration','step02'))
end
filename.refModel = fullfile(dataFolder,'registration','step02','ref_shape.stl');
filename.refDistmap_abs = fullfile(dataFolder,'registration','step02','ref_shape_abs.nii.gz');
filename.refDistmap_signed = fullfile(dataFolder,'registration','step02','ref_shape_signed.nii.gz');
copyfile(fullfile(dataFolder,'registration','step01','mean_shape.stl'),...
filename.refModel)
% Use the same spatial dimensions as the dimension of all the original
% distance maps.
dims = load(fullfile(dataFolder,'distmap','dimensions.mat'));
stl2distmap(...
filename.refModel,...
'distmap_signed',filename.refDistmap_signed,...
'distmap_abs',filename.refDistmap_abs,...
'padding',dims.padding,...
'spacex',dims.voxelsize(1),...
'spacey',dims.voxelsize(2),...
'spacez',dims.voxelsize(3));
% Register distance maps to the mean shape of step 1
results02 = register_distance_maps(...
fullfile(dataFolder,'distmap','abs'),...
filename.refDistmap_abs,...
filename.refModel,...
parfile,...
'transformFolder',fullfile(dataFolder,'registration','step02','transform'),...
'resultFolder',fullfile(dataFolder,'registration','step02'));
%% Visually inspect the result of registration
% Load the registration results
n_surf = size(results02.X,3);
% NOTE: The transformed reference surface should closely match the original
% model.
% The residual error is the average distance, across vertices, of the
% transformed reference surface to the original surface (calculated by
% interpolating the absolute distance map of the original surface at the
% vertices of the transformed reference.)
figure('Color','w','Name','Original and transformed reference')
for nr = 1 : n_surf
subplot(1,6,nr)
% Original (aligned) surface model
ho = plotSurface(fullfile(dataFolder,'surface','aligned',sprintf('%02d.stl',nr)),...
'EdgeColor','r','EdgeAlpha',0.2);
% Reference model transformed to target
TR = triangulation(results02.mean_shape.ConnectivityList, results02.X(:,:,nr));
hr = plotSurface(TR,'FaceColor','none','EdgeColor','k');
legend([ho hr],{'original','transformed reference'})
title(sprintf('Residual error: %.2f mm',mean(results02.res_dist(:,nr))))
end
% Set some axes properties
set(findobj(gcf,'Type','Axes'),'Clipping','off')
axis(findobj(gcf,'Type','Axes'),'off')
% Link the viewing angle and axis limits of all subplots.
hlink = linkprop(findobj(gcf,'Type','Axes'),{'XLim','YLim','ZLim','View'});
view(0,15)
%% Show the mean shape and the individual shapes
figure('Color','w');hold on
hm = plotSurface(results02.mean_shape,...
'FaceColor','none',...
'EdgeColor','k',...
'LineWidth',2);
view(0,15)
colors = linspecer(n_surf);
h = zeros(1,n_surf);
for nr = 1 : n_surf
TR = triangulation(results02.mean_shape.ConnectivityList,results02.X(:,:,nr));
h(nr) = plotSurface(TR,...
'EdgeColor',colors(nr,:),'FaceColor',colors(nr,:),...
'EdgeAlpha',0.5,'FaceAlpha',0.2);
txt{nr} = sprintf('shape%02d',nr);
end
legend([hm h],[{'mean shape'} txt]);
set(gca,'CLipping','off')
view(0,15)