Skip to content

Commit c8414c6

Browse files
committed
Added code for doing line-shifting in-place
1 parent 3945e2f commit c8414c6

40 files changed

+607
-140
lines changed

.gitignore

+1
Original file line numberDiff line numberDiff line change
@@ -8,4 +8,5 @@ parfor_progress.txt
88
*.stdouterr
99
foreground-classifiers/data/
1010
foreground-classifiers/2021-03-17/
11+
RESTORED/
1112

build_tile_index_dirwalk_callback.m

+2-1
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,8 @@
1313
[~,file_base_name] = fileparts2(current_folder_relative_path) ;
1414
acquisition_file_name = [file_base_name '-ngc.acquisition'] ;
1515
acquisition_file_path = fullfile(root_folder_absolute_path, current_folder_relative_path, acquisition_file_name) ;
16-
if exist(acquisition_file_path, 'file') ,
16+
if ismember(acquisition_file_name, file_names) ,
17+
%if exist(acquisition_file_path, 'file') , % don't hit the FS again
1718
[ijk1, xyz] = read_lattice_position_from_acquisition_file(acquisition_file_path) ;
1819
updated_state = struct() ;
1920
updated_state.ijk1_from_tile_index = [state.ijk1_from_tile_index ; ijk1] ;
File renamed without changes.
File renamed without changes.
File renamed without changes.

copy_raw_tile_neighborhood.m

+4-4
Original file line numberDiff line numberDiff line change
@@ -27,19 +27,19 @@ function copy_raw_tile_neighborhood(hood_raw_tile_root_folder_path, ...
2727
center_tile_relative_path) ; %#ok<ASGLU>
2828
tile_count = length(relative_path_from_tile_index) ;
2929

30-
% For each relative path, make a link
30+
% For each relative path, make a copy
3131
for tile_index = 1 : tile_count ,
3232
tile_relative_path = relative_path_from_tile_index{tile_index} ;
3333
destination_path = fullfile(hood_raw_tile_root_folder_path, tile_relative_path) ;
3434
source_path = fullfile(raw_tile_root_folder_path, tile_relative_path) ;
3535
destination_parent_path = fileparts(destination_path) ;
3636
ensure_folder_exists(destination_parent_path) ;
37-
system_from_list_with_error_handling({'ln', '-s', '-T', source_path, destination_path}) ;
37+
system_from_list_with_error_handling({'cp', '-R', '-T', source_path, destination_path}) ;
3838
end
3939

40-
% Make links for the metadata, etc
40+
% Make copies for the metadata, etc
4141
file_name = 'sample-metadata.txt' ;
4242
source_path = fullfile(raw_tile_root_folder_path, file_name) ;
4343
destination_path = fullfile(hood_raw_tile_root_folder_path, file_name) ;
44-
system_from_list_with_error_handling({'ln', '-s', '-T', source_path, destination_path}) ;
44+
system_from_list_with_error_handling({'cp', '-T', source_path, destination_path}) ;
4545
end

copy_raw_tile_neighborhood_as_links.m

+45
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
function copy_raw_tile_neighborhood_as_links(hood_raw_tile_root_folder_path, ...
2+
raw_tile_root_folder_path, ...
3+
center_tile_relative_path, ...
4+
sample_memo_folder_path, ...
5+
do_force_computation)
6+
% Build the full tile index
7+
raw_tile_index = compute_or_read_from_memo(sample_memo_folder_path, ...
8+
'raw-tile-index', ...
9+
@()(build_raw_tile_index(raw_tile_root_folder_path)), ...
10+
do_force_computation) ;
11+
full_tile_index_from_full_tile_ijk1 = raw_tile_index.tile_index_from_tile_ijk1 ;
12+
full_tile_ijk1_from_full_tile_index = raw_tile_index.ijk1_from_tile_index ;
13+
xyz_from_full_tile_index = raw_tile_index.xyz_from_tile_index ;
14+
relative_path_from_full_tile_index = raw_tile_index.relative_path_from_tile_index ;
15+
%full_tile_grid_shape = size(full_tile_index_from_full_tile_ijk1)
16+
%full_tile_count = length(relative_path_from_full_tile_index)
17+
18+
% Get out just the tiles around the center_tile
19+
[tile_index_from_tile_ijk1, ...
20+
tile_ijk1_from_tile_index, ...
21+
xyz_from_tile_index, ...
22+
relative_path_from_tile_index] = ...
23+
extract_tiles_near_tile(full_tile_index_from_full_tile_ijk1, ...
24+
full_tile_ijk1_from_full_tile_index, ...
25+
xyz_from_full_tile_index, ...
26+
relative_path_from_full_tile_index, ...
27+
center_tile_relative_path) ; %#ok<ASGLU>
28+
tile_count = length(relative_path_from_tile_index) ;
29+
30+
% For each relative path, make a link
31+
for tile_index = 1 : tile_count ,
32+
tile_relative_path = relative_path_from_tile_index{tile_index} ;
33+
destination_path = fullfile(hood_raw_tile_root_folder_path, tile_relative_path) ;
34+
source_path = fullfile(raw_tile_root_folder_path, tile_relative_path) ;
35+
destination_parent_path = fileparts(destination_path) ;
36+
ensure_folder_exists(destination_parent_path) ;
37+
system_from_list_with_error_handling({'ln', '-s', '-T', source_path, destination_path}) ;
38+
end
39+
40+
% Make links for the metadata, etc
41+
file_name = 'sample-metadata.txt' ;
42+
source_path = fullfile(raw_tile_root_folder_path, file_name) ;
43+
destination_path = fullfile(hood_raw_tile_root_folder_path, file_name) ;
44+
system_from_list_with_error_handling({'ln', '-s', '-T', source_path, destination_path}) ;
45+
end

