OSCR

Protocol for simultaneous in vivo two-photon imaging and locomotion quantification during olfactory stimulation and pharmacology in walking Drosophila.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches
  1. [1] § Step-by-step method details › Data analysis ↔ cia_utils.py, lines 743–848 · score 0.67 · ball radius, fictrac dat, frame rate, mm
  2. [2] § Step-by-step method details › Data analysis ↔ main.py, lines 124–184 · score 0.53 · ball radius, frame rate, dat, fictrac

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,617 lines · 61 KB · no license · 1 match

  1. # -*- coding: utf-8 -*-
  2. #-*- coding : utf-8-*-
  3. # coding:unicode_escape
  5. import os
  6. import cv2
  7. import json
  8. import pandas as pd
  9. import numpy as np
  10. import matplotlib.pyplot as plt
  11. import matplotlib.patches as patches
  12. from matplotlib.colors import NoNorm
  13. from glob import glob
  14. from scipy.stats import circmean
  16. from cia_consts import *
  17. from cia_photodiode import get_stim_range, get_all_sti_pos, find_nearest
  18. import imageio
  19. from matplotlib.animation import FuncAnimation
  20. COLORS = ["k", "r", "g", "b", "y", "c", "m", "gray", "pink", "springgreen", "deepskyblue", "yellow",]
  23. def load_memmap(filename, mode='r'):
  24. file_to_load = filename
  25. filename = os.path.split(filename)[-1]
  26. fpart = filename.split('_')[1:-1] # The filename encodes the structure of the map
  27. d1, d2, d3, T, order = int(fpart[-9]), int(fpart[-7]), int(fpart[-5]), int(fpart[-1]), fpart[-3]
  28. Yr = np.memmap(file_to_load, mode=mode, shape=(d1 * d2 * d3, T), dtype=np.float32, order=order)
  29. if d3 == 1:
  30. dims = (d1, d2)
  31. else:
  32. dims = (d1, d2, d3)
  33. return np.reshape(Yr.T, [T] + list(dims), order='F')
  35. def calc_avg_frame(m, parent):
  36. avg_frame = np.mean(m, axis=0)
  37. eq = norm_img(avg_frame)
  38. cv2.imwrite(parent + "/i_avg.png", eq)
  39. std_frame = np.std(m, axis=0)
  40. eq2 = norm_img(std_frame)
  41. cv2.imwrite(parent + "/i_std.png", eq2)
  42. max_frame = np.max(m, axis=0)
  43. eq3 = norm_img(max_frame)
  44. cv2.imwrite(parent + "/i_max.png", eq3)
  45. eqc = np.transpose([eq3, eq2, eq], (1, 2, 0))
  46. cv2.imwrite(parent + "/i_all.png", eqc)
  48. # cluster_roi(m, parent)
  49. # show_cluster_roi(m, parent)
  50. # th2, res = cv2.threshold(eq3.astype(np.uint8), 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
  51. # cv2.imwrite(parent + "/i_test.png", res)
  53. def cluster_roi(m, parent):
  54. eq3 = norm_img(np.max(m, axis=0))*2 # NOTE: m: 1400*128*128
  55. th2, img = cv2.threshold(eq3.astype(np.uint8), 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
  56. #img = cv2.morphologyEx(img, cv2.MORPH_OPEN, cv2.getStructuringElement(cv2.MORPH_RECT, (1, 1)))
  57. p_l = np.nonzero(img)
  58. n = len(p_l[0])
  59. cluster_m = np.zeros(img.shape, dtype=np.uint8)
  60. # for i in range(n):
  61. # print(i)
  62. # for j in range(n):
  63. # cor_map[i, j] = cor(m[:, p_l[0][i], p_l[1][i]], m[:, p_l[0][j], p_l[1][j]])
  64. X = m[:, p_l[0], p_l[1]].T # NOTE: n*1400
  66. from pyclustering.cluster.kmeans import kmeans
  67. from pyclustering.cluster.center_initializer import kmeans_plusplus_initializer
  68. from pyclustering.utils.metric import type_metric, distance_metric
  69. metric = distance_metric(type_metric.USER_DEFINED, func=lambda a, b: 1-cor(a, b))
  70. initial_centers = kmeans_plusplus_initializer(X, 10).initialize()
  71. kmeans_instance = kmeans(X, initial_centers, metric=metric)
  72. kmeans_instance.process()
  73. clusters = kmeans_instance.get_clusters()
  74. pred = np.zeros((X.shape[0]))
  75. cf = []
  76. for i, c in enumerate(clusters):
  77. pred[c] = i + 1
  78. cf.extend(c)
  80. plot_fast_cor_map(X, cf)
  81. # plt.show()
  82. plt.savefig(parent + "/i_cor.png")
  84. # from sklearn.cluster import KMeans, DBSCAN
  85. # estimator = KMeans(n_clusters=20, max_iter=500)
  86. # estimator = DBSCAN(eps=0.8, min_samples=50, metric=lambda a, b: 1-cor(a, b))
  88. # estimator.fit(X)
  89. # pred = estimator.labels_
  91. cluster_m[p_l[0], p_l[1]] = pred
  92. cv2.imwrite(parent + "/i_test.png", norm_img(cluster_m))
  94. for i in range(int(np.max(pred))):
  95. show_largest_for_roi(m, cluster_m, i)
  96. plt.savefig(parent + "/i_largest_%d.png" % i)
  97. # plt.imshow(cluster_m)
  99. def plot_cor_map(X, cf):
  100. n = len(cf)
  101. cor_map = np.zeros((n, n))
  102. for i, cf1 in enumerate(cf):
  103. print(i, "/", n)
  104. for j, cf2 in enumerate(cf):
  105. cor_map[i][j] = cor(X[cf1], X[cf2])
  106. plt.imshow(cor_map, cmap="bwr")
  107. plt.colorbar()
  109. def plot_fast_cor_map(X, cf):
  110. from scipy.stats import zscore
  111. Xz = zscore(X, axis=1)
  112. n = len(cf)
  113. cor_map = np.zeros((n, n))
  114. for i, cf1 in enumerate(cf):
  115. print(i, "/", n)
  116. for j, cf2 in enumerate(cf):
  117. cor_map[i][j] = np.dot(Xz[cf1], Xz[cf2])/n
  118. plt.imshow(cor_map, cmap="bwr")
  119. plt.colorbar()
  121. def show_cluster_roi(m, parent):
  122. cluster_m = cv2.imread(parent + "/i_test.png")
  123. for k, i in enumerate(np.unique(cluster_m)):
  124. show_largest_for_roi(m, cluster_m, i)
  125. plt.savefig(parent + "/i_largest_%d.png" % k)
  127. def show_largest_for_roi(m, cluster_m, label):
  128. p_l = np.nonzero(cluster_m == label)
  129. X = m[:, p_l[0], p_l[1]]
  130. n = X.shape[1]
  131. # if n < 1600:
  132. # cor_map = np.zeros((n, n))
  133. # for i in range(n):
  134. # print(i, "/", n)
  135. # for j in range(n):
  136. # cor_map[i][j] = cor(X[:, i], X[:, j])
  137. # plt.imshow(cor_map)
  138. # plt.colorbar()
  139. # plt.show()
  140. f = np.argmax(np.sum(X, axis=1))
  141. plt.figure(figsize=(6, 3))
  142. plt.title(str(f))
  143. fig, axs = plt.subplots(1, 2)
  144. axs[0].imshow(m[f])
  145. axs[0].scatter(p_l[1], p_l[0], alpha=1, color="r", s=1)
  146. axs[1].imshow(m[f])
  148. def norm_img(img):
  149. pmax = np.percentile(img, 99)
  150. pmin = np.percentile(img, 1)
  151. r_norm = 255.0 * (img - pmin) / (pmax - pmin)
  152. np.clip(r_norm, 0, 255, r_norm)
  153. return r_norm.astype(np.uint8)
  155. def load_roi(roi_file, shape):
  156. if roi_file.endswith(".zip"):
  157. return load_roi_zip(roi_file, shape)
  158. elif roi_file.endswith(".npy"):
  159. roi = np.load(roi_file, allow_pickle=True)
  160. names = np.arange(len(roi))
  161. xy = [r[:, 0, :] for r in roi]
  162. return names, roi_contours_to_points(xy, shape), xy
  164. def load_roi_zip(roizip, shape):
  165. from read_roi import read_roi_zip# read_roi_file
  166. rois = read_roi_zip(roizip)
  168. xy = []
  169. names = []
  170. for d in rois.values():
  171. if d["type"] == "polygon":
  172. xy.append(np.array(tuple(zip(d["x"], d["y"]))))
  173. names.append(d["name"])
  174. elif d["type"] == "rectangle":
  175. l, t, w, h = d["left"], d["top"], d["width"], d["height"]
  176. xy.append(np.array([[l, t], [l + w, t], [l + w, t + h], [l, t + h]]))
  177. names.append(d["name"])
  179. return names, roi_contours_to_points(xy, shape), xy
  181. def roi_contours_to_points(contours_xy, shape):
  182. xy = []
  183. for c in contours_xy:
  184. temp = np.zeros(shape)
  185. cv2.drawContours(temp, [c.astype(int)], 0, color=1, thickness=-1)
  186. xy.append(temp.nonzero())
  187. # plt.imshow(temp, cmap="Greys_r")
  188. # plt.show()
  189. return xy
  191. def calc_all_roi_F(roizip, m, parent):
  192. print("Do calc F...")
  193. names, points, contours = load_roi(roizip, m[0].shape)
  194. ret = [[] for i in range(len(names))]
  195. for i, mi in enumerate(m):
  196. for j, xyj in enumerate(points):
  197. ret[j].append(float(mi[xyj].mean()))
  198. # plot_lines(ret, names)
  199. plot_rois(contours, m)
  200. plt.savefig(parent + "/roi.png")
  202. F = pd.DataFrame(np.array(ret).T, columns=names)
  203. F.to_csv(parent + "/F.csv", index=False)
  204. plot_lines([F[n] for n in names], names)
  205. plt.savefig(parent + "/F.png")
  207. zs_csv = parent + "/zscore.csv"
  208. zs = get_zscore(F)
  209. zs.to_csv(zs_csv, index=False)
  210. write_roi_dFF_video(zs, contours, m, zs_csv)
  212. dFF_csv = parent + "/dFF.csv"
  213. dFF = get_dFF(F)
  214. dFF.to_csv(dFF_csv, index=False)
  215. write_roi_dFF_video(dFF, contours, m, dFF_csv)
  216. print(dFF_csv)
  217. return dFF_csv, zs_csv
  219. def write_roi_dFF_video(dFF, contours, m, name):
  220. return
  221. # fig = plt.figure(figsize=(6, 6))
  222. # ax = plt.gca()
  223. # ax.set_position([0, 0, 1, 1], which="both")
  224. # fig.canvas.draw()
  225. # w, h = fig.canvas.get_width_height()
  226. if os.path.exists(name + "_roi.avi"):
  227. return
  228. h, w = m[0].shape
  229. output_video = cv2.VideoWriter(name + "_roi.avi", cv2.VideoWriter_fourcc(*"DIVX"), 30, (w, h))
  230. cmap = plt.cm.get_cmap("Greys_r")
  231. dFFn = norm_img(np.array(dFF))
  232. for i, mi in enumerate(m):
  233. print(i)
  234. # ax.cla()
  235. # ax.imshow(mi, cmap="Greys_r")
  236. # for j, xys in enumerate(contours):
  237. # ax.add_patch(plt.Polygon(xys, linewidth=0, color=cmap(dFFn.loc[i, j])))
  238. # write_fig_to_video(output_video, fig, w, h)
  239. img = np.zeros(mi.shape, dtype=np.uint8)#np.array(mi, dtype=np.uint8)
  240. for j, xys in enumerate(contours):
  241. cv2.drawContours(img, [contours[j].astype(int)], 0, color=cmap(dFFn[i, j])[0]*255, thickness=-1)
  242. img_bgr = cv2.cvtColor(img, cv2.COLOR_GRAY2BGR)
  243. output_video.write(img_bgr)
  244. output_video.release()
  246. def write_zscore_video(m, name, pv_idx=None, contours=None):
  247. # n, h, w = m.shape
  248. cmap = plt.cm.get_cmap("plasma")
  249. mg = np.array([cv2.GaussianBlur(mi, (9, 9), 0) for mi in m])
  250. me = np.mean(mg, axis=0)
  251. m1 = mg - me
  252. m1[m1 < 0] = 0
  253. m1 = norm_img(m1)
  255. m2 = []
  256. for i, mi in enumerate(m1):
  257. img = mi.astype(np.uint8)#cv2.GaussianBlur(mi.astype(np.uint8), (3, 3), 0)
  258. if pv_idx is not None and i < len(pv_idx) and pv_idx[i] >= 0:
  259. # cv2.putText(img, str(pva[i]), (6, 40), cv2.FONT_HERSHEY_PLAIN, 1, (255, 255, 128), 1)
  260. cv2.drawContours(img, [contours[pv_idx[i]].astype(int)], 0, color=(255, 255, 255), thickness=1)
  261. m2.append(img)
  262. write_video(name, (cmap(m2)[:,:,:,:3] * 255).astype(np.uint8), cv2.COLOR_RGB2BGR, FICTRAC_RATE)
  264. # m_mean = np.mean(m2, axis=0)
  265. # m_std = np.std(m2, axis=0)
  266. # mz = []
  267. # for mi in m2:
  268. # zs = (mi - m_mean) / m_std
  269. # zs = np.nan_to_num(zs)
  270. # zs[zs<0]=0
  271. # mz.append(zs)
  272. # mz = norm_img(mz)
  273. # write_video(name + "_zscore.avi", (cmap(mz)[:,:,:,:3] * 255).astype(np.uint8), cv2.COLOR_RGB2BGR, FICTRAC_RATE)
  275. # def write_fig_to_video(output_video, fig, w, h, desc=None, save_img_path=None):
  276. # fig.canvas.draw()
  277. # img = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8)
  278. # img.shape = (h*2, w*2, 3)
  279. # img_bgr = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)
  280. # desc and cv2.putText(img_bgr, desc, (6, 20), cv2.FONT_HERSHEY_PLAIN, 1, (255, 255, 128), 1)
  281. # output_video.write(img_bgr)
  282. # save_img_path and cv2.imwrite(save_img_path, img_bgr)
  283. # return img_bgr
  285. def plot_lines(lines, names, ylim=None, xlim_r=1):
  286. n = len(lines)
  287. fig, axes = plt.subplots(n, 1, sharex=True, figsize=(30, 12), dpi=300)
  288. plt.subplots_adjust(left=0.05, right=0.99, top=0.95, bottom=0.05, hspace=0)
  289. for i, r in enumerate(lines):
  290. ax = axes[n - i - 1]
  291. ax.plot(r, c=COLORS[i%12])
  292. ax.set_ylabel(names[i], rotation=0, fontsize=6)
  293. if ylim is not None:
  294. ax.set_ylim(ylim)
  295. if xlim_r is not None:
  296. ax.set_xlim((0, xlim_r * len(r)))
  298. def plot_lines_with_err(py_l, name, USE_SEM=False):
  299. py = np.nanmean(py_l, axis=0)
  300. pe = np.nanstd(py_l, axis=0)
  301. if USE_SEM:
  302. pe = pe / np.sqrt(np.count_nonzero(~np.isnan(py_l), axis=0))
  303. plt.figure()
  304. plt.plot(np.transpose(py_l), alpha=0.1)
  305. plt.plot(py, "--")
  306. plt.fill_between(range(len(py_l[0])), py - pe, py + pe, alpha=0.1)
  307. plt.savefig(name)
  308. plt.close()
  310. def plot_lines_in_one(lines, names):
  311. plt.figure(figsize=(15, 10), dpi=300)
  312. for i, n in enumerate(lines):
  313. plt.plot(n, linewidth=1, c=COLORS[i%12], alpha=0.5)
  315. def plot_pva(ax, dff, unwrap=False, c="k", offset=0):
  316. # max_idx = np.argmax(dff.T, axis=0) + offset
  317. pv_dir, pv_len = calc_pva(dff)
  318. if unwrap:
  319. pv_dir = unwrap_dir(pv_dir)
  320. pi = dff.shape[1]/2
  321. plot_angle(ax, (np.array(pv_dir)+np.pi)*pi/np.pi+offset, c, pi)
  323. def plot_hot(df, save_name, img_rate, is_PB=False):
  324. plt.figure(figsize=(24, 2), dpi=300)
  325. plt.subplots_adjust(left=0.02, right=0.99)
  326. dff = df.to_numpy()
  327. plt.pcolor(dff.T, cmap="jet")
  328. if is_PB:
  329. plot_pva(plt.gca(), dff[:, :9])
  330. plot_pva(plt.gca(), dff[:, 9:], offset=9)
  331. else:
  332. plot_pva(plt.gca(), dff)
  334. frames = len(dff)
  335. seconds = int(frames / img_rate)
  336. labels = np.linspace(0, seconds, 9)
  337. plt.xticks(labels * img_rate, labels)
  338. # plt.colorbar()
  339. plt.savefig(save_name + ".png")
  341. plt.figure(figsize=(24, 2), dpi=300)
  342. plt.subplots_adjust(left=0.02, right=0.99)
  343. if is_PB:
  344. plot_pva(plt.gca(), dff[:, :9], unwrap=True)
  345. plot_pva(plt.gca(), dff[:, 9:], unwrap=True, c="gray", offset=9)
  346. else:
  347. plot_pva(plt.gca(), dff, unwrap=True)
  348. plt.xticks(labels * img_rate, labels)
  349. plt.xlim(0, seconds * img_rate)
  350. plt.savefig(save_name + "_unwrap.png")
  352. def lim_dir(dir1, pi=np.pi):
  353. if dir1 > pi:
  354. dir1 -= 2*pi
  355. elif dir1 < -pi:
  356. dir1 += 2*pi
  357. return dir1
  359. def lim_dir_l(dir_l, pi=np.pi):
  360. dir1 = dir_l.copy()
  361. dir1[dir1 > pi] -= 2*pi
  362. dir1[dir1 < -pi] += 2*pi
  363. return dir1
  365. def unwrap_dir(v, pi=np.pi):
  366. ret = []
  367. li = v[0]
  368. offset = 0
  369. for i in v:
  370. i += offset
  371. d = i - li
  372. if d > pi:
  373. offset -= 2*pi
  374. i -= 2*pi
  375. elif d < -pi:
  376. offset += 2*pi
  377. i += 2*pi
  378. ret.append(i)
  379. li = i
  380. return ret
  382. def unwrap_dir_win(v, rate, sec=10):
  383. win = int(sec*rate)
  384. ret = []
  385. for i in range(0, len(v), win):
  386. u = unwrap_dir(v[i:i + win])
  387. ret.extend([d-u[0] for d in u])
  388. return ret
  390. def plot_angle(ax, al, c, pi=np.pi, xs=None):
  391. # al = unwrap_dir(al) #[lim_dir(a-2) for a in al]
  392. if xs is None:
  393. xs = range(len(al))
  394. last = 0
  395. start = 0
  396. for i, a in enumerate(al):
  397. if abs(a - last) > pi:
  398. ax.plot(xs[start:i], al[start:i], c=c, lw=0.6)
  399. start = i
  400. last = a
  401. ax.plot(xs[start:len(al)], al[start:], c=c, lw=0.6)
  403. def plot_scatter(ax, al, c, xs=None):
  404. if xs is None:
  405. xs = range(len(al))
  406. ax.scatter(xs, al, c=c, s=0.1)
  408. def rotate_unit_vec(rad):
  409. return np.array([np.cos(rad), np.sin(rad)])
  411. def calc_pva(f):
  412. # pv_len = np.max(f, axis=1)
  413. # pv_dir = (np.argmax(f, axis=1) + 0.5) / 8 * np.pi - np.pi
  415. vs = np.array([rotate_unit_vec(2*np.pi*r/f.shape[1]-np.pi) for r in range(f.shape[1])])
  416. pv = f.dot(vs)
  417. pv_len = np.sqrt(np.sum(pv ** 2, axis=1)) / f.shape[1]
  418. pv_dir = np.arctan2(pv[:, 1], pv[:, 0])
  419. return pv_dir, pv_len
  421. def bump_amplitude(f):
  422. bump_a = np.max(f, axis=1) - np.min(f, axis=1)
  423. nor_bump_a = []
  424. max = np.mean(bump_a[np.argpartition(bump_a,-50)[-50:]])
  425. min = np.min(bump_a)
  426. for x in bump_a:
  427. x = float(x - min) / (max - min)
  428. nor_bump_a.append(x)
  429. return nor_bump_a
  431. def bin_bump_speed(f,ft_speed,ft_vx,ft_vy,ft_vz): #for now, ft_speed-mm/s, ft_vx`ft_vy`ft_vz-rad/s
  432. nor_bump_a = bump_amplitude(f)
  433. all_speed = {'0-':[], '0~1':[] , '1~2':[],'2~3':[],'3~4':[],'4~5':[],'5~6':[],'6~7':[],'7~8':[],'8+':[]}
  434. forward_speed = {'0-':[], '0~1':[] , '1~2':[],'2~3':[],'3~4':[],'4~5':[],'5~6':[],'6~7':[],'7~8':[],'8~9':[],'9~10':[],'10+':[]}
  435. lateral_speed = {'-4-':[], '-4~-3':[] , '-3~-2':[],'-2~-1':[],'-1~0':[],'0~1':[],'1~2':[],'2~3':[],'3~4':[],'4+':[]}
  436. # rotation_speed = {}
  437. for i in range(len(ft_speed) - 10):
  438. mean_ft_speed = np.mean(ft_speed[i:i+10])
  439. mean_ft_vx = np.mean(ft_vx[i:i+10] * 8) # forward speed
  440. mean_ft_vy = np.mean(ft_vy[i:i+10] * 8) # lateral speed
  442. if mean_ft_speed < 0:
  443. all_speed['0-'].append(np.mean(nor_bump_a[i:i+10]))
  444. elif mean_ft_speed > 0 and mean_ft_speed <= 1:
  445. all_speed['0~1'].append(np.mean(nor_bump_a[i:i+10]))
  446. elif mean_ft_speed > 1 and mean_ft_speed <= 2:
  447. all_speed['1~2'].append(np.mean(nor_bump_a[i:i+10]))
  448. elif mean_ft_speed > 2 and mean_ft_speed <= 3:
  449. all_speed['2~3'].append(np.mean(nor_bump_a[i:i+10]))
  450. elif mean_ft_speed > 3 and mean_ft_speed <= 4:
  451. all_speed['3~4'].append(np.mean(nor_bump_a[i:i+10]))
  452. elif mean_ft_speed > 4 and mean_ft_speed <= 5:
  453. all_speed['4~5'].append(np.mean(nor_bump_a[i:i+10]))
  454. elif mean_ft_speed > 5 and mean_ft_speed <= 6:
  455. all_speed['5~6'].append(np.mean(nor_bump_a[i:i+10]))
  456. elif mean_ft_speed > 6 and mean_ft_speed <= 7:
  457. all_speed['6~7'].append(np.mean(nor_bump_a[i:i+10]))
  458. elif mean_ft_speed > 8:
  459. all_speed['8+'].append(np.mean(nor_bump_a[i:i+10]))
  461. if mean_ft_vx < 0:
  462. forward_speed['0-'].append(np.mean(nor_bump_a[i:i+10]))
  463. elif mean_ft_vx > 0 and mean_ft_vx <= 1:
  464. forward_speed['0~1'].append(np.mean(nor_bump_a[i:i + 10]))
  465. elif mean_ft_vx > 1 and mean_ft_vx <= 2:
  466. forward_speed['1~2'].append(np.mean(nor_bump_a[i:i + 10]))
  467. elif mean_ft_vx > 2 and mean_ft_vx <= 3:
  468. forward_speed['2~3'].append(np.mean(nor_bump_a[i:i + 10]))
  469. elif mean_ft_vx > 3 and mean_ft_vx <= 4:
  470. forward_speed['3~4'].append(np.mean(nor_bump_a[i:i + 10]))
  471. elif mean_ft_vx > 4 and mean_ft_vx <= 5:
  472. forward_speed['4~5'].append(np.mean(nor_bump_a[i:i + 10]))
  473. elif mean_ft_vx > 5 and mean_ft_vx <= 6:
  474. forward_speed['5~6'].append(np.mean(nor_bump_a[i:i + 10]))
  475. elif mean_ft_vx > 6 and mean_ft_vx <= 7:
  476. forward_speed['6~7'].append(np.mean(nor_bump_a[i:i + 10]))
  477. elif mean_ft_vx > 7 and mean_ft_vx <= 8:
  478. forward_speed['7~8'].append(np.mean(nor_bump_a[i:i + 10]))
  479. elif mean_ft_vx > 8 and mean_ft_vx <= 9:
  480. forward_speed['8~9'].append(np.mean(nor_bump_a[i:i + 10]))
  481. elif mean_ft_vx > 10:
  482. forward_speed['10+'].append(np.mean(nor_bump_a[i:i + 10]))
  484. if mean_ft_vy < -4:
  485. lateral_speed['-4-'].append(np.mean(nor_bump_a[i:i+10]))
  486. elif mean_ft_vy > -4 and mean_ft_vy < -3:
  487. lateral_speed['-4~-3'].append(np.mean(nor_bump_a[i:i+10]))
  488. elif mean_ft_vy > -3 and mean_ft_vy < -2:
  489. lateral_speed['-3~-2'].append(np.mean(nor_bump_a[i:i+10]))
  490. elif mean_ft_vy > -2 and mean_ft_vy < -1:
  491. lateral_speed['-2~-1'].append(np.mean(nor_bump_a[i:i+10]))
  492. elif mean_ft_vy > -1 and mean_ft_vy < 0:
  493. lateral_speed['-1~0'].append(np.mean(nor_bump_a[i:i+10]))
  494. elif mean_ft_vy > 0 and mean_ft_vy < 1:
  495. lateral_speed['0~1'].append(np.mean(nor_bump_a[i:i+10]))
  496. elif mean_ft_vy > 1 and mean_ft_vy < 2:
  497. lateral_speed['1~2'].append(np.mean(nor_bump_a[i:i+10]))
  498. elif mean_ft_vy > 2 and mean_ft_vy < 3:
  499. lateral_speed['2~3'].append(np.mean(nor_bump_a[i:i+10]))
  500. elif mean_ft_vy > 3 and mean_ft_vy < 4:
  501. lateral_speed['3~4'].append(np.mean(nor_bump_a[i:i+10]))
  502. elif mean_ft_vy > 4:
  503. lateral_speed['4+'].append(np.mean(nor_bump_a[i:i+10]))
  505. return all_speed, forward_speed, lateral_speed
  507. def unify_sample(f, ts, fps, n): # no interp
  508. idx = []
  509. inter = 1.0 / fps
  510. i = 0
  511. for j in range(n):
  512. t = j * inter
  513. if ts[i+1] < t:
  514. while i < len(ts) - 1 and ts[i+1] < t:
  515. i += 1
  516. if i+1 >= len(f):
  517. break
  518. if t < ts[i]:
  519. idx.append(i)
  520. continue
  521. if i >= len(ts) - 1:
  522. break
  523. if t - ts[i] > ts[i+1] - t:
  524. idx.append(i+1)
  525. else:
  526. idx.append(i)
  527. return f[idx], idx
  529. def calc_offset(ft, pv):
  530. n = min(len(ft), len(pv))
  531. o = ft[:n] - pv[:n]
  532. # o[o < 0] += 2*np.pi
  533. o[o > np.pi] -= 2*np.pi
  534. o[o < -np.pi] += 2*np.pi
  535. return o
  537. def plot_slide_cor(ax, ft, pv, rate):
  538. for win_sec in [10]:#4, 16, 64
  539. win = int(win_sec * rate)
  540. ax.plot(*slide_cor(ft, pv, win, rate))
  542. def circular_std(o):
  543. o = o[~np.isnan(o)]
  544. # s = np.sin(o)
  545. # c = np.cos(o)
  546. # return np.sqrt(-2 * np.log(np.sqrt(np.sum(s) ** 2 + np.sum(c) ** 2) / np.linalg.norm(o)))
  547. from scipy.stats import circstd
  548. return circstd(o)
  550. def format_time(t):
  551. if t < 60:
  552. return "%.2f" % t
  553. return "%d:%.2f" % (t//60, t%60)
  555. def non_nan(s):
  556. return s[~np.isnan(s)]
  558. def plot_dFF_MB_grating(dFF_name, exp_info, unify_rate=None, use_fictrac = False):
  559. parent = os.path.dirname(dFF_name)
  560. if use_fictrac:
  561. fictrac_name = get_ft_dat(parent)
  562. ft = pd.read_csv(fictrac_name, header=None).to_numpy() # (15079frame, 23field) 50Hz
  563. bar_name = get_ft_bar(parent)
  564. pd_name = get_pd_h5(parent)
  566. dFF = pd.read_csv(dFF_name).to_numpy()[:, 1:] # (2200frame, 16roi) 6.736Hz
  567. import h5py
  568. ni = h5py.File(pd_name, "r")
  569. pd_info_raw = ni['AI']["photodiode"][:, 0] # (1640499frame,) 5000Hz
  570. frame_counter = ni['CI']["FrameCounter"][:, 0]
  571. frame_out = ni['DI']["FrameOut"][:, 0]
  572. img_rate = exp_info["frameRate"] # only one z-step in this experiment
  573. if not unify_rate:
  574. unify_rate = img_rate
  575. use_img_rate = True
  576. seconds = len(pd_info_raw) / PD_RATE # int(len(dFF) / img_rate)
  577. n = int(seconds * unify_rate)
  578. seconds = n / unify_rate
  580. frame_start = np.nonzero(np.diff(frame_counter) > 0)[0]
  581. volume_frame = 1
  582. frame_time = frame_start[int(exp_info["steps"] / 2)::volume_frame] / PD_RATE
  584. # smooth filter for PD data
  585. n = 100
  586. pd_info_raw = np.convolve(pd_info_raw, np.ones((n,)) / n, mode='same')
  587. temp = get_all_sti_pos(pd_info_raw > 0.02)
  588. sti_position_duration = []
  589. for cc in temp:
  590. if cc[1] > PD_RATE:
  591. sti_position_duration.append(cc)
  593. sti_ima_frame = []
  594. for cc in sti_position_duration:
  595. sti_ima_frame.append([find_nearest(frame_start,cc[0]),find_nearest(frame_start,cc[0]+cc[1])])
  597. fig, axs = plt.subplots(np.shape(dFF)[1]+1, 1, figsize=(20, 10), dpi=300, sharex=True)
  598. axs[0].set_title(parent.split('\\')[-2] + ' ' + parent.split('\\')[-1])
  599. ts = np.linspace(1,len(dFF),len(dFF)) / img_rate
  600. for i in range(np.shape(dFF)[1]):
  601. axs[i].plot(ts, dFF[:,i])
  602. axs[i].set_ylim( np.min(dFF[:,i]), np.max(dFF[:,i]) )
  603. for kk in sti_ima_frame:
  604. axs[i].add_patch(
  605. patches.Rectangle(
  606. (kk[0] / img_rate, np.min(dFF[:,i])),
  607. (kk[1]-kk[0]) / img_rate,
  608. np.max(dFF[:,i]),
  609. edgecolor=None,
  610. facecolor='blue',
  611. alpha = 0.3
  612. ))
  613. axs[i].set_ylabel('ROI' + str(i + 1) + ' _deltaF/F')
  614. axs[i].set_xlabel('t/s')
  615. plt.savefig(os.path.join(parent,'all_ROI.png'))
  616. plt.savefig(os.path.join(parent, 'all_ROI.pdf'))
  618. each_sti = []
  619. for kk in sti_ima_frame:
  620. each_sti.append((dFF[kk[0] - 100 : kk[1] + 100,:]).tolist())
  621. pass
  622. json.dump(each_sti, open(parent + "/each_sti.txt", "w"))
  627. ###########################################################################
  628. # m = load_memmap(glob(os.path.join(parent, "*.mmap"))[0])
  629. #
  630. # m2 = []
  631. # cmap = plt.cm.get_cmap("gray")
  632. # for i, mi in enumerate(m):
  633. # img = mi.astype(np.uint8)#cv2.GaussianBlur(mi.astype(np.uint8), (3, 3), 0)
  634. # m2.append(img)
  635. # write_video(parent + "/mc_sti.avi",(cmap(m2)[:,:,:,:3] * 255).astype(np.uint8),need_time=False,sti_img_frame=sti_ima_frame)
  640. # try:
  641. # whole_sti_start_end = (dFF[sti_ima_frame[0][0] - 150 : sti_ima_frame[-1][-1] + 150]).tolist()
  642. # whole_sti_pos = sti_ima_frame - sti_ima_frame[0][0] + 150
  643. # except:
  644. # whole_sti_start_end = []
  645. # whole_sti_pos = np.nan
  646. # json.dump( [whole_sti_start_end, whole_sti_pos.tolist()]
  647. # , open(parent + "/whole_sti_start_end.txt", "w"))
  649. pass
  651. def plot_dFF_MB_UVorLaser(dFF_name, exp_info, unify_rate=None, use_fictrac = False):
  652. parent = os.path.dirname(dFF_name)
  653. if use_fictrac:
  654. fictrac_name = get_fictrac_dat(parent)
  655. ft = pd.read_csv(fictrac_name, header=None).to_numpy() # (15079frame, 23field) 50Hz
  656. bar_name = get_ft_bar(parent)
  657. pd_name = get_pd_h5(parent)
  659. dFF = pd.read_csv(dFF_name).to_numpy()[:, 1:] # (2200frame, 16roi) 6.736Hz
  660. import h5py
  661. ni = h5py.File(pd_name, "r")
  662. pd_info_raw = ni['AI']["UVLED"][:, 0] # the DAQ input of the matlab (UV or laser signal & fictrac trigger signal)
  663. frame_counter = ni['CI']["FrameCounter"][:, 0]
  664. frame_out = ni['DI']["FrameOut"][:, 0]
  665. img_rate = exp_info["frameRate"] # only one z-step in this experiment
  666. if not unify_rate:
  667. unify_rate = img_rate
  668. use_img_rate = True
  669. seconds = len(pd_info_raw) / PD_RATE # int(len(dFF) / img_rate)
  670. n = int(seconds * unify_rate)
  671. seconds = n / unify_rate
  673. frame_start = np.nonzero(np.diff(frame_counter) > 0)[0]
  675. if exp_info.get("zFastEnable"):
  676. img_rate = exp_info["frameRate"] / (exp_info["steps"] + exp_info["flybackFrames"])
  677. volume_frame = exp_info["steps"] + exp_info["flybackFrames"]
  678. frame_time = frame_start[int(exp_info["steps"] / 2)::volume_frame] / PD_RATE
  679. frame_start=frame_start[int(exp_info["steps"] / 2)::volume_frame]
  680. else:
  681. img_rate = exp_info["frameRate"]
  682. volume_frame = 1
  683. frame_time = frame_start / PD_RATE
  684. # volume_frame = 1
  685. # frame_time = frame_start[int(exp_info["steps"] / 2)::volume_frame] / PD_RATE
  687. #a = pd_info_raw > 0.04
  688. temp = get_all_sti_pos(pd_info_raw > 0.8)
  689. sti_position_duration = []
  690. for cc in temp:
  691. if cc[1] > PD_RATE/2:
  692. sti_position_duration.append(cc)
  694. sti_ima_frame = []
  695. for cc in sti_position_duration:
  696. sti_ima_frame.append([find_nearest(frame_start,cc[0]),find_nearest(frame_start,cc[0]+cc[1])])
  698. fig, axs = plt.subplots(np.shape(dFF)[1]+1, 1, figsize=(20, 10), dpi=300, sharex=True)
  699. axs[0].set_title(parent.split('\\')[-2] + ' ' + parent.split('\\')[-1])
  700. ts = np.linspace(1,len(dFF),len(dFF)) / img_rate
  701. for i in range(np.shape(dFF)[1]):
  702. axs[i].plot(ts, dFF[:,i])
  703. axs[i].set_ylim( np.min(dFF[:,i]), np.max(dFF[:,i]) )
  704. for kk in sti_ima_frame:
  705. axs[i].add_patch(
  706. patches.Rectangle(
  707. (kk[0] / img_rate, np.min(dFF[:,i])),
  708. (kk[1]-kk[0]) / img_rate,
  709. np.max(dFF[:,i]),
  710. edgecolor=None,
  711. facecolor='blue',
  712. alpha = 0.3
  713. ))
  714. axs[i].set_ylabel('ROI' + str(i + 1))
  715. axs[i].set_xlabel('t/s')
  716. plt.savefig(os.path.join(parent,'all_ROI.png'))
  717. # plt.savefig(os.path.join(parent, 'all_ROI.pdf'))
  719. each_sti = []
  720. for kk in sti_ima_frame:
  721. each_sti.append((dFF[kk[0] - int(img_rate * 2) : kk[1] + int(img_rate * 2),:]).tolist())
  722. pass
  723. json.dump(each_sti, open(parent + "/each_sti.txt", "w"))
  724. try:
  725. whole_sti_start_end = (dFF[sti_ima_frame[0][0] - 150 : sti_ima_frame[-1][-1] + 150]).tolist()
  726. whole_sti_pos = sti_ima_frame - sti_ima_frame[0][0] + 150
  727. except:
  728. whole_sti_start_end = whole_sti_pos = []
  729. json.dump( [whole_sti_start_end, whole_sti_pos.tolist()]
  730. , open(parent + "/whole_sti_start_end.txt", "w"))
  732. # ########################################################################### write video ######################################
  733. # m = load_memmap(glob(os.path.join(parent, "*.mmap"))[0])
  734. #
  735. # m2 = []
  736. # cmap = plt.cm.get_cmap("gray")
  737. # for i, mi in enumerate(m):
  738. # img = mi.astype(np.uint8)#cv2.GaussianBlur(mi.astype(np.uint8), (3, 3), 0)
  739. # m2.append(img)
  740. # write_video(parent + "/mc_sti.avi",(cmap(m2)[:,:,:,:3] * 255).astype(np.uint8),need_time=False,sti_img_frame=sti_ima_frame)
  743. def plot_dFF_MB_fictrac(dFF_name, exp_info, unify_rate=None, use_fictrac = True):
  744. parent = os.path.dirname(dFF_name)
  745. if use_fictrac:
  746. fictrac_name = get_fictrac_dat(parent)
  747. ft = pd.read_csv(fictrac_name, header=None).to_numpy() # (15079frame, 23field) 50Hz
  748. #bar_name = get_ft_bar(parent)
  749. pd_name = get_pd_h5(parent)
  750. dFF = pd.read_csv(dFF_name).to_numpy()[:, 1:] # (2200frame, 16roi) 6.736Hz
  751. import h5py
  752. ni = h5py.File(pd_name, "r")
  753. pd_info_raw = ni['AI']["side_camera"][:, 0] # the DAQ input of the matlab (UV or laser signal & fictrac trigger signal)
  754. frame_counter = ni['CI']["FrameCounter"][:, 0]
  755. frame_out = ni['DI']["FrameOut"][:, 0]
  757. frame_start = np.nonzero(np.diff(frame_counter) > 0)[0]
  758. if exp_info.get("zFastEnable"):
  759. img_rate = exp_info["frameRate"] / (exp_info["steps"] + exp_info["flybackFrames"])
  760. volume_frame = exp_info["steps"] + exp_info["flybackFrames"]
  761. frame_time = frame_start[int(exp_info["steps"] / 2)::volume_frame] / PD_RATE
  762. frame_start=frame_start[int(exp_info["steps"] / 2)::volume_frame]
  763. else:
  764. img_rate = exp_info["frameRate"]
  765. volume_frame = 1
  766. frame_time = frame_start / PD_RATE
  768. unify_rate = img_rate
  771. # "int(exp_info["steps"]/2)" get the middle frame of each volumn
  773. #fictrac_start = np.nonzero(np.diff(pd_info_raw) > 1)[0]
  775. temp = get_all_sti_pos(pd_info_raw > 1)
  776. sti_position_duration = []
  777. for cc in temp:
  778. if cc[1] > 100:
  779. sti_position_duration.append(cc)
  781. fictrac_start = []
  782. for cc in sti_position_duration:
  783. fictrac_start.append(find_nearest(frame_start,cc[0]))
  785. fictrac_time = np.array(sti_position_duration)[:,0] / PD_RATE
  787. sti_ima_frame = [find_nearest(frame_start, sti_position_duration[0][0]), find_nearest(frame_start,sti_position_duration[-1][0])]
  789. dFF_unifyToFictrac = dFF[fictrac_start][0:-1]
  790. #Fictrac_data = ft[1:-2][:] # discard the first and last frame
  792. if len(dFF_unifyToFictrac) > len(ft):
  793. dFF_unifyToFictrac = dFF_unifyToFictrac[:len(ft)]
  794. else:
  795. ft = ft[:len(dFF_unifyToFictrac)]
  796. Fictrac_data = ft
  797. ft_speed = Fictrac_data[:, 18] * FICTRAC_RATE * BALL_RADIUS # mm/s all_speed
  798. ft_vy = -Fictrac_data[:, 5] * FICTRAC_RATE * BALL_RADIUS # rad/s lateral_speed
  799. ft_vx = Fictrac_data[:, 6] * FICTRAC_RATE * BALL_RADIUS # rad/s forward_speed
  800. ft_vz = -Fictrac_data[:, 7] * FICTRAC_RATE # rad/s rotation_speed
  803. for roi_index in range(np.shape(dFF)[1]):
  804. window = np.ones(int(10)) / float(10)
  805. dFF_unifyToFictrac_avg = np.convolve(dFF_unifyToFictrac[:,roi_index], window, 'same')
  806. ft_speed_avg = np.convolve(ft_speed, window, 'same')
  807. ft_vx_avg = np.convolve(ft_vx, window, 'same')
  808. ft_vy_avg = np.convolve(ft_vy, window, 'same')
  809. ft_vz_avg = np.convolve(ft_vz, window, 'same')
  810. ts = np.linspace(1, len(ft_speed_avg), len(ft_speed_avg)) / FICTRAC_RATE
  812. fig, axs = plt.subplots(4, 1, figsize=(20, 10), dpi=300, sharex=True)
  814. axs[0].set_title(os.path.basename(os.path.dirname(parent)) + os.path.basename(parent))
  815. axs[0].plot(ts, dFF_unifyToFictrac_avg, 'black')
  816. axs[0].set_ylabel('$\Delta$F/F')
  817. twin0 = axs[0].twinx()
  818. twin0.plot(ts, ft_speed_avg, 'red')
  819. twin0.set_ylabel('all speed mm/s', color='r', fontsize=8)
  821. axs[1].plot(ts, dFF_unifyToFictrac_avg, 'black')
  822. axs[1].set_ylabel('$\Delta$F/F')
  823. twin0 = axs[1].twinx()
  824. twin0.plot(ts, ft_vx_avg, 'g')
  825. twin0.set_ylabel('forward speed mm/s', color='g', fontsize=8)
  827. axs[2].plot(ts, dFF_unifyToFictrac_avg, 'black')
  828. axs[2].set_ylabel('$\Delta$F/F')
  829. twin0 = axs[2].twinx()
  830. twin0.plot(ts, ft_vy_avg, 'b')
  831. twin0.set_ylabel('lateral speed mm/s', color='b', fontsize=8)
  833. axs[3].plot(ts, dFF_unifyToFictrac_avg, 'black')
  834. axs[3].set_ylabel('$\Delta$F/F')
  835. twin0 = axs[3].twinx()
  836. twin0.plot(ts, ft_vz_avg, 'orange')
  837. twin0.set_ylabel('rotation speed rad/s', color='orange', fontsize=8)
  838. axs[3].set_xlabel('t/s')
  839. plt.savefig(os.path.join(os.path.dirname(dFF_name), dFF_name +str(roi_index+1)+ "_fictrac.png"))
  841. json.dump([dFF_unifyToFictrac.tolist(), ft_speed.tolist(), ft_vx.tolist(), ft_vy.tolist(), ft_vz.tolist()], open(parent + "/dFF_speed.txt", "w"))
  843. # seconds = len(pd_info_raw) / PD_RATE # int(len(dFF) / img_rate)
  844. # n = int(seconds * unify_rate)
  845. # #dFF_info, dff_unify_idx = unify_sample(dFF, frame_time, unify_rate, n)
  846. # ft_info, ft_unify_idx = unify_sample(ft, fictrac_time - fictrac_time[0], unify_rate, n)
  848. pass
  851. def plot_dFF_fictrac(dFF_name, exp_info, unify_rate=None):
  852. parent = os.path.dirname(dFF_name)
  853. fictrac_name = get_ft_dat(parent)
  854. bar_name = get_ft_bar(parent)
  855. pd_name = get_pd_h5(parent)
  856. ft = pd.read_csv(fictrac_name, header=None).to_numpy() # (15079frame, 23field) 50Hz
  857. dFF = pd.read_csv(dFF_name).to_numpy()[:, 1:] # (2200frame, 16roi) 6.736Hz
  858. import h5py
  859. ni = h5py.File(pd_name, "r")
  860. pd_info_raw = ni['AI']["photodiode"][:, 0] # (1640499frame,) 5000Hz
  861. frame_counter = ni['CI']["FrameCounter"][:, 0]
  862. frame_out = ni['DI']["FrameOut"][:, 0]
  864. use_img_rate = False
  865. img_rate = exp_info["frameRate"]/(exp_info["steps"] + exp_info["flybackFrames"])
  866. if not unify_rate:
  867. unify_rate = img_rate
  868. use_img_rate = True
  870. seconds = len(pd_info_raw) / PD_RATE #int(len(dFF) / img_rate)
  871. n = int(seconds * unify_rate)
  872. seconds = n / unify_rate
  873. frame_start = np.nonzero(np.diff(frame_counter) > 0)[0]
  874. volume_frame = exp_info["steps"] + exp_info["flybackFrames"]
  875. frame_time = frame_start[int(exp_info["steps"]/2)::volume_frame] / PD_RATE
  876. # frame_time = np.arange(frame_start[0]/PD_RATE, seconds, 1.0/img_rate)
  878. pd_info = down_sample(pd_info_raw[:int(seconds * PD_RATE)], n)
  879. pd_info_50 = down_sample(pd_info_raw[:int(seconds * PD_RATE)], int(seconds * 50))
  880. # frame_info = down_sample(frame_out[:int(seconds * PD_RATE)], n)/10.0
  881. dFF_info, dff_unify_idx = unify_sample(dFF, frame_time, unify_rate, n)
  882. ts = ft[:, -2]/1000
  883. if max(ts - ts[0]) <= 0:
  884. ts = np.arange(0, len(ts)) / FICTRAC_RATE
  885. ft_info, ft_unify_idx = unify_sample(ft, ts - ts[0], unify_rate, n)
  886. ft_bar_info = load_ft_bar(bar_name)
  887. if ft_bar_info is not None:
  888. ft_bar, ft_unify_idx = unify_sample(ft_bar_info, ts - ts[0], unify_rate, n)
  890. pd_thresh = 0.14
  891. ft_range = get_stim_range(pd_info_50 > pd_thresh)
  892. ft_range_t = np.array(ft_range[0]) / 50
  893. ft_range_u = (ft_range_t * unify_rate).astype(int)
  894. dff_range_t = [frame_time[0], frame_time[-1]]
  895. ft_unify_idx = np.concatenate([[0]*ft_range_u[0], ft_unify_idx]).astype(int)
  896. json.dump([ft_range_t.tolist(), dff_range_t, ft_unify_idx.tolist(), dff_unify_idx, ft_bar_info.tolist() if ft_bar_info is not None else []], open(dFF_name + "_time_info.txt", "w"))
  897. # exit(0)
  899. # ft_xs = np.arange(len(ft_info)) + ft_range_u[0]
  900. ft_heading = np.concatenate([[np.nan]*ft_range_u[0], ft_info[:, 16]])
  901. # ft_heading = lim_dir_l(ft_heading)
  902. if ft_bar_info is not None:
  903. ft_bar = np.concatenate([[np.nan]*ft_range_u[0], 2*np.pi - ft_bar])
  904. # ft_bar = lim_dir_l(ft_bar)
  905. else:
  906. ft_bar = ft_heading
  907. ft_speed = np.concatenate([[np.nan]*ft_range_u[0], ft_info[:, 18]]) * FICTRAC_RATE * BALL_RADIUS # mm/s
  908. ft_vy = np.concatenate([[np.nan]*ft_range_u[0], -ft_info[:, 5]]) * FICTRAC_RATE # rad/s
  909. ft_vx = np.concatenate([[np.nan]*ft_range_u[0], ft_info[:, 6]]) * FICTRAC_RATE
  910. ft_vz = np.concatenate([[np.nan]*ft_range_u[0], -ft_info[:, 7]]) * FICTRAC_RATE
  911. #########################################################################################################################################################
  912. all_speed, forward_speed, lateral_speed = bin_bump_speed(dFF_info,ft_speed,ft_vx,ft_vy,ft_vz)
  913. json.dump([all_speed, forward_speed, lateral_speed], open(parent + "/bin_bump_speed.txt", "w"))
  914. #return
  915. ###############################################################################################################################################################
  916. roi_n = dFF_info.shape[1]
  917. if use_img_rate:
  918. fig, axs = plt.subplots(roi_n, 1, figsize=(10, 30), dpi=300, sharex=True)
  919. pv_dir, pv_len = calc_pva(dFF_info)
  920. cor_s, cor_e = int(ft_range_u[0]+5*unify_rate), min(ft_range_u[1], len(ft_bar), len(pv_dir))
  922. ft_bar_cor = ft_bar[cor_s:cor_e]
  923. for i in range(roi_n):
  924. dffi = dFF_info[cor_s:cor_e, i]
  925. axs[i].scatter(ft_bar_cor, dffi, s=1)
  926. plt.savefig(parent + "/ft_pv_tune")
  927. #return
  929. fig, axs = plt.subplots(10, 1, figsize=(20, 10), dpi=300, sharex=True)
  930. plt.subplots_adjust(left=0.02, right=0.99)
  932. axs[0].pcolor(dFF_info.T, cmap="jet")
  933. is_pb = dFF_name.find("-PB-") > 0
  934. if is_pb:
  935. pv_dir, pv_len = calc_pva(dFF_info[:, :9])
  936. # pv_dir = lim_dir_l(np.array(pv_dir))
  937. pv_dir2, pv_len2 = calc_pva(dFF_info[:, 9:])
  938. # pv_dir2 = lim_dir_l(np.array(pv_dir2))
  939. # max_idx = np.argmax(dFF_info.T[:9], axis=0)
  940. # pv_dir = (max_idx + 0.5)*np.pi/4.5
  941. plot_angle(axs[0], (pv_dir+np.pi) * 4.5 / np.pi + 9, "k", 4.5)
  942. else:
  943. pv_dir, pv_len = calc_pva(dFF_info)
  944. # pv_dir = lim_dir_l(np.array(pv_dir))
  945. pv_dir += np.pi
  946. pv_dir2, pv_len2 = pv_dir, pv_len
  947. # max_idx = np.argmax(dFF_info.T, axis=0)
  948. # pv_dir = (max_idx + 0.5)*np.pi/8
  949. pi = roi_n / 2
  950. # plot_angle(axs[0], (np.array(pv_dir) + np.pi) * pi / np.pi, "k", pi)
  951. plot_angle(axs[0], pv_dir * pi / np.pi, "k", pi)
  953. # axs[1].plot(frame_info, alpha=0.2, c="k")
  954. axs[1].plot(pd_info, c="b")
  956. plot_angle(axs[2], pv_dir, c="k")
  957. # if is_pb:
  958. # ft_bar = -ft_bar
  959. # plot_angle(axs[2], pv_dir2, c="gray")
  960. plot_angle(axs[2], ft_bar, "b") # bar
  961. # plot_angle(axs[2], ft_heading, "g") # heading
  962. axs[2].set_ylim(0, 2*np.pi)
  963. # axs[2].set_ylim(-np.pi, np.pi)
  964. # axs[2].add_patch(plt.Rectangle((0, -np.pi), n, np.pi / 4, alpha=0.1, color="k", linewidth=0))
  965. # axs[2].add_patch(plt.Rectangle((0, 3*np.pi/4), n, np.pi / 4, alpha=0.1, color="k", linewidth=0))
  966. axs[2].add_patch(plt.Rectangle((0, 0), n, np.pi / 4, alpha=0.1, color="k", linewidth=0))
  967. axs[2].add_patch(plt.Rectangle((0, 7*np.pi/4), n, np.pi / 4, alpha=0.1, color="k", linewidth=0))
  969. offset = calc_offset(ft_bar, pv_dir)
  970. axs[3].plot(offset)
  972. # plot_angle(axs[4], lim_dir_l(pv_dir), c="k")
  973. ft_bar_correct = lim_dir_l(ft_bar - circmean(non_nan(offset)))
  974. ft_bar_correct[ft_bar_correct < 0] += 2 * np.pi
  975. plot_angle(axs[0], ft_bar_correct / (2*np.pi) * roi_n, "r") # heading - offset
  976. # axs[4].set_ylim(-np.pi, np.pi)
  978. plot_scatter(axs[5], unwrap_dir_win(ft_bar, unify_rate), "b")
  979. plot_scatter(axs[5], unwrap_dir_win(pv_dir, unify_rate), "k")
  980. axs[5].set_ylim(-6, 6)
  982. axs[6].plot(unwrap_dir(ft_bar), c="b", lw=0.5)
  983. axs[6].plot(unwrap_dir(ft_heading), c="g", lw=0.5)
  984. axs[6].plot(unwrap_dir(smooth_angle(pv_dir, 10)), "k--", lw=0.5)
  985. axs[6].plot(unwrap_dir(pv_dir), c="k", lw=0.5)
  987. plot_slide_cor(axs[-3], ft_bar, pv_dir, unify_rate)
  988. axs[-3].set_ylim(-1, 1)
  990. axs[-2].plot(ft_speed, c="r") # speed
  991. axs[-1].plot(zscore1(ft_speed), c="r", alpha=0.5)
  992. axs[-1].plot(zscore1(pv_len), c="k", alpha=0.5)
  994. # plt.show()
  995. ticks = np.arange(0, len(pd_info), 5*unify_rate)
  996. axs[-1].set_xticks(ticks)
  997. axs[-1].set_xticklabels(["%.1f"%l for l in ticks/unify_rate])
  999. cor_s, cor_e = int(ft_range_u[0]+5*unify_rate), min(ft_range_u[1], len(ft_bar), len(pv_dir))
  1000. if cor_s > cor_e:
  1001. print("error: ft range not found!")
  1002. exit(0)
  1003. ft, pv = ft_bar[cor_s:cor_e], pv_dir[cor_s:cor_e]
  1004. cor_unwrap = cor(unwrap_dir(ft), unwrap_dir(pv))
  1005. cor_cir = cir_cor(ft, pv)
  1006. axs[0].set_title("%s fps: %.2f, speed: %.2f, dFF (%s, %s), fictrac (%s, %s), dFF~heading: %.2f (unwrap) %.2f (circ)" % (
  1007. os.path.basename(parent), unify_rate, np.nanmean(ft_speed),
  1008. format_time(dff_range_t[0]), format_time(dff_range_t[1]), format_time(ft_range_t[0]), format_time(ft_range_t[1]),
  1009. cor_unwrap, cor_cir
  1010. ))
  1011. plt.savefig(os.path.join(os.path.dirname(dFF_name), dFF_name + "_fictrac.png"))
  1012. json.dump({"cor_unwrap": cor_unwrap, "cor_cir": cor_cir}, open(parent + "/info.txt", "w"))
  1013. ft_pv = [
  1014. [cor_s/unify_rate, cor_e/unify_rate], ft.tolist(), pv.tolist(), pv_len[cor_s:cor_e].tolist(), ft_speed[cor_s:cor_e].tolist(),
  1015. ft_vx[cor_s:cor_e].tolist(), ft_vy[cor_s:cor_e].tolist(), ft_vz[cor_s:cor_e].tolist()
  1016. ]
  1017. if pv_dir2 is not None:
  1018. ft_pv.append(pv_dir2[cor_s:cor_e].tolist())
  1019. ft_pv.append(pv_len2[cor_s:cor_e].tolist())
  1020. json.dump(ft_pv, open(parent + "/ft_pv.txt", "w"))
  1021. json.dump(offset.tolist(),open(parent + "/offset.txt", "w"))
  1024. plt.figure()
  1025. plt.hist(offset, bins=50)
  1026. plt.title("%s fps: %.2f, mean: %.2f SD: %.2f cSD: %.2f" % (
  1027. os.path.basename(parent), unify_rate,
  1028. np.nanmean(offset), np.nanstd(offset), circular_std(offset)
  1029. ))
  1030. plt.savefig(os.path.join(os.path.dirname(dFF_name), dFF_name + "_offset.png"))
  1032. # ft_win = window * FICTRAC_RATE
  1033. # dFF_win = window * IMG_RATE
  1034. # s = int(min(len(ft) / ft_win, len(dFF) / dFF_win))
  1035. # cor_dir, cor_speed = [], []
  1036. # for i in range(s):
  1037. # ft_s = ft[int(i * ft_win): int((i+1) * ft_win)]
  1038. # dFF_s = dFF[int(i * dFF_win): int((i+1) * dFF_win)]
  1039. # pv_dir, pv_len = calc_pva(dFF_s)
  1040. # # plt.figure("pv_dir");plt.plot(unwrap_dir(pv_dir))
  1041. # # plt.figure("heading");plt.plot(unwrap_dir(ft_s[:, 16]))
  1042. # # plt.figure("pv_len");plt.plot(pv_len)
  1043. # # plt.figure("speed");plt.plot(ft_s[:, 18])
  1044. # # plt.show()
  1045. # cor_dir.append(cor_equal_len(unwrap_dir(pv_dir), unwrap_dir(ft_s[:, 16])))
  1046. # cor_speed.append(cor_equal_len(pv_len, ft_s[:, 18]))
  1047. # plt.plot(cor_dir)
  1048. # plt.plot(cor_speed)
  1049. # plt.show()
  1051. def plot_rois(rois, m):
  1052. plt.figure(figsize=(3, 3), dpi=300)
  1053. plt.axis("off")
  1054. ax = plt.gca()
  1055. ax.set_position([0, 0, 1, 1], which="both")
  1056. ax.imshow(norm_img(np.std(m, axis=0)), cmap=plt.cm.gray)
  1057. for i, xys in enumerate(rois):
  1058. ax.add_patch(plt.Polygon(xys, alpha=0.6, fill=False, linewidth=1, color=COLORS[i%12]))
  1059. x, y = xys.min(axis=0)/2 + xys.max(axis=0)/2
  1060. ax.text(x, y, str(i), color="r")
  1062. def get_part(fname):
  1063. if fname.find("-PB-") >= 0:
  1064. return "PB"
  1065. elif fname.find("-FB-") >= 0:
  1066. return "FB"
  1067. elif fname.find("MB") >= 0:
  1068. return "manual draw"
  1069. elif fname.find("57C10") >= 0:
  1070. return "manual draw"
  1071. elif fname.find("D161") >= 0:
  1072. return "manual draw"
  1073. elif fname.find("BU") >= 0:
  1074. return "manual draw"
  1075. elif fname.find("AOTU") >= 0:
  1076. return "manual draw"
  1077. elif fname.find("EB") >= 0:
  1078. # return "manual draw"
  1079. return "EB"
  1080. else:
  1081. return "manual draw"
  1083. def plot_fictrac(fname, parent, win=10, step=1, speed_min=0.5):
  1084. df = pd.read_csv(fname, header=None).to_numpy().T
  1085. heading = df[16]
  1086. speed = df[18]
  1087. inv = speed*FICTRAC_RATE < speed_min
  1088. heading[inv] = np.nan
  1089. stand = (np.count_nonzero(inv) / len(heading))
  1090. print("stand still %.2f" % stand)
  1091. vec_x, vec_y = np.cos(df[16]), np.sin(df[16])
  1092. mx, my = [], []
  1093. for i in range(0, len(heading) - win * FICTRAC_RATE, step * FICTRAC_RATE):
  1094. mx.append(np.nanmean(vec_x[i: i+win * FICTRAC_RATE]))
  1095. my.append(np.nanmean(vec_y[i: i+win * FICTRAC_RATE]))
  1096. mean_vec = np.nanmean(np.sqrt(np.array(mx)**2+np.array(my)**2))
  1097. plt.figure()
  1098. plt.title("walk %.2f, mean %.2f" % (1-stand, mean_vec))
  1099. plt.axis("equal")
  1100. plt.hist2d(mx, my, bins=20, range=((-1, 1), (-1, 1)))
  1101. plt.savefig(parent + "_mean_vec.png")
  1102. plt.figure()
  1103. plt.plot(speed*FICTRAC_RATE)
  1104. plt.savefig(fname + "_speed.png")
  1105. return
  1107. def get_pd_h5(parent):
  1108. pp = os.path.dirname(parent)
  1109. #if (not parent[-1].isdigit()) or parent[-1] == 'D' or parent[-1] == 'g' or parent[-1] == 'r' or parent[-1] == 'V' or parent[-6:] == 'MBON11':n
  1110. if (not parent[-3:].isdigit()):
  1111. ppn = os.path.basename(parent).rstrip("+") + '_000'
  1112. else:
  1113. ppn = os.path.basename(parent).rstrip("-").rstrip("+").replace("_0", "_0")
  1114. return os.path.join(pp, "PD_" + ppn + "/Episode001.h5")
  1116. def get_fictrac_dat(parent):
  1117. file_folder = os.path.basename(parent)
  1118. temp = file_folder.split('_')
  1119. if len(temp) == 5:
  1120. data_folder = temp[2][0:2] + temp[3][0] + '-0'
  1121. elif len(temp) == 6:
  1122. data_folder = temp[2][0:2] + temp[3][0] + '-' + str(int(temp[-1]))
  1124. if temp[0] == 'D161':
  1125. data_folder = data_folder.replace('L','')
  1126. if temp[0] == 'VC102':
  1127. # if len(temp) == 5:
  1128. # data_folder = temp[2][0:2] + temp[3] + '-0'
  1129. # elif len(temp) == 6:
  1130. # data_folder = temp[2][0:2] + temp[3] + '-' + str(int(temp[-1]))
  1131. if 'R' in data_folder:
  1132. data_folder = data_folder.replace('R', '')
  1133. elif 'L' in data_folder:
  1134. data_folder = data_folder.replace('L', '')
  1135. return glob(os.path.join(os.path.dirname(parent), data_folder)+ "/*.dat")[0]
  1136. pass
  1138. def get_ft_dat(parent):
  1139. pp = os.path.dirname(parent)
  1140. ppn = os.path.basename(parent).rstrip("-").rstrip("+")
  1141. #ppn = os.path.basename(parent)
  1142. tt = ppn.split("-")
  1143. if len(tt) > 2:
  1144. return glob(os.path.join(pp, tt[2] + "-" + tt[-1]) + "/*.dat")[0]
  1145. else:
  1146. t3 = tt[0].split("_F")
  1147. t4 = t3[1].split("_")
  1148. t5 = t4[-1].split("-")
  1149. if len(tt) == 2:
  1150. n = tt[-1]
  1151. else:
  1152. n = t5[-1]
  1153. if n == 'CLOSED':
  1154. n = '000'
  1155. return glob(os.path.join(pp, "F" + t4[0] + "-" + str(int(n))) + "/*.dat")[0]
  1158. def get_ft_bar(parent):
  1159. return os.path.join(os.path.dirname(get_ft_dat(parent)), "bar_position.mat")
  1161. def write_fictrac_stim_dFF_video(parent):
  1162. ft_name = get_ft_dat(parent)
  1163. ft_folder = os.path.dirname(ft_name)
  1164. ft_range_t, dff_range, ft_idx, dff_idx, ft_bar_info = json.load(open(parent + "/zscore.csv_time_info.txt", "r"))
  1165. m = load_memmap(glob(os.path.join(parent, "*.mmap"))[0])
  1166. m_color = m[dff_idx]
  1168. c = json.load(open(parent + "/ft_pv.txt", "r")) #50Hz, fictrac start +5s
  1169. t_range, heading, pv, pvl, speed, vx, vy, vx, pv2, pvl2 = c
  1170. pva = np.concatenate([[np.nan] * int((ft_range_t[0]+5)*50), pv]) # after proc_csv
  1171. names, points, contours = load_roi(parent + "/roi.npy", m[0].shape)
  1172. roi_n = len(contours) - 1 # 16 or 18
  1173. pv_contour_idx = ((pva*roi_n/(2*np.pi)) - 0.5).astype(int) + 1 # (0~16)->(0-15)+1
  1174. write_zscore_video(m_color, parent + "/v_color.avi", pv_contour_idx, contours)
  1175. # return
  1179. if not ft_bar_info:
  1180. df = pd.read_csv(ft_name, header=None).to_numpy().T
  1181. heading = 2*np.pi - df[16]
  1182. else:
  1183. heading = ft_bar_info
  1184. m_stim = []
  1185. w, h = 320, 160
  1186. barh = 5
  1187. for d in heading:
  1188. img = np.zeros((h, w), dtype=np.uint8)
  1189. mid = int(d * w / (2*np.pi))
  1190. img[:, mid - barh:mid + barh] = 255
  1191. if mid < barh:
  1192. img[:, 0 :mid + barh] = 255
  1193. img[:, w-(barh-mid):w] = 255
  1194. if mid + barh >= w:
  1195. img[:, mid:w] = 255
  1196. img[:,0:mid+barh-w] = 255
  1197. m_stim.append(img)
  1198. write_video(parent + "/v_stim.avi", np.array(m_stim)[ft_idx], fps=FICTRAC_RATE)
  1199. #embed_video(ft_folder + "/fictrac-debug.avi", parent + "/v_stim.avi", parent + "/v_color.avi", ft_range_t[0], ft_range_t[0])
  1200. embed_video(os.path.join(ft_folder,"fictrac-debug.avi"), os.path.join(parent,'v_stim.avi') , os.path.join(parent,'v_color.avi'), ft_range_t[0],
  1201. ft_range_t[0])
  1203. def realtime_heading_PVA(frames, heading, pv,parent):
  1204. filenames = []
  1205. for i in range(0,len(heading),200):
  1207. plt.figure(figsize=(30, 3))
  1208. if i > 100:
  1209. plt.plot(heading[i-100:i])
  1210. plt.plot(pv[i-100:i])
  1211. else:
  1212. plt.plot(heading[0:i])
  1213. plt.plot(pv[0:i])
  1214. filename = f'{i}.png'
  1215. filenames.append(filename)
  1216. plt.savefig(filename)
  1217. plt.close()
  1219. with imageio.get_writer('mygif.gif', mode='I') as writer:
  1220. for filename in filenames:
  1221. image = imageio.imread(filename)
  1222. writer.append_data(image)
  1224. fourcc = cv2.VideoWriter_fourcc('D', 'I', 'V', 'X')
  1225. video = cv2.VideoWriter(os.path.join(parent,'heading_pva.avi'), fourcc, 10, (3000,300))
  1226. for item in filenames:
  1227. img = cv2.imread(item)
  1228. video.write(img)
  1229. video.release()
  1231. for filename in set(filenames):
  1232. os.remove(filename)
  1233. pass
  1236. def plot_dFF(dff_csv, fr):
  1237. dFF = pd.read_csv(dff_csv)
  1238. names = dFF.columns
  1239. print(names)
  1241. plot_lines([dFF[n] for n in names], names)
  1242. plt.savefig(dff_csv + ".png")
  1243. plot_hot(dFF[dFF.keys()[1:]], dff_csv + "_hot", fr, is_PB=(dff_csv.find("-PB-") > 0))
  1244. # plot_lines([auto_cor(dFF[n][200:]) for n in names], names, ylim=(-0.2, 0.5), xlim_r=0.505)
  1245. # plt.savefig(parent + "/dFF_cor.png")
  1247. def get_dFF_by_baseline_range(F, baseline_range, bg_idx=0):
  1248. names = F.columns
  1249. bg = F[names[bg_idx]].copy()
  1250. for g in names:
  1251. F[g] = F[g] - bg
  1252. baseline = np.mean(F[g][baseline_range[0]:baseline_range[1]])
  1253. F[g] = (F[g] - baseline) / baseline
  1254. return F
  1256. def get_dFF(F, bg_correction = True):
  1257. names = F.columns
  1258. bg = F[names[0]]
  1259. for g in names[1:]:
  1260. if np.mean(bg) < np.mean(F[g]) + 5:
  1261. if bg_correction:
  1262. f = F[g] - bg
  1263. else:
  1264. f = F[g] - np.mean(bg)
  1265. else:
  1266. f = F[g]
  1267. f = f.mask(f < 0, np.nan)
  1268. baseline = np.nanmean(f[f <= np.nanpercentile(f, 5)]) # NOTE: F as the mean of the lower 5% (MaimonG_Nat17)
  1269. if baseline <= 1:
  1270. baseline = 1
  1271. print(g, "baseline", baseline)
  1272. F[g] = (f - baseline) / baseline
  1273. return F
  1275. def get_zscore(F):
  1276. from scipy.stats import zscore
  1277. zs = zscore(F, axis=0)
  1278. return pd.DataFrame(zs, columns=F.columns)
  1280. def zscore1(xs):
  1281. m, s = np.nanmean(xs), np.nanstd(xs)
  1282. return (xs - m) / s
  1284. def cov1(xs, ys):
  1285. n = len(xs)
  1286. mx, my = np.mean(xs), np.mean(ys)
  1287. return (xs - mx).dot(ys - my) / (n - 1)
  1289. def cor1(xs, ys):
  1290. n = len(xs)
  1291. return zscore1(xs).dot(zscore1(ys)) / (n - 1)
  1293. def cor2(xs, ys):
  1294. return cov1(xs, ys) / np.std(xs) / np.std(ys)
  1296. def cir_cor(xs, ys): # circstat circ_corrcc
  1297. mx, my = circmean(xs), circmean(ys)
  1298. num = np.sin(xs - mx).dot(np.sin(ys - my))
  1299. den = np.sqrt(np.sum(np.sin(xs - mx) ** 2) * np.sum(np.sin(ys - my) ** 2))
  1300. return num / den
  1302. def auto_cor(xs):
  1303. xm = xs - xs.mean()
  1304. xn = np.sum(xm ** 2)
  1305. return np.correlate(xm, xm, "same") / xn
  1307. def cor(x, y):
  1308. # return np.dot(x, y)/np.count_nonzero(y>0.5)
  1309. return np.corrcoef(x, y)[0][1]
  1311. def cor_equal_len(x, y):
  1312. if len(x) > len(y):
  1313. return cor(down_sample(x, len(y)), y)
  1314. else:
  1315. return cor(x, down_sample(y, len(x)))
  1317. def slide_cor(x, y, win, rate):
  1318. n = min(len(x), len(y))
  1319. rx, ry = [], []
  1320. for i in range(0, n - win, int(rate)):
  1321. c = cor(unwrap_dir(x[i:i + win]), unwrap_dir(y[i:i + win]))
  1322. rx.append(i + win/2.0)
  1323. ry.append(c)
  1324. return rx, ry
  1326. def down_sample(x, n):
  1327. x = np.array(x, dtype=np.float)
  1328. index_arr = np.linspace(0, len(x) - 1, num=n, dtype=np.float)
  1329. index_floor = np.array(index_arr, dtype=np.int)
  1330. index_ceil = index_floor + 1
  1331. index_rem = index_arr - index_floor
  1333. val1 = x[index_floor]
  1334. val2 = x[index_ceil % len(x)]
  1335. interp = val1 * (1.0 - index_rem) + val2 * index_rem
  1336. assert (len(interp) == n)
  1337. return interp
  1339. def view_img_seq(imgs, shifts_rig=None):
  1340. global g_total_frame, g_frame, g_is_input_begin, g_input_int
  1341. g_total_frame, h, w = imgs.shape
  1342. g_frame = 0
  1343. g_is_input_begin = False
  1344. g_input_int = 0
  1346. def plot_one_frame(f):
  1347. global g_frame
  1348. g_ax.cla()
  1349. g_ax.imshow(imgs[f].astype(int), cmap=plt.cm.gray, norm=NoNorm())
  1350. # g_ax.set_xlabel("%02d:%02.2f" % (t_sec / 60, t_sec % 60))
  1351. g_ax.set_xlabel(str(f))
  1352. if shifts_rig is not None:
  1353. g_ax.set_title(str(shifts_rig[g_frame]))
  1354. g_ax.grid(True)
  1355. g_ax.set_xticks(np.linspace(0, imgs[f].shape[1], 13))
  1356. g_frame = f
  1357. def on_slider(val):
  1358. plot_one_frame(int(val))
  1359. def onkey(event):
  1360. print(event.key)
  1361. global g_frame, g_is_input_begin, g_input_int
  1362. if event.key == "left":
  1363. g_frame -= 1
  1364. elif event.key == "right":
  1365. g_frame += 1
  1366. elif event.key == "enter":
  1367. g_frame = g_input_int
  1368. g_input_int = 0
  1369. elif event.key in list([*"1234567890"]):
  1370. g_input_int = int(event.key) + g_input_int * 10
  1371. print("input: %d" % g_input_int)
  1372. return
  1373. else:
  1374. g_input_int = 0
  1375. return
  1376. if g_frame >= g_total_frame:
  1377. g_frame = g_total_frame - 1
  1378. if g_frame < 0:
  1379. g_frame = 0
  1380. g_slider.set_val(g_frame)
  1381. event.canvas.draw()
  1383. from matplotlib.widgets import Slider
  1384. fig, g_ax = plt.subplots(figsize=(w/25, h/25))
  1385. plt.subplots_adjust(top=0.95, bottom=0.1)
  1386. g_slider = Slider(plt.axes([0.1, 0.03, 0.8, 0.03]), "", valmin=0, valmax=g_total_frame - 1, valfmt="%d", valinit=0)
  1387. g_slider.on_changed(on_slider)
  1388. plot_one_frame(0)
  1389. fig.canvas.mpl_connect('key_press_event', onkey)
  1391. def load_exp_xml(xml):
  1392. if not os.path.exists(xml):
  1393. return None
  1394. def get_value(ss, key):
  1395. p = ss.find(key + "=")
  1396. return ss[p:].split("\"")[1] if p >= 0 else None
  1397. s = open(xml, "r").readlines()
  1398. ret = {}
  1399. for ss in s:
  1400. if ss.find("ThorZPiezo") >= 0:
  1401. ret["steps"] = int(get_value(ss, "steps"))
  1402. ret["stepSizeUM"] = float(get_value(ss, "stepSizeUM"))
  1403. ret["startPos"] = float(get_value(ss, "startPos"))
  1404. elif ss.find("<Streaming") >= 0:
  1405. ret["frames"] = int(get_value(ss, "frames"))
  1406. ret["zFastEnable"] = int(get_value(ss, "zFastEnable"))
  1407. ret["flybackFrames"] = int(get_value(ss, "flybackFrames"))
  1408. elif ss.find("<LSM") >= 0:
  1409. ret["name"] = get_value(ss, "name")
  1410. ret["pixelX"] = int(get_value(ss, "pixelX"))
  1411. ret["pixelY"] = int(get_value(ss, "pixelY"))
  1412. ret["pixelSizeUM"] = float(get_value(ss, "pixelSizeUM"))
  1413. ret["frameRate"] = float(get_value(ss, "frameRate")) #167.569
  1414. return ret
  1416. def load_ft_bar(bar_name):
  1417. if not os.path.exists(bar_name):
  1418. return None
  1419. import scipy.io as sio
  1420. bar_pos = sio.loadmat(bar_name)["bar_position"][0]
  1421. p = np.min(np.nonzero(np.isnan(bar_pos)))
  1422. return (bar_pos[:p] - 64) / 1648 * 2 * np.pi # 64~1712
  1424. def write_video(path, m, cvt=cv2.COLOR_GRAY2BGR, fps=30, need_time=True, sti_img_frame = None):
  1425. h, w = m[0].shape[:2]
  1426. output_video = cv2.VideoWriter(path, cv2.VideoWriter_fourcc(*"DIVX"), fps, (w, h))
  1427. for i, mi in enumerate(m):
  1428. #img_bgr = cv2.cvtColor(mi, cvt)
  1429. img_bgr = mi
  1430. if need_time:
  1431. cv2.putText(img_bgr, format_time(i/fps), (6, 20), cv2.FONT_HERSHEY_PLAIN, 1, (255, 255, 0), 1)
  1432. if sti_img_frame is not None:
  1433. for j in range(len(sti_img_frame)):
  1434. if i>sti_img_frame[j][0] and i < sti_img_frame[j][1]:
  1435. img = np.ones((20, 20,3), dtype=np.uint8) * 255
  1436. img_bgr[0:20, 0:20, :] = img
  1437. pass
  1438. output_video.write(img_bgr)
  1439. output_video.release()
  1441. def embed_video(v1, v2, v3, t2, t3, offset2=(0, 320), offset3=(160, 320)):
  1442. cap1 = cv2.VideoCapture(v1)
  1443. cap2 = cv2.VideoCapture(v2)
  1444. cap3 = cv2.VideoCapture(v3)
  1445. n1 = int(cap3.get(cv2.CAP_PROP_FRAME_COUNT))
  1446. cap2.set(cv2.CAP_PROP_POS_FRAMES, t2 * FICTRAC_RATE)
  1447. cap3.set(cv2.CAP_PROP_POS_FRAMES, t3 * FICTRAC_RATE)
  1448. w, h = int(cap1.get(cv2.CAP_PROP_FRAME_WIDTH)), int(cap1.get(cv2.CAP_PROP_FRAME_HEIGHT))
  1449. output_video = cv2.VideoWriter(v2 + "_combine.avi", cv2.VideoWriter_fourcc(*"DIVX"), FICTRAC_RATE, (w, h))
  1450. for seq in range(0, n1):
  1451. ret, img1 = cap1.read()
  1452. if not ret:
  1453. break
  1454. ret, img2 = cap2.read()
  1455. if not ret:
  1456. break
  1457. ret, img3 = cap3.read()
  1458. img3 = resize_image(img3,320,320)
  1459. if not ret:
  1460. break
  1461. h, w, c = img2.shape
  1462. img1[offset2[0]:offset2[0]+h, offset2[1]:offset2[1]+w, :] = img2
  1463. h, w, c = img3.shape
  1464. img1[offset3[0]:offset3[0]+h, offset3[1]:offset3[1]+w, :] = img3
  1465. output_video.write(img1)
  1466. output_video.release()
  1468. def resize_image(image, width, height,COLOUR=[0,0,0]):
  1469. h, w, layers = image.shape
  1470. if h > height:
  1471. ratio = height/h
  1472. image = cv2.resize(image,(int(image.shape[1]*ratio),int(image.shape[0]*ratio)))
  1473. h, w, layers = image.shape
  1474. if w > width:
  1475. ratio = width/w
  1476. image = cv2.resize(image,(int(image.shape[1]*ratio),int(image.shape[0]*ratio)))
  1477. h, w, layers = image.shape
  1478. if h < height and w < width:
  1479. hless = height/h
  1480. wless = width/w
  1481. if(hless < wless):
  1482. image = cv2.resize(image, (int(image.shape[1] * hless), int(image.shape[0] * hless)))
  1483. else:
  1484. image = cv2.resize(image, (int(image.shape[1] * wless), int(image.shape[0] * wless)))
  1485. h, w, layers = image.shape
  1486. if h < height:
  1487. df = height - h
  1488. df /= 2
  1489. image = cv2.copyMakeBorder(image, int(df), int(df), 0, 0, cv2.BORDER_CONSTANT, value=COLOUR)
  1490. if w < width:
  1491. df = width - w
  1492. df /= 2
  1493. image = cv2.copyMakeBorder(image, 0, 0, int(df), int(df), cv2.BORDER_CONSTANT, value=COLOUR)
  1494. image = cv2.resize(image,(width,height),interpolation=cv2.INTER_AREA)
  1495. return image
  1498. def merge_z_slices(fname, n, method=None):
  1499. import tifffile
  1500. with tifffile.TiffFile(fname) as tffl:
  1501. input_arr = tffl.asarray() # 25200*128*128
  1502. # write_video(fname + ".avi", norm_img(input_arr))
  1503. # plt.figure()
  1504. # a = input_arr[0].flatten()
  1505. # plt.hist(a[a>0.01], bins=50)
  1506. # plt.savefig(fname + "_hist.png")
  1507. chs = []
  1508. tn = int(len(input_arr) / n)
  1509. #fname_ch_list = []
  1510. for i in range(n):
  1511. ch0 = input_arr[i:tn*n:n]
  1512. # avg_frame = np.mean(ch0, axis=0)
  1513. # w = 255.0 / np.max(avg_frame)
  1514. # print(w)
  1515. eq = ch0 #* w*w
  1516. # eq[eq > 0] = 255
  1517. chs.append(eq)
  1518. #tifffile.imwrite(fname + "_ch%d.tif" % i, eq)
  1519. #fname_ch_list.append(fname + "_ch%d.tif" % i)
  1520. # write_video(fname + "_ch%d.avi"%i, norm_img(eq))
  1521. cv2.imwrite(fname + "_avg%d.png" % i, norm_img(np.mean(eq, axis=0)))
  1522. volu_arr = input_arr.reshape((np.shape(ch0)[0],np.shape(ch0)[1],np.shape(ch0)[2],n))
  1523. tifffile.imwrite(fname + "_ch_all.tif", volu_arr)
  1524. if method == "max":
  1525. avg_tif = np.max(chs, axis=0)
  1526. else:
  1527. avg_tif = np.mean(chs, axis=0)
  1528. tifffile.imwrite(fname + "_avg.tif", avg_tif)
  1529. return fname + "_avg.tif", fname + "_ch_all.tif"
  1531. def merge_info(date_folder):
  1532. info_l = []
  1533. for f in glob(date_folder + "/*/info.txt"):
  1534. c = json.load(open(f, "r"))
  1535. print(f, c)
  1536. exp_name = os.path.basename(os.path.dirname(f))
  1537. t = exp_name.split("-")
  1538. pair = "-".join(t[1:-1])
  1539. trail = t[-1]
  1540. stim = t[-2]
  1541. part = t[-3]
  1542. fly = t[2]
  1543. info_l.append([pair, fly, trail, part, stim, c["cor_unwrap"], c["cor_cir"]])
  1544. df = pd.DataFrame.from_records(info_l, columns=["pair", "fly", "trail", "part", "stim", "cor_unwrap", "cor_cir"])
  1545. import seaborn as sns
  1546. for x in ["cor_unwrap", "cor_cir"]:
  1547. plt.figure()
  1548. sns.set_theme(style="darkgrid")
  1549. sns.catplot(x=x, y="stim", hue="fly", col="part", data=df)
  1550. plt.xlim(-1, 1)
  1551. plt.tight_layout()
  1552. plt.savefig(date_folder + "/" + x)
  1554. def smooth_angle(s, win):
  1555. hw = int(win/2)
  1556. s1 = np.concatenate([[s[0]] * hw, s, [s[-1]] * hw])
  1557. ret = []
  1558. for i in range(len(s)):
  1559. m = circmean(s1[i:i+win])
  1560. if m > np.pi:
  1561. m -= np.pi * 2
  1562. ret.append(m)
  1563. return ret
  1565. def proc_ft_pv(fname):
  1566. c = json.load(open(fname, "r"))
  1567. t_range, heading, pv, pvl, speed, vx, vy, vx, pv2, pvl2 = c
  1568. n = len(heading)
  1569. ts = np.linspace(t_range[0], t_range[1], n)
  1570. # plot_lines(c[:7], "heading pv pvl speed vx vy vz".split())
  1571. pvs = smooth_angle(pv, 10)
  1572. pvs2 = smooth_angle(pv2, 10)
  1574. fig, axs = plt.subplots(3, 1, figsize=(8, 4), sharex=True, dpi=300)
  1575. # plt.subplots_adjust(left=0.02, right=0.99)
  1576. axs[0].plot(ts, unwrap_dir(np.array(heading)), c="m", lw=0.5)
  1577. axs[0].plot(ts, unwrap_dir(pvs), "k--", lw=0.5)
  1578. axs[0].plot(ts, unwrap_dir(pv), c="k", lw=0.5)
  1579. # axs[0].plot(ts, unwrap_dir(pvs2), c="gray", lw=0.5)
  1580. print(cor(unwrap_dir(np.array(heading)), unwrap_dir(pvs)))
  1581. # axs[0].plot(ts, unwrap_dir(pv2), c="gray", lw=0.5)
  1582. # axs[1].plot(ts, -np.array(heading), c="m", lw=0.5)
  1583. axs[1].plot(ts, pv, c="k", lw=0.5)
  1584. # axs[1].plot(ts, pvs, c="g", lw=0.5)
  1585. axs[1].plot(ts, pv2, c="gray", lw=0.5)
  1586. axs[2].plot(ts, zscore1(speed), c="m", lw=0.5)
  1587. axs[2].plot(ts, zscore1(pvl), c="k", lw=0.5)
  1588. axs[2].plot(ts, zscore1(pvl2), c="gray", lw=0.5)
  1589. plt.tight_layout()
  1590. # plt.show()
  1591. plt.savefig(fname + ".png")
  1592. """
  1593. // frame_count
  1594. ss << _cnt << ", ";
  1595. // rel_vec_cam[3] | error
  1596. ss << _dr_cam[0] << ", " << _dr_cam[1] << ", " << _dr_cam[2] << ", " << _err << ", ";
  1597. // rel_vec_world[3]
  1598. ss << _dr_lab[0] << ", " << _dr_lab[1] << ", " << _dr_lab[2] << ", ";
  1599. // abs_vec_cam[3]
  1600. ss << _r_cam[0] << ", " << _r_cam[1] << ", " << _r_cam[2] << ", ";
  1601. // abs_vec_world[3]
  1602. ss << _r_lab[0] << ", " << _r_lab[1] << ", " << _r_lab[2] << ", ";
  1603. // integrated xpos | integrated ypos | integrated heading
  1604. ss << _posx << ", " << _posy << ", " << _heading << ", ";
  1605. // direction (radians) | speed (radians/frame)
  1606. ss << _step_dir << ", " << _step_mag << ", ";
  1607. // integrated x movement | integrated y movement (mouse output equivalent)
  1608. ss << _intx << ", " << _inty << ", ";
  1609. // timestamp | sequence number
  1610. ss << _ts << ", " << _seq << std::endl;
  1611. """
  1613. if __name__ == '__main__':
  1614. ft, pv = json.load(open(r"D:\exp_2p\data\EPG\210823\CX1001-7F-F6-PB-reverseCLOSED-7+\ft_pv.txt", "r"))
  1615. print(cor(ft, pv))
  1616. print(cor(unwrap_dir(ft), unwrap_dir(pv)))
  1617. print(cir_cor(ft, pv))

cia_utils.py at commit 9c55e7a, no license · at the source

Overview

Authors: Qiantao Lv1,2,3,4, Yunming Gao2,3,4, Jing Ning2,3,4, Yi Sun1,2,3,4
ORCID iDs: Yi Sun
  1. Zhejiang University, Hangzhou, Zhejiang 310058, China
  2. Westlake Laboratory of Life Sciences and Biomedicine, Hangzhou, Zhejiang 310024, China
  3. Key Laboratory of Growth Regulation and Translational Research of Zhejiang Province, School of Life Sciences, Westlake University, Hangzhou, Zhejiang 310030, China
  4. Institute of Basic Medical Sciences, Westlake Institute for Advanced Study, Hangzhou, Zhejiang 310024, China
Institutions: Westlake University (China); Zhejiang University (China)
Journal: STAR protocols, volume 7, issue 3, article 104659
Dates: published online 25 June 2026; in print June 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1016/j.xpro.2026.104659 · PMID 42348421 · PMCID PMC13333336 · OpenAlex W7165931131
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: histology / microscopy (modality), optical imaging (calcium, voltage, 2-photon) (modality), drosophila (organism), clinical / translational (subfield)
Methods: Spectral & time-frequency, Connectivity, Statistics, fMRI & imaging, Single-unit activity, calcium imaging, Smoothing, state filtering, decompositions
Keywords: Microscopy, Model Organisms, Neuroscience
Topic: Neurobiology and Insect Physiology Research (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Funding: Westlake University (WU2022C011); Westlake Laboratory of Life Sciences and Biomedicine (W101386022101); Westlake Education Foundation
Citations: not cited yet (Europe PMC); 15 references in the paper
Research resources: Python RRID:SCR_008394, Illustrator RRID:SCR_010279

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

SunLabWestlake/cia

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 9c55e7ac2ebd9ea6b11ccb8a013fa9d76119f856, 4 June 2026
Languages: Python (7)
Size: 10 files, 7 scripts
Software Heritage: not archived
Found in: “Data and code availability”
Holds: README, environment (requirements.txt)
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: Matplotlib (5 files), NumPy (5 files), pandas (5 files), h5py (3 files), OpenCV (3 files), SciPy (3 files), CaImAn (2 files), imageio (1 file), Pillow (1 file), seaborn (1 file), tifffile (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
9 files
At the source: github.com/SunLabWestlake/cia

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 7 scripts, each with its path and the digest of its content;
  • 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1016/j.xpro.2026.104659.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 3 keywords, 3 funders, 15 references, 2 RRIDs.

Cite

This paper

Lv, Q., Gao, Y., Ning, J., & Sun, Y. (2026). Protocol for simultaneous in vivo two-photon imaging and locomotion quantification during olfactory stimulation and pharmacology in walking Drosophila. STAR protocols, 7(3), 104659. https://doi.org/10.1016/j.xpro.2026.104659

BibTeX

@article{lv2026protocol,
author = {Lv, Qiantao and Gao, Yunming and Ning, Jing and Sun, Yi},
title = {{Protocol for simultaneous in vivo two-photon imaging and locomotion quantification during olfactory stimulation and pharmacology in walking Drosophila}},
journal = {STAR protocols},
year = {2026},
month = jun,
volume = {7},
number = {3},
pages = {104659},
publisher = {Elsevier},
issn = {2666-1667},
doi = {10.1016/j.xpro.2026.104659},
url = {https://doi.org/10.1016/j.xpro.2026.104659},
pmid = {42348421},
pmcid = {PMC13333336}
}

RIS

TY - JOUR
AU - Lv, Qiantao
AU - Gao, Yunming
AU - Ning, Jing
AU - Sun, Yi
TI - Protocol for simultaneous in vivo two-photon imaging and locomotion quantification during olfactory stimulation and pharmacology in walking Drosophila
T2 - STAR protocols
J2 - STAR Protoc
PY - 2026
DA - 2026/06/25
VL - 7
IS - 3
SP - 104659
SN - 2666-1667
PB - Elsevier
DO - 10.1016/j.xpro.2026.104659
UR - https://doi.org/10.1016/j.xpro.2026.104659
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.xpro.2026.104659",
"type": "article-journal",
"title": "Protocol for simultaneous in vivo two-photon imaging and locomotion quantification during olfactory stimulation and pharmacology in walking Drosophila",
"container-title": "STAR protocols",
"author": [
{
"family": "Lv",
"given": "Qiantao"
},
{
"family": "Gao",
"given": "Yunming"
},
{
"family": "Ning",
"given": "Jing"
},
{
"family": "Sun",
"given": "Yi"
}
],
"container-title-short": "STAR Protoc",
"volume": "7",
"issue": "3",
"page": "104659",
"DOI": "10.1016/j.xpro.2026.104659",
"PMID": "42348421",
"PMCID": "PMC13333336",
"ISSN": "2666-1667",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.xpro.2026.104659",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
25
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-72437-1 [code]
High-speed whole-brain imaging in Drosophila.
Journal: Nature communications
In common: OpenCV, h5py, Pillow, 4 other tools, drosophila, 6 references
[2] doi:10.1016/j.isci.2026.116206 [code]
Gut distension evokes rapid neural dynamics in vagal and hindbrain populations of larval zebrafish.
Journal: iScience
In common: CaImAn, tifffile, OpenCV, 6 other tools, optical imaging (calcium, voltage, 2-photon), 2 references
[3] doi:10.1126/sciadv.adv3770 [code]
Evolution of a central dopamine circuit underlies adaptation of a light-evoked sensorimotor response in the blind cavefish.
Journal: Science advances
In common: CaImAn, tifffile, OpenCV, 7 other tools, 1 reference
[4] doi: [code]
Real-time closed-loop feedback system for mouse mesoscale cortical signal and movement control
Journal: eLife
In common: imageio, tifffile, OpenCV, 7 other tools, optical imaging (calcium, voltage, 2-photon)
[5] doi:10.1038/s41467-026-72709-w [code]
An epifluorescence microscope design for naturalistic behavior and cellular activity in freely moving Caenorhabditis elegans.
Journal: Nature communications
In common: imageio, tifffile, OpenCV, 7 other tools, histology / microscopy
[6] doi:10.1364/boe.600665 [code]
NeuroSeg-MF: robust neuron segmentation in two-photon Ca&lt;sup&gt;2+&lt;/sup&gt; imaging using multi-feature fusion and detection-guided SAM.
Journal: Biomedical optics express
In common: tifffile, OpenCV, Pillow, 5 other tools, optical imaging (calcium, voltage, 2-photon), 2 references
[7] doi:10.1364/boe.605322 [code]
Generalized plaque digitization framework for multi-dimensional mesoscopic images.
Journal: Biomedical optics express
In common: imageio, tifffile, OpenCV, 7 other tools
[8] doi:10.7554/elife.109717 [code]
Retrosplenial cortex enables context-dependent goal-directed sensorimotor transformation.
Journal: eLife
In common: imageio, tifffile, OpenCV, 7 other tools
[9] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: imageio, tifffile, OpenCV, 7 other tools
[10] doi:10.1038/s41598-026-57519-w [code]
Automated segmentation of neurons and spinal cord structures in immunofluorescence images using SpineDL.
Journal: Scientific reports
In common: imageio, tifffile, OpenCV, 6 other tools, histology / microscopy

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.