Skip to content

Commit 4670b57

Browse files
epnevepnev
committed
3d initialization
1 parent 4ead56b commit 4670b57

File tree

3 files changed

+68
-66
lines changed

3 files changed

+68
-66
lines changed

‎initialize_components.m‎

Lines changed: 44 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -51,21 +51,40 @@
5151
if ~isfield(options, 'tsub'), options.tsub = 1; end; tsub = options.tsub;
5252
if tsub == 1; fprintf('No temporal downsampling is performed. Consider temporal downsampling if the recording is very long. \n'); end
5353

54-
[d1,d2,T] = size(Y);
55-
d1s = ceil(d1/ssub); %size of downsampled image
56-
d2s = ceil(d2/ssub);
54+
ndimsY = ndims(Y)-1;
55+
sY = size(Y);
56+
d = sY(1:ndimsY);
57+
T = sY(end);
58+
59+
ds = d;
60+
ds(1:2) = ceil(d(1:2)/ssub); % do not subsample along z axis
61+
%d1s = ceil(d1/ssub); %size of downsampled image
62+
%d2s = ceil(d2/ssub);
5763
Ts = floor(T/tsub); %reduced number of frames
5864
% spatial downsampling
59-
if ssub~=1, Y_ds = imresize(Y, [d1s, d2s], 'box'); else Y_ds = Y; end
65+
fprintf('starting resampling \n')
66+
if ssub~=1;
67+
if ndimsY == 2; Y_ds = imresize(Y, [ds(1), ds(2)], 'box'); end
68+
if ndimsY == 3;
69+
Y_ds = zeros([ds(1:2),T,ds(end)]);
70+
for z = 1:ds(3)
71+
Y_ds(:,:,:,z) = imresize(squeeze(Y(:,:,z,:)), [ds(1), ds(2)], 'box');
72+
end
73+
Y_ds = permute(Y_ds,[1,2,4,3]);
74+
end
75+
else
76+
Y_ds = Y;
77+
end
6078
% temporal downsampling
6179
if tsub~=1
62-
Y_ds = squeeze(mean(reshape(Y_ds(:, :, 1:(Ts*tsub)),d1s, d2s, tsub, Ts), 3));
80+
if ndimsY == 2; Y_ds = squeeze(mean(reshape(Y_ds(:, :, 1:(Ts*tsub)),ds(1), ds(2), tsub, Ts), 3)); end
81+
if ndimsY == 3; Y_ds = squeeze(mean(reshape(Y_ds(:, :, :, 1:(Ts*tsub)),ds(1), ds(2), ds(3), tsub, Ts), 4)); end
6382
end
6483

6584
if strcmpi(options.init_method,'greedy')
6685
% run greedy method
6786
fprintf('Initializing components with greedy method \n');
68-
[Ain, Cin, bin, fin] = greedyROI2d(Y_ds, K, options);
87+
[Ain, Cin, bin, fin] = greedyROI(Y_ds, K, options);
6988
elseif strcmpi(options.init_method,'sparse_NMF')
7089
% run sparse_NMF method
7190
fprintf('Initializing components with sparse NMF \n');
@@ -75,19 +94,26 @@
7594
end
7695

