-
Notifications
You must be signed in to change notification settings - Fork 1
/
clean_up_skeleton_component.m
116 lines (95 loc) · 5.37 KB
/
clean_up_skeleton_component.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
function [A_decimated, xyz_from_decimated_node_index] = ...
clean_up_skeleton_component(A_input, ...
xyz_from_input_node_index, ...
size_threshold, ...
smoothing_filter_width, ...
length_threshold, ...
do_visualize, ...
sampling_interval)
% Plot the input
if do_visualize ,
f = figure('Color', 'w') ;
ax = axes(f) ;
l0 = gplot3(ax, A_input, xyz_from_input_node_index, 'b') ; % blue
hold(ax, 'on') ;
drawnow
end
% Get a spanning tree, and the deleted edges
input_node_count = size(A_input, 1) ;
input_node_index_from_input_node_index = (1:input_node_count)' ;
degree = full(sum(A_input)) ;
degree_modified = degree ;
degree_modified(degree==0) = +inf ;
[~, min_degree_node_index] = min(degree_modified) ; % min degree is hopefully one, but sometimes not possible (e.g. if graph is just a loop)
dA_spanning = rooted_tree_from_connected_graph(A_input, min_degree_node_index) ;
A_spanning = max(dA_spanning, dA_spanning') ;
A_leftovers = A_input - A_spanning ; % All the edges that were deleted from A_input to make it into a tree
if do_visualize
l1 = gplot3(ax, A_spanning, xyz_from_input_node_index, 'Color', [0 0.5 0]) ; % green
drawnow
end
%%
[A_pruned, xyz_from_pruned_node_index, input_node_index_from_pruned_node_index] = ...
prune_tree_with_id_tracking(A_spanning, xyz_from_input_node_index, input_node_index_from_input_node_index, length_threshold, do_visualize) ;
if do_visualize
l2 = gplot3(ax, A_pruned, xyz_from_pruned_node_index, 'r') ; % red
drawnow
end
%%
pruned_node_count = length(input_node_index_from_pruned_node_index) ;
if pruned_node_count < size_threshold ,
% fprintf('Component with id %d fails to meet the size threshold (%d) after pruning, so discarding. (It contains %d nodes.)\n', ...
% component_id, ...
% size_threshold, ...
% pruned_node_count) ;
A_decimated = sparse([]) ;
xyz_from_decimated_node_index = zeros(0, 3) ;
return
end
% Decompose tree into chains
pruned_as_chains = chains_from_tree(A_pruned, input_node_index_from_pruned_node_index) ;
% Smooth the z coords of chain nodes
smoothed_xyz_from_pruned_node_index = smooth_chains(pruned_as_chains, xyz_from_pruned_node_index, smoothing_filter_width) ;
%xyz_smoothed = xyz_pruned ;
% Visualize the smoothed tree
if do_visualize
l3 = gplot3(ax, A_pruned, smoothed_xyz_from_pruned_node_index, 'Color', [1 2/3 0]) ; % orange
drawnow
end
% Decimate chains
G_decimated_as_chains_using_pruned_node_indices = decimate_chains(pruned_as_chains, smoothed_xyz_from_pruned_node_index, sampling_interval) ;
% Convert chains to edges
G_decimated_as_edges_using_pruned_node_indices = edges_from_chains(G_decimated_as_chains_using_pruned_node_indices) ;
% Defragment the node ids
[G_decimated_as_edges_using_decimated_node_indices, xyz_from_decimated_node_index, pruned_node_index_from_decimated_node_index] = ...
defragment_node_ids_in_edges(G_decimated_as_edges_using_pruned_node_indices, smoothed_xyz_from_pruned_node_index) ;
input_node_index_from_decimated_node_index = input_node_index_from_pruned_node_index(pruned_node_index_from_decimated_node_index) ;
% Convert edges to (sparse) adjacency
decimated_node_count = length(pruned_node_index_from_decimated_node_index) ;
A_decimated = undirected_adjacency_from_edges(G_decimated_as_edges_using_decimated_node_indices, decimated_node_count) ;
% Visualize the decimated tree
if do_visualize
l4 = gplot3(ax, A_decimated, xyz_from_decimated_node_index, 'Color', [0.6 0 0.8]) ; % purple
legend([l0 l1 l2 l3 l4], {'original', 'tree', 'pruned', 'smoothed', 'decimated'}) ;
axis(ax, 'vis3d') ;
%title(ax, sprintf('Component %d', component_id)) ;
xlabel(ax, 'x (um)') ;
ylabel(ax, 'y (um)') ;
zlabel(ax, 'z (um)') ;
drawnow
end
% Compute
is_in_decimated_from_input_node_index = false(input_node_count, 1) ;
is_in_decimated_from_input_node_index(input_node_index_from_decimated_node_index) = true ;
A_leftovers_from_decimated_node_index = A_leftovers(is_in_decimated_from_input_node_index, is_in_decimated_from_input_node_index) ;
% % Convert to a named tree
% digits_needed_for_index = floor(log10(max_component_id)) + 1 ;
% tree_name_template = sprintf('auto-cc-%%0%dd', digits_needed_for_index) ; % e.g. 'tree-%04d'
% tree_name = sprintf(tree_name_template, component_id) ;
% A_result = named_tree_from_undirected_graph(A_decimated, xyz_result, tree_name) ;
% % Compute the output file path
% tree_mat_file_name = sprintf('%s.mat', tree_name) ;
% tree_mat_file_path = fullfile(output_folder_path, tree_mat_file_name);
% Write full tree as a .mat file
%save_named_tree_as_mat(tree_mat_file_path, named_tree) ;
end