-
Notifications
You must be signed in to change notification settings - Fork 14
/
spTop.m
86 lines (70 loc) · 2.21 KB
/
spTop.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
function [svLst,dReconSp,riseX] = spTop(dat,dF,lmLoc,evtSpatialMask,opts,ff)
% in case starting, or ending signals are not detected
ext = 2;
[H,W,T] = size(dat);
backGround = mean(dat-dF,3);
dat2 = zeros(H,W,T+2*ext,'single');
dat2(:,:,ext+1:end-ext) = dat;
for i = 1:ext
dat2(:,:,i) = backGround;
dat2(:,:,end-i+1) = backGround;
end
dat = dat2; clear dat2;
dF2 = zeros(H,W,T+2*ext,'single');
dF2(:,:,ext+1:end-ext) = dF;
dF = dF2; clear dF2;
T = T + 2*ext;
lmLoc = lmLoc + H*W*ext;
dReconSp = [];
if isempty(evtSpatialMask)
evtSpatialMask = ones(H,W);
end
if ~isfield(opts,'thrSvSig')
opts.thrSvSig = 4;
end
% grow seeds
resCell = cell(numel(lmLoc),1);
lblMap = zeros(H,W,T,'uint32');
opts1 = opts;
opts1.maxStp = 1;
lmAll = zeros(H,W,T,'logical');
lmAll(lmLoc) = true;
[resCell,lblMap] = burst.growSeed(dF,dF,resCell,lblMap,lmLoc,lmAll,opts1,0);
for pp=1:40
fprintf('Grow %d\n',pp);
if exist('ff','var')
waitbar(0.1+pp/80,ff);
end
[resCell,lblMap] = burst.growSeed(dF,dF,resCell,lblMap,lmLoc,lmAll,opts1,pp);
end
% clean super voxels
fprintf('Cleaning super voxels by size ...\n')
lblMap = lblMap.*uint32(evtSpatialMask);
lblMap = burst.filterAndFillSp(lblMap);
if exist('ff','var')
waitbar(0.7,ff);
end
fprintf('Cleaning super voxels by z score...\n')
zVec1 = stat.getSpZ(dat,lblMap,opts.varEst);
spLst = label2idx(lblMap);
spLst = spLst(zVec1>opts.thrSvSig);
lblMap = zeros(size(lblMap),'uint32');
for nn=1:numel(spLst)
lblMap(spLst{nn}) = nn;
end
if exist('ff','var')
waitbar(0.8,ff);
end
% remove extension
lblMap = lblMap(:,:,ext+1:end-ext);
dat = dat(:,:,ext+1:end-ext);
% Extend and re-fit each patch, estimate delay, reconstruct signal
fprintf('Extending super voxels ...\n')
[lblMap,riseX] = burst.getSpDelay(dat,lblMap,opts);
if exist('ff','var')
waitbar(1,ff);
end
% poolobj = gcp('nocreate');
% delete(poolobj);
svLst = label2idx(lblMap);
end