get_tile_flip_state.m

+26
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
function flip_metadata = get_tile_flip_state(tile_path, sample_metadata)
2+
if sample_metadata.sample_metadata_version<=1 ,
3+
flip_metadata = struct() ;
4+
flip_metadata.is_x_flipped = sample_metadata.is_x_flipped ;
5+
flip_metadata.is_y_flipped = sample_metadata.is_y_flipped ;
6+
elseif sample_metadata.sample_metadata_version==2 ,
7+
xlineshift_file_path = fullfile(tile_path, 'Xlineshift.txt') ;
8+
if exist(xlineshift_file_path, 'file') ,
9+
% Tile is already line-shifted, which means it's also flipped on both axes
10+
flip_metadata = struct() ;
11+
flip_metadata.is_x_flipped = true ;
12+
flip_metadata.is_y_flipped = true ;
13+
else
14+
% Tile has not been line-shifted, so use the value from the sample metadata
15+
flip_metadata = struct() ;
16+
flip_metadata.is_x_flipped = sample_metadata.is_x_flipped ;
17+
flip_metadata.is_y_flipped = sample_metadata.is_y_flipped ;
18+
end
19+
else
20+
tile_metadata = read_tile_metadata(tile_path) ;
21+
flip_metadata = struct() ;
22+
flip_metadata.is_x_flipped = tile_metadata.is_x_flipped ;
23+
flip_metadata.is_y_flipped = tile_metadata.is_y_flipped ;
24+
end
25+
end
26+
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
function advance_state_for_fix_line_shift_in_place(input_root_folder, tile_relative_path, do_run_in_debug_mode, ...
2+
min_shift, max_shift, ...
3+
tif_file_name_from_channel_index, original_tif_file_name_from_channel_index, shifted_tif_file_name_from_channel_index, ...
4+
neuron_channel_index0, original_is_x_flipped, original_is_y_flipped, ...
5+
state_index)
6+
input_folder_path = fullfile(input_root_folder, tile_relative_path) ;
7+
8+
function rename(old_file_name, new_file_name)
9+
tif_file_path = fullfile(input_folder_path, old_file_name) ;
10+
original_tif_file_path = fullfile(input_folder_path, new_file_name) ;
11+
system_from_list_with_error_handling({'mv', '-T', tif_file_path, original_tif_file_path}) ;
12+
end
13+
14+
function delete_file(file_name)
15+
file_path = fullfile(input_folder_path, file_name) ;
16+
system_from_list_with_error_handling({'rm', '-f', file_path}) ;
17+
end
18+
19+
if state_index==1 ,
20+
% Rename ddddd-ngc.0.tif to dddd-ngc.0.original.tif
21+
cellfun(@rename, tif_file_name_from_channel_index, original_tif_file_name_from_channel_index) ;
22+
elseif state_index==2 ,
23+
% Run the line-shift code and output Xlineshift.txt
24+
neuron_channel_index = neuron_channel_index0 + 1 ;
25+
neurons_channel_tif_file_name = original_tif_file_name_from_channel_index{neuron_channel_index} ;
26+
write_xlineshift_txt_in_place(input_folder_path, neurons_channel_tif_file_name, min_shift, max_shift, do_run_in_debug_mode) ;
27+
elseif state_index==3 ,
28+
% Generate the line-shifted and (possibly) flipped stacks, and the
29+
% tile-metadata.txt file
30+
write_line_shifted_stacks_and_tile_metadata_in_place(input_folder_path, ...
31+
original_tif_file_name_from_channel_index, ...
32+
shifted_tif_file_name_from_channel_index, ...
33+
original_is_x_flipped, ...
34+
original_is_y_flipped)
35+
elseif state_index==4 ,
36+
% Rename ddddd-ngc.0.shifted.tif to dddd-ngc.0.tif
37+
cellfun(@rename, shifted_tif_file_name_from_channel_index, tif_file_name_from_channel_index) ;
38+
elseif state_index==5 ,
39+
% Delete the original stacks
40+
cellfun(@delete_file, original_tif_file_name_from_channel_index) ;
41+
% Delete Xlineshift.txt
42+
delete_file('Xlineshift.txt') ;
43+
elseif state_index==6 ,
44+
% this is the final state, nothing to do
45+
else
46+
error('Illegal state index: %g', state_index) ;
47+
end
48+
end
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
function [result, state_count, tif_file_names, original_tif_file_names, shifted_tif_file_names] = ...
2+
determine_fix_line_shift_tile_state(input_root_folder, tile_relative_path, channel_count)
3+
% Construct absolute paths to input folder for this tile
4+
input_folder_path = fullfile(input_root_folder, tile_relative_path) ;
5+
[~, tile_base_name] = fileparts2(tile_relative_path) ;
6+
7+
% Generate the names of the tif files
8+
zero_based_channel_index_from_channel_index = (1:channel_count)-1 ;
9+
tif_file_names = ...
10+
arrayfun(@(index)([tile_base_name '-ngc.' num2str(index) '.tif']), ...
11+
zero_based_channel_index_from_channel_index, ...
12+
'UniformOutput', false) ;
13+
original_tif_file_names = ...
14+
arrayfun(@(index)([tile_base_name '-ngc.' num2str(index) '.original.tif']), ...
15+
zero_based_channel_index_from_channel_index, ...
16+
'UniformOutput', false) ;
17+
shifted_tif_file_names = ...
18+
arrayfun(@(index)([tile_base_name '-ngc.' num2str(index) '.shifted.tif']), ...
19+
zero_based_channel_index_from_channel_index, ...
20+
'UniformOutput', false) ;
21+
22+
% Generate a list of all the files we need to worry about
23+
name_from_file_index = [ tif_file_names, original_tif_file_names, {'Xlineshift.txt'}, shifted_tif_file_names, {'tile-metadata.txt'} ] ;
24+
25+
% Need to build a state_count x file_count matrix, with
26+
% +1 for all files that must be present at that state
27+
% -1 for all files that must be absent at that state
28+
% 0 for all files that might be present or might be absent at that state
29+
file_count = length(name_from_file_index) ;
30+
state_count = 6 ;
31+
file_sign_from_state_index_from_file_index = zeros(state_count, file_count) ;
32+
% The initial state, with ddddd-ngc.0.tif, ddddd-ngc.1.tif, and nothing else
33+
file_sign_from_state_index_from_file_index(1,:) = ...
34+
[ repmat(+1, [1 channel_count]) repmat(-1, [1 channel_count]) -1 repmat(-1, [1 channel_count]) -1 ] ; %#ok<REPMAT>
35+
file_sign_from_state_index_from_file_index(2,:) = ...
36+
[ repmat(-1, [1 channel_count]) repmat(+1, [1 channel_count]) -1 repmat(-1, [1 channel_count]) -1 ] ; %#ok<REPMAT>
37+
file_sign_from_state_index_from_file_index(3,:) = ...
38+
[ repmat(-1, [1 channel_count]) repmat(+1, [1 channel_count]) +1 repmat(-1, [1 channel_count]) -1 ] ; %#ok<REPMAT>
39+
file_sign_from_state_index_from_file_index(4,:) = ...
40+
[ repmat(-1, [1 channel_count]) repmat(+1, [1 channel_count]) +1 repmat(+1, [1 channel_count]) +1 ] ; %#ok<REPMAT>
41+
file_sign_from_state_index_from_file_index(5,:) = ...
42+
[ repmat(+1, [1 channel_count]) repmat(+1, [1 channel_count]) +1 repmat(-1, [1 channel_count]) +1 ] ; %#ok<REPMAT>
43+
file_sign_from_state_index_from_file_index(6,:) = ...
44+
[ repmat(+1, [1 channel_count]) repmat(-1, [1 channel_count]) -1 repmat(-1, [1 channel_count]) +1 ] ; %#ok<REPMAT>
45+
46+
% Derive the matrices that indicate which files must exist vs cant exist
47+
must_exist_from_state_index_from_file_index = ( file_sign_from_state_index_from_file_index > 0 ) ;
48+
cant_exist_from_state_index_from_file_index = ( file_sign_from_state_index_from_file_index < 0 ) ;
49+
50+
% Determine which files exist
51+
does_exist_from_file_index = cellfun(@(file_name)(exist(fullfile(input_folder_path, file_name), 'file')), name_from_file_index) ;
52+
is_missing_from_file_index = ~does_exist_from_file_index ;
53+
54+
% Determine what states we could be in
55+
is_possible_from_state_index = ...
56+
all(~must_exist_from_state_index_from_file_index | does_exist_from_file_index,2) & ...
57+
all(~cant_exist_from_state_index_from_file_index | is_missing_from_file_index,2) ;
58+
59+
% If can't determine a unique state, just error out---don't try to be clever
60+
possible_state_count = sum(is_possible_from_state_index) ;
61+
if possible_state_count == 0 ,
62+
error('Can''t determine state for tile %s. All possible states ruled out.', tile_relative_path) ;
63+
elseif possible_state_count == 1 ,
64+
% no nothing, this is what we want
65+
else
66+
state_index_from_possible_state_index = find(is_possible_from_state_index) ;
67+
error('Can''t determine unique state for tile %s. Possible states are: %s', tile_relative_path, mat2str(state_index_from_possible_state_index)) ;
68+
end
69+
70+
% Extract the current state index
71+
result = find(is_possible_from_state_index) ;
72+
end

line-fix/determine_line_shift.m

+11
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
function shift = determine_line_shift(raw_stack, min_shift, max_shift, do_run_in_debug_mode)
2+
% binarize it to eliminate spatial non-uniformity bias
3+
serial_stack = double(raw_stack(:)) ;
4+
middle_value = mean(quantile(serial_stack,[0.0001 0.9999])) ;
5+
stack = (raw_stack>middle_value) ;
6+
[shift, shift_float] = determine_line_shift_core(stack, min_shift, max_shift, false, do_run_in_debug_mode) ;
7+
% check if shift is closer to halfway. 0.4<|shift-round(shift)|<0.6
8+
if abs(abs(round(shift_float,2)-round(shift_float,0))-0.5) < 0.1 ,
9+
[shift, shift_float] = determine_line_shift_core(stack, min_shift, max_shift, true, do_run_in_debug_mode) ; %#ok<ASGLU>
10+
end
11+
end

find_shift.m line-fix/determine_line_shift_core.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function [rounded_shift, shift] = find_shift(stack, min_shift, max_shift, do_3d, do_run_in_debug_mode)
1+
function [rounded_shift, shift] = determine_line_shift_core(stack, min_shift, max_shift, do_3d, do_run_in_debug_mode)
22
if do_3d ,
33
im1 = double(stack(1:2:end,:,:)) ;
44
im2 = double(stack(2:2:end,:,:)) ;