7796
% refine with HALS
78-
fprintf('Refining initial estimates with HALS \n');
79-
[Ain, Cin, bin, fin] = HALS_2d(Y_ds, full(Ain), Cin, bin, fin, options);
80-
97+
fprintf('Refining initial estimates with HALS...');
98+
%[Ain1, Cin1, bin1, fin1] = HALS_2d(Y_ds, full(Ain), Cin, bin, fin, options);
99+
[Ain, Cin, bin, fin] = HALS(Y_ds, full(Ain), Cin, bin, fin, options);
100+
fprintf(' done \n');
81101
%% upsample Ain, Cin, bin, fin
82-
center = ssub*com(Ain,d1s,d2s);
83-
Ain = imresize(reshape(Ain, d1s, d2s, sum(K)), [d1, d2]);
84-
Ain = reshape(Ain, d1*d2, []);
85-
bin_temp = reshape(bin, d1s, d2s, options.nb);
86-
bin = zeros(d1,d2,options.nb);
87-
for i = 1:options.nb
88-
bin(:,:,i) = imresize(bin_temp(:,:,i), [d1, d2]);
89-
end
90-
bin = reshape(bin, [], options.nb);
102+
if ndimsY == 2; center = ssub*com(Ain,ds(1),ds(2)); else center = ssub*com(Ain,ds(1),ds(2),ds(3)); end
103+
%Ain = imresize(reshape(Ain, ds(1), ds(2), sum(K)), d);
104+
%Ain = imresize(reshape(full(Ain), [ds, sum(K)]), d);
105+
Ain = imresize(reshape(full(Ain), [ds(1),ds(2), sum(K)*prod(ds)/ds(1)/ds(2)]),[d(1),d(2)]); %,prod(d)/d(1)/d(2)*sum(K)]);
106+
Ain = sparse(reshape(Ain, prod(d), []));
107+
%bin_temp = reshape(bin, ds(1), ds(2), options.nb);
108+
%bin = zeros(d(1),d(2),options.nb);
109+
bin = imresize(reshape(bin,[ds(1),ds(2), options.nb*prod(ds)/ds(1)/ds(2)]),[d(1),d(2)]);
110+
bin = reshape(bin,prod(d),[]);
111+
% bin_temp = reshape(bin, [ds, options.nb]);
112+
% bin = zeros([d,options.nb]);
113+
% for i = 1:options.nb
114+
% bin(:,:,i) = imresize(bin_temp(:,:,i), d);
115+
% end
116+
%bin = reshape(bin, [], options.nb);
91117
Cin = imresize(Cin, [sum(K), Ts*tsub]);
92118
fin = imresize(fin, [options.nb, Ts*tsub]);
93119
if T ~= Ts*tsub

‎preprocess_data.m‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,7 @@
5353
fprintf('Estimating the noise power for each pixel from a simple PSD estimate...');
5454
[sn,psdx] = get_noise_fft(Y,options);
5555
P.sn = sn(:);
56+
%P.psdx = psdx;
5657
fprintf(' done \n');
5758

5859
% cluster pixels based on PSD

‎utilities/initialize_components_3d.m‎ renamed to ‎utilities/initialize_components_2d.m‎

Lines changed: 23 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function [Ain, Cin, bin, fin, center] = initialize_components_3d(Y, K, tau, options)
1+
function [Ain, Cin, bin, fin, center] = initialize_components_2d(Y, K, tau, options)
22

33
% Initalize components using a greedy approach followed by hierarchical
44
% alternative least squares (HALS) NMF. Optional use of spatio-temporal
@@ -51,70 +51,45 @@
5151
if ~isfield(options, 'tsub'), options.tsub = 1; end; tsub = options.tsub;
5252
if tsub == 1; fprintf('No temporal downsampling is performed. Consider temporal downsampling if the recording is very long. \n'); end
5353

54-
ndimsY = ndims(Y)-1;
55-
sY = size(Y);
56-
d = sY(1:ndimsY);
57-
T = sY(end);
58-
59-
ds = d;
60-
ds(1:2) = ceil(d(1:2)/ssub); % do not subsample along z axis
61-
%d1s = ceil(d1/ssub); %size of downsampled image
62-
%d2s = ceil(d2/ssub);
54+
[d1,d2,T] = size(Y);
55+
d1s = ceil(d1/ssub); %size of downsampled image
56+
d2s = ceil(d2/ssub);
6357
Ts = floor(T/tsub); %reduced number of frames
6458
% spatial downsampling
65-
fprintf('starting resampling \n')
66-
if ssub~=1;
67-
if ndimsY == 2; Y_ds = imresize(Y, [ds(1), ds(2)], 'box'); end
68-
if ndimsY == 3;
69-
Y_ds = zeros([ds(1:2),T,ds(end)]);
70-
for z = 1:ds(3)
71-
Y_ds(:,:,:,z) = imresize(squeeze(Y(:,:,z,:)), [ds(1), ds(2)], 'box');
72-
end
73-
Y_ds = permute(Y_ds,[1,2,4,3]);
74-
end
75-
else
76-
Y_ds = Y;
77-
end
59+
if ssub~=1, Y_ds = imresize(Y, [d1s, d2s], 'box'); else Y_ds = Y; end
7860
% temporal downsampling
7961
if tsub~=1
80-
if ndimsY == 2; Y_ds = squeeze(mean(reshape(Y_ds(:, :, 1:(Ts*tsub)),ds(1), ds(2), tsub, Ts), 3)); end
81-
if ndimsY == 3; Y_ds = squeeze(mean(reshape(Y_ds(:, :, :, 1:(Ts*tsub)),ds(1), ds(2), ds(3), tsub, Ts), 4)); end
62+
Y_ds = squeeze(mean(reshape(Y_ds(:, :, 1:(Ts*tsub)),d1s, d2s, tsub, Ts), 3));
8263
end
83-
64+
options_ds = options;
65+
options_ds.d1 = d1s;
66+
options_ds.d2 = d2s;
8467
if strcmpi(options.init_method,'greedy')
8568
% run greedy method
8669
fprintf('Initializing components with greedy method \n');
87-
if ndimsY == 2; [Ain, Cin, bin, fin] = greedyROI(Y_ds, K, options); end
88-
if ndimsY == 3; [Ain, Cin, bin, fin] = greedyROI(Y_ds, K, options); end %, gSig, gSiz); end
70+
[Ain, Cin, bin, fin] = greedyROI2d(Y_ds, K, options_ds);
8971
elseif strcmpi(options.init_method,'sparse_NMF')
9072
% run sparse_NMF method
9173
fprintf('Initializing components with sparse NMF \n');
92-
[Ain,Cin,bin,fin] = sparse_NMF_initialization(Y_ds,K,options);
74+
[Ain,Cin,bin,fin] = sparse_NMF_initialization(Y_ds,K,options_ds);
9375
else
9476
error('Unknown initialization method')
9577
end
9678

