-
Notifications
You must be signed in to change notification settings - Fork 0
/
buildNewContactArrayMulti.m
157 lines (113 loc) · 7.71 KB
/
buildNewContactArrayMulti.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
function [contacts, params] = buildNewContactArray(array, params)
% Version 0.5.0 SAH 110827
if nargin == 1;
disp('Build new contact array using default analysis parameters')
% create default parameters
contacts=cell(1,length(array.trials));
params=struct;
params.touchThresh = [.1 .1 .3 .3]; %Touch threshold for go (protraction, retraction), no-go (protraction,retraction). Check with Parameter Estimation cell
params.goProThresh = 0; % Mean curvature above this value indicates probable go protraction, below it, a go retraction trial.
params.nogoProThresh = 0; % Mean curvature above this value indicates probable nogo protraction, below it, a nogo retraction trial.
params.poleOffset = .535; % Time where pole becomes accessible
params.poleEndOffset = .13; % Time between start of pole exit and inaccessiblity
params.spikeSynapticOffset = .02; % Estimated synaptic delay for assigning spikes to contact epochs
params.tid=0; % Trajectory id
params.framesUsed=1:length(array.trials{find(array.whiskerTrialInds,1)}.whiskerTrial.time{1});
params.curveMultiplier=1.5;
params.trialcolors = {'g','r','k','b'};
params.baselineCurve = [0 .02]; % To subtract from the baseline curve for contact detection
params.savedir = 'Z:\users\Andrew\Whisker Project\SingleUnit\'
disp(params)
else
disp('Using custom analysis parameters')
end
for tidx = 1:numel(array.trials{5}.whiskerTrial.time);
% Prefetch all relevant data to prevent numerous slow calls to the trial
% array classes
params.tidx = tidx;
whiskerTIN = find(array.whiskerTrialInds);
goTrialNums = cat(1,array.missTrialNums(:),array.hitTrialNums(:));
trialNums = array.trialNums;
pinDescentOnsetTime(whiskerTIN) = cellfun(@(x)x.behavTrial.pinDescentOnsetTime,array.trials(whiskerTIN),'UniformOutput',0);
pinAscentOnsetTime(whiskerTIN) = cellfun(@(x)x.behavTrial.pinAscentOnsetTime,array.trials(whiskerTIN),'UniformOutput',0);
time(whiskerTIN) = cellfun(@(x)x.whiskerTrial.time{tidx}, array.trials(whiskerTIN),'UniformOutput',0);
distanceToPoleCenter(whiskerTIN) = cellfun(@(x)x.whiskerTrial.distanceToPoleCenter{tidx}, array.trials(whiskerTIN),'UniformOutput',0);
kappa(whiskerTIN) = cellfun(@(x)x.whiskerTrial.deltaKappa{tidx}, array.trials(whiskerTIN),'UniformOutput',0);
thetaAtBase(whiskerTIN) = cellfun(@(x)x.whiskerTrial.thetaAtBase{tidx}, array.trials(whiskerTIN),'UniformOutput',0);
thetaAtContact(whiskerTIN) = cellfun(@(x)x.whiskerTrial.thetaAtContact{tidx}, array.trials(whiskerTIN),'UniformOutput',0);
meanContactCurve = zeros(max(whiskerTIN),1);
trialContactType = zeros(max(whiskerTIN),1);
% Make initial guesses for contact periods
disp('Determining contact periods')
for k = whiskerTIN
if max(goTrialNums==trialNums(k)) %determine trial type
% Finds indexes of time periods of contacts by distance to pole - curvature
touchThreshi=max(params.touchThresh(1:2)); % go trials
contacts{k}.contactInds{tidx} = time{k}(distanceToPoleCenter{k} - params.curveMultiplier*abs(kappa{k}-params.baselineCurve(1)) - ...
(params.curveMultiplier*5*(kappa{k}-params.baselineCurve(1))).^2 < touchThreshi);
else
% Finds indexes of time periods of contacts by distance to pole - curvature
touchThreshi=max(params.touchThresh(3:4)); % nogo trials
contacts{k}.contactInds{tidx} = time{k}(distanceToPoleCenter{k} - params.curveMultiplier*abs(kappa{k}-params.baselineCurve(2)) - ...
(params.curveMultiplier*5*(kappa{k}-params.baselineCurve(2))).^2 < touchThreshi);
end
% Crops indexs to only pole available times
contacts{k}.contactInds{tidx} = contacts{k}.contactInds{tidx}(contacts{k}.contactInds{tidx} >...
params.poleOffset+pinDescentOnsetTime{k} & contacts{k}.contactInds{tidx} < params.poleEndOffset + pinAscentOnsetTime{k} );
[~,~,contacts{k}.contactInds{tidx}]=intersect(contacts{k}.contactInds{tidx}, time{k});
% Calculate the mean curvature during contact period
meanContactCurve(k)=nanmean(kappa{k}(contacts{k}.contactInds{tidx}));
end
for j=1:2
% Determine if contacts are protraction or retraction and refine contact
% periods
disp('Determining contact directions and refining contact periods')
for k = whiskerTIN
if max(goTrialNums==trialNums(k)) %determine trial type
if meanContactCurve(k) < params.goProThresh
touchThreshi=params.touchThresh(1); % use go / protraction threshold
trialContactType(k)=1;
% Finds indexes of time periods of contacts by distance to pole - curvature
contacts{k}.contactInds{tidx} = time{k}(distanceToPoleCenter{k}...
- params.curveMultiplier*abs(kappa{k}-params.baselineCurve(1))...
- (params.curveMultiplier*5*(kappa{k}-params.baselineCurve(1))).^2 <touchThreshi);
elseif meanContactCurve(k) > params.goProThresh
touchThreshi=params.touchThresh(2); % use go / retraction threshold
trialContactType(k)=2;
% Finds indexes of time periods of contacts by distance to pole - curvature
contacts{k}.contactInds{tidx} = time{k}(distanceToPoleCenter{k}...
- params.curveMultiplier*abs(kappa{k}-params.baselineCurve(1))...
- (params.curveMultiplier*5*(kappa{k}-params.baselineCurve(1))).^2 <touchThreshi);
else
trialContactType(k)=0;
end
else
if meanContactCurve(k) < params.nogoProThresh
touchThreshi=params.touchThresh(3); % use no-go / protraction threshold
trialContactType(k)=3;
% Finds indexes of time periods of contacts by distance to pole - curvature
contacts{k}.contactInds{tidx} = time{k}(distanceToPoleCenter{k}...
- params.curveMultiplier*abs(kappa{k}-params.baselineCurve(2))...
- (params.curveMultiplier*5*(kappa{k}-params.baselineCurve(2))).^2 <touchThreshi);
elseif meanContactCurve(k) > params.nogoProThresh
touchThreshi=params.touchThresh(4); % use no-go / retraction threshold
trialContactType(k)=4;
% Finds indexes of time periods of contacts by distance to pole - curvature
contacts{k}.contactInds{tidx} = time{k}(distanceToPoleCenter{k}...
- params.curveMultiplier*abs(kappa{k}-params.baselineCurve(2))...
- (params.curveMultiplier*5*(kappa{k}-params.baselineCurve(2))).^2 <touchThreshi);
else
trialContactType(k)=0;
end
end
% Crops contact indexes to only pole available times
contacts{k}.contactInds{tidx} = ...
contacts{k}.contactInds{tidx}(contacts{k}.contactInds{tidx} > params.poleOffset + pinDescentOnsetTime{k} ...
& contacts{k}.contactInds{tidx} < params.poleEndOffset+ pinAscentOnsetTime{k} );
[~,~,contacts{k}.contactInds{tidx}]=intersect(contacts{k}.contactInds{tidx},time{k}); % c, ia are unused
% Recalculate mean contact curvature with refined contacts
meanContactCurve(k)=nanmean(kappa{k}(contacts{k}.contactInds{tidx}));
contacts{k}.trialContactType=trialContactType(k);
end
end
end