Folded in code to label centerpoints as machine vs human

.gitignore

@@ -13,5 +13,4 @@ purge.lock
 *.stdout
 *.stderr
 *.stdouterr
-2018-08-01
-2018-10-01
+20??-??-??

cclf.m

@@ -0,0 +1,4 @@
+close all ;
+clear ;
+
+

cellfilt.m

@@ -0,0 +1,4 @@
+function result = cellfilt(predicate, cell_array)
+    is_predicate_true = cellfun(predicate, cell_array) ;
+    result = cell_array(is_predicate_true) ;    
+end

collect_consensus_neurons.m

@@ -0,0 +1,6 @@
+function result = collect_consensus_neurons(consensus_swcs_folder_path)    
+    [consensus_swc_file_paths, consensus_neuron_names] = collect_consensus_swc_file_paths(consensus_swcs_folder_path) ;
+    consensus_neurons = cellfun(@load_swc, consensus_swc_file_paths, 'UniformOutput', false) ;
+    result = struct('consensus_neurons', {consensus_neurons}, ...
+                    'consensus_neuron_names', {consensus_neuron_names}) ;        
+end

collect_consensus_swc_file_paths.m

@@ -0,0 +1,34 @@
+function [consensus_swc_paths, neuron_names_with_consensus_swc_paths]  = collect_consensus_swc_file_paths(consensus_swcs_folder_path)    
+    % E.g. one consensus .swc file path within the consensus_swcs_folder_path
+    %
+    % /groups/mousebrainmicro/mousebrainmicro/shared_tracing/Finished_Neurons/2018-10-01
+    %
+    % is
+    % 
+    % /groups/mousebrainmicro/mousebrainmicro/shared_tracing/Finished_Neurons/2018-10-01/G-014/consensus/2018-10-01_G-014_consensus.swc
+
+    date_as_string = leaf_file_name_from_path(consensus_swcs_folder_path) ;
+    neuron_folder_path_template = fullfile(consensus_swcs_folder_path, 'G-*') ;
+    neuron_names = simple_dir(neuron_folder_path_template) ;
+
+    function result = consensus_swc_path_from_neuron_name(neuron_name)
+        lower_case_file_name = sprintf('%s_%s_consensus.swc', date_as_string, neuron_name) ;
+        lower_case_file_path = fullfile(consensus_swcs_folder_path, neuron_name, 'consensus', lower_case_file_name) ;
+        if exist(lower_case_file_path, 'file') ,
+            result = lower_case_file_path ;
+        else
+            upper_case_file_name = sprintf('%s_%s_Consensus.swc', date_as_string, neuron_name) ;
+            upper_case_file_path = fullfile(consensus_swcs_folder_path, neuron_name, 'Consensus', upper_case_file_name) ;
+            if exist(upper_case_file_path, 'file') ,
+                result = upper_case_file_path ;
+            else
+                result = '' ;
+            end
+        end
+    end
+
+    consensus_swc_paths_but_with_empties = cellfun(@consensus_swc_path_from_neuron_name, neuron_names, 'UniformOutput', false) ;
+    has_consensus_swc = cellfun(@is_nonempty, consensus_swc_paths_but_with_empties) ;
+    consensus_swc_paths = consensus_swc_paths_but_with_empties(has_consensus_swc) ;
+    neuron_names_with_consensus_swc_paths = neuron_names(has_consensus_swc) ;
+end

compute_or_read_from_memo.m

@@ -0,0 +1,21 @@
+function result = compute_or_read_from_memo(cache_folder_path, ...
+                                            memo_base_name, ...
+                                            compute_function) 
+    % A function to help with the common pattern of wanting to cache expensive-to-compute quantities on disk                                    
+    memo_file_name = horzcat(memo_base_name, '.mat') ;
+    memo_file_path = fullfile(cache_folder_path, memo_file_name) ;
+    if exist(cache_folder_path, 'file') ,
+        if exist(memo_file_path, 'file') ,        
+            load(memo_file_path, 'result') ;
+        else
+            result = feval(compute_function) ;
+            save(memo_file_path, '-mat', '-v7.3', 'result') ;
+        end
+    else
+        mkdir(cache_folder_path) ;
+        % We know memo_file_name doesn't exist, b/c we just created its
+        % parent folder
+        result = feval(compute_function) ;
+        save(memo_file_path, '-mat', '-v7.3', 'result') ;        
+    end
+end

configparser.m

@@ -0,0 +1,80 @@
+function [out] = configparser(configfile)
+%CONFIGPARSER reads a config file and return a structure array based on the
+%fields of the configfile
+%
+% [OUTPUTARGS] = CONFIGPARSER(configfile) 
+%
+% Inputs:
+%
+% Outputs:
+%
+% Examples:
+%
+% See also: 
+
+% $Author: base $	$Date: 2016/01/14 15:55:07 $	$Revision: 0.1 $
+% Copyright: HHMI 2016
+
+if nargin<1
+    configfile = 'myparam'
+end
+
+fid = fopen(configfile);
+out = [];
+tline = fgetl(fid);
+while ischar(tline)
+    if isempty(tline)
+    else
+        out = checkcase(tline,out);
+    end
+    tline = fgetl(fid);
+end
+fclose(fid);
+
+gt=fieldnames(out);
+for ii=1:length(gt)
+    fieldtxt = strtrim(gt{ii});
+    txt = strtrim(out.(gt{ii}));
+    if strcmp(gt{ii},'HEADER')
+    elseif strcmp(gt{ii},'Tform')
+        out.(fieldtxt) = cellfun(@str2num,txt);
+    else
+        if iscell(txt)
+            out.(fieldtxt) = cellfun(@str2num,txt);
+        elseif txt(1)=='''' % string, most likely path
+            out.(fieldtxt) = eval(txt);
+        else
+            out.(fieldtxt) = str2num(txt);
+        end
+    end
+end
+
+end
+
+function out = checkcase(tline,out)
+if tline(1)=='#' % header skip
+    if isfield(out,'HEADER')
+        out.HEADER = sprintf('%s%s\n',out.HEADER,tline);
+    else
+        out.HEADER = sprintf('%s\n',tline);
+    end
+else
+    fd = strsplit(tline,'=');
+    if length(fd)<2
+        % check for : as splitter
+        fd = strsplit(tline,':');
+    end
+    if length(fd)<2
+        error('Couldnt find any argument with "=" or ":"')
+    end
+    fd_sub = strsplit(strtrim(fd{2}),' ');
+    if length(fd_sub)>2
+        out.(deblank(fd{1})) = fd_sub;
+    else
+        out.(deblank(fd{1})) = fd{2};
+    end
+end
+
+end
+
+

investigate_export_issues.m

@@ -0,0 +1,21 @@
+original_neuron = load_swc('test-tree/auto111_cc-0231.swc')
+exported_neuron = load_swc('auto111_cc-0231-as-exported.swc')
+exported_neuron - original_neuron  % all off by same amount
+mean(original_neuron)
+diff(original_neuron,1)  
+diff(exported_neuron,1)  % same (up to floating-point error) as last value
+[origin, spacing] = load_transform_txt('/nrs/mouselight/SAMPLES/2018-10-01/transform.txt')
+origin = [69445.01944245792, 12917.29486937314, 30198.96941474185] ;  % origin from Workstation
+spacing(3)/spacing(1)  % about 4
+spacing(3)/spacing(2)  % about 4
+original_neuron_xyz = original_neuron(:,3:5)
+exported_neuron_xyz = exported_neuron(:,3:5)
+diff(original_neuron_xyz)
+dr_original = diff(original_neuron_xyz)
+bsxfun(@times, original_neuron_xyz-origin, 1./spacing)  % half-integers
+bsxfun(@times, exported_neuron_xyz-origin, 1./spacing)  % not half-integers
+exported_centroid = mean(exported_neuron_xyz)
+bsxfun(@times, exported_centroid-origin, 1./spacing)  % not half-integers, or integers...
+exported_neuron_xyz - original_neuron_xyz  % all columns the same
+original_centroid = mean(original_neuron_xyz)
+exported_centroid-original_centroid

is_nonempty.m

@@ -0,0 +1,3 @@
+function result = is_nonempty(x)
+    result = ~isempty(x) ;
+end

label_machine_centerpoints.m

@@ -0,0 +1,36 @@
+function consensus_neurons_with_machine_centerpoints_labelled = ...
+        label_machine_centerpoints(output_folder_path, ...
+                                   cache_folder_path, ...
+                                   fragments_folder_path, ...
+                                   consensus_swcs_folder_path, ...
+                                   sample_folder_path)
+    all_fragment_centerpoints = compute_or_read_from_memo(cache_folder_path, ...
+                                                          'all_fragment_centerpoints',  ...
+                                                          @()(read_all_fragment_centerpoints(fragments_folder_path))) ;
+    consensus_neurons_and_names = compute_or_read_from_memo(cache_folder_path, ...
+                                                  'consensus_neurons',  ...
+                                                  @()(collect_consensus_neurons(consensus_swcs_folder_path))) ;    
+    consensus_neurons =  consensus_neurons_and_names.consensus_neurons ;
+    consensus_neuron_names = consensus_neurons_and_names.consensus_neuron_names ;
+
+    [~, spacing] = load_transform_txt(fullfile(sample_folder_path, 'transform.txt')) ;
+    consensus_neurons_with_machine_centerpoints_labelled = ...
+        compute_or_read_from_memo(output_folder_path, ...
+                                  'consensus_neurons_with_machine_centerpoints_labelled',  ...
+                                  @()(label_machine_centerpoints_in_neurons(consensus_neurons, all_fragment_centerpoints, spacing))) ;                                  
+
+    % Write the swc's that don't exist already                          
+    consensus_neuron_count = length(consensus_neurons) ;
+    swc_folder_path = fullfile(output_folder_path, 'consensus-neurons-with-machine-centerpoints-labelled-as-swcs') ;
+    if consensus_neuron_count>0 && ~exist(swc_folder_path, 'file') ,
+        mkdir(swc_folder_path) ;
+    end
+    for i = 1:consensus_neuron_count ,
+        consensus_neuron_name = consensus_neuron_names{i} ;
+        swc_file_name = sprintf('%s.swc', consensus_neuron_name) ;
+        swc_file_path = fullfile(swc_folder_path, swc_file_name) ;
+        if ~exist(swc_file_path, 'file') ,
+            save_swc(swc_file_path, consensus_neurons_with_machine_centerpoints_labelled{i}, consensus_neuron_name) ;
+        end
+    end
+end

label_machine_centerpoints_2018_07_02.m

@@ -0,0 +1,11 @@
+sample_date = '2018-07-02' ;
+this_file_path = mfilename('fullpath') ;
+this_folder_path = fileparts(this_file_path) ;
+output_folder_path = fullfile(this_folder_path, sample_date) ;
+cache_folder_path = fullfile(this_folder_path, sample_date) ;
+fragments_folder_path = sprintf('/groups/mousebrainmicro/mousebrainmicro/cluster/Reconstructions/%s/prob0_swcs/frags', sample_date) ;
+consensus_swcs_folder_path = sprintf('/groups/mousebrainmicro/mousebrainmicro/shared_tracing/Finished_Neurons/%s', sample_date) ;
+sample_folder_path = sprintf('/nrs/mouselight/SAMPLES/%s', sample_date) ;
+
+consensus_neurons_with_machine_centerpoints_labelled = ...
+    label_machine_centerpoints(output_folder_path, cache_folder_path, fragments_folder_path, consensus_swcs_folder_path, sample_folder_path) ;

label_machine_centerpoints_2018_08_01.m

@@ -0,0 +1,11 @@
+sample_date = '2018-08-01' ;
+this_file_path = mfilename('fullpath') ;
+this_folder_path = fileparts(this_file_path) ;
+output_folder_path = fullfile(this_folder_path, sample_date) ;
+cache_folder_path = fullfile(this_folder_path, sample_date) ;
+fragments_folder_path = sprintf('/groups/mousebrainmicro/mousebrainmicro/cluster/Reconstructions/%s/prob0_swcs/frags', sample_date) ;
+consensus_swcs_folder_path = sprintf('/groups/mousebrainmicro/mousebrainmicro/shared_tracing/Finished_Neurons/%s', sample_date) ;
+sample_folder_path = sprintf('/nrs/mouselight/SAMPLES/%s', sample_date) ;
+
+consensus_neurons_with_machine_centerpoints_labelled = ...
+    label_machine_centerpoints(output_folder_path, cache_folder_path, fragments_folder_path, consensus_swcs_folder_path, sample_folder_path) ;

label_machine_centerpoints_2018_10_01.m

@@ -0,0 +1,8 @@
+output_folder_path = '/home/taylora/cache/gt/2018-10-01' ;
+cache_folder_path = '/home/taylora/cache/gt/2018-10-01' ;
+fragments_folder_path = '/groups/mousebrainmicro/mousebrainmicro/cluster/Reconstructions/2018-10-01/prob0_swcs/frags' ;
+consensus_swcs_folder_path = '/groups/mousebrainmicro/mousebrainmicro/shared_tracing/Finished_Neurons/2018-10-01' ;
+sample_folder_path = '/nrs/mouselight/SAMPLES/2018-10-01' ;
+
+consensus_neurons_with_machine_centerpoints_labelled = ...
+    label_machine_centerpoints(output_folder_path, cache_folder_path, fragments_folder_path, consensus_swcs_folder_path, sample_folder_path) ;

label_machine_centerpoints_in_neuron.m

@@ -0,0 +1,12 @@
+function result = label_machine_centerpoints_in_neuron(neuron, kd_tree, machine_centerpoints, spacing)
+    neuron_centerpoints = neuron(:,3:5) ;
+    index_of_nearest_machine_point = knnsearch(kd_tree, neuron_centerpoints) ;
+    nearest_machine_point = machine_centerpoints(index_of_nearest_machine_point, :) ;
+    offset_to_nearest_machine_point = nearest_machine_point - neuron_centerpoints ;    
+    offset_threshold = spacing/2 ;  % um 
+    is_offset_below_threshold = bsxfun(@lt, offset_to_nearest_machine_point, offset_threshold) ;
+    is_a_machine_point = all(is_offset_below_threshold, 2) ;  
+    new_structure_identifier = 42 + is_a_machine_point ;
+    result = neuron ;
+    result(:,2) = new_structure_identifier ;
+end

label_machine_centerpoints_in_neurons.m

@@ -0,0 +1,12 @@
+function result = label_machine_centerpoints_in_neurons(neurons, machine_centerpoints, spacing)
+    % neurons is a cell array, each element an centerpoint_count x 7 double array
+    % Each row is an SWC-style record: id, type, x, y, z, something i forget,
+    % and parent id.  The root has a parent id of -1, and parents come before
+    % children.
+    %
+    % machine_centerpoints is a set of machine-annotated points, each row an x
+    % y z triple
+
+    kd_tree  = KDTreeSearcher(machine_centerpoints) ;        
+    result = cellfun(@(neuron)(label_machine_centerpoints_in_neuron(neuron, kd_tree, machine_centerpoints, spacing)), neurons, 'UniformOutput', false) ;
+end

leaf_file_name_from_path.m

@@ -0,0 +1,4 @@
+function result = leaf_file_name_from_path(path) 
+    [~,base,ext] = fileparts(path) ;
+    result = horzcat(base, ext) ;
+end

load_swc.m

@@ -0,0 +1,73 @@
+function swc_data = load_swc(swc_file_name)
+    %LOADSWC reads a JW formatted swc file and return color/header and data
+    %fields
+    %
+    % [OUTPUTARGS] = LOADSWC(INPUTARGS) Explain usage here
+    %
+    % Examples:
+    %     [swcData,offset,color, header] = loadSWC(swcfile);
+    %     swcData(:,3:5) = swcData(:,3:5) + ones(size(swcData,1),1)*offset;
+    %     swcData(:,3:5) = swcData(:,3:5)*scale;
+    %
+    % Provide sample usage code here
+    %
+    % See also: List related files here
+
+    % $Author: base $	$Date: 2015/08/14 12:09:52 $	$Revision: 0.1 $
+    % Copyright: HHMI 2015
+
+    %% parse a swc file
+    offset = [0 0 0] ;
+    header = cell(1);
+
+    fid = fopen(swc_file_name);
+    this_line = fgets(fid);
+    skipline = 0;
+    header_line_count = 1;
+    centerpoint_count = 0 ;
+    while ischar(this_line)
+        % assumes that all header info is at the top of the swc file
+        if strcmp(this_line(1), '#')
+            skipline = skipline + 1;
+            header{header_line_count} = this_line;
+            header_line_count = header_line_count+1;
+        else
+            centerpoint_count = centerpoint_count + 1 ;
+            this_line = fgets(fid);            
+            continue;
+            %break
+        end
+        % if get here, tline contains a header line
+        if contains(this_line, 'OFFSET') ,
+            tokens = strsplit(deblank(this_line)) ;
+            tokens_as_double = cellfun(@str2double, tokens) ;
+            numeric_tokens = tokens_as_double(~isnan(tokens_as_double)) ;
+            if isequal(size(numeric_tokens), [1 3]) ,
+                offset = numeric_tokens ;
+            else
+                error('Unable to parse OFFSET header line: "%s"', this_line) ;
+            end
+        elseif strcmp(this_line(1:8),'# COLOR ')
+            color =cellfun(@str2double,strsplit(deblank(this_line(9:end)),','));  %#ok<NASGU>
+        end
+        this_line = fgets(fid);
+    end
+    fclose(fid);
+
+    fid = fopen(swc_file_name);
+    for i=1:skipline
+        this_line = fgets(fid);  %#ok<NASGU>
+    end
+    swc_data = zeros(centerpoint_count,7) ;
+    this_line = fgets(fid);
+    tl = 1;
+    while ischar(this_line)
+        swc_data(tl,:) = str2num(this_line);  %#ok<ST2NM>
+        tl = tl+1;
+        this_line = fgets(fid);
+    end
+    fclose(fid);
+    raw_centerpoints = swc_data(:,3:5) ;
+    centerpoints = bsxfun(@plus, raw_centerpoints, offset) ;    
+    swc_data(:,3:5) = centerpoints ;
+end

load_transform_txt.m

@@ -0,0 +1,11 @@
+function [origin, spacing] = load_transform_txt(file_name)
+    transform_as_raw_struct = configparser(file_name) ;
+        % has fields, ox, oy, oz, sx, sy, sz, all in nm
+        % also nl, the number of levels in the octree
+    origin = [ transform_as_raw_struct.ox transform_as_raw_struct.oy transform_as_raw_struct.oz]/1e3 ;  % nm -> um
+    spacing = [ transform_as_raw_struct.sx transform_as_raw_struct.sy transform_as_raw_struct.sz]/1e3/2^(transform_as_raw_struct.nl-1) ;  % nm -> um
+        % determine the spacing of the voxels at the bottom level
+    % Note that this is not *exactly* the origin used by Jaws as of this writing (2019-07-16).
+    % But it's intended to be the location of the voxel *center* of the
+    % 'lowest-index' voxel in the stack.
+end

makecellfun.m

@@ -0,0 +1,9 @@
+function makecellfun(lambda, cell_array, file_name_template)
+    n = length(cell_array) ;
+    for i = 1:n ,
+        file_name = sprintf(file_name_template, i) ;
+        if ~exist(file_name, 'file') ,
+            feval(lambda, cell_array{i}, file_name) ;
+        end
+    end
+end

one_based_ijk_from_xyz.m

@@ -0,0 +1,7 @@
+function result = one_based_ijk_from_xyz(xyz, origin, spacing)
+    % xyz is n x 3, in um
+    % offset is 1 x 3, in um
+    % spacing is 1 x 3, in um
+    ijk_zero_based = zero_based_ijk_from_xyz(xyz, origin, spacing) ;
+    result = ijk_zero_based + 1 ;
+end