<jupyter_start><jupyter_text>Harringtonine CropArray Example ---  Notebook summary - Load a microscope image of video- Tracking spots on the image and generate a pandas dataframe with the spots locations- Creating a croparray with the image and dataframe- Signal quantification and plotting- Visualization of croparray with Napari----  Importing libraries----<jupyter_code># To manipulate arrays
import numpy as np 
from skimage.io import imread
import matplotlib.pyplot as plt 
from matplotlib.path import Path
import pylab as pyl
import seaborn as sns; sns.set()  
import pathlib # for working with windows paths
import sys
import cv2
import trackpy as tp
#!pip install shapely
from shapely.geometry import Polygon 
from shapely.geometry import Point
current_dir = pathlib.Path().absolute()
croparray_dir = current_dir.parents[0].joinpath('croparray')
sys.path.append(str(croparray_dir))
import crop_array_tools as ca
# %matplotlib inline 
plt.style.use('dark_background')
# Napari 
%gui qt5 
import napari
from napari.utils import nbscreenshot
# Magicgui
from magicgui import magicgui
import datetime
import pathlib
import pandas as pd<jupyter_output><empty_output><jupyter_text>Functions ___<jupyter_code>def TrackWithMasks(img_max, particle_diameter ,min_m, mask_include, mask_exclude):
    f = tp.batch(img_max, diameter=particle_diameter,minmass=min_m)
    f_list = []
    for i in np.arange(len(f['frame'].unique())):
        f0 = f[f['frame']==i]
        f1 = f0.copy()

        # If no masks, include everything and exclude nothing
        if len(mask_include) == 0 :
            mask_include1 = [[0,0],[0,10000000],[10000000,10000000],[10000000,0]]
        else:
            mask_include1 = mask_include
        if len(mask_exclude) == 0:
            mask_exclude1 = [[0,0],[0,-10],[-10,-10]]
        else:
            mask_exclude1 = mask_exclude

        mask_in = Polygon(mask_include1)  # This is a polygon that defines the mask
        mypts = np.transpose([f1.y,f1.x]) # These are the points detected by trackpy, notice the x/y inversion for napari
        f1['Include']=[mask_in.contains(Point(mypts[i])) for i in np.arange(len(mypts))] # Check if pts are on/in polygon mask  

        # # Label points in nucleus if polygon mask exists
        # if polygon != None:
        mask_out = Polygon(mask_exclude1)  # This is a polygon that defines the mask
        f1['Exclude']=[mask_out.contains(Point(mypts[i])) for i in np.arange(len(mypts))] # Check if pts are on/in polygon mask  
        
        f_list.append(f1[(f1['Include']==True) & (f1['Exclude']==False)])
    f_all = pd.concat(f_list)
    return f_all<jupyter_output><empty_output><jupyter_text>Loading data (and track file, if already available)----<jupyter_code>#Video directory
img_4D =  imread(os.path.join(dir,img_4D_filename))
img_4D.shape
# Converting the video to Croparray format
img_croparray = np.expand_dims(img_4D,axis=0) # 
img_croparray.shape # dimensions MUST be (fov, f , z, y, x, ch)
img_croparray.shape