9779
% refine with HALS
98-
fprintf('Refining initial estimates with HALS...');
99-
%[Ain1, Cin1, bin1, fin1] = HALS_2d(Y_ds, full(Ain), Cin, bin, fin, options);
100-
[Ain, Cin, bin, fin] = HALS(Y_ds, full(Ain), Cin, bin, fin, options);
101-
fprintf(' done \n');
80+
fprintf('Refining initial estimates with HALS \n');
81+
[Ain, Cin, bin, fin] = HALS_2d(Y_ds, full(Ain), Cin, bin, fin, options);
82+
10283
%% upsample Ain, Cin, bin, fin
103-
if ndimsY == 2; center = ssub*com(Ain,ds(1),ds(2)); else center = ssub*com(Ain,ds(1),ds(2),ds(3)); end
104-
%Ain = imresize(reshape(Ain, ds(1), ds(2), sum(K)), d);
105-
%Ain = imresize(reshape(full(Ain), [ds, sum(K)]), d);
106-
Ain = imresize(reshape(full(Ain), [ds(1),ds(2), sum(K)*prod(ds)/ds(1)/ds(2)]),[d(1),d(2)]); %,prod(d)/d(1)/d(2)*sum(K)]);
107-
Ain = sparse(reshape(Ain, prod(d), []));
108-
%bin_temp = reshape(bin, ds(1), ds(2), options.nb);
109-
%bin = zeros(d(1),d(2),options.nb);
110-
bin = imresize(reshape(bin,[ds(1),ds(2), options.nb*prod(ds)/ds(1)/ds(2)]),[d(1),d(2)]);
111-
bin = reshape(bin,prod(d),[]);
112-
% bin_temp = reshape(bin, [ds, options.nb]);
113-
% bin = zeros([d,options.nb]);
114-
% for i = 1:options.nb
115-
% bin(:,:,i) = imresize(bin_temp(:,:,i), d);
116-
% end
117-
%bin = reshape(bin, [], options.nb);
84+
center = ssub*com(Ain,d1s,d2s);
85+
Ain = imresize(reshape(Ain, d1s, d2s, sum(K)), [d1, d2]);
86+
Ain = reshape(Ain, d1*d2, []);
87+
bin_temp = reshape(bin, d1s, d2s, options.nb);
88+
bin = zeros(d1,d2,options.nb);
89+
for i = 1:options.nb
90+
bin(:,:,i) = imresize(bin_temp(:,:,i), [d1, d2]);
91+
end
92+
bin = reshape(bin, [], options.nb);
11893
Cin = imresize(Cin, [sum(K), Ts*tsub]);
11994
fin = imresize(fin, [options.nb, Ts*tsub]);
12095
if T ~= Ts*tsub

0 commit comments

Comments
 (0)