-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbsoid_mt.m
126 lines (123 loc) · 9.45 KB
/
bsoid_mt.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
function [g_num,perc_unk,trans_p] = bsoid_mt(data,fps,pix_cm,smth_hstry,smth_futr,min_trans) % data{m}, pixel/cm
%BSOID_MT Classifies 3D pose estimation (DeepLabCut) output for rodent open field behavior using manual thresholding criteria.
%
% [G_LABEL,G_NUM,PERC_UNK] = BSOID_MT(DATA,PIX_CM) outputs classified behaviors based on DeepLabCut analysis
%
% INPUTS:
% DATA A cell-array containing either one or multiple rodents' estimated positions of the 6-body parts (snout, 4 paws, and base of tail)
% over time recorded from the bottom-up perspective. Rows represent chronological frame numbers. Columns represent x and y coordinates of
% individual body parts as follows:
% Columns 1 & 2 tracks snout; columns 3 to 6 tracks the two front paws; columns 7 to 10 tracks the two hind paws;
% columns 11 & 12 tracks the base of the tail. Tested on tracking data generated by DeepLabCut 2.0.
% FPS Rounded frame rate, can use VideoReader/ffmpeg(linux command) to automatically detect the input video fps.
% PIX_CM Pixel-to-cm conversion to adapt our manual thresholds to the user's set-up.
% SMTH_HSTRY BOXCAR smoothing using number of frames from before. Default ~40ms before.
% SMTH_FUTR BOXCAR smoothing using number of frames from after. Default ~40ms after.
% MIN_TRANS Minimum number of frames for transition to occur. Default ~100ms.
%
% G_NUM Classified behavioral state, numbered.
% PERC_UNK Percentage of total number of frames that are undefined.
% TRANS_P Transition probability matrix with Markov assumption (rows as current, columns as next).
%
% Created by Alexander Hsu, Date: 051519
% Contact [email protected]
if nargin < 3
error('Please input dataset, frame rate, & pixel-to-cm conversion!')
end
if nargin < 4
smth_hstry = round(0.05/(1/fps));
smth_futr = round(0.05/(1/fps));
end
if nargin < 6
min_trans = round(0.1/(1/fps))-1;
end
fprintf('Obtaining features from dataset... \n');
for m = 1:length(data) % For each csv file you uploaded.
data{m} = data{m}*24/pix_cm; % With our zoom and resolution, our threshold was set according to the ~24 pixel/cm set-up.
%% Obtain features
clear fpd_norm cfp_pt_norm sn_pt_norm sn_pt_ang sn_disp pt_disp hpR_disp hpL_disp fpR_disp fpL_disp;
cfp = [mean([data{m}(:,3),data{m}(:,5)],2),mean([data{m}(:,4),data{m}(:,6)],2)]; % Center of front paws x,y
cfp_pt = [cfp(:,1) - data{m}(:,11), cfp(:,2) - data{m}(:,12)]; % Center of front paw to proximal tail x,y
sn_pt = [data{m}(:,1) - data{m}(:,11), data{m}(:,2) - data{m}(:,12)]; % Snout to proximal tail x,y
for i = 1:length(data{m}) % Euclidean distance of x,y, since position means nothing
fpd_norm(i) = norm(data{m}(i,3:4)-data{m}(i,5:6)); % Front paw R to L euclidean distance
cfp_pt_norm(i) = norm(cfp_pt(i,:)); % Center of front paws to proximal tail euclidean distance
sn_pt_norm(i) = norm(sn_pt(i,:)); % Snout to proximal tail euclidean distance, i.e. body length
end
fpd_norm_smth{m} = movmean(fpd_norm,[smth_hstry,smth_futr]); % Reduce label noise
sn_cfp_norm_smth{m} = movmean(sn_pt_norm-cfp_pt_norm,[smth_hstry,smth_futr]); % Reduce label noise
sn_pt_norm_smth{m} = movmean(sn_pt_norm,[smth_hstry,smth_futr]); % Reduce label noise
for k = 1:length(data{m})-1 % Velocity and angle over time
b_3d = [sn_pt(k+1,:),0]; a_3d = [sn_pt(k,:),0]; c = cross(b_3d,a_3d);
sn_pt_ang(k) = sign(c(3))*180/pi*atan2(norm(c),dot(sn_pt(k,:),sn_pt(k+1,:))); % Body angle, arctan between body
sn_disp(k) = norm(data{m}(k+1,1:2)-data{m}(k,1:2)); % Snout displacement over time
pt_disp(k) = norm(data{m}(k+1,11:12)-data{m}(k,11:12)); % Proximal tail displacement over time
hpR_disp(k) = norm(data{m}(k+1,7:8)-data{m}(k,7:8)); % Right hind paw displacement over time
hpL_disp(k) = norm(data{m}(k+1,9:10)-data{m}(k,9:10)); % Left hind paw displacement over time
fpR_disp(k) = norm(data{m}(k+1,3:4)-data{m}(k,3:4)); % Right forepaw displacement over time
fpL_disp(k) = norm(data{m}(k+1,5:6)-data{m}(k,5:6)); % Left forepaw displacement over time
end
sn_pt_ang_smth{m} = movmean(sn_pt_ang,[smth_hstry,smth_futr]); % Reduce label noise
sn_disp_smth{m} = movmean(sn_disp,[smth_hstry,smth_futr]); % Reduce label noise
pt_disp_smth{m} = movmean(pt_disp,[smth_hstry,smth_futr]); % Reduce label noise
hpR_disp_smth{m} = movmean(hpR_disp,[smth_hstry,smth_futr]); % Reduce label noise
hpL_disp_smth{m} = movmean(hpL_disp,[smth_hstry,smth_futr]); % Reduce label noise
fpR_disp_smth{m} = movmean(fpR_disp,[smth_hstry,smth_futr]); % Reduce label noise
fpL_disp_smth{m} = movmean(fpL_disp,[smth_hstry,smth_futr]); % Reduce label noise
%% Classify action based on parameters specified above
%%% Rest/Pause is defined as minimal change in body angle, and minimal displacement of all labels.
idx_rest{m} = find(abs(sn_pt_ang_smth{m}) < 1 & sn_disp_smth{m} > 0 & sn_disp_smth{m} < 2 & pt_disp_smth{m} < 2 & ...
fpR_disp_smth{m} < 2 & fpL_disp_smth{m} < 2 & sn_cfp_norm_smth{m}(:,2:end) > 0);
%%% Rear is defined as either snout is gone, front paws are gone (reared up), or front paws coming together (rearing up).
idx_rear{m} = find(sn_disp_smth{m} == 0 & abs(sn_pt_ang_smth{m}) < 1 | ...
fpR_disp_smth{m} == 0 & abs(sn_pt_ang_smth{m}) < 1 | ...
fpL_disp_smth{m} == 0 & abs(sn_pt_ang_smth{m}) < 1 | ...
sn_disp_smth{m} > 2 & abs(fpR_disp_smth{m}-fpL_disp_smth{m}) < 5 & sn_cfp_norm_smth{m}(:,2:end) > 0 & ...
fpd_norm_smth{m}(:,2:end) < 20 & pt_disp_smth{m} >= 2 & fpR_disp_smth{m} >= 2 & fpL_disp_smth{m} >= 2);
%%% Groom is defined as either snout invisible, front paws visible and close, and no change in angle or tail position
%%%% snout is visible, front paws in/visible and far, hands are in front of snout, and no change in angle or tail position
%%%% or snout visibly moving, both front paws visible, one of which moves, and no change in angle or tail position
idx_groom{m} = find(sn_disp_smth{m} == 0 & fpd_norm_smth{m}(:,2:end) < 10 & fpR_disp_smth{m} >= 2 & abs(sn_pt_ang_smth{m}) < 1 & ...
pt_disp_smth{m} < 2 | sn_disp_smth{m} == 0 & fpd_norm_smth{m}(:,2:end) < 10 & fpL_disp_smth{m} >= 2 & ...
abs(sn_pt_ang_smth{m}) < 1 & pt_disp_smth{m} < 2 | ...
sn_disp_smth{m} > 0 & fpR_disp_smth{m} >= 2 & abs(sn_pt_ang_smth{m}) < 1 & hpR_disp_smth{m} < 2 & hpL_disp_smth{m} < 2 & pt_disp_smth{m} < 2 | ...
sn_disp_smth{m} > 0 & fpL_disp_smth{m} >= 2 & abs(sn_pt_ang_smth{m}) < 1 & hpR_disp_smth{m} < 2 & hpL_disp_smth{m} < 2 & pt_disp_smth{m} < 2 | ...
sn_disp_smth{m} > 0 & abs(sn_pt_ang_smth{m}) < 1 & sn_cfp_norm_smth{m}(:,2:end) < 0 & pt_disp_smth{m} < 2 | ...
sn_disp_smth{m} >= 2 & fpR_disp_smth{m} > 0 & fpL_disp_smth{m} > 0 & abs(sn_pt_ang_smth{m}) < 1 & pt_disp_smth{m} < 2);
%%% Orientation of body is defined as non-locomoting body angle change
idx_ori_right{m} = find(pt_disp_smth{m} < 2 & sn_pt_ang_smth{m} >= 1);
idx_ori_left{m} = find(pt_disp_smth{m} < 2 & sn_pt_ang_smth{m} <= -1);
%%% Locomotion is defined as body displacement with either front paw
idx_loc{m} = find(pt_disp_smth{m} >= 2 & fpR_disp_smth{m} >= 2 | pt_disp_smth{m} >= 2 & fpL_disp_smth{m} >= 2 | ...
pt_disp_smth{m} >= 2 & hpR_disp_smth{m} >= 2 | pt_disp_smth{m} >= 2 & hpL_disp_smth{m} >= 2);
%%% Head movement is defined as only snout movement with stationary body, and no change in angle or tail position
idx_headmov{m} = find(abs(sn_pt_ang_smth{m}) < 1 & sn_disp_smth{m} >= 2 & hpR_disp_smth{m} < 2 & hpL_disp_smth{m} < 2);
%% Low pass filter of unreasonable fast state changes
MC = []; MC(:) = 0; MC(idx_rest{m},1) = 1; MC(idx_loc{m},1) = 5; MC(idx_ori_left{m},1) = 6; MC(idx_ori_right{m},1) = 7;
MC(idx_rear{m},1) = 2; MC(idx_groom{m},1) = 3; MC(idx_headmov{m},1) = 4;
stch = find(diff(MC)~=0); % Remove short transitions
for i = 1:length(stch)-1
if stch(i+1) - stch(i) <= min_trans
MC(stch(i):stch(i+1)) = MC(stch(i));
end
end
%% Reassign sandwiched unidentified behavior to prior and posterior, if the same
for k = 1:length(stch)-2
if MC(stch(k+1)) == 0 && MC(stch(k)) == MC(stch(k+2))
MC(stch(k):stch(k+2)) = MC(stch(k));
end
end
g_num{m} = MC; % Store behavior per .csv
%% Pull out transitional probability
perc_unk{m} = length(find(g_num{m}==0))/(length(data{m})-1); % Undefined percentage
binEdges = 1:length(MC); [~,states] = histc(nonzeros(MC),binEdges); % Histc just asks which bin each pose belongs to
binTicks = binEdges(1:end-1); % Don't need the last bin edge
ActMat = zeros(length(unique(MC))-1); fromBinNos = states(1:end-1); nextBinNos = states(2:end);
for i = 1:length(fromBinNos)
ActMat(fromBinNos(i),nextBinNos(i)) = ActMat(fromBinNos(i),nextBinNos(i)) + 1;
end
for j = 1:length(unique(MC))-1
trans_p{m}(j,:) = ActMat(j,:)/sum(ActMat(j,:));
end
end
return