print("croparray format shape [fov, f , z, y, x, ch] = ", img_croparray.shape)<jupyter_output>croparray format shape [fov, f , z, y, x, ch] =  (1, 65, 13, 512, 512, 3)<jupyter_text>Optional tracking: Max projection, masking, and tracking---- Just view video to determine what are the best z planes<jupyter_code>viewer1 = napari.view_image(img_croparray[0,:,:,:,:,1])<jupyter_output><empty_output><jupyter_text>Now do max-projection on best-z planes<jupyter_code>best_zs = [1,10]
img_max = np.max(img_croparray[0,:,best_zs[0]:best_zs[1],:,:,1],axis=1)
img_max.shape<jupyter_output><empty_output><jupyter_text>Now you can create a mask to exclude regions: Shapes will be excluded and Shapes [1] will be included in tracks<jupyter_code>viewer2 = napari.view_image(np.max(img_max,axis=0))  # Max of max, should be a single frame
# Use this if you don't have a mask:
# f_all = TrackWithMasks(img_max,7,1000,[],[])
particle_diameter = 7
min_mass = 1000
viewer3 = napari.view_image(img_max)
f_all = TrackWithMasks(img_max,particle_diameter,min_mass,viewer.layers['Shapes [1]'].data[0][:,-2:],viewer.layers['Shapes'].data[0][:,-2:])
data = f_all[['frame','y','x']].values
layer_name = 'Spots all frames'
viewer3.add_points(data, size = 7, edge_color = 'yellow', symbol='ring', name = layer_name)
f_all
f_all.to_csv(os.path.join(dir,img_4D_filename[:-4]+'.csv'))<jupyter_output><empty_output><jupyter_text>Convert f to crop_array format<jupyter_code>#only if you didn't track:
spots = f_all.copy() # Nice to copy; seems it can cause to overwrite otherwise
spots['id']=spots.index
spots.rename(columns={'x': 'xc','y': 'yc', 'frame': 'f','signal':'signal_tp'}, 
inplace=True, errors='raise')
spots['fov']=0
spots.rename(columns={'particle':'id'})
spots = spots[['fov','id','f','yc','xc','signal_tp','Include','Exclude']] # keeping signal out of curiousity... want to compare to disk-donut measurements  
spots.head()<jupyter_output><empty_output><jupyter_text>Create Crop Array____ Create a crop array from 4D movie<jupyter_code>my_ca = ca.create_crop_array(img_croparray,spots,xy_pad=5, dxy=130, dz=500, dt=1, units=['nm','min'], name = os.path.join(dir,img_4D_filename))
my_ca<jupyter_output>Original video dimensions:  (1, 65, 13, 512, 512, 3)
Padded video dimensions:  (1, 65, 13, 524, 524, 3)
Max # of spots per frame:  36
Shape of numpy array to hold all crop intensity data:  (1, 36, 65, 13, 11, 11, 3)
Shape of xc and yc numpy arrays:  (1, 36, 65, 3)
Shape of extra my_layers numpy array:  (4, 1, 36, 65)<jupyter_text>Save the crop array____<jupyter_code>my_ca.to_netcdf(os.path.join(dir,img_4D_filename[:-4]+'.nc') )<jupyter_output><empty_output><jupyter_text>Quantify signal intensity through time____<jupyter_code># Measure signals and plot average signal through time, creating 'best_z' layer and 'signal' layer
ca.measure_signal(my_ca, ref_ch=1, disk_r=3, roll_n=3)
my_ca.best_z.mean('n').sel(fov=0,ch=1).rolling(t=3,min_periods=1).mean().plot.imshow(col='t',col_wrap=10,robust=True,xticks=[],yticks=[],size=1.5,cmap='gray', vmin=0, vmax =500)<jupyter_output><empty_output><jupyter_text>Let's compare our disk-donut 'signal' layer (acquired from 3D image) to trackpy's (acquired from max-projection):<jupyter_code># Let's compare our intensity numbers to those from trackpy:
my_ca.where(my_ca.signal>0).plot.scatter(x='signal',y='signal_tp',col='ch',hue='ch',colors=['red','limegreen','blue'],levels=[0,1,2,3])<jupyter_output><empty_output><jupyter_text>Let's look at average signal vs time<jupyter_code># Let's look at average signal vs time
start_sig = my_ca.signal.mean('n').sel(t=slice(0,4)).mean('t')
end_sig = 0# my_ca.signal.mean('n').sel(t=slice(15,20)).mean('t')
norm_sig = (my_ca.signal.mean('n') - end_sig)/(start_sig - end_sig)
sns.set_palette(['limegreen','limegreen','blue'])
norm_sig.sel(fov=0,ch=1).plot.line(x='t',hue='ch')<jupyter_output><empty_output><jupyter_text>Now let's just use trackpy's values:<jupyter_code># Let's look at average signal vs time
start_sig = my_ca.signal_tp.mean('n').sel(t=slice(0,4)).mean('t')
end_sig = 0# my_ca.signal_tp.mean('n').sel(t=slice(15,20)).mean('t')
norm_sig = (my_ca.signal_tp.mean('n') - end_sig)/(start_sig - end_sig)
sns.set_palette(['limegreen','limegreen','blue'])
norm_sig.sel(fov=0).plot.line(x='t',hue='ch')<jupyter_output><empty_output><jupyter_text>I guess trackpy and the disk donut method do a very good job at getting the intensities of spots. Although note that trackpy got the values from the max-intensity projection. Interesting.  Visualize crop array montage with Napari___ Now let's see a montage of the selected spots' best-z planes:<jupyter_code># view the action of montage showing an n x n crop_array through time
viewer = napari.view_image(ca.montage(my_ca.sel(fov=0,ch=0).best_z,row='n',col='t'),contrast_limits=[60,800])<jupyter_output><empty_output><jupyter_text>Optional: Create Track Array___<jupyter_code>f_all
# Note, if you actually wanted to track, you could use the following:
# link tracks
max_distance_movement = 15
track_skip_frames = 3
min_trajectory_length = 10
t = tp.link(f_all, max_distance_movement, memory=track_skip_frames)
t1 = tp.filter_stubs(t, min_trajectory_length)
t1['particle'] = t1['particle']+1 # VERY IMPORTANT NOT TO HAVE TRACK IDs WITH VALUES = 0 WHEN MAKING CROP ARRAYS AS ZERO IS DEFAULT EMPTY VALUE
# Compare the number of particles in the unfiltered and filtered data.
print('Before:', t['particle'].nunique())
print('After:', t1['particle'].nunique()) 
# only if you tracked:
spots = t1.copy()
spots.rename(columns={'x': 'xc','y': 'yc', 'frame': 'f','signal':'signal_tp','particle':'id'}, 
inplace=True, errors='raise')


