text: from __future__ import print_function, unicode_literals, absolute_import, division
  from Functions.tool_modules import *
  
  import numpy as np
  from glob import glob
  from natsort import natsorted
  from tifffile import imread,imwrite
  import pandas as pd
  pd.options.mode.chained_assignment = None
  
  from lmfit import Model,Parameter,Parameters
  from sklearn.metrics import r2_score
  from scipy.interpolate import interp1d
  from tensorflow.keras.models import load_model
  import tensorflow as tf
  
  import matplotlib.pyplot as plt
  import matplotlib
  import matplotlib.animation as animation
  matplotlib.rcParams["image.interpolation"] = None
  from matplotlib.ticker import MultipleLocator
  import matplotlib as mpl
  
  import os
  os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' 
  import shutil
  import time
  from tqdm import tqdm
  import sys
  import gc
  import multiprocessing
  
  from csbdeep.utils import Path, normalize
  from csbdeep.io import save_tiff_imagej_compatible
  from stardist import random_label_cmap, _draw_polygons, export_imagej_rois
  from stardist.models import StarDist2D
  
  import btrack
  import napari
  from btrack.constants import BayesianUpdates
  
  from pympler.tracker import SummaryTracker
  #tracker_memory = SummaryTracker()
  
  t0_chrono = time.perf_counter()
  
  # MOVIE PARAMETERS
  PxToUm = float(ConfigSectionMap("../../../config.ini","MovieSettings")["pxtoum"])
  FrameToMin = float(ConfigSectionMap("../../../config.ini","MovieSettings")["frametomin"])
  len_movie = int(ConfigSectionMap("../../../config.ini","MovieSettings")["len_movie"])
  shape_x = int(ConfigSectionMap("../../../config.ini","MovieSettings")["shape_x"])
  shape_y = int(ConfigSectionMap("../../../config.ini","MovieSettings")["shape_y"])
  
  #SEARCH RADIUS
  search_radius_nk = int(ConfigSectionMap("../../../config.ini","SearchRadii")["search_radius_nk"])
  
  
  np.random.seed(6)
  lbl_cmap = random_label_cmap()
  
  bf_channel = 0
  red_channel = 1
  green_channel = 2
  blue_channel = 3
  
  tracklength_threshold = 1
  
  do_segmentation_and_measure = True
  do_classification = True
  do_neighbour_search = True
  	
  files = glob("../aligned/*.tif")
  filename = os.path.split(files[0])[-1][:-8]
  
  stack = np.array(imread(files[0]))
  
  blue = np.array(stack[:,blue_channel,:,:])
  green = np.array(stack[:,green_channel,:,:])
  red = np.array(stack[:,red_channel,:,:])
  
  stack = np.moveaxis(stack,1,-1)
  print(np.shape(stack))
  
  #del stack
  #gc.collect()
  
  if do_segmentation_and_measure:
  
  	gpus = tf.config.experimental.list_physical_devices("GPU")
  	if gpus:
  		try:
  			for gpu in gpus:
  				tf.config.experimental.set_memory_growth(gpu,True)
  		except RuntimeError as e:
  			print(e)
  
  	print("##########################################")
  	print("############# STARDIST 2D ################")
  	print("##########################################")
  	print("Segmentation of the green channel...")
  
  	# Segmentation model
  	model =  StarDist2D(None, name='stardist_nk_feb', basedir='/home/limozin/Desktop/segment_nks/models/') #StarDist2D.from_pretrained('2D_versatile_fluo')
  
  	axis_norm = (0,1)   # normalize channels independently
  
  	X_pred = [normalize(x,1,99.8,axis=axis_norm) for x in tqdm(stack)]
  
  	rows = []
  	patch = 10
  	mask = create_circular_mask(2*patch,2*patch,((2*patch)//2,(2*patch)//2),patch)
  	green = np.array([np.pad(f,[patch+1,patch+1],mode="constant") for f in green])
  	red = np.array([np.pad(f,[patch+1,patch+1],mode="constant") for f in red])
  
  	recreate_folder("../output/nk_intensity/")
  	recreate_folder("../output/tables/")
  	recreate_folder("../labels_nk/")
  
  	nk_index = 0
  	rows = []
  	for i in tqdm(range(len(X_pred))):
  		# Segment CFSE for each frame
  		Y_pred,details = model.predict_instances(X_pred[i], n_tiles=model._guess_n_tiles(X_pred[i]), show_tile_progress=False)
  		save_tiff_imagej_compatible("../labels_nk/"+f"{str(i).zfill(4)}.tif", Y_pred, axes='YX')
  		coords = details["coord"]
  		igreens = []; ireds = [];
  		for co in coords:
  			# Measure PI & CFSE signal around each NK in disk
  			x = np.mean(co[1,:])
  			y = np.mean(co[0,:])
  			xmin = int(x); xmax = int(x)+2*int(patch);
  			ymin = int(y); ymax = int(y)+2*int(patch);
  			local_green = np.multiply(green[i,ymin:ymax,xmin:xmax],mask)
  			green_measurement = np.mean(local_green[local_green!=0.])
  			local_red = np.multiply(red[i,ymin:ymax,xmin:xmax],mask)
  			red_measurement = np.mean(local_red[local_red!=0.])
  			
  			igreens.append(green_measurement)
  			ireds.append(red_measurement)
  			row_labels = ["NK_ID","X","Y","T","IG","IR"]
  			row = [nk_index,x,y,float(i),green_measurement,red_measurement]
  			rows.append(row)
  			
  			nk_index+=1
  		
  		plt.scatter(igreens,ireds,alpha=0.1)
  		plt.ylim(np.amin(red),np.amax(red))
  		plt.xlim(np.amin(green),np.amax(green))
  		plt.xlabel("intensity green [a.u.]")
  		plt.ylabel("intensity red [a.u.]")
  		plt.title(f"Frame {i}")
  		plt.savefig(f"../output/nk_intensity/nk_intensity_t{i}.png")
  		plt.close()
  
  	df_nk = pd.DataFrame(rows,columns=row_labels)
  	df_nk.to_csv("../output/tables/table_nks.csv")
  			
  
  	del X_pred
  	del Y_pred
  	del details
  	del model
  	gc.collect()
  
  if do_classification:
  
  	print("##########################################")
  	print("############# CLASSIFICATION #############")
  	print("##########################################")
  	print("Preparing for the classification of the NKs as dead or alive based on their PI & CFSE intensity values...")	
  
  	df = pd.read_csv("../output/tables/table_nks.csv")
  	df.sort_values(by="T")
  
  	times = np.unique(df["T"].to_numpy())
  
  	t0 = times[0]; tf = times[-1];
  	df0 = df[(df["T"]==t0)]
  	dff = df[(df["T"]==tf)]
  
  	Ig0 = df0.IG.to_numpy()
  	Ir0 = df0.IR.to_numpy()
  
  	Igf = dff.IG.to_numpy()
  	Irf = dff.IR.to_numpy()
  
  
  	def onclick(event):
  	    global ix1, iy1, ix2, iy2, center_coords
  	    ix, iy = event.xdata, event.ydata
  	    print('x = %d, y = %d'%(ix, iy))
  	    center_coords.append([ix, iy])
  
  	    return(center_coords)
  
  	validate=False
  	while validate==False:
  		center_coords = []
  
  		fig = plt.figure()
  		ax = fig.add_subplot(111)
  		ax.scatter(Ig0,Ir0,alpha=0.1,color="tab:blue",label="t0")
  		ax.scatter(Igf,Irf,alpha=0.1,color="tab:red",label='tf')
  		ax.set_xlabel("intensity green [a.u.]")
  		ax.set_ylabel("intensity red [a.u.]")
  		ax.legend()
  		cid = fig.canvas.mpl_connect('button_press_event', onclick)
  		plt.show()
  
  		center_coords = np.array(center_coords)
  		xpos = center_coords[:,0]
  		ypos = center_coords[:,1]
  		a = (ypos[1] - ypos[0]) / (xpos[1] - xpos[0])
  		b = ypos[0] - a*xpos[0]
  		print(f"y = {a} * x + {b}")
  
  		xplot = np.linspace(np.amin(df.IG.to_numpy()),np.amax(df.IG.to_numpy()),1000)
  		plt.plot(xplot,a*xplot+b,"--")
  		plt.scatter(Ig0,Ir0,alpha=0.1,color="tab:blue",label="t0")
  		plt.scatter(Igf,Irf,alpha=0.1,color="tab:red",label='tf')
  		plt.xlabel("intensity green [a.u.]")
  		plt.ylabel("intensity red [a.u.]")
  		plt.legend()
  		plt.pause(2)
  		plt.close()
  
  		val = str(input("Validate? [y/n]"))
  		if val=='y':
  			validate=True
  
  	test = np.array([ir <= (a*ig + b) for ir,ig in zip(df.IR.to_numpy(),df.IG.to_numpy())],dtype=np.int8)
  
  	df["STATUS"] = test
  
  	df.to_csv("../output/tables/table_nks_w_death.csv")
  
  	plt.scatter(df.IG.to_numpy(),df.IR.to_numpy(),c=df.STATUS.to_numpy(),alpha=0.01)
  	plt.xlabel("intensity green [a.u.]")
  	plt.ylabel("intensity red [a.u.]")
  	plt.savefig("../output/nk_intensity/final.png",dpi=300)
  	plt.show()
  	
  if do_neighbour_search:
  
  	print("##########################################")
  	print("########### NEIGHBOUR SEARCH #############")
  	print("##########################################")
  	print("Attributing each TC to its NK cell neighbours...")
  
  	df_tc = pd.read_csv("../visual_table_checked.csv")
  	df_nk = pd.read_csv("../output/tables/table_nks_w_death.csv")
  
  	option_movie = False
  
  	unique_tc = np.unique(df_tc.TID.to_numpy())
  	unique_times = np.unique(df_tc["T"].to_numpy())
  
  
  	mask = 	create_circular_mask(2*search_radius_nk,2*search_radius_nk,((2*search_radius_nk)//2,(2*search_radius_nk)//2),search_radius_nk)
  
  	rows = []
  	for tid in tqdm(unique_tc):
  
  		df_tc_at_track = df_tc[(df_tc["TID"]==tid)]
  		indices = df_tc_at_track.index
  
  		x_array = df_tc_at_track.X.to_numpy()/PxToUm
  		y_array = df_tc_at_track.Y.to_numpy()/PxToUm
  		cell_class = df_tc_at_track.CLASS.to_numpy()
  		death_time = df_tc_at_track.X0.to_numpy()
  
  		border_cdt = np.all([np.all(x_array>search_radius_nk),np.all(y_array>search_radius_nk),np.all(x_array<(blue[0].shape[1] - search_radius_nk)),np.all(y_array<(blue[0].shape[0] - search_radius_nk))])
  		if border_cdt:
  			for t,time in enumerate(unique_times):
  				df_nk_at_t = df_nk[(df_nk["T"]==time)]
  				x_nks = df_nk_at_t.X.to_numpy()
  				y_nks = df_nk_at_t.Y.to_numpy()
  				nk_indices = df_nk_at_t.NK_ID.to_numpy()
  				nk_status = df_nk_at_t.STATUS.to_numpy()
  				radius_test = ((x_nks - x_array[t])**2 + (y_nks - y_array[t])**2 < search_radius_nk**2)
  				neighbours = np.array(nk_indices[np.where(radius_test)[0]])
  				neighbours_status = np.array(nk_status[np.where(radius_test)[0]])
  				neighbours_alive = neighbours[(neighbours_status==1)]
  				neighbours_dead = neighbours[(neighbours_status==0)]
  
  
  				xmin = int(x_array[t]) - search_radius_nk; xmax = int(x_array[t]) + search_radius_nk;
  				ymin = int(y_array[t]) - search_radius_nk; ymax = int(y_array[t]) + search_radius_nk;
  
  
  				if option_movie:
  					recreate_folder("../output/movies_cell_centered/")
  					local_blue = np.multiply(blue[t,ymin:ymax,xmin:xmax],mask)
  					plt.imshow(local_blue,cmap="gray")
  					plt.scatter(x_nks[np.where(radius_test)[0]] - (x_array[t]-search_radius_nk),y_nks[np.where(radius_test)[0]] - (y_array[t]-search_radius_nk),c="green",marker="x")
  					plt.xticks([])
  					plt.yticks([])
  					plt.savefig(f"../output/movies_cell_centered/TC_{tid}_time_{t}.png",dpi=300,bbox_inches="tight")
  					plt.close()
  
  				rel_nk_x = x_nks[np.where(radius_test)[0]] - (x_array[t])
  				rel_nk_y = y_nks[np.where(radius_test)[0]] - (y_array[t])
  				relative_dist = [[a,b] for a,b in zip(rel_nk_x,rel_nk_y)]
  				relative_radii = [np.sqrt(a**2 + b**2) for a,b in zip(rel_nk_x,rel_nk_y)]
  
  				row_labels = ["TC_INDEX","TID","X","Y","T","CLASS","T0","NK_NEIGHS","NK_NEIGH_STATUS","NK_DISTANCES","NK_RADII","NK_NBR_NEIGH","NK_NBR_NEIGH_ALIVE","NK_NBR_NEIGH_DEAD"]
  				row = [indices[t],tid,x_array[t],y_array[t],time,cell_class[t],death_time[t],neighbours,neighbours_status,relative_dist,relative_radii,len(neighbours),len(neighbours_alive),len(neighbours_dead)]
  				rows.append(row)
  
  	df = pd.DataFrame(rows,columns=row_labels)
  	df.sort_values(by=["TID","T"])
  	df.to_csv("../output/tables/tc_w_nk_neighbours.csv")
  
  
  	
  