# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %matplotlib notebook
from colicoords import load, CellPlot, CellListPlot, IterCellPlot, iter_subplots, save
import matplotlib.pyplot as plt
import numpy as np
import os
from collections import namedtuple
from tqdm.auto import tqdm

from scipy.signal import medfilt
from skimage.feature import peak_local_max

c41_02 = load('c41_02_binary_opt.hdf5')
c41_03 = load('c41_03_binary_opt.hdf5')

# +
clp_c41_02 = CellListPlot(c41_02)
clp_c41_03 = CellListPlot(c41_03)

fig, ax = plt.subplots()
clp_c41_02.hist_intensity(ax=ax, data_name='g500', linewidth=0, label='c41_2')
clp_c41_03.hist_intensity(ax=ax, data_name='g500', linewidth=0, label='c41_3', alpha=0.75)

plt.savefig('intensity comparison.png')
plt.legend()
# -

storm_dtype = [('x', float), ('y', float), ('intensity', float), ('frame', int)]
def add_peakfind(cell, med=9, thd=7500, min_dst=5):
    img = cell.data.data_dict['g500']
    mf = medfilt(img, 9)
    img_bg = img - mf
    
    cell.data.add_data(img_bg, 'fluorescence', 'flu_mf')
    peaks = peak_local_max(img_bg, min_distance=min_dst, threshold_abs=thd)
    y, x = peaks.T
    
    new_storm = np.empty(len(x), dtype=storm_dtype)
    new_storm['x'] = x
    new_storm['y'] = y
    new_storm['intensity'] = np.ones_like(x)
    new_storm['frame'] = np.ones_like(x)
    
    cell.data.add_data(new_storm, 'storm', 'storm_thd_{}'.format(thd))


len(c41_02), len(c41_03)

# +
c41_02_new = c41_02.copy()
[add_peakfind(c) for c in tqdm(c41_02_new)]

c41_03_new = c41_03.copy()
[add_peakfind(c) for c in tqdm(c41_03_new)]


''

# +
icp = IterCellPlot(c41_02_new)
fig, axes = iter_subplots(2, 1, figsize=(8,6))

icp.imshow('g500', ax=axes[0])
icp.plot_storm(data_name='storm_thd_7500', ax=axes[0])

icp.imshow('flu_mf', ax=axes[1])
icp.plot_storm(data_name='storm_thd_7500', ax=axes[1])

plt.tight_layout()
fig.display()

# +
icp = IterCellPlot(c41_03_new)
fig, axes = iter_subplots(2, 1, figsize=(8,6))

icp.imshow('g500', ax=axes[0])
icp.plot_storm(data_name='storm_thd_7500', ax=axes[0])

icp.imshow('flu_mf', ax=axes[1])
icp.plot_storm(data_name='storm_thd_7500', ax=axes[1])

plt.tight_layout()
fig.display()

# +
labels = ['c41_02', 'c41_03']
fig, ax = plt.subplots()
nums = []
for cells in [c41_02_new, c41_03_new]:
    num = [len(c.data.data_dict['storm_thd_7500']) for c in cells]
    nums.append(num)

ax.hist(nums, bins = np.arange(15), label=labels, density=True)

ax.legend()
    
#fig.text(0.04, 0.5, 'Number of spots', va='center', rotation='vertical')
plt.ylabel('Fraction of cells')
plt.xlabel('Number of spots')
plt.savefig('spots per cell_c41 epec escc.png')
# -

save('c41_02_with_spots.hdf5', c41_02_new)
save('c41_03_with_spots.hdf5', c41_03_new)

StarcoderdataJupyterScriptsDedupFiltered