Skip to content

Commit 9983437

Browse files
adamltayloradamltaylor
committed
Added option to stitch() to do a cold stitch
1 parent 43821f0 commit 9983437

File tree

4 files changed

+348
-280
lines changed

4 files changed

+348
-280
lines changed

‎stitching/+util/getcontolpixlocations.m

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,9 @@
1-
function [st,ed] = getcontolpixlocations(scopeloc, params, scopeparams)
1+
function [st,ed] = getcontolpixlocations(scopeloc, params)
22
gridix = scopeloc.gridix ;
33
loc = scopeloc.loc ;
44
dims = params.imagesize ;
5-
scopeparams1 = scopeparams(1) ;
6-
imsize_um = scopeparams1.imsize_um ;
7-
scopeparams1dims = scopeparams1.dims ;
5+
imsize_um = params.imsize_um ;
6+
scopeparams1dims = params.imagesize ;
87
targetoverlap_um = 4*[5 5 5]; %in um
98
N = params.Ndivs;
109
gridix_xyz = gridix(:,1:3) ;

‎stitching/compute_targets_from_z_matches.m

Lines changed: 214 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,214 @@
1+
function [targets_from_tile_index, cpg_k0_values_from_tile_index, tile_k_from_run_layer_index] = ...
2+
compute_targets_from_z_matches(baseline_targets_from_tile_index, ...
3+
do_cold_stitch, ...
4+
regpts, ...
5+
cpg_ij1s, ...
6+
default_cpg_k0_values, ...
7+
scopeloc, ...
8+
scopeparams, ...
9+
has_k_plus_1_tile_from_tile_index, ...
10+
tile_ijk_from_tile_index, ...
11+
curvemodel, ...
12+
baseline_affine_transform_from_tile_index, ...
13+
tile_ij1s, ...
14+
tileneighbors, ...
15+
params, ...
16+
do_apply_field_correction)
17+
18+
% Get some dimensions
19+
tile_count = length(scopeloc.filepath) ;
20+
cpg_i_count = params.Ndivs+1 ; % Number of i/x values in the per-tile control point grid (traditionally 5)
21+
cpg_j_count = params.Ndivs+1 ; % Number of j/y values in the per-tile control point grid (traditionally 5)
22+
cpg_ij_count = cpg_i_count * cpg_j_count ; % Number of points in each k/z layer in the per-tile control point grid (traditionally 25)
23+
24+
% The first pass of the z planes of the control points, for each tile
25+
% If a cold stitch, these will also be the final pass values.
26+
targets_from_tile_index = baseline_targets_from_tile_index ;
27+
baseline_cpg_k0_values_from_tile_index = repmat(default_cpg_k0_values', [1 tile_count]) ;
28+
tile_k_from_tile_index = tile_ijk_from_tile_index(:,3) ; % column
29+
tile_k_from_layer_index = unique(tile_k_from_tile_index) ; % layer of tiles, that is
30+
tile_k_from_run_layer_index = tile_k_from_layer_index(1:end-1)' ; % a "run layer" is a layer that will actually be run
31+
cpg_k0_values_from_tile_index = baseline_cpg_k0_values_from_tile_index ;
32+
33+
% If a cold stitch, exit early
34+
if do_cold_stitch ,
35+
return
36+
end
37+
38+
% Collect some information about the point correspondences for each tile
39+
has_zero_z_face_matches_from_tile_index = false(tile_count,1) ;
40+
match_statistics = nan(tile_count, 8) ;
41+
for tile_index = 1 : tile_count ,
42+
% pix stats
43+
this_tile_regpts = regpts{tile_index} ;
44+
match_coords = this_tile_regpts.X ; % match_count x 3, the coordinates of matched landmarks in this tile
45+
if isempty(match_coords) ,
46+
has_zero_z_face_matches_from_tile_index(tile_index) = true ;
47+
continue
48+
end
49+
match_coords_in_neighbor = this_tile_regpts.Y ; % match_count x 3, the coordinates of matched landmarks in the z+1 tile, in same order as layer
50+
match_z_in_both_tiles = [ match_coords(:,3) match_coords_in_neighbor(:,3) ] ; % match_count x 2
51+
median_match_z_in_both_tiles = round(median(match_z_in_both_tiles,1)) ; % 1x2
52+
min_match_z_in_both_tiles = round(min(match_z_in_both_tiles,[],1)) ; % 1x2
53+
max_match_z_in_both_tiles = round(max(match_z_in_both_tiles,[],1)) ; % 1x2
54+
tile_index_of_z_plus_1_tile = this_tile_regpts.neigs(4) ;
55+
match_statistics(tile_index,:) = ...
56+
[ tile_index tile_index_of_z_plus_1_tile ...
57+
median_match_z_in_both_tiles min_match_z_in_both_tiles max_match_z_in_both_tiles ] ;
58+
end
59+
60+
% Sort out which we'll run with which thing
61+
% The tiles we run with optimpertile are "anchors".
62+
% The tiles we run with nomatchoptim are "floaters".
63+
is_floater_from_tile_index = has_k_plus_1_tile_from_tile_index & has_zero_z_face_matches_from_tile_index ;
64+
is_anchor_from_tile_index = ~is_floater_from_tile_index ;
65+
66+
% Make an interpolator for each z layer, use it to shift the control point
67+
% targets for the tiles in the layer, and the z+1 layer
68+
run_layer_count = length(tile_k_from_run_layer_index) ;
69+
for run_layer_index = 1 : run_layer_count ,
70+
% Sort out which tiles are in this layer
71+
tile_k = tile_k_from_run_layer_index(run_layer_index) ;
72+
fprintf(' Layer %d of %d, tile k/z = %d\n', run_layer_index, run_layer_count, tile_k);
73+
is_in_this_layer_from_tile_index = (tile_k_from_tile_index'==tile_k) ;
74+
tile_index_from_layer_tile_index = find(is_in_this_layer_from_tile_index);
75+
if isempty(tile_index_from_layer_tile_index) ,
76+
fprintf('No tiles found in layer with tile k/z = %d!!\n', tile_k) ;
77+
continue
78+
end
79+
80+
% get interpolants based on paired descriptors
81+
[Fx_layer, Fy_layer, Fz_layer, Fx_next_layer, Fy_next_layer, Fz_next_layer, XYZ_original, XYZ_neighbor_original, outliers] =...
82+
util.getInterpolants(tile_index_from_layer_tile_index, ...
83+
regpts, ...
84+
baseline_affine_transform_from_tile_index, ...
85+
tile_ij1s, ...
86+
params, ...
87+
curvemodel, ...
88+
do_apply_field_correction) ;
89+
90+
% Show some debugging output if called for
91+
if params.debug ,
92+
vector_field_3d_debug_script(scopeloc, scopeparams, params, XYZ_original, XYZ_neighbor_original, outliers) ;
93+
end
94+
95+
% If too few matched landmarks, don't proceed with this layer
96+
if isempty(Fx_layer) || size(Fx_layer.Points,1) < 10 ,
97+
fprintf(' Layer with k/z = %d has too few matches to proceed.\n', tile_k) ;
98+
continue
99+
end
100+
101+
% Print the number of matches in this layer
102+
layer_used_match_count = size(Fx_layer.Points,1) ;
103+
fprintf(' Layer with k/z = %d total used matches: %d\n', tile_k, layer_used_match_count) ;
104+
105+
% Go through all the tiles which either have no z+1 tile or for which we have
106+
% nonzero z-face matches, and optimize the CPG for each. Then, compute the
107+
% target values for the moved z-planes of the CPG, using the interpolators to
108+
% try to make the matched pairs each go to the same point.
109+
is_anchor_from_layer_tile_index = is_anchor_from_tile_index(tile_index_from_layer_tile_index) ;
110+
tile_index_from_layer_anchor_index = tile_index_from_layer_tile_index(is_anchor_from_layer_tile_index) ;
111+
for tile_index = tile_index_from_layer_anchor_index ,
112+
neighbor_tile_index = tileneighbors(tile_index, 7) ; % the z+1 tile
113+
114+
if do_apply_field_correction ,
115+
this_tile_curve_model = curvemodel(:,:,tile_index) ;
116+
field_corrected_cpg_ij1s = util.fcshift(this_tile_curve_model, order, tile_ij1s, tile_shape_ijk, cpg_ij1s) ;
117+
else
118+
field_corrected_cpg_ij1s = cpg_ij1s ;
119+
end
120+
field_corrected_cpg_ij0s = field_corrected_cpg_ij1s - 1 ;
121+
122+
[targets_at_cpg_k_indices_3_and_4, ...
123+
other_tile_targets_at_cpg_k_indices_1_and_2, ...
124+
new_cpg_k0_values, ...
125+
k_plus_1_tile_cpg_k0_values] = ...
126+
optimpertile(tile_index, ...
127+
params, ...
128+
neighbor_tile_index, ...
129+
baseline_affine_transform_from_tile_index, ...
130+
match_statistics, ...
131+
cpg_k0_values_from_tile_index, ...
132+
field_corrected_cpg_ij0s, ...
133+
Fx_layer, ...
134+
Fy_layer, ...
135+
Fz_layer, ...
136+
Fx_next_layer, ...
137+
Fy_next_layer, ...
138+
Fz_next_layer, ...
139+
default_cpg_k0_values) ;
140+
cpg_k0_values_from_tile_index(:,tile_index) = new_cpg_k0_values ;
141+
if ~isnan(neighbor_tile_index) ,
142+
cpg_k0_values_from_tile_index(:,neighbor_tile_index) = k_plus_1_tile_cpg_k0_values ;
143+
end
144+
targets_from_tile_index(2*cpg_ij_count+1:end,:,tile_index) = targets_at_cpg_k_indices_3_and_4 ;
145+
if ~isempty(other_tile_targets_at_cpg_k_indices_1_and_2) ,
146+
targets_from_tile_index(1:2*cpg_ij_count,:,neighbor_tile_index) = other_tile_targets_at_cpg_k_indices_1_and_2 ;
147+
end
148+
end % for tile_index = ...
149+
150+
% Get the tile ijk indices of the anchors and the floaters
151+
is_floater_from_tile_within_layer_index = is_floater_from_tile_index(tile_index_from_layer_tile_index) ;
152+
tile_index_from_layer_floater_index = tile_index_from_layer_tile_index(is_floater_from_tile_within_layer_index) ;
153+
layer_floater_count = length(tile_index_from_layer_floater_index) ;
154+
tile_ijk_from_layer_anchor_index = tile_ijk_from_tile_index(tile_index_from_layer_anchor_index,:) ;
155+
tile_ijk_from_layer_floater_index = tile_ijk_from_tile_index(tile_index_from_layer_floater_index,:) ;
156+
157+
% Build an index of the closest anchor tile to each floater tile
158+
nearest_layer_anchor_index_from_layer_floater_index = knnsearch(tile_ijk_from_layer_anchor_index,tile_ijk_from_layer_floater_index,'K',1) ;
159+
160+
% Build an index of all the anchor tiles in the "18-neighborhood" of each
161+
% floater tile. If that set is empty for any floater tile, augment with
162+
% whatver the nearest anchor tile is.
163+
nearby_layer_anchor_indices_from_layer_floater_index = rangesearch(tile_ijk_from_layer_anchor_index,tile_ijk_from_layer_floater_index,sqrt(2)) ;
164+
% "nearby" here seems to mean not the 6-neighborhood nor the 26-neighborhood, but
165+
% the 18-neighborhood: All the neighbors that differ by 1 in at most 2
166+
% dimensions. So the whole rubick's cube (3^3==27), minus self (-1), minus the
167+
% corners (-8). So 27-1-8==18, the '18-neighborhood'.
168+
for layer_floater_index = 1 : layer_floater_count ,
169+
layer_anchor_index_from_nearby_tile_index = nearby_layer_anchor_indices_from_layer_floater_index{layer_floater_index} ;
170+
if isempty(layer_anchor_index_from_nearby_tile_index) ,
171+
nearby_layer_anchor_indices_from_layer_floater_index{layer_floater_index} = ...
172+
nearest_layer_anchor_index_from_layer_floater_index(layer_floater_index);
173+
end
174+
end
175+
176+
% For each floater tile, optimize its CPG and recompute targets, if needed
177+
for layer_floater_index = 1 : layer_floater_count ,
178+
% Get the tile index, and the z+1 tile index
179+
tile_index = tile_index_from_layer_floater_index(layer_floater_index) ;
180+
neighbor_tile_index = tileneighbors(tile_index,7) ; % the z/k+1 tile
181+
182+
% Get a fallback value for the CPG k0/z levels, by taking the median of the the
183+
% CPG k0 values for the nearby anchor tiles.
184+
layer_anchor_index_from_nearby_tile_index = nearby_layer_anchor_indices_from_layer_floater_index{layer_floater_index} ;
185+
tile_index_from_nearby_tile_index = tile_index_from_layer_anchor_index(layer_anchor_index_from_nearby_tile_index) ;
186+
cpg_k0_values_from_nearby_tile_index = cpg_k0_values_from_tile_index(:,tile_index_from_nearby_tile_index) ;
187+
local_default_cpg_k0_values = round(median(cpg_k0_values_from_nearby_tile_index,2))' ;
188+
189+
[targets_at_cpg_k_indices_3_and_4, ...
190+
other_tile_targets_at_cpg_k_indices_1_and_2, ...
191+
new_cpg_k0_values, ...
192+
k_plus_1_tile_cpg_k0_values] = ...
193+
nomatchoptim(tile_index, ...
194+
params, ...
195+
neighbor_tile_index, ...
196+
baseline_affine_transform_from_tile_index, ...
197+
match_statistics, ...
198+
cpg_k0_values_from_tile_index(:,tile_index), ...
199+
cpg_k0_values_from_tile_index(:,neighbor_tile_index), ...
200+
field_corrected_cpg_ij0s, ...
201+
Fx_layer, ...
202+
Fy_layer, ...
203+
Fz_layer, ...
204+
Fx_next_layer, ...
205+
Fy_next_layer, ...
206+
Fz_next_layer, ...
207+
local_default_cpg_k0_values) ;
208+
cpg_k0_values_from_tile_index(:,tile_index) = new_cpg_k0_values ;
209+
cpg_k0_values_from_tile_index(:,neighbor_tile_index) = k_plus_1_tile_cpg_k0_values ;
210+
targets_from_tile_index(2*cpg_ij_count+1:end,:,tile_index) = targets_at_cpg_k_indices_3_and_4 ;
211+
targets_from_tile_index(1:2*cpg_ij_count,:,neighbor_tile_index) = other_tile_targets_at_cpg_k_indices_1_and_2 ;
212+
end
213+
end % for loop over layers
214+
end

0 commit comments

Comments
 (0)