fix_line_shift.m line-fix/fix_line_shift.m

+2-9
Original file line numberDiff line numberDiff line change
@@ -40,15 +40,8 @@ function fix_line_shift(input_root_folder, tile_relative_path, output_root_folde
4040
neurons_channel_tif_file_path = fullfile(input_folder_path, neurons_channel_tif_file_name) ;
4141
raw_stack = read_16bit_grayscale_tif(neurons_channel_tif_file_path) ;
4242

43-
% binarize it to eliminate spatial non-uniformity bias
44-
serial_stack = double(raw_stack(:)) ;
45-
middle_value = mean(quantile(serial_stack,[0.0001 0.9999])) ;
46-
stack = (raw_stack>middle_value) ;
47-
[shift, shift_float] = find_shift(stack, min_shift, max_shift, false, do_run_in_debug_mode) ;
48-
% check if shift is closer to halfway. 0.4<|shift-round(shift)|<0.6
49-
if abs(abs(round(shift_float,2)-round(shift_float,0))-0.5) < 0.1 ,
50-
[shift, shift_float] = find_shift(stack, min_shift, max_shift, true, do_run_in_debug_mode) ; %#ok<ASGLU>
51-
end
43+
% Determine the optimal line-shift
44+
shift = determine_line_shift(raw_stack, min_shift, max_shift, do_run_in_debug_mode) ;
5245

5346
% Make sure the output folder exists
5447
ensure_folder_exists(output_folder_path) ;

line-fix/fix_line_shift_in_place.m

+31
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
function fix_line_shift_in_place(input_root_folder, tile_relative_path, original_sample_metadata, do_run_in_debug_mode)
2+
% Deal with args
3+
if ~exist('do_run_in_debug_mode', 'var') || isempty(do_run_in_debug_mode) ,
4+
do_run_in_debug_mode = false ;
5+
end
6+
7+
% Set some internal parmeters
8+
channel_count = 2 ;
9+
min_shift = -9 ;
10+
max_shift = +9 ;
11+
12+
% Figure out which channel we'll use to compute the line shift
13+
neuron_channel_index0 = original_sample_metadata.neuron_channel_index ;
14+
original_is_x_flipped = original_sample_metadata.is_x_flipped ;
15+
original_is_y_flipped = original_sample_metadata.is_y_flipped ;
16+
17+
% Determine the current state by probing the file system
18+
[initial_state_index, state_count, tif_file_names, original_tif_file_names, shifted_tif_file_names] = ...
19+
determine_fix_line_shift_tile_state(input_root_folder, tile_relative_path, channel_count) ;
20+
21+
% Try to advance from the current state to the done state
22+
for state_index = initial_state_index : (state_count-1) ,
23+
advance_state_for_fix_line_shift_in_place(input_root_folder, tile_relative_path, do_run_in_debug_mode, ...
24+
min_shift, max_shift, ...
25+
tif_file_names, original_tif_file_names, shifted_tif_file_names, ...
26+
neuron_channel_index0, original_is_x_flipped, original_is_y_flipped, ...
27+
state_index) ;
28+
end
29+
% If get here without error, that must have worked
30+
end
31+
File renamed without changes.
File renamed without changes.
File renamed without changes.
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
function result = post_line_fix_sample_metadata_from_original_sample_metadata(original_sample_metadata)
2+
% Drop the is_*_flipped fields, because now those are specified for each tile.
3+
result = rmfield(original_sample_metadata, {'is_x_flipped', 'is_y_flipped'}) ;
4+
end
File renamed without changes.
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
function write_line_shifted_stacks_and_tile_metadata_in_place(tile_folder_path, ...
2+
original_tif_file_name_from_channel_index, ...
3+
shifted_tif_file_name_from_channel_index, ...
4+
original_is_x_flipped, ...
5+
original_is_y_flipped)
6+
% Generate the line-shifted and (possibly) flipped stacks, and the
7+
% tile-metadata.txt file.
8+
xlineshift_file_path = fullfile(tile_folder_path, 'Xlineshift.txt') ;
9+
shift = read_xlineshift_file(xlineshift_file_path) ;
10+
11+
% Write the line-shifted tif stacks
12+
channel_count = length(original_tif_file_name_from_channel_index) ;
13+
for channel_index = 1 : channel_count ,
14+
input_tif_file_name = original_tif_file_name_from_channel_index{channel_index} ;
15+
input_tif_file_path = fullfile(tile_folder_path, input_tif_file_name) ;
16+
output_tif_file_name = shifted_tif_file_name_from_channel_index{channel_index} ;
17+
output_tif_file_path = fullfile(tile_folder_path, output_tif_file_name) ;
18+
stack = read_16bit_grayscale_tif(input_tif_file_path) ;
19+
stack(2:2:end,:,:) = circshift(stack(2:2:end,:,:), shift, 2) ;
20+
% shift every other y level, starting with the second, by shift voxels, in x
21+
% Apply any needed flips---want the output to be flipped in x and y, since
22+
% that's what the rest of the pipeline expects
23+
if original_is_x_flipped ,
24+
% do nothing, this is how we want it
25+
else
26+
stack = fliplr(stack) ;
27+
end
28+
if original_is_y_flipped ,
29+
% do nothing, this is how we want it
30+
else
31+
stack = flipud(stack) ;
32+
end
33+
write_16bit_grayscale_tif(output_tif_file_path, stack) ;
34+
end
35+
36+
% Write the tile metadate file
37+
tile_metadata_file_path = fullfile(tile_folder_path, 'tile-metadata.txt') ;
38+
is_x_flipped = true ;
39+
is_y_flipped = true ;
40+
write_tile_metadata(tile_metadata_file_path, shift, is_x_flipped, is_y_flipped) ;
41+
end
File renamed without changes.

0 commit comments

Comments
 (0)