spots['fov']=0
spots.rename(columns={'particle':'id'})
spots = spots[['fov','id','f','yc','xc','signal_tp']] # keeping signal out of curiousity... want to compare to disk-donut measurements  
spots.head()
my_ca2 = ca.create_crop_array(img_croparray,spots,xy_pad=particle_diameter//2+1, dxy=130, dz=500, dt=1, units=['nm','min'])
# Measure signals and plot average signal through time, creating 'best_z' layer and 'signal' layer
ca.measure_signal(my_ca2, ref_ch=1, disk_r=3, roll_n=3)
import xarray as xr
import pandas as pd
# Since ids correspond to tracks, we can organize tracks in rows  
my_ids = np.unique(my_ca2.id) # Find all unique ids
my_ids = my_ids[1:] # remove the '0' ID used as filler in Crop Arrays
my_ids
# Get a list of xarrays for each unique id 
my_das = []
for i in np.arange(len(my_ids)):
    temp = my_ca2.groupby('id')[my_ids[i]].reset_index('stacked_fov_n_t').reset_coords('n',drop=True).reset_coords('fov',drop=True).swap_dims({'stacked_fov_n_t':'t'})
    my_das.append(temp)
    del temp
# Concatenate the xarrays together to make a new xarray dataset in a track array format (each track on separate row). Here 'n' is replaced by 'tracks'
my_taz = xr.concat(my_das, dim=pd.Index(my_ids, name='track_id'), fill_value=255) # fill_value=0 so keep int instead of moving to floats with NaNs
my_taz = my_taz.transpose('track_id','fov','n','t','z','y','x','ch', missing_dims='ignore') # reorder for napari
my_taz
# view the action of montage showing an n x n crop_array through time
viewer = napari.view_image(ca.montage(my_taz.sel(ch=1).best_z,row='track_id',col='t'),contrast_limits=[60,800])
my_taz.isel(track_id=[12,16,20]).sel(ch=1).signal.plot.line(x='t',col='track_id',col_wrap=3)<jupyter_output><empty_output><jupyter_text>Filenames----<jupyter_code># Data filename and directory
dir = r'X:\Tim'
dir = r'Z:\galindo\1_Imaging_Data\20220210_metabolites\PEP_10mM'
#img_4D_max_filename = r'MAX_Chamber02_HT_Cell01.tif'
img_4D_filename = r'Chamber02_HT_Cell02.tif'
img_4D_filename = r'Cell02.tif'
img_4D_filename = r'Cell04.tif'<jupyter_output><empty_output><jupyter_text>Trying to track in x, y, and z: Measure signals and plot average signal through time, creating 'best_z' layer and 'signal' layer<jupyter_code>max_distance_movement = 5
track_skip_frames = 3
my_frame = 15
min_trajectory_length=2
my_list = []
for my_frame in np.arange(len(img_croparray[0,:])): #np.arange(15,16,1):
    f = tp.batch(img_croparray[0,my_frame,:,:,:,1], diameter=7,minmass=1000)
    t1 = tp.link(f, max_distance_movement, memory=track_skip_frames)
    t = tp.filter_stubs(t1, min_trajectory_length)
    sort_t = t.sort_values(['particle']).groupby('particle').aggregate('max').reset_index().rename(columns={'frame':'z'})
    sort_t['t']=my_frame
    my_list.append(sort_t)
my_df = pd.concat(my_list)
my_df.reset_index()
import pandas as pd
my_df = pd.concat(my_list)
my_df.reset_index()
my_df['raw_mass'].hist()
plt.figure(figsize=(20,10))
tp.annotate(my_df[(my_df['t']==28)&(my_df['raw_mass']<40000)], img_4D_max_real[28]);
np.array([len(t[t['particle']==i].x) for i in np.arange(len(t.particle))])
f[f['frame']==6]<jupyter_output><empty_output><jupyter_text>GUI ___<jupyter_code>from napari.layers import Image, Shapes
from napari.types import LabelsData, ImageData, ShapesData
import napari.types
import trackpy as tp
#viewer.dims.current_step[0] This is how you can access the slider slice in napari

# @magicgui(call_button = 'Load File')
# def open_file(
#     number_of_channels: int,
#     filename = pathlib.Path('/some/path.ext')
# ) -> napari.types.LayerDataTuple:
#     img = imread(filename)
#     my_channel_axis = np.where(np.array(img.shape) == number_of_channels)[0][0]
#     img = np.moveaxis(img_4D,my_channel_axis,0)
#     return [(img[i],{'colormap':'gray','name':'Channel'+str(i)}) for i in np.arange(len(img))]
@magicgui(call_button = 'Make Tracking Channel Layer')
def get_channel(
    image: Image,
    channel_axis: int,
    channel_num: int,
) -> napari.types.LayerDataTuple:
    return (np.expand_dims(np.moveaxis(image.data,channel_axis,0)[channel_num],axis=channel_axis),{'name':image.name+'Ch '+str(channel_num)})

@magicgui(call_button = 'Make Max Projection')
def max_project(
    image1: Image,
    projection_axis: int,
    start_slice: int,
    stop_slice: int,
) -> napari.types.LayerDataTuple:
    return (np.expand_dims(np.max(image1.data, axis=projection_axis),axis=projection_axis),{'name':'Max of '+image1.name})

@magicgui(call_button = 'Detect Spots in Max Projection')
def detect_spots(
    image: Image,
    mask_include: Shapes,
    mask_exclude: Shapes,
    test_frame: int,
    size = 9,
    min_m = 1000,
    detect_spots_in_all_frames = False,
) -> napari.types.LayerDataTuple:
    if detect_spots_in_all_frames:
        f_all = TrackWithMasks(np.squeeze(image.data),size,min_m,mask_include.data[0][:,-2:],mask_exclude.data[0][:,-2:])
        print(f_all)
        print(np.squeeze(image.data).shape)
        data = f_all[['frame','y','x']].values
        layer_name = 'Spots all frames'
    else:
        f = tp.locate(np.squeeze(image.data[test_frame]),diameter=size,minmass=min_m)
        print(f)
        data = f[['y','x']].values
        layer_name = 'Spots frame '+str(test_frame)
    #save_croparray.dir.value = 'test'
    #save_croparray.show()
    return [(data,{'size':size,'edge_color':'yellow','opacity':0.4,'symbol':'ring','name':layer_name}, 'points'),
            (np.max(image,axis=0), {'name': 'Max Projection', 'blending': 'additive'})]

viewer = napari.Viewer()
#viewer = napari.view_image(img_max, name="My Image")
#viewer.window.add_dock_widget(detect_spots)
#viewer.window.add_dock_widget(max_project)
#viewer.window.add_dock_widget(open_file, name='Step 1',area='right')
viewer.window.add_dock_widget(get_channel, name='Get Video',area='right')
dw2 = viewer.window.add_dock_widget(max_project, name='Max Project',area='right')
dw3 = viewer.window.add_dock_widget(detect_spots, name ='Tracking',area='right')
viewer.window._qt_window.tabifyDockWidget(dw2, dw3)
#viewer.layers.events.changed.connect(detect_spots.reset_choices)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14829849130640688, 'fraction_numerical': 0.02449090048034639, 'mean_word_length': 5.943165805542508, 'pattern_counts': {'":': 0, '<': 72, '<?xml version=': 0, '>': 76, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '309496', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5738, 'n_tokens_neox': 5430, 'n_words': 1176, 'path': 'notebooks/ToolDevelopment_Harringtonine_Morisaki.ipynb', 'repo_name': 'Colorado-State-University-Stasevich-Lab/croparray', 'size': '14781'}