-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathget_footprint.py
188 lines (148 loc) · 5.31 KB
/
get_footprint.py
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
178
179
180
181
182
183
184
185
186
187
188
from NPUltra_analysis.waveform_analysis import get_amps
import numpy as np
import scipy
def find_intercept(p,xlim,ylim,deg):
'''Pretty clunky way of finding the edge intercept coordinate of a straight line from the max amplitude channel, given a desired input angle.
p = tuple of max amp channel (x,y) in pixel space
xlim = x limit of the probe in pixel space (7)
ylim = y limit of the probe in pixel space (47)
deg = defined angle'''
if deg==0:
m = 0
b = p[1] - m*p[0]
x = xlim[1]
y = p[1]
elif 0<deg<90:
m = np.tan(np.deg2rad(deg))
b = p[1] - m*p[0]
dx = abs(p[0] - xlim[1])
dy = abs(p[1] - ylim[1])
if dx<dy:
x = xlim[1]
y = m*x+b
if y>ylim[1]:
y = ylim[1]
x = (y-b)/m
else:
y = ylim[1]
x = (y-b)/m
if x>xlim[1]:
x = xlim[1]
y = m*x+b
elif deg == 90:
m = 0
b = p[1] - m*p[0]
x = p[0]
y = ylim[1]
elif 90<deg<180:
m = np.tan(np.deg2rad(deg))
b = p[1] - m*p[0]
dx = abs(p[0] - xlim[0])
dy = abs(p[1] - ylim[1])
if dx<dy:
x = xlim[0]
y = m*x+b
if y>ylim[1]:
y = ylim[1]
x = (y-b)/m
else:
y = ylim[1]
x = (y-b)/m
if x<xlim[0]:
x = xlim[0]
y = m*x+b
elif deg == 180:
m = 0
b = p[1] - m*p[0]
x = xlim[0]
y = p[1]
elif 180<deg<270:
m = np.tan(np.deg2rad(deg))
b = p[1] - m*p[0]
dx = abs(p[0] - xlim[0])
dy = abs(p[1] - ylim[0])
if dx<dy:
x = xlim[0]
y = m*x+b
if y<ylim[0]:
y = ylim[0]
x = (y-b)/m
else:
y = ylim[0]
x = (y-b)/m
if x<xlim[0]:
x = xlim[0]
y = m*x+b
elif deg == 270:
m = 0
b = p[1] - m*p[0]
x = p[0]
y = ylim[0]
elif 270<deg<360:
m = np.tan(np.deg2rad(deg))
b = p[1] - m*p[0]
dx = abs(p[0] - xlim[1])
dy = abs(p[1] - ylim[0])
if dx<dy:
x = xlim[1]
y = m*x+b
if y<ylim[0]:
y = ylim[0]
x = (y-b)/m
else:
y = ylim[0]
x = (y-b)/m
if x>xlim[1]:
x = xlim[1]
y = m*x+b
elif deg == 360:
m = 0
b = p[1] - m*p[0]
x = xlim[1]
y = p[1]
return x,y
def tolerant_mean(arrs):
'''Allows for averaging across axis=1 in a 2D ragged array. Averaging does not include arrays with len < current idx'''
lens = [len(i) for i in arrs]
arr = np.ma.empty((np.max(lens),len(arrs)))
arr.mask = True
for idx, l in enumerate(arrs):
arr[:len(l),idx] = l
return arr.mean(axis = -1), arr.std(axis=-1)
def get_footprint_radius(unit, threshold=30,scale=6,unit_shape=(384,82),n_rows=48,n_col=8):
''' Given a unit waveform of shape 384 channels x n time samples, return distance metric calculated as the vector length from the maximum amplitude channel about which
the average amplitude is less than or equal to threshold.
unit = matrix of 384 channels x t samples. Is reshaped into an NP Ultra configuration for distance purposes.
threshold = arbitrary voltage threshold used to find the footprint.'''
if unit_shape[0]//n_col==n_rows:
# get the maximum amplitude per channel.
amps = get_amps(unit).reshape(n_rows,n_col)
# subtract noise by averaging the amplitudes from the 10 lowest amplitude channels and subtracting from all channels.
amps = amps-np.mean(np.sort(amps.reshape(n_rows*n_col))[:10])
# work-around for oddly shaped mean unit waveform arrays
else:
if unit_shape[0]%n_col==0:
amps = get_amps(unit).reshape(unit_shape[0]//n_col,n_col)
amps = amps-np.mean(np.sort(amps.reshape((unit_shape[0]//n_col)*n_col))[:10])
elif unit_shape[0]%n_col==1:
unit = unit[1:,:]
amps = get_amps(unit).reshape(unit.shape[0]//n_col,n_col)
amps = amps-np.mean(np.sort(amps.reshape((unit.shape[0]//n_col)*n_col))[:10])
# get channel index for the max. amp. channel in both unwrapped and grid coords.
max_chan = np.where(amps == np.max(amps))
row,col = max_chan[1][0],max_chan[0][0]
vectors = []
for deg in np.arange(0,360):
x,y = find_intercept((row,col),(0,(amps.shape[1])-1),(0,(amps.shape[0])-1),deg)
hyp = int(round(np.sqrt((x-row)**2 + (y-col)**2)*scale))
c,d = np.linspace(row, x, int(hyp)), np.linspace(col, y, int(hyp))
zi = scipy.ndimage.map_coordinates(amps, np.vstack((d,c)))
vectors.append(zi)
mean_vector = tolerant_mean(vectors)[0]
try:
footprint = np.where(mean_vector<=threshold)[0][0]
except:
footprint=scale
if footprint<scale:
footprint=scale
return footprint