diff --git a/CNMFSetParms.m b/CNMFSetParms.m index 706987a..f15fc1c 100644 --- a/CNMFSetParms.m +++ b/CNMFSetParms.m @@ -10,6 +10,7 @@ % dataset info 'd1 ' % number of rows 'd2 ' % number of cols + 'd3 ' % number of planes (for 3d imaging, default: 1) % INITIALIZATION (initialize_components.m) 'ssub ' % spatial downsampling factor (default: 1) 'tsub ' % temporal downsampling factor (default: 1) @@ -164,6 +165,7 @@ % dataset info {[]} {[]} + {1} % INITIALIZATION (initialize_components.m) {1} {1} diff --git a/utilities/HALS.m b/utilities/HALS.m new file mode 100644 index 0000000..b77984a --- /dev/null +++ b/utilities/HALS.m @@ -0,0 +1,86 @@ +function [A, C, b, f] = HALS(Y, A, C, b, f, params) +%% Hierarchical alternating least square method for solving NMF problem +% Y = A*C + b*f + +%input: +% Y: d1 X d2 X T, raw data. It will be reshaped to (d1*d2) X T in this +% function +% A: (d1*d2) X K, initial value of spatial components +% C: K X T, initial value of temporal components +% b: (d1*d2) X 1, initial value of background spatial component +% f: 1 X T, initial value of background temporal component +% params: parameters used in this function. +% bSiz: blur size. A box kernel (bSiz X bSiz) will be convolved +% with each neuron's initial spatial component, then all nonzero +% pixels will be picked as pixels to be updated, and the rest will be +% forced to be 0. +% maxIter: maximum iteration of iterating HALS. + +% Author: Pengcheng Zhou, Columbia University, based on a python +% implementation from Johannes Friedrich, Columbia University, 2015. + +%[d1, d2, ~] = size(Y); +dimY = ndims(Y)-1; +sizY = size(Y); +d = sizY(1:dimY); + +if isfield(params, 'bSiz'), bSiz = params.bSiz; else bSiz=3; end +if isfield(params, 'maxIter'), maxIter = params.maxIter; else maxIter=5; end + +blur_kernel = ones(bSiz); +if dimY == 3; blur_kernel(:,:,3:end) = []; end +blur_kernel = blur_kernel/numel(blur_kernel); +if ismatrix(A) + A = reshape(A,[d,size(A,2)]); + ind_A = (imfilter(A,blur_kernel)>0); +end +A = sparse(reshape(A,prod(d),[])); +ind_A = sparse(reshape(ind_A,prod(d),[])); +K = size(A,2); +% if ndims(A)==3 +% ind_A = (imfilter(A, blur_kernel)>0); +% A = reshape(A, d1*d2, []); +% ind_A = reshape(ind_A, d1*d2, []); +% else +% ind_A = (imfilter(reshape(A, d1,d2,[]), blur_kernel)>0); +% ind_A = reshape(ind_A, d1*d2, []); +% end +% ind_A = sparse(ind_A); %indicator of nonnero pixels +% K = size(A, 2); %number of neurons + +%% update spatial and temporal components neuron by neurons +%Yres = reshape(Y, d1*d2, []) - A*C - b*f; +Yres = reshape(Y, prod(d), []) - A*C - b*f; + +for miter=1:maxIter + for mcell = 1:K + ind_pixels = ind_A(:, mcell); + tmp_Yres = Yres(ind_pixels, :); + + % update temporal component of the neuron + c0 = C(mcell, :); + a0 = A(ind_pixels, mcell); + norm_a2 = norm(a0, 2)^2; + C(mcell, :) = max(0, c0 + a0'*tmp_Yres/norm_a2); + tmp_Yres = tmp_Yres + a0*(c0-C(mcell,:)); + + % update spatial component of the neuron + norm_c2 = norm(C(mcell,:),2)^2; + tmp_a = max(0, a0 + tmp_Yres*C(mcell, :)'/norm_c2); + A(ind_pixels, mcell) = tmp_a; + Yres(ind_pixels,:) = tmp_Yres + (a0-tmp_a)*C(mcell,:); + end + + % update temporal component of the background + f0 = f; + b0 = b; + norm_b2 = norm(b0,2)^2; + f = max(0, f0 + b0'*Yres/norm_b2); + Yres = Yres + b0*(f0-f); + % update spatial component of the background + norm_f2 = norm(f, 2)^2; + b = max(0, b0 + Yres*f'/norm_f2); + Yres = Yres + (b0-b)*f; + + %fprintf('Iteration %d, the norm of residual is %.4f\n', miter, norm(Yres, 'fro')); +end \ No newline at end of file diff --git a/utilities/com.m b/utilities/com.m index 8009063..e7484c7 100644 --- a/utilities/com.m +++ b/utilities/com.m @@ -11,10 +11,15 @@ if nargin < 4 d3 = 1; end +if d3 == 1 + ndim = 2; +else + ndim = 3; +end nr = size(A,2); Coor.x = kron(ones(d2*d3,1),(1:d1)'); Coor.y = kron(ones(d3,1),kron((1:d2)',ones(d1,1))); Coor.z = kron((1:d3)',ones(d2*d1,1)); cm = [Coor.x, Coor.y, Coor.z]'*A/spdiags(sum(A)',0,nr,nr); -cm = cm(1:nargin-1,:)'; \ No newline at end of file +cm = cm(1:ndim,:)'; \ No newline at end of file diff --git a/utilities/determine_search_location.m b/utilities/determine_search_location.m index 6ae2dd4..4fe17ba 100644 --- a/utilities/determine_search_location.m +++ b/utilities/determine_search_location.m @@ -28,12 +28,18 @@ end [d,nr] = size(A); + if ~isfield(params,'d1'); d1 = sqrt(d); params.d1 = d1; else d1 = params.d1; end % # of rows if ~isfield(params,'d2'); d2 = sqrt(d); params.d2 = d2; else d2 = params.d2; end % # of columns +if ~isfield(params,'d3'); d3 = 1; params.d3 = d3; else d3 = params.d3; end if ~isfield(params,'min_size') || isempty(params.min_size); min_size = 3; else min_size = params.min_size; end % minimum size of ellipse axis if ~isfield(params,'max_size') || isempty(params.max_size); max_size = 8; else max_size = params.max_size; end % maximum size of ellipse axis if ~isfield(params,'dist') || isempty(params.dist); dist = 3; else dist = params.dist; end % expansion factor of ellipse -if ~isfield(params,'se') || isempty(params.se); expandCore = strel('disk',4,0); else expandCore = params.se; end % morphological element (for 'dilate') +if ~isfield(params,'se') || isempty(params.se); % morphological element (for 'dilate') + if d3 == 1; expandCore = strel('disk',4,0); else expandCore = strel(ones(4,4,2)); end +else + expandCore = params.se; +end if strcmpi(method,'ellipse'); method = 'ellipse'; elseif strcmpi(method,'dilate'); method = 'dilate'; @@ -43,21 +49,37 @@ IND = false(d,nr); switch method case 'ellipse' - Coor.x = kron(ones(d2,1),(1:d1)'); - Coor.y = kron((1:d2)',ones(d1,1)); + Coor.x = kron(ones(d2*d3,1),(1:d1)'); + Coor.y = kron(ones(d3,1),kron((1:d2)',ones(d1,1))); + Coor.z = kron((1:d3)',ones(d1*d2,1)); if ~(dist==Inf) % determine search area for each neuron - cm = zeros(nr,2); % vector for center of mass + %cm = zeros(nr,2); % vector for center of mass + cm = com(A,d1,d2,d3); + if d3 == 1 + cm = [cm,ones(nr,1)]; + end Vr = cell(nr,1); IND = zeros(d,nr); % indicator for distance - cm(:,1) = Coor.x'*A(:,1:nr)./sum(A(:,1:nr)); - cm(:,2) = Coor.y'*A(:,1:nr)./sum(A(:,1:nr)); % center of mass for each components + %cm(:,1) = Coor.x'*A(:,1:nr)./sum(A(:,1:nr)); + %cm(:,2) = Coor.y'*A(:,1:nr)./sum(A(:,1:nr)); % center of mass for each components for i = 1:nr % calculation of variance for each component and construction of ellipses - Vr{i} = ([Coor.x - cm(i,1), Coor.y - cm(i,2)]'*spdiags(A(:,i),0,d,d)*[Coor.x - cm(i,1), Coor.y - cm(i,2)])/sum(A(:,i)); - [V,D] = eig(Vr{i}); - cor = [Coor.x - cm(i,1),Coor.y - cm(i,2)]; + if d3 == 1 + Vr{i} = ([Coor.x - cm(i,1), Coor.y - cm(i,2)]'*spdiags(A(:,i),0,d,d)*[Coor.x - cm(i,1), Coor.y - cm(i,2)])/sum(A(:,i)); + [V,D] = eig(Vr{i}); + cor = [Coor.x - cm(i,1),Coor.y - cm(i,2)]; + else + Vr{i} = ([Coor.x - cm(i,1), Coor.y - cm(i,2), Coor.z - cm(i,3)]'*spdiags(A(:,i),0,d,d)*[Coor.x - cm(i,1), Coor.y - cm(i,2), Coor.z - cm(i,3)])/sum(A(:,i)); + [V,D] = eig(Vr{i}); + cor = [Coor.x - cm(i,1),Coor.y - cm(i,2),Coor.z - cm(i,3)]; + end d11 = min(max_size^2,max(min_size^2,D(1,1))); - d22 = min(max_size^2,max(min_size^2,D(2,2))); - IND(:,i) = sqrt((cor*V(:,1)).^2/d11 + (cor*V(:,2)).^2/d22)<=dist; % search indeces for each component + d22 = min(max_size^2,max(min_size^2,D(2,2))); + if d3 == 1 + IND(:,i) = sqrt((cor*V(:,1)).^2/d11 + (cor*V(:,2)).^2/d22)<=dist; + else + d33 = min((max_size/2)^2,max((min_size/2)^2,D(3,3))); + IND(:,i) = sqrt((cor*V(:,1)).^2/d11 + (cor*V(:,2)).^2/d22 + (cor*V(:,3)).^2/d33)<=dist; % search indeces for each component + end end else IND = true(d,nr); diff --git a/utilities/determine_search_location_2d.m b/utilities/determine_search_location_2d.m new file mode 100644 index 0000000..9d40dfa --- /dev/null +++ b/utilities/determine_search_location_2d.m @@ -0,0 +1,73 @@ +function IND = determine_search_location_2d(A,method,params) + +% determine the search location for updating each spatial component +% INPUTS: +% A: current estimate of spatial component ( d x nr sparse matrix) +% method: method to be used of determining search locations +% available methods: +% 'ellipse': draw an ellipse centered at the center of mass with +% axes the first two principal components and expand it to +% determine the search location (default) +% 'dilate' : grow the spatial footprint by dilating the current +% footprint +% if a different argument is passed then search location covers the whole field of view +% params: hyper-parameter struct for the various methods (see default settings below for description) +% +% OUTPUT: +% IND: binary d x nr matrix. IND(i,j) = 1 if pixel i is included in the search location of component j +% +% Written by Eftychios A. Pnevmatikakis, Simons Foundation +% with input from Weijian Yang, Columbia University + +if nargin < 3 + params = []; +end + +if nargin < 2 || isempty(method) + method = 'ellipse'; % default method +end + +[d,nr] = size(A); +if ~isfield(params,'d1'); d1 = sqrt(d); params.d1 = d1; else d1 = params.d1; end % # of rows +if ~isfield(params,'d2'); d2 = sqrt(d); params.d2 = d2; else d2 = params.d2; end % # of columns +if ~isfield(params,'min_size') || isempty(params.min_size); min_size = 3; else min_size = params.min_size; end % minimum size of ellipse axis +if ~isfield(params,'max_size') || isempty(params.max_size); max_size = 8; else max_size = params.max_size; end % maximum size of ellipse axis +if ~isfield(params,'dist') || isempty(params.dist); dist = 3; else dist = params.dist; end % expansion factor of ellipse +if ~isfield(params,'se') || isempty(params.se); expandCore = strel('disk',4,0); else expandCore = params.se; end % morphological element (for 'dilate') + +if strcmpi(method,'ellipse'); method = 'ellipse'; +elseif strcmpi(method,'dilate'); method = 'dilate'; +else fprintf('Method not recongnized. Search location equals the entire field of view. \n'); +end + +IND = false(d,nr); +switch method + case 'ellipse' + Coor.x = kron(ones(d2,1),(1:d1)'); + Coor.y = kron((1:d2)',ones(d1,1)); + if ~(dist==Inf) % determine search area for each neuron + cm = zeros(nr,2); % vector for center of mass + Vr = cell(nr,1); + IND = zeros(d,nr); % indicator for distance + cm(:,1) = Coor.x'*A(:,1:nr)./sum(A(:,1:nr)); + cm(:,2) = Coor.y'*A(:,1:nr)./sum(A(:,1:nr)); % center of mass for each components + for i = 1:nr % calculation of variance for each component and construction of ellipses + Vr{i} = ([Coor.x - cm(i,1), Coor.y - cm(i,2)]'*spdiags(A(:,i),0,d,d)*[Coor.x - cm(i,1), Coor.y - cm(i,2)])/sum(A(:,i)); + [V,D] = eig(Vr{i}); + cor = [Coor.x - cm(i,1),Coor.y - cm(i,2)]; + d11 = min(max_size^2,max(min_size^2,D(1,1))); + d22 = min(max_size^2,max(min_size^2,D(2,2))); + IND(:,i) = sqrt((cor*V(:,1)).^2/d11 + (cor*V(:,2)).^2/d22)<=dist; % search indeces for each component + end + else + IND = true(d,nr); + end + case 'dilate' + A = threshold_components(A,params); + for i = 1:nr + A_temp = imdilate(reshape(full(A(:,i)),d1,d2),expandCore); + IND(:,i) = A_temp(:)>0; + end + otherwise + IND = true(d,nr); +end diff --git a/utilities/threshold_components.m b/utilities/threshold_components.m index 698030b..8e79565 100644 --- a/utilities/threshold_components.m +++ b/utilities/threshold_components.m @@ -10,27 +10,32 @@ % Written by: % Eftychios A. Pnevmatikakis, Simons Foundation, 2015 - defoptions.nrgthr = 0.9999; % energy threshold + defoptions.nrgthr = 0.9999; % energy threshold defoptions.clos_op = strel('square',3); % morphological operator for closing - defoptions.medw = [3,3]; % size of median filter + defoptions.medw = [3,3]; % size of median filter if ~isfield(options,'nrgthr') || isempty(options.nrgthr); options.nrgthr = defoptions.nrgthr; end if ~isfield(options,'clos_op') || isempty(options.clos_op); options.clos_op = defoptions.clos_op; end if ~isfield(options,'medw') || isempty(options.medw); options.medw = defoptions.medw; end - + if ~isfield(options,'d3') || isempty(options.d3); options.d3 = 1; end + [d,nr] = size(A); Ath = spalloc(d,nr,nnz(A)); for i = 1:nr - A_temp = full(reshape(A(:,i),options.d1,options.d2)); - A_temp = medfilt2(A_temp,options.medw); + A_temp = reshape(full(A(:,i)),options.d1,options.d2,options.d3); + for z = 1:options.d3 + A_temp(:,:,z) = medfilt2(A_temp(:,:,z),options.medw); + end A_temp = A_temp(:); [temp,ind] = sort(A_temp(:).^2,'ascend'); temp = cumsum(temp); ff = find(temp > (1-options.nrgthr)*temp(end),1,'first'); - BW = zeros(options.d1,options.d2); + BW = zeros(options.d1,options.d2,options.d3); BW(ind(ff:d)) = 1; - BW = imclose(BW,options.clos_op); - [L,NUM] = bwlabel(BW,8); + for z = 1:options.d3 + BW(:,:,z) = imclose(BW(:,:,z),options.clos_op); + end + [L,NUM] = bwlabeln(BW,8*(options.d3==1) + 6*(options.d3~=1)); if NUM > 0 nrg = zeros(NUM,1); for l = 1:NUM diff --git a/utilities/threshold_components_2d.m b/utilities/threshold_components_2d.m new file mode 100644 index 0000000..0f9a565 --- /dev/null +++ b/utilities/threshold_components_2d.m @@ -0,0 +1,45 @@ +function Ath = threshold_components_2d(A,options) + +% post processing of spatial components +% for each component perform the following: +% (i) perform median filtering +% (ii) keep only pixels that contibute up to a level of total energy +% (iii) perform morphological closing +% (iv) extract largest connected component + +% Written by: +% Eftychios A. Pnevmatikakis, Simons Foundation, 2015 + + defoptions.nrgthr = 0.9999; % energy threshold + defoptions.clos_op = strel('square',3); % morphological operator for closing + defoptions.medw = [3,3]; % size of median filter + + if ~isfield(options,'nrgthr') || isempty(options.nrgthr); options.nrgthr = defoptions.nrgthr; end + if ~isfield(options,'clos_op') || isempty(options.clos_op); options.clos_op = defoptions.clos_op; end + if ~isfield(options,'medw') || isempty(options.medw); options.medw = defoptions.medw; end + + [d,nr] = size(A); + Ath = spalloc(d,nr,nnz(A)); + for i = 1:nr + A_temp = full(reshape(A(:,i),options.d1,options.d2)); + A_temp = medfilt2(A_temp,options.medw); + A_temp = A_temp(:); + [temp,ind] = sort(A_temp(:).^2,'ascend'); + temp = cumsum(temp); + ff = find(temp > (1-options.nrgthr)*temp(end),1,'first'); + BW = zeros(options.d1,options.d2); + BW(ind(ff:d)) = 1; + BW = imclose(BW,options.clos_op); + [L,NUM] = bwlabel(BW,8); + if NUM > 0 + nrg = zeros(NUM,1); + for l = 1:NUM + ff = find(L==l); + nrg(l) = sum(A(ff,i).^2); + end + [~,indm] = max(nrg); + ff = find(L==indm); + Ath(ff,i) = A(ff,i); + end + end +end \ No newline at end of file