Compute the optical crosstalk in two 48-pixel SPAD arrays from the 48-spot smFRET-PAX setup.
fname = 'data/2017-10-16_00_DCR.hdf5'
fname
'data/2017-10-16_00_DCR.hdf5'
from pathlib import Path
assert Path(fname).is_file(), 'File not found.'
mlabel = Path(fname).stem
mlabel
'2017-10-16_00_DCR'
from itertools import combinations
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['font.sans-serif'].insert(0, 'Arial')
plt.rcParams['font.size'] = 14
from tqdm import tnrange, tqdm_notebook
from fretbursts import *
- Optimized (cython) burst search loaded. - Optimized (cython) photon counting loaded. -------------------------------------------------------------- You are running FRETBursts (version 0.6.5). If you use this software please cite the following paper: FRETBursts: An Open Source Toolkit for Analysis of Freely-Diffusing Single-Molecule FRET Ingargiola et al. (2016). http://dx.doi.org/10.1371/journal.pone.0160716 --------------------------------------------------------------
import time
time.ctime()
'Mon Oct 30 16:50:38 2017'
def coincidence_py(timestamps1, timestamps2):
coinc = 0
i1, i2 = 0, 0
while (i1 < timestamps1.size) and (i2 < timestamps2.size):
if timestamps1[i1] == timestamps2[i2]:
coinc += 1
i1 += 1
i2 += 1
elif timestamps1[i1] > timestamps2[i2]:
i2 += 1
elif timestamps1[i1] < timestamps2[i2]:
i1 += 1
return coinc
%load_ext Cython
%%cython
cimport numpy as np
def coincidence(np.int64_t[:] timestamps1, np.int64_t[:] timestamps2):
cdef np.int64_t coinc, i1, i2, size1, size2
size1 = timestamps1.size
size2 = timestamps2.size
coinc = 0
i1, i2 = 0, 0
while i1 < size1 and i2 < size2:
if timestamps1[i1] == timestamps2[i2]:
coinc += 1
i1 += 1
i2 += 1
elif timestamps1[i1] > timestamps2[i2]:
i2 += 1
elif timestamps1[i1] < timestamps2[i2]:
i1 += 1
return coinc
def crosstalk_probability(t1, t2, tol=1e-6, max_iter=100):
"""Estimate crosstalk probability between two pixels in a SPAD array.
Given two input arrays of timestamps `t1` and `t2`, estimate
the crosstalk probability using Poisson statistics without
approximations.
Arguments:
t1, t2 (array of ints): arrays of timestamps from DCR measurements
for the two pixels to be measured. Timestamps need to be
integers and coincidences are detected when values in the two
arrays are equal. These arrays need to be rescaled, if
coincidence need to be computed on a delta t larger than \
a single timestamp unit.
tol (float): tollerance for iterative equasion solution
max_iter (int): max iterations used to solve the equation
Returns:
A 3-element tuple:
- crosstalk probability
- crosstalk probability standard deviation
- number of iterations used for the estimation.
"""
T = (max((t1.max(), t2.max())) - min((t1.min(), t2.min())))
C = coincidence(t1, t2)
# Find C_c by solving eq. (1) iteratively
C_c, C_u_prev = 0, 0
for i in range(max_iter):
C_u = ((1 - np.exp(-(t1.size - C_c)/T)) *
(1 - np.exp(-(t2.size - C_c)/T)) * T)
C_c = C - C_u
if np.abs(C_u - C_u_prev) < tol:
break
C_u_prev = C_u
P_c = C_c / (t1.size + t2.size - C_c)
if C_c <= 0:
sigma = np.nan
else:
sigma = np.sqrt(C_c) / (t1.size + t2.size - C_c)
return P_c, sigma, i
def crosstalk(dx, spot1, spot2, divide=1, det=0):
"""
Calculate the crosstalk between two pixels in a SPAD array for data object, dx.
Arguments:
spot1, spot2 (ints): pixel pair for
combinations possible in 48 spots. The pixel pair
tuple needs to be intergers. The optical crosstalk
across pixels is computed for a pair of pixels.
divide (int): integer division of timestamps by `divide` are the
timestamps used for coincidence counting.
`divide` rescales the timestamps so that adjacent timestamps
become the same timestamp after division.
det (int): 0, 1 for D, A SPADs respectively.
Returns:
A 3-element tuple:
- A float corresponding to optical crosstalk between the
selected pixels in the D/A SPAD array.
- A float corresponding to the standard deviation of optical crosstalk
between the selected pixels in the D/A SPAD array.
- A float indicating the probability of optical crosstalk between
selected pixels.
"""
ph1 = dx.ph_times_m[spot1][dx.detectors[spot1] == det] // divide
ph2 = dx.ph_times_m[spot2][dx.detectors[spot2] == det] // divide
delta_t = d.clk_p*divide
if ph1.size == 0 or ph2.size == 0:
return np.nan, np.nan, 0
ct, sigma, i = crosstalk_probability(ph1, ph2)
return ct, sigma, i
def dist(manta_shape, ich1, ich2):
"""Compute distance between two pixels on the SPAD array.
"""
row1, col1 = np.where(manta_shape == ich1)
row2, col2 = np.where(manta_shape == ich2)
del_x = abs(col2 - col1)
del_y = abs(row2 - row1)
dist = float(np.sqrt((del_x)**2 + (del_y)**2)) #multipy by 500 to convert to microns
return(dist)
def savefig(fname, **kwargs):
"""Save figure with default arguments."""
import os
basename = os.path.basename(fname)
dir_ = f'figures/{mlabel}_'
kwargs_ = dict(dpi=180, bbox_inches='tight',
frameon=True, facecolor='white', transparent=False)
kwargs_.update(kwargs)
plt.savefig(dir_ + basename, **kwargs_)
print('Saved: ', dir_ + basename)
d = loader.photon_hdf5(fname)
spot = 31
det = 0
ph = d.ph_times_m[spot][d.detectors[spot] == det] * d.clk_p
counts, bins = np.histogram(ph, bins=np.arange(0, d.time_max, 5))
fig, ax = plt.subplots(figsize=(24, 4))
ax.plot(bins[1:], counts, label='Det %d, spot %d' % (det, spot))
plt.legend()
display(fig)
plt.close(fig)
d.detectors[0].shape
(6846745,)
d.ph_times_m[0].shape
(6846745,)
manta_shape = np.arange(48).reshape(4, 12)[::-1].T[::-1]
manta_shape
array([[47, 35, 23, 11], [46, 34, 22, 10], [45, 33, 21, 9], [44, 32, 20, 8], [43, 31, 19, 7], [42, 30, 18, 6], [41, 29, 17, 5], [40, 28, 16, 4], [39, 27, 15, 3], [38, 26, 14, 2], [37, 25, 13, 1], [36, 24, 12, 0]])
NOTE: With
combinations(range(48), 2)
we generate all the unique combinations (i.e. no permutations) of 48 pixel IDs. See itertools.combinations for more info.
We start by computing all paris of pixels (pix_pairs
) and all the distances for all the pairs of pixels (distances
):
pix_pairs = []
distances = []
for ich1, ich2 in combinations(range(48), 2):
pix_pairs.append((ich1, ich2))
distances.append(dist(manta_shape, ich1, ich2))
Then we compute the crosstalk for all pairs of pixes on both arrays.
Results are stored in two DataFrame (i.e. tables) called pairdata0
and pairdata1
for donor and acceptor detector respectively.
def compute_crosstalk_detector(det):
xtalk_err_list = []
for ich1, ich2 in tqdm_notebook(combinations(range(48), 2), total=len(distances),
desc=f'Det. {det}', leave=True):
res = crosstalk(d, ich1, ich2, divide=4, det=det)
xtalk_err_list.append(res)
if res[2] > 50:
print('i', ich1, ich2, flush=True)
if np.isnan(res[1]):
print(f' W: Crosstalk for detector {det} pair {ich1}, {ich2} is NaN.', flush=True)
xtalk_err_array = np.array(xtalk_err_list)
pairdata = pd.DataFrame(columns=['crosstalk', 'distance', 'error'])
pairdata['crosstalk'] = xtalk_err_array[:, 0]
pairdata['error'] = xtalk_err_array[:, 1]
pairdata['distance'] = distances
pairdata['pix1'] = np.array(pix_pairs)[:,0]
pairdata['pix2'] = np.array(pix_pairs)[:,1]
return pairdata
def compute_crosstalk_both_dectectors(recompute=False):
res = []
for det in (0, 1):
fname_crosstalk = f'results/{mlabel}_crosstalk_data_detector{det}.csv'
if Path(fname_crosstalk).exists():
print(f'- Loading crosstalk for detector {det} from cache', flush=True)
pairdata = pd.read_csv(fname_crosstalk, index_col=0)
else:
print(f'- Computing crosstalk for detector {det}', flush=True)
pairdata = compute_crosstalk_detector(det)
pairdata.to_csv(fname_crosstalk)
res.append(pairdata)
return res
pairdata0, pairdata1 = compute_crosstalk_both_dectectors(recompute=False)
pairdata = pairdata1
- Loading crosstalk for detector 0 from cache - Loading crosstalk for detector 1 from cache
pairdata.head()
crosstalk | distance | error | pix1 | pix2 | |
---|---|---|---|---|---|
0 | 0.000772 | 1.0 | 1.481646e-05 | 0 | 1 |
1 | 0.000014 | 2.0 | 2.193197e-06 | 0 | 2 |
2 | 0.000007 | 3.0 | 2.729624e-06 | 0 | 3 |
3 | 0.000002 | 4.0 | 8.092910e-07 | 0 | 4 |
4 | 0.000010 | 5.0 | 3.699850e-06 | 0 | 5 |
pairdata.plot(x='distance', y='crosstalk', kind='scatter', alpha=0.3, marker='H', s=30,
logy=True, ylim=(2e-8, 10e-3));
plt.ylabel(u'crosstalk prob.')
plt.xlabel(u'Distance (unit = 500 μm)')
plt.title('Crosstalk vs distance')
Text(0.5,1,'Crosstalk vs distance')
Optical crosstalk between pixels is maximal at shorter distances.
pairdata0.crosstalk.max(), pairdata1.crosstalk.max()
(0.0011385712427402512, 0.001042788265769392)
Optical crosstalk along the x,y axes is determined by selecting a pair of neighboring pixels.
dist = 1
mask = pairdata['distance'] == dist
plt.errorbar(x=np.arange(mask.sum()), y=pairdata['crosstalk'][mask]*1e3, yerr=3*pairdata['error'][mask]*1e3, fmt='H')
plt.title('Crosstalk ($\pm 3\sigma$), pair distance = %.1f pixels' % dist);
plt.ylabel(r'Probability $\times \; 10^{-3}$')
Text(0,0.5,'Probability $\\times \\; 10^{-3}$')
Optical crosstalk along the diagonal is determined by selecting a pair of cater-cornered pixels and calculating the distance along the diagonal with the Pythagorian Theorem.
dist = np.sqrt(2)
mask = pairdata['distance'] == dist
plt.errorbar(x=np.arange(mask.sum()), y=pairdata['crosstalk'][mask]*1e3, yerr=3*pairdata['error'][mask]*1e3, fmt='H')
plt.title('Crosstalk ($\pm 3\sigma$), pair distance = %.1f pixels' % dist);
plt.ylabel(r'Probability $\times \; 10^{-3}$');
Here, optical crosstalk along the x,y axes is calculated for a pair of pixels separated by a distance equal to two pixels.
dist = 2
mask = pairdata['distance'] == dist
plt.errorbar(x=np.arange(mask.sum()), y=pairdata['crosstalk'][mask]*1e3, yerr=3*pairdata['error'][mask]*1e3, fmt='H')
plt.title('Crosstalk ($\pm 3\sigma$), pair distance = %.1f pixels' % dist);
plt.ylabel(r'Probability $\times \; 10^{-3}$');
def pixel_rowcol_to_ch(d, row, col):
detectors = d.setup['detectors']
i = np.where((detectors['position'] == (row, col)).sum(axis=1) == 2)[0][0]
ich = detectors['spot'][i]
return ich
def pixel_ch_to_rowcol(d, ich):
detectors = d.setup['detectors']
ich = np.where(detectors['spot'] == ich)
pix_pos = detectors['position'][ich][0]
return pix_pos
def pair_coord(pix1, pix2):
"""
Returns (i, j) position in the expanded `a` array used to generate
heatmap of optical crosstalk for a pair of pixels
`pix1` and `pix2` (pixel ids, [0, 47]).
"""
v1 = pixel_ch_to_rowcol(d, pix1)
v2 = pixel_ch_to_rowcol(d, pix2)
return (v1[0] + v2[0]), (v1[1] + v2[1])
nshape = len(manta_shape[0,:])*2 - 1, len(manta_shape[:,1])*2 - 1
nshape
(7, 23)
a = np.zeros(nshape)
a[1::2, 1::2] = np.nan
a[0::2, 0::2] = np.nan
plt.imshow(a, interpolation='none', cmap='YlGnBu')
plt.grid(False)
dist = 1
mask = pairdata['distance'] == dist
pairdata_d1 = pairdata.loc[mask]
pairdata_d1.head()
crosstalk | distance | error | pix1 | pix2 | |
---|---|---|---|---|---|
0 | 0.000772 | 1.0 | 0.000015 | 0 | 1 |
11 | 0.000893 | 1.0 | 0.000036 | 0 | 12 |
47 | 0.000892 | 1.0 | 0.000013 | 1 | 2 |
58 | 0.000913 | 1.0 | 0.000016 | 1 | 13 |
93 | 0.000958 | 1.0 | 0.000018 | 2 | 3 |
for p, pair in pairdata_d1.iterrows():
i, j = pair_coord(pair.pix1, pair.pix2)
a[j, i] = pair.crosstalk
#print(i, j, pair.pix1, pair.pix2, pair.crosstalk)
dist = np.sqrt(2)
mask = pairdata['distance'] == dist
pairdata_diag = pairdata.loc[mask]
pairdata_diag
crosstalk | distance | error | pix1 | pix2 | |
---|---|---|---|---|---|
12 | 0.000140 | 1.414214 | 0.000014 | 0 | 13 |
57 | 0.000138 | 1.414214 | 0.000006 | 1 | 12 |
59 | 0.000146 | 1.414214 | 0.000006 | 1 | 14 |
103 | 0.000139 | 1.414214 | 0.000007 | 2 | 13 |
105 | 0.000134 | 1.414214 | 0.000003 | 2 | 15 |
... | ... | ... | ... | ... | ... |
1033 | 0.000151 | 1.414214 | 0.000004 | 33 | 44 |
1035 | 0.000195 | 1.414214 | 0.000011 | 33 | 46 |
1047 | 0.000148 | 1.414214 | 0.000002 | 34 | 45 |
1049 | 0.000139 | 1.414214 | 0.000001 | 34 | 47 |
1060 | 0.000163 | 1.414214 | 0.000008 | 35 | 46 |
66 rows × 5 columns
for p, pair in pairdata_diag.iterrows():
i = pair_coord(pair.pix1, pair.pix2)[0]
j = pair_coord(pair.pix1, pair.pix2)[1]
a[j, i] = pair.crosstalk
#print(i, j, pair.pix1, pair.pix2, pair.crosstalk)
m = manta_shape.T[::-1]
m
array([[11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0], [23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12], [35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24], [47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36]])
fig, ax = plt.subplots(figsize=(16, 3.5))
xs = np.arange(12)
ys = np.repeat(np.arange(4)[np.newaxis,:].T, xs.size, axis=1)
plt.plot(xs, ys.T, 'o', ms=20, mew=1, mec='k', color='white')
im = plt.imshow(a*1e2, interpolation='none', cmap='Blues', vmin=0, vmax=0.15,
extent=(-0.25, 11.25, -0.25, 3.25))
for row in range(4):
for col in range(12):
ax.text(col+0.01, row-0.01, str(m[row,col]), va='center', ha='center', fontsize=12)
im.cmap.set_under(alpha=0)
plt.colorbar()
plt.xticks(range(12))
plt.yticks(range(4))
plt.grid(False)
def reflection_distance(dist, thickness=1):
return np.sqrt((dist/2)**2 + thickness*2)*2
fig, AX = plt.subplots(2, 3, figsize=(20, 10), sharex='col', sharey='col')
ax1a = AX[0, 0]
ax2a = AX[0, 1]
ax3a = AX[0, 2]
ax1d = AX[1, 0]
ax2d = AX[1, 1]
ax3d = AX[1, 2]
plt.subplots_adjust(hspace=0.1, wspace=0.3)
xtext, ytext = 0.95, 0.95
font_kw = dict(fontsize=26, va='top', ha='right')
pairdata = pairdata1
ax1a.scatter('distance', 'crosstalk', data=pairdata.loc[pairdata.distance == 1], alpha=0.3, s=30, c='C0')
ax1a.scatter('distance', 'crosstalk', data=pairdata.loc[pairdata.distance == np.sqrt(2)], alpha=0.3, s=30, c='C1')
ax1a.scatter('distance', 'crosstalk', data=pairdata.loc[pairdata.distance > np.sqrt(2)], alpha=0.2, s=30, c='grey')
ax1a.set_yscale('log')
ax1a.set_ylim(2e-7, 2e-3)
ax1a.set_ylabel(r'Crosstalk Probability')
#ax1a.set_xlabel(u'Distance (unit = 500 μm)')
ax1a.set_title('Crosstalk vs distance')
ax1a.text(xtext, ytext, 'A', transform=ax1a.transAxes, **font_kw)
x = np.unique(pairdata.distance)
A = pairdata.loc[pairdata.distance == 1, 'crosstalk'].mean()
y = A / x**2
ax1a.plot(x, y, ls='--', color='k')
ax1a.text(0.5, 0.58, r'$\propto R^{-2}$', transform=ax1a.transAxes, fontsize=24)
# A = pairdata.loc[pairdata.distance == np.sqrt(2), 'crosstalk'].mean() * 2
# y = A / x**2
# ax1a.plot(x, y, ls='--', color='k')
# A = 3e-3
# y = A / reflection_distance(x)**2
# ax1a.plot(x, y, ls='-', color='k')
dist = 1
mask = pairdata['distance'] == dist
ax2a.errorbar(x=np.arange(mask.sum()), y=pairdata['crosstalk'][mask]*1e3, yerr=3*pairdata['error'][mask]*1e3, fmt='H', color='C0')
ax2a.set_title('Horizontal/vertical pairs');
ax2a.set_ylabel(r'Crosstalk Probability $\times \; 10^{-3}$')
#ax2a.set_xlabel('Pixel pair')
ax2a.text(xtext, ytext, 'B', transform=ax2a.transAxes, **font_kw)
dist = np.sqrt(2)
mask = pairdata['distance'] == dist
ax3a.errorbar(x=np.arange(mask.sum()), y=pairdata['crosstalk'][mask]*1e3, yerr=3*pairdata['error'][mask]*1e3, fmt='H', color='C1')
ax3a.set_title('Diagonal pairs');
ax3a.set_ylabel(r'Crosstalk Probability $\times \; 10^{-3}$')
#ax3a.set_xlabel('Pixel pair')
ax3a.text(xtext, ytext, 'C', transform=ax3a.transAxes, **font_kw)
pairdata = pairdata0
ax1d.scatter('distance', 'crosstalk', data=pairdata.loc[pairdata.distance == 1], alpha=0.3, s=30, c='C0')
ax1d.scatter('distance', 'crosstalk', data=pairdata.loc[pairdata.distance == np.sqrt(2)], alpha=0.3, s=30, c='C1')
ax1d.scatter('distance', 'crosstalk', data=pairdata.loc[pairdata.distance > np.sqrt(2)], alpha=0.2, s=30, c='grey')
ax1d.set_yscale('log')
ax1d.set_ylim(2e-7, 2e-3)
ax1d.set_ylabel(r'Crosstalk Probability')
ax1d.set_xlabel(u'Distance (unit = 500 μm)')
#ax1d.set_title('Crosstalk vs distance')
ax1d.text(xtext, ytext, 'D', transform=ax1d.transAxes, **font_kw)
x = np.unique(pairdata.distance)
A = pairdata.loc[pairdata.distance == 1, 'crosstalk'].mean()
y = A / x**2
ax1d.plot(x, y, ls='--', color='k')
ax1d.text(0.5, 0.58, r'$\propto R^{-2}$', transform=ax1d.transAxes, fontsize=24)
dist = 1
mask = pairdata['distance'] == dist
ax2d.errorbar(x=np.arange(mask.sum()), y=pairdata['crosstalk'][mask]*1e3, yerr=3*pairdata['error'][mask]*1e3, fmt='H', color='C0')
#ax2d.set_title('Horizontal/vertical pairs');
ax2d.set_ylabel(r'Crosstalk Probability $\times \; 10^{-3}$')
ax2d.set_xlabel('Pixel pair')
ax2d.text(xtext, ytext, 'E', transform=ax2d.transAxes, **font_kw)
dist = np.sqrt(2)
mask = pairdata['distance'] == dist
ax3d.errorbar(x=np.arange(mask.sum()), y=pairdata['crosstalk'][mask]*1e3, yerr=3*pairdata['error'][mask]*1e3, fmt='H', color='C1')
#ax3d.set_title('Diagonal pairs');
ax3d.set_ylabel(r'Crosstalk Probability $\times \; 10^{-3}$')
ax3d.set_xlabel('Pixel pair')
ax3d.text(xtext, ytext, 'F', transform=ax3d.transAxes, **font_kw);
savefig('crosstalk_both_detectors')
Saved: figures/2017-10-16_00_DCR_crosstalk_both_detectors
def compute_heatmap(pairdata):
a = np.zeros(nshape)*np.nan
dist = 1
mask = pairdata['distance'] == dist
pairdata_d1 = pairdata.loc[mask]
for p, pair in pairdata_d1.iterrows():
i, j = pair_coord(pair.pix1, pair.pix2)
a[j, i] = pair.crosstalk
#print(i, j, pair.pix1, pair.pix2, pair.crosstalk)
dist = np.sqrt(2)
mask = pairdata['distance'] == dist
pairdata_diag = pairdata.loc[mask]
for p, pair in pairdata_diag.iterrows():
i, j = pair_coord(pair.pix1, pair.pix2)
a[j, i] = pair.crosstalk
#print(i, j, pair.pix1, pair.pix2, pair.crosstalk)
return a
def plot_crosstalk_heatmap(a, ax, title='Spatial distribution of crosstalk', cmap='Blues', vmin=0, vmax=1.5):
m = manta_shape.T[::-1]
xs = np.arange(12)
ys = np.repeat(np.arange(4)[np.newaxis,:].T, xs.size, axis=1)
ax.plot(xs, ys.T, 'o', ms=20, mew=1, mec='k', color='white')
im = ax.imshow(a*1e3, interpolation='none', cmap=cmap, vmin=vmin, vmax=vmax,
extent=(-0.25, 11.25, -0.25, 3.25))
for row in range(4):
for col in range(12):
ax.text(col+0.01, row-0.01, str(m[row,col]), va='center', ha='center', fontsize=12)
ax.set_xticks(range(12))
ax.set_yticks(range(4));
ax.set_title(title)
ax.text(1.14, 0.5, r'Crosstalk Probability $\times \; 10^{-3}$', rotation=90, va='center', transform=ax.transAxes)
return im
a0 = compute_heatmap(pairdata0)
a1 = compute_heatmap(pairdata1)
fig, (ax0, ax1) = plt.subplots(2, 1, figsize=(14, 8))
im1 = plot_crosstalk_heatmap(a1, ax1, title='Acceptor SPAD array', cmap='Oranges')
fig.colorbar(im1, ax=ax1)
im0 = plot_crosstalk_heatmap(a0, ax0, title='Donor SPAD array', cmap='Blues')
fig.colorbar(im0, ax=ax0)
savefig('crosstalk_heatmaps_both_detectors')
Saved: figures/2017-10-16_00_DCR_crosstalk_heatmaps_both_detectors