-
Notifications
You must be signed in to change notification settings - Fork 0
/
groupingLAP.m
177 lines (152 loc) · 6.15 KB
/
groupingLAP.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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
function group = groupingLAP(fileName,dirName)
% pass the TRACKS struct as input
nbFrames = 120;
LifeTime = 4; %4
alpha = 2; % default 2, shrinkage speed factor->(1,2]
MaxTimeSpan = 30;%25; %15
delta_t_max = sqrt(MaxTimeSpan); % 4 or more?
coneAngPos = cos(pi/4);%to 15 from 45 (/12 from /4) /4=default
coneAngNeg = cos(pi/18);%10 - /18=default
trackAng = cos(pi/3);%to 15 from 60 (/12 from /3) /3=default
constVel = 1; %debug100
duConst = 1000000;
if nargin == 0
[fileName,dirName] = uigetfile('*.mat','Choose a .mat file');
load([dirName,filesep,fileName]);
else
load([dirName,filesep,fileName]);
end
%check wether the "groups" subdirectory exists or not
[success, msg, msgID] = mkdir(dirName(1:end-12), 'groups1'); % default 'groups'
if (success ~= 1)
error(msgID, msg);
elseif (~isempty(msg))
fprintf('Directory "groups" already exists.\n');
else
fprintf('Directory "groups" has been created.\n');
end
indx = find( [tracks.len] >= LifeTime);
TimeSpan = MaxTimeSpan;%min(max([tracks.len]),MaxTimeSpan);
traj = tracks(indx);
leIndx = length(indx);
% calculate trajectory information
for i = 1:leIndx
traj(i).endID = traj(i).startID + traj(i).len - 1;
dY = traj(i).points(end,1) - traj(i).points(1,1);
dX = traj(i).points(end,2) - traj(i).points(1,2);
traj(i).vec = [dX; dY];
traj(i).vel = sqrt(dY^2+dX^2)/(traj(i).len-1);
end
% exclude tracks with NO motion
traj = traj(find([traj.vel]>0)); % SOME DONT MOVE!?
leIndx = length(traj);%fiNa = 'drop30a.xlsx';xlswrite(fiNa, [traj.vel])
traj_vel = sort([traj.vel]);
percentileVel = 0.95;+6
v_max = traj_vel(round(length([traj.vel])*percentileVel)); % pixels
v_med = median(traj_vel);
distance_max = v_med * MaxTimeSpan;
% get start & end points of all tracks for distance matrix calculation
for i = 1:leIndx
E(i,:)=traj(i).points(end,1:2);
B(i,:)=traj(i).points(1,1:2);
end
[listStartID,indxStartID] = sort([traj.startID]);
M = sparse(leIndx,leIndx); % angle multiplication matrix
C = sparse(leIndx,leIndx); % distance cut-off matrix
progressText(0,'First FOR-loop') % Create text
for i = 1:leIndx
if TimeSpan > nbFrames-traj(i).endID
T = nbFrames-traj(i).endID;
else
T = TimeSpan;
end
if (traj(i).endID < (nbFrames - LifeTime))
aux1 = find(listStartID>(traj(i).endID+1));% skip next frame
if ~isempty(aux1)
firstIndx = aux1(1);
aux2 = find(listStartID>(traj(i).endID+T+1));
if ~isempty(aux2)
finalIndx = aux2(1)-1;
else
finalIndx = length(listStartID);
end
a = indxStartID(firstIndx:finalIndx);
magTrajVec = sqrt(sum(traj(i).vec.^2,1));
magTrajA = sqrt(sum([traj(a).vec]'.^2,2));
cos_trackAng = ([traj(a).vec]'*traj(i).vec)./(magTrajA*magTrajVec);
listA = find(cos_trackAng>trackAng);
aa = a(listA);
cos_trA = cos_trackAng(listA);
if ~isempty(aa)
dxCo = []; dyCo = [];
for j = 1:length(aa)
dyCo(j) = traj(aa(j)).points(1,1) - traj(i).points(end,1); % get it out of the loop - to the other loop up
dxCo(j) = traj(aa(j)).points(1,2) - traj(i).points(end,2);
end
aVecs = ([dxCo; dyCo])';
magAvecs = sqrt(sum(aVecs.^2,2));
magAvecs(find(magAvecs==0)) = 0.001; % some new/t+4 traj begin where an old/t one ended
cos_coneAng = (aVecs*traj(i).vec)./(magAvecs*magTrajVec);
posC = find(cos_coneAng>coneAngPos);
negC = find(cos_coneAng<-coneAngNeg);
if ~isempty(posC)
% C(i,aa(posC)) = traj(i).vel * sqrt([traj(aa(posC)).startID]-traj(i).endID); %* constVel;
delta_t_gr = [traj(aa(posC)).startID]-traj(i).endID;
C(i,aa(posC)) = v_max * min( delta_t_max , delta_t_gr ); % delta_t_max = sqrt(MaxTimeSpan)
M(i,aa(posC)) = abs(cos_trA(posC) - cos_coneAng(posC));
end
if ~isempty(negC)
% C(i,aa(negC)) = traj(i).vel * sqrt([traj(aa(negC)).startID]-traj(i).endID); %* constVel;
delta_t_sh = [traj(aa(negC)).startID]-traj(i).endID;
C(i,aa(negC)) = min( distance_max , v_max * alpha * delta_t_sh ); % distance_max = v_med * MaxTimeSpan
M(i,aa(negC)) = abs(cos_trA(negC) - abs(cos_coneAng(negC)));
end
end
end
end
progressText(i/leIndx);
end
R = max(C(:));
% D=createSparseDistanceMatrix(E,B,R); % dimentions - leIndx
D=createDistanceMatrix(E,B); % dimentions - leIndx
radIndx = find(D<=C);
clear C;
C = sparse(leIndx,leIndx); % distance cut-off matrix
C(radIndx) = D(radIndx).*M(radIndx); % overwrite old cut-off matrix to generate cost matrix
[links12, links21] = lap(C,-1,0,1);
lnk = links12(1:leIndx);
lnkIndx = find(lnk<=leIndx);
leLnkIndx = length(lnkIndx);
group=struct('list',[]);
assocTracks = zeros(leIndx,1);
grNb = 0;
progressText(0,'Second FOR-loop') % Create text
for i = 1:leLnkIndx
k = lnkIndx(i);
if assocTracks(k) == 0
if assocTracks(lnk(k)) == 0
grNb = grNb + 1;
assocTracks(k) = grNb;
assocTracks(lnk(k)) = grNb;
group(grNb).list = [k lnk(k)];
else
assocTracks(k) = assocTracks(lnk(k));
group(assocTracks(lnk(k))).list = [k,group(assocTracks(lnk(k))).list];
end
else
if assocTracks(lnk(k)) == 0
assocTracks(lnk(k)) = assocTracks(k);
group(assocTracks(k)).list = [group(assocTracks(k)).list,lnk(k)];
else
% there are two groups to be linked together
group(assocTracks(k)).list = [group(assocTracks(k)).list,group(assocTracks(lnk(k))).list];
% update all members of group(assocTracks(link(k))).list in traj
group(assocTracks(lnk(k))).list = [];
end
end
progressText(i/leLnkIndx);
end
% s = 2;
% strg=sprintf('%%.%dd',s);
% indxStr=sprintf(strg,10*alpha);
save([dirName(1:end-12),filesep,'groups1',filesep,'group'],'group') % -13 % default 'groups'