-
Notifications
You must be signed in to change notification settings - Fork 146
/
demo_script.m
126 lines (100 loc) · 5.13 KB
/
demo_script.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
117
118
119
120
121
122
123
124
125
126
clear;
%% load file
gcp; % start cluster
addpath(genpath('utilities'));
nam = 'demoMovie.tif'; % insert path to tiff stack here
sframe=1; % user input: first frame to read (optional, default 1)
num2read=2000; % user input: how many frames to read (optional, default until the end)
Y = bigread2(nam,sframe,num2read);
%Y = Y - min(Y(:));
if ~isa(Y,'double'); Y = double(Y); end % convert to single
[d1,d2,T] = size(Y); % dimensions of dataset
d = d1*d2; % total number of pixels
%% Set parameters
K = 40; % number of components to be found
tau = 4; % std of gaussian kernel (size of neuron)
p = 2; % order of autoregressive system (p = 0 no dynamics, p=1 just decay, p = 2, both rise and decay)
merge_thr = 0.8; % merging threshold
options = CNMFSetParms(...
'd1',d1,'d2',d2,... % dimensions of datasets
'search_method','dilate','dist',3,... % search locations when updating spatial components
'deconv_method','constrained_foopsi',... % activity deconvolution method
'temporal_iter',2,... % number of block-coordinate descent steps
'fudge_factor',0.98,... % bias correction for AR coefficients
'merge_thr',merge_thr,... % merging threshold
'gSig',tau...
);
%% Data pre-processing
[P,Y] = preprocess_data(Y,p);
%% fast initialization of spatial components using greedyROI and HALS
[Ain,Cin,bin,fin,center] = initialize_components(Y,K,tau,options,P); % initialize
% display centers of found components
Cn = correlation_image(Y); %reshape(P.sn,d1,d2); %max(Y,[],3); %std(Y,[],3); % image statistic (only for display purposes)
figure;imagesc(Cn);
axis equal; axis tight; hold all;
scatter(center(:,2),center(:,1),'mo');
title('Center of ROIs found from initialization algorithm');
drawnow;
%% manually refine components (optional)
refine_components = false; % flag for manual refinement
if refine_components
[Ain,Cin,center] = manually_refine_components(Y,Ain,Cin,center,Cn,tau,options);
end
%% update spatial components
Yr = reshape(Y,d,T);
[A,b,Cin] = update_spatial_components(Yr,Cin,fin,[Ain,bin],P,options);
%% update temporal components
P.p = 0; % set AR temporarily to zero for speed
[C,f,P,S,YrA] = update_temporal_components(Yr,A,b,Cin,fin,P,options);
%% classify components
[ROIvars.rval_space,ROIvars.rval_time,ROIvars.max_pr,ROIvars.sizeA,keep] = classify_components(Y,A,C,b,f,YrA,options);
%% further classification with cnn_classifier
try % matlab 2017b or later is needed
[ind,value] = cnn_classifier(A,[d1,d2],'cnn_model',0.2);
catch
ind = true(size(A,2),1);
end
%% display kept and discarded components
A_keep = A(:,(keep & ind));
C_keep = C((keep & ind),:);
figure;
subplot(121); montage(extract_patch(A(:,(keep & ind)),[d1,d2],[30,30]),'DisplayRange',[0,0.15]);
title('Kept Components');
subplot(122); montage(extract_patch(A(:,~(keep & ind)),[d1,d2],[30,30]),'DisplayRange',[0,0.15])
title('Discarded Components');
%% merge found components
[Am,Cm,K_m,merged_ROIs,Pm,Sm] = merge_components(Yr,A_keep,b,C_keep,f,P,S,options);
%%
display_merging = 1; % flag for displaying merging example
if and(display_merging, ~isempty(merged_ROIs))
i = 1; %randi(length(merged_ROIs));
ln = length(merged_ROIs{i});
figure;
set(gcf,'Position',[300,300,(ln+2)*300,300]);
for j = 1:ln
subplot(1,ln+2,j); imagesc(reshape(A_keep(:,merged_ROIs{i}(j)),d1,d2));
title(sprintf('Component %i',j),'fontsize',16,'fontweight','bold'); axis equal; axis tight;
end
subplot(1,ln+2,ln+1); imagesc(reshape(Am(:,K_m-length(merged_ROIs)+i),d1,d2));
title('Merged Component','fontsize',16,'fontweight','bold');axis equal; axis tight;
subplot(1,ln+2,ln+2);
plot(1:T,(diag(max(C_keep(merged_ROIs{i},:),[],2))\C_keep(merged_ROIs{i},:))');
hold all; plot(1:T,Cm(K_m-length(merged_ROIs)+i,:)/max(Cm(K_m-length(merged_ROIs)+i,:)),'--k')
title('Temporal Components','fontsize',16,'fontweight','bold')
drawnow;
end
%% refine estimates excluding rejected components
Pm.p = p; % restore AR value
[A2,b2,C2] = update_spatial_components(Yr,Cm,f,[Am,b],Pm,options);
[C2,f2,P2,S2,YrA2] = update_temporal_components(Yr,A2,b2,C2,f,Pm,options);
%% do some plotting
[A_or,C_or,S_or,P_or] = order_ROIs(A2,C2,S2,P2); % order components
K_m = size(C_or,1);
[C_df,~] = extract_DF_F(Yr,A_or,C_or,P_or,options); % extract DF/F values (optional)
figure;
[Coor,json_file] = plot_contours(A_or,Cn,options,1); % contour plot of spatial footprints
%savejson('jmesh',json_file,'filename'); % optional save json file with component coordinates (requires matlab json library)
%% display components
plot_components_GUI(Yr,A_or,C_or,b2,f2,Cn,options)
%% make movie
make_patch_video(A_or,C_or,b2,f2,Yr,Coor,options)