-
Notifications
You must be signed in to change notification settings - Fork 148
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
7 changed files
with
258 additions
and
20 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |