Skip to content

Commit

Permalink
update demo scripts and replace the demo data with a new one
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhoupc
zhoupc committed Oct 23, 2017
1 parent e0cf489 commit 69f7a7f
Show file tree
Hide file tree
Showing 7 changed files with 220 additions and 165 deletions.
2 changes: 1 addition & 1 deletion .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,4 +17,4 @@ deconvolveCa/example*
deconvolveCa/constrained-foopsi*
deconvolveCa/MCMC*

demos/compare/
demos/temp
14 changes: 12 additions & 2 deletions ca_source_extraction/@Sources2D/Sources2D.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -448,8 +448,18 @@ function updateTemporal_nb(obj, Y)
[merged_ROIs, newIDs] = MergeNeighbors(obj, dmin, method)

%------------------------------------------------------------------EXPORT---%
function save_neurons(obj)
folder_nm = [obj.P.log_folder, 'neurons'];
function save_neurons(obj, folder_nm)
if ~exist('folder_nm', 'var') || isempty(folder_nm)
folder_nm = [obj.P.log_folder, 'neurons'];
end
if exist(folder_nm, 'dir')
temp = cd();
cd(folder_nm);
delete *;
cd(temp);
else
mkdir(folder_nm);
end
obj.viewNeurons([], obj.C_raw, folder_nm);
end

Expand Down
Binary file added demos/data_1p.tif
View file
Binary file not shown.
Binary file removed demos/data_endoscope.tif
View file
Binary file not shown.
93 changes: 52 additions & 41 deletions demos/demo_batch_1p.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,107 +3,118 @@

%% choose multiple datasets or just one
neuron = Sources2D();
nams = {'./data_endoscope.tif'}; % you can put all file names into a cell array; when it's empty, manually select files
nams = {'./data_1p.tif'}; % you can put all file names into a cell array; when it's empty, manually select files
nams = neuron.select_multiple_files(nams); %if nam is [], then select data interactively

%% parameters
% ------------------------- COMPUTATION ------------------------- %
pars_envs = struct('memory_size_to_use', 8, ... % GB, memory space you allow to use in MATLAB
'memory_size_per_patch', 0.3, ... % GB, space for loading data within one patch
'memory_size_per_patch', 0.5, ... % GB, space for loading data within one patch
'patch_dims', [64, 64],... %GB, patch size
'batch_frames', 500); % number of frames per batch

% ------------------------- SPATIAL ------------------------- %
'batch_frames', 1000); % number of frames per batch
% ------------------------- SPATIAL ------------------------- %
gSig = 3; % pixel, gaussian width of a gaussian kernel for filtering the data. 0 means no filtering
gSiz = 11; % pixel, neuron diameter
gSiz = 13; % pixel, neuron diameter
ssub = 1; % spatial downsampling factor
with_dendrites = false; % with dendrites or not
with_dendrites = true; % with dendrites or not
if with_dendrites
% determine the search locations by dilating the current neuron shapes
updateA_search_method = 'dilate'; %#ok<UNRCH>
updateA_bSiz = 20;
updateA_dist = neuron.options.dist;
updateA_bSiz = 5;
updateA_dist = neuron.options.dist;
else
% determine the search locations by selecting a round area
updateA_search_method = 'ellipse'; %#ok<UNRCH>
updateA_dist = 5;
updateA_bSiz = neuron.options.dist;
end
spatial_constraints = struct('connected', true, 'circular', false); % you can include following constraints: 'circular'
spatial_algorithm = 'hals'; % algorithm for updating spatial components
spatial_algorithm = 'hals';

% ------------------------- TEMPORAL ------------------------- %
Fs = 6; % frame rate
Fs = 10; % frame rate
tsub = 1; % temporal downsampling factor
deconv_options = struct('type', 'ar1', ... % model of the calcium traces. {'ar1', 'ar2'}
'method', 'foopsi', ... % method for running deconvolution {'foopsi', 'constrained', 'thresholded'}
'smin', -3, ... % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
'smin', -5, ... % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
'optimize_pars', true, ... % optimize AR coefficients
'optimize_b', true);% optimize the baseline);
nk = 5; % detrending the slow fluctuation. usually 1 is fine (no detrending)
% when changed, try some integers smaller than total_frame/(Fs*30)
detrend_method = 'spline'; % compute the local minimum as an estimation of trend.
'optimize_b', true, ...% optimize the baseline);
'max_tau', 100); % maximum decay time (unit: frame);

nk = 3; % detrending the slow fluctuation. usually 1 is fine (no detrending)
% when changed, try some integers smaller than total_frame/(Fs*30)
detrend_method = 'spline'; % compute the local minimum as an estimation of trend.

% ------------------------- BACKGROUND ------------------------- %
bg_model = 'ring'; % model of the background {'ring', 'svd'(default), 'nmf'}
nb = 1; % number of background sources for each patch (only be used in SVD and NMF model)
bg_neuron_factor = 1.5;
ring_radius = round(bg_neuron_factor * gSiz); % when the ring model used, it is the radius of the ring used in the background model.
%otherwise, it's just the width of the overlapping area
bg_neuron_factor = 1.4;
ring_radius = round(bg_neuron_factor * gSiz); % when the ring model used, it is the radius of the ring used in the background model.
%otherwise, it's just the width of the overlapping area
num_neighbors = 50; % number of neighbors for each neuron

% ------------------------- MERGING ------------------------- %
show_merge = false; % if true, manually verify the merging step
merge_thr = 0.85; % thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
show_merge = false; % if true, manually verify the merging step
merge_thr = 0.65; % thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
method_dist = 'max'; % method for computing neuron distances {'mean', 'max'}
dmin = 5; % minimum distances between two neurons. it is used together with merge_thr
dmin_only = 1; % merge neurons if their distances are smaller than dmin_only.
dmin_only = 2; % merge neurons if their distances are smaller than dmin_only.
merge_thr_spatial = [0.8, 0.4, -inf]; % merge components with highly correlated spatial shapes (corr=0.8) and small temporal correlations (corr=0.1)

% ------------------------- INITIALIZATION ------------------------- %
K = []; % maximum number of neurons per patch. when K=[], take as many as possible.
K = []; % maximum number of neurons per patch. when K=[], take as many as possible.
min_corr = 0.8; % minimum local correlation for a seeding pixel
min_pnr = 10; % minimum peak-to-noise ratio for a seeding pixel
min_pixel = 4^2; % minimum number of nonzero pixels for each neuron
min_pnr = 8; % minimum peak-to-noise ratio for a seeding pixel
min_pixel = gSig^2; % minimum number of nonzero pixels for each neuron
bd = 0; % number of rows/columns to be ignored in the boundary (mainly for motion corrected data)
frame_range = []; % when [], uses all frames
save_initialization = false; % save the initialization procedure as a video.
use_parallel = true; % use parallel computation for parallel computing
show_init = true; % show initialization results
choose_params = false; % manually choose parameters
center_psf = true; % set the value as true when the background fluctuation is large (usually 1p data)
% set the value as false when the background fluctuation is small (2p)
frame_range = []; % when [], uses all frames
save_initialization = false; % save the initialization procedure as a video.
use_parallel = true; % use parallel computation for parallel computing
show_init = true; % show initialization results
choose_params = true; % manually choose parameters
center_psf = true; % set the value as true when the background fluctuation is large (usually 1p data)
% set the value as false when the background fluctuation is small (2p)

% ------------------------- Residual ------------------------- %
min_corr_res = 0.7;
min_pnr_res = 6;
seed_method_res = 'auto'; % method for initializing neurons from the residual
update_sn = true;

% ---------------------- WITH MANUAL INTERVENTION -------------------- %
with_manual_intervention = true;
with_manual_intervention = false;

% ------------------------- FINAL RESULTS ------------------------- %
save_demixed = true; % save the demixed file or not
kt = 1; % frame intervals
save_demixed = true; % save the demixed file or not
kt = 3; % frame intervals

% ------------------------- UPDATE ALL ------------------------- %
neuron.updateParams('gSig', gSig, ... % -------- spatial --------
neuron.updateParams('gSig', gSig, ... % -------- spatial --------
'gSiz', gSiz, ...
'ring_radius', ring_radius, ...
'ssub', ssub, ...
'search_method', updateA_search_method, ...
'bSiz', updateA_bSiz, ...
'dist', updateA_bSiz, ...
'spatial_constraints', spatial_constraints, ...
'tsub', tsub, ... % -------- temporal --------
'deconv_options', deconv_options, ...
'spatial_algorithm', spatial_algorithm, ...
'tsub', tsub, ... % -------- temporal --------
'deconv_options', deconv_options, ...
'nk', nk, ...
'detrend_method', detrend_method, ...
'background_model', bg_model, ... % -------- background --------
'background_model', bg_model, ... % -------- background --------
'nb', nb, ...
'ring_radius', ring_radius, ...
'num_neighbors', num_neighbors, ...
'merge_thr', merge_thr, ... % -------- merging ---------
'dmin', dmin, ...
'method_dist', method_dist, ...
'min_corr', min_corr, ... % ----- initialization -----
'min_corr', min_corr, ... % ----- initialization -----
'min_pnr', min_pnr, ...
'min_pixel', min_pixel, ...
'bd', bd, ...
'center_psf', center_psf);
neuron.Fs = Fs;
neuron.Fs = Fs;

%% distribute data and be ready to run source extraction
neuron.getReady_batch(pars_envs);
Expand Down
41 changes: 17 additions & 24 deletions demos/demo_endoscope.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,15 +4,15 @@
neuron_full Ybg_weights; %#ok<NUSED> % global variables, don't change them manually

%% select data and map it to the RAM
nam = './data_endoscope.tif';
nam = './data_1p.tif';
cnmfe_choose_data;

%% create Source2D class object for storing results and parameters
Fs = 6; % frame rate
Fs = 10; % frame rate
ssub = 1; % spatial downsampling factor
tsub = 1; % temporal downsampling factor
gSig = 3; % width of the gaussian kernel, which can approximates the average neuron shape
gSiz = 11; % maximum diameter of neurons in the image plane. larger values are preferred.
gSiz = 13; % maximum diameter of neurons in the image plane. larger values are preferred.
neuron_full = Sources2D('d1',d1,'d2',d2, ... % dimensions of datasets
'ssub', ssub, 'tsub', tsub, ... % downsampleing
'gSig', gSig,... % sigma of the 2D gaussian that approximates cell bodies
Expand All @@ -34,11 +34,12 @@
merge_thr = [1e-1, 0.85, 0]; % thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
dmin = 1;
%% options for running deconvolution
neuron_full.options.deconv_options = struct('type', 'ar1', ... % model of the calcium traces. {'ar1', 'ar2'}
neuron.options.deconv_options = struct('type', 'ar1', ... % model of the calcium traces. {'ar1', 'ar2'}
'method', 'foopsi', ... % method for running deconvolution {'foopsi', 'constrained', 'thresholded'}
'smin', -3, ... % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
'smin', -5, ... % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
'optimize_pars', true, ... % optimize AR coefficients
'optimize_b', true);% optimize the baseline);
'optimize_b', true, ...% optimize the baseline);
'max_tau', 100); % maximum decay time (unit: frame);

%% downsample data for fast and better initialization
sframe=1; % user input: first frame to read (optional, default:1)
Expand All @@ -62,8 +63,8 @@
K = []; % maximum number of neurons to search within each patch. you can use [] to search the number automatically

min_corr = 0.8; % minimum local correlation for a seeding pixel
min_pnr = 10; % minimum peak-to-noise ratio for a seeding pixel
min_pixel = 3^2; % minimum number of nonzero pixels for each neuron
min_pnr = 8; % minimum peak-to-noise ratio for a seeding pixel
min_pixel = gSig^2; % minimum number of nonzero pixels for each neuron
bd = 1; % number of rows/columns to be ignored in the boundary (mainly for motion corrected data)
neuron.updateParams('min_corr', min_corr, 'min_pnr', min_pnr, ...
'min_pixel', min_pixel, 'bd', bd);
Expand Down Expand Up @@ -189,41 +190,33 @@

b0 = mean(Y,2)-neuron.A*mean(neuron.C,2);
Ybg = bsxfun(@plus, Ybg, b0-mean(Ybg, 2));
neuron.orderROIs('snr');

%% display neurons
dir_neurons = sprintf('%s%s%s_neurons%s', dir_nm, filesep, file_nm, filesep);
if exist('dir_neurons', 'dir')
temp = cd();
cd(dir_neurons);
delete *;
cd(temp);
else
mkdir(dir_neurons);
end
neuron.viewNeurons([], neuron.C_raw, dir_neurons);
close(gcf);
neuron.save_neurons(dir_neurons);

%% display contours of the neurons
neuron.Coor = neuron.get_contours(0.8); % energy within the contour is 80% of the total
figure;
Cnn = correlation_image(neuron.reshape(Ysignal(:, 1:5:end), 2), 4);
neuron.Coor = plot_contours(neuron.A, Cnn, 0.8, 0, [], neuron.Coor, 2);
colormap winter;
neuron.show_contours(0.6);
colormap gray;
axis equal; axis off;
title('contours of estimated neurons');

% plot contours with IDs
% [Cn, pnr] = neuron.correlation_pnr(Y(:, round(linspace(1, T, min(T, 1000)))));
figure;
Cn = imresize(Cn, [d1, d2]);
plot_contours(neuron.A, Cn, 0.8, 0, [], neuron.Coor, 2);
colormap winter;
neuron.show_contours(0.6);
colormap gray;
title('contours of estimated neurons');

%% check spatial and temporal components by playing movies
save_avi = false;
avi_name = 'play_movie.avi';
neuron.Cn = Cn;
neuron.runMovie(Ysignal, [0, 50], save_avi, avi_name);
neuron.runMovie(Ysignal, [0, 150], save_avi, avi_name);

%% save video
kt = 3; % play one frame in every kt frames
Expand Down
Loading

0 comments on commit 69f7a7f

Please sign in to comment.