Update AnalyzeSequenceV2.ipynb (f0581cc6) · Commits · Nicolas Elie / REDPOL-open

AnalyzeSequenceV2.ipynb

+1 −1

Changes for AnalyzeSequenceV2.ipynb: 1 added line, 1 removed line.
Original line number	Diff line number	Diff line
		%% Cell type:markdown id: tags:

		![REDPOL_LOGO](https://git.unicaen.fr/nicolas.elie/redpol-open/-/raw/master/media/REDPOL_LOGO.jpg)
		## _Script to analyze our raw image sequences were taken with a colour camera mounted on a stereomicroscope._

		Anaïd Gouveneaux (1,2),
		Mamoudou Sano (3),
		Nicolas Elie (3),
		Apolline Chabenat (1,2),
		Céline Thomasse (2),
		Christelle Jozet-Alves (2),
		Anne-Sophie Darmaillacq (2),
		Ludovic Dickel(2),
		Thomas Knigge (1),
		Cécile Bellanger (2),

		- 1 Normandie Univ, UNILEHAVRE, UMR-I02, Environmental stress and biomonitoring of aquatic environments (SEBIO), 76600 Le Havre, France
		- 2 Normandie Univ, UNICAEN, UMR 6552, Rennes 1-CNRS (EthoS) - Cephalopod cognitive neuroethology (NECC), 14000 Caen, France
		- 3 Normandie Univ, UNICAEN, CMAbio3 - Centre de Microscopie Appliquée à la biologie

		> The file format is VSI from Olympus

		%% Cell type:markdown id: tags:

		## ----------------------------------------------------------------------------------

		%% Cell type:markdown id: tags:

		## Import the required modules

		%% Cell type:code id: tags:

		``` python
		import numpy as np
		import os
		from csbdeep.utils import normalize
		from stardist.models import StarDist2D
		import cv2
		import javabridge
		import bioformats
		import pandas as pd
		from shapely.geometry import Polygon
		from shutil import make_archive,rmtree
		from zipfile import ZipFile
		axis_norm = (0,1) # normalize channels independently
		from PIL import Image
		import math
		from skimage import io,filters,morphology,measure
		```

		%% Cell type:markdown id: tags:

		## Load Model

		%% Cell type:code id: tags:

		``` python
		model = StarDist2D(None, name='stardist_multiclass', basedir='models')
		```

		%% Output

		Loading network weights from 'weights_best.h5'.
		Loading thresholds from 'thresholds.json'.
		Using default values: prob_thresh=0.65, nms_thresh=0.75.

		%% Cell type:markdown id: tags:

		## Call java and using bioformat to read vsi files.

		%% Cell type:code id: tags:

		``` python
		def _init_logger():
		"""This is so that Javabridge doesn't spill out a lot of DEBUG messages
		during runtime.
		From CellProfiler/python-bioformats.
		"""
		rootLoggerName = javabridge.get_static_field("org/slf4j/Logger",
		"ROOT_LOGGER_NAME",
		"Ljava/lang/String;")

		rootLogger = javabridge.static_call("org/slf4j/LoggerFactory",
		"getLogger",
		"(Ljava/lang/String;)Lorg/slf4j/Logger;",
		rootLoggerName)

		logLevel = javabridge.get_static_field("ch/qos/logback/classic/Level",
		"WARN",
		"Lch/qos/logback/classic/Level;")

		javabridge.call(rootLogger,
		"setLevel",
		"(Lch/qos/logback/classic/Level;)V",
		logLevel)
		```

		%% Cell type:markdown id: tags:

		## Read metadata of vsi files

		%% Cell type:code id: tags:

		``` python
		def get_metadata(filename):
		"""Read the meta data and return the metadata object.
		"""
		meta = bioformats.get_omexml_metadata(filename)
		metadata = bioformats.omexml.OMEXML(meta)
		return metadata
		```

		%% Cell type:markdown id: tags:

		## Javabridge: running and interacting with the JVM from Python

		%% Cell type:code id: tags:

		``` python
		javabridge.start_vm(class_path=bioformats.JARS,run_headless=False)
		logger = _init_logger()
		```

		%% Cell type:code id: tags:

		``` python
		def Ruifrok(ImageColor,color):
		"""
		stack : Colour image RGB
		stain : "H E" "H DAB" "H E DAB" "H AEC", choose among these staining
		If you indicate "US", for User Values, fill : R1,G1,B1,R2,G2,B2,R3,G3,B3

		N.Elie at Centre de Morphologie Mathematique, France (in python)

		// G.Landini at bham ac uk
		// This plugin implements stain separation using the colour deconvolution
		// method described in:
		//
		// Ruifrok AC, Johnston DA. Quantification of histochemical
		// staining by color deconvolution. Analytical & Quantitative
		// Cytology & Histology 2001; 23: 291-299.
		//
		// The code is based on "Color separation-30", a macro for NIH Image kindly provided
		// by A.C. Ruifrok. Thanks Arnout! AND Colour_Deconvolution in IMAGEJ
		//
		// The plugin assumes images generated by color subtraction (i.e. light-absorbing dyes
		// such as those used in bright field histology or ink on printed paper) but the dyes
		// should not be neutral grey.
		//
		// I strongly suggest to read the paper reference above to understand how to determine
		// new vectors and how the whole procedure works.
		//
		// The plugin works correctly when the background is neutral (white to light grey),
		// so background subtraction and colour correction must be applied to the images before
		// processing.
		//
		// The plugin provides a number of "built in" stain vectors some of which were determined
		// experimentally in our lab (marked GL), but you may have to determine your own vectors to
		// provide a more accurate stain separation, depending on the stains and methods you use.
		// Ideally, vector determination should be done on slides stained with only one colour
		// at a time (using the "From ROI" interactive option).
		//
		// The plugin takes an RGB image and returns three 8-bit images. If the specimen is
		// stained with a 2 colour scheme (such as H & E) the 3rd image represents the
		// complimentary of the first two colours (i.e. green).
		//
		// Please be very careful about how to interpret the results of colour deconvolution
		// when analysing histological images.
		// Most staining methods are not stochiometric and so optical density of the chromogen
		// may not correlate well with the quantity of the reactants.
		// This means that optical density of the colour may not be a good indicator of
		// the amount of material stained.
		//
		// Read the paper!

		"""

		if color=="H DAB":Kernel=np.mat('-0.042843347519395891, 2.0723332642033467, -1.5074071436349348; -1.4160021820814035, 0.83451202938550129, 1.1640033205088089; 1.9795963018710074, -2.2092659159211636, 0.93911051955365066')

		imageR, imageV, imageB = ImageColor.transpose( 2,0,1 )


		imageR=np.float64(imageR)
		imageV=np.float64(imageV)
		imageB=np.float64(imageB)

		log255=math.log(255.0)
		imageRLog=-((255.0*np.log((imageR+1)/255.0))/log255)
		imageVLog=-((255.0*np.log((imageV+1)/255.0))/log255)
		imageBLog=-((255.0*np.log((imageB+1)/255.0))/log255)

		# Faire une boucle de 3 passages, chacun calcul le nouveau canal
		for i in range(3):
		Rscaled=imageRLog*Kernel[i,0]
		Vscaled=imageVLog*Kernel[i,1]
		Bscaled=imageBLog*Kernel[i,2]
		output = np.exp(-((Rscaled + Vscaled + Bscaled) - 255.0) * log255 / 255.0)
		tab=np.where(output>255)
		output[tab]=255
		if i==0:NstackR=np.floor(output+0.5)
		if i==1:NstackG=np.floor(output+0.5)
		if i==2:NstackB=np.floor(output+0.5)

		NstackR=np.uint8(NstackR)
		NstackG=np.uint8(NstackG)
		NstackB=np.uint8(NstackB)
		return NstackR,NstackG,NstackB
		```

		%% Cell type:markdown id: tags:

		## Main function: application of the Stardist model on each image of the vsi file, detection and classification of chromatophores, saving of results and creation of optional additional result files.

		%% Cell type:code id: tags:

		``` python
		def AI_OnePile(number,pathFile, output_path, start=None,End=None, OptionS="9",):
		"""
		number: number corresponding to vsi file
		pathFile: path of vsi file
		output_path : folder to save data
		start : Index of the first image in stack
		End : Index of the last image in stack

		OptionS : 9 , just file result
		0, add video file
		1, add images files : BINARY
		2, add images files : LABEL
		3, add files pickle to keep indivudial dataframe.

		Example : "0-2" --> add video file and Label images.
		"""

		# Create result file
		fileResults=open(output_path + number+'_'+start+'-'+str(End)+"_resultat.csv","w")
		fileResults.write("Sequence;Image area;Number of Light chromatophores;Number of Dark chromatophores;Total area of Light chromatophores;Total area of Dark chromatophores;Average size of Light chromatophores;Average size of Dark chromatophores;Median size of Light chromatophores; Median size of Dark chromatophores;Standard deviation of the size of Light chromatophores;Standard deviation of the size of Dark chromatophores\n")

		# path of image
		fileEXT=pathFile+"/Process_"+number+".vsi"

		# Bioformat reading
		ImageVSI=bioformats.ImageReader(path=fileEXT)

		# Get metadata of vsi file
		metadata = get_metadata(fileEXT)

		# Get number of images of vsi file
		if End is None:
		End=metadata.image(0).Pixels.SizeT
		if start is None:
		start=0

		if "0" in OptionS:
		# Criteria to create video capture *****************************************************************************************
		#FIXME: Dimension of vsi in our example
		unnoyau=[0,0,1544,1038]
		height=unnoyau[3]
		width=unnoyau[2]
		# Define the codec and create VideoWriter object
		frame_size=(width,height)
		fourcc =cv2.VideoWriter_fourcc('M','J','P','G')
		# path of video output
		out2 = cv2.VideoWriter(output_path+number+'_'+start+'-'+str(End)+'.avi', fourcc, 5,frame_size)
		# ***************************************************************************************************************************

		Couleur=[]
		Labels=[]

		# OPTION: To initialise sequence to save images
		if ("1" in OptionS) or ("2" in OptionS):
		multilist1=[]
		multilist2=[]

		if ("2" in OptionS):
		multilist3=[]
		multilist4=[]

		if ("3" in OptionS):
		multilist3=[]
		multilist4=[]
		try:
		os.mkdir(os.path.join(output_path,str(number)))
		except OSError:
		pass

		# Browse images in vsi
		for stack in range(int(start),int(End)):

		image=ImageVSI.read(c=None,t=stack,series=0,rescale=False)
		X=cv2.cvtColor(image, cv2.COLOR_RGB2BGR)
		X1 = normalize(X, 1,99.8, axis=axis_norm)
		label, details = model.predict_instances(X1)

		"Use it Stardist classifer"
		#df = pd.DataFrame()
		#df['class']=details['class_id']
		#df['coord']=[Polygon(np.squeeze(x)) for x in np.transpose(details['coord'],(0,2,1))]
		#df['area']=[x.area for x in df['coord']]

		# Compute all data for chromatophores
		#stats = df.groupby('class')['area'].agg(['sum','mean', 'std','median','size'])

		"Use it a thresold on Ruifrok image deconvolution"
		ruif0,ruif1,ruif2=Ruifrok(image,"H DAB")
		filt=filters.median(ruif0,morphology.disk(1))

		props = measure.regionprops_table(label, filt, properties=['label', 'mean_intensity','area'])
		PourClass=pd.DataFrame(props)
		# dark <= 128 [1] and light >128 [2]
		Threshold_Value=220
		Threshold_Value=245#220
		PourClass['class']= np.where(PourClass["mean_intensity"] <= Threshold_Value, 1, 2)
		stats = PourClass.groupby('class')['area'].agg(['sum','mean', 'std','median','size'])

		"""
		Class 1 : dark chromatophore
		Class 2 : light chromatophore
		"""
		if ("0" in OptionS) or ("1" in OptionS) or ("2" in OptionS):


		# Criteria to create video capture *****************************************************************************************
		# with stardist classifer : threshold = details['class_id']==1
		#threshold= np.where(PourClass["intensity_mean"] <= Threshold_Value, True, False)
		#arr = np.insert(threshold, 0, False, axis=None)
		arr1=np.full(np.max(label)+1, False, dtype=bool)
		arr2=np.full(np.max(label)+1, False, dtype=bool)
		for row in PourClass.itertuples():
		if row._4==1:
		arr1[row.label]=True
		if row._4==2:
		arr2[row.label]=True

		Class1 = arr1[label]



		# with stardist classifer : threshold = details['class_id']==2
		#threshold= np.where(PourClass["mean_intensity"] > Threshold_Value, True, False)
		#arr = np.insert(threshold, 0, False, axis=None)
		Class2 = arr2[label]


		LabelClass2=label.copy()
		LabelClass2[Class1] = 0
		LabelClass1=label.copy()
		LabelClass1[Class2] = 0

		if ("0" in OptionS):

		Rouge, Vert, Bleu= image.transpose(2, 0, 1)

		Tab1 = np.where(cv2.morphologyEx(np.float32(LabelClass1),cv2.MORPH_GRADIENT,np.ones((3, 3), np.uint8)) >0)
		Rouge[Tab1] = 255
		Vert[Tab1] = 0
		Bleu[Tab1] = 0
		Tab1 = np.where(cv2.morphologyEx(np.float32(LabelClass2),cv2.MORPH_GRADIENT,np.ones((3, 3), np.uint8)) >0)
		Rouge[Tab1] = 0
		Vert[Tab1] = 0
		Bleu[Tab1] = 255
		Icolor = np.ones((height,width,3), 'uint8')
		Icolor[..., 0] = Bleu
		Icolor[..., 1] = Vert
		Icolor[..., 2] = Rouge
		out2.write(Icolor)

		# OPTION: To append binary image, mask of classes
		if ("1" in OptionS):
		multilist1.append(Image.fromarray(Class1))
		multilist2.append(Image.fromarray(Class2))

		# OPTION: To append label classes
		if ("2" in OptionS):
		multilist3.append(Image.fromarray(LabelClass1))
		multilist4.append(Image.fromarray(LabelClass2))


		# ***************************************************************************************************************************

		text=str(stack)+";"+str(1544*1038)+";"+str(stats['size'][2])+";"+str(stats['size'][1])+";"+str(stats['sum'][2])+";"+str(stats['sum'][1])+";"\
		+str(stats['mean'][2])+";"+str(stats['mean'][1])+";"+str(stats['median'][2])+";"+str(stats['median'][1])+";"\
		+str(stats['std'][2])+";"+str(stats['std'][1])\
		+"\n"
		fileResults.write(text)

		if ("3" in OptionS):
		# OPTION: Export data by image 1 image = 1 zip file
		df.to_pickle(os.path.join(output_path,str(number),number+'_'+str(stack)+'.zip')) # to read unpickled_df = pd.read_pickle(file.zip)



		# OPTION: to save image mask
		if ("1" in OptionS) :
		multilist1[0].save(output_path + number+'_'+start+'-'+str(End)+"_dark_binary.tif",compression="tiff_deflate", save_all=True, append_images=multilist1[1:])
		multilist2[0].save(output_path + number+'_'+start+'-'+str(End)+"_light_binary.tif",compression="tiff_deflate", save_all=True, append_images=multilist2[1:])

		if ("2" in OptionS):
		multilist3[0].save(output_path + number+'_'+start+'-'+str(End)+"_dark_label.tif",compression="tiff_deflate", save_all=True, append_images=multilist3[1:])
		multilist4[0].save(output_path + number+'_'+start+'-'+str(End)+"_light_label.tif",compression="tiff_deflate", save_all=True, append_images=multilist4[1:])



		if ("0" in OptionS):
		out2.release()
		fileResults.close()

		if ("3" in OptionS):
		file = "IndividualDataFrame" # zip file name
		make_archive(os.path.join(output_path,file), "zip", os.path.join(output_path,str(number))) # zipping the directory
		rmtree(os.path.join(output_path,str(number)))
		```

		%% Cell type:code id: tags:

		``` python
		# Browse and find all vsi images inside folders and sub-folders
		def listdirectory(path):
		fichier=[]
		for root, dirs, files in os.walk(path):
		for i in files:
		if ".vsi" in i:
		path,ext=os.path.splitext(i)
		if (os.path.isfile(os.path.join(root, path+".vsi"))):
		fichier.append(os.path.join(root, i))
		return sorted(fichier)
		```

		%% Cell type:code id: tags:

		``` python
		# path where images are stored
		pathchemin="Examples"

		# path where files results are recorded
		output_path="Examples/"

		# call function to find vsi files
		listfile=listdirectory(pathchemin)

		for fileVSI in listfile:
		# get number of vsi files
		# all files in our example are this format : number_Process.vsi
		number=fileVSI.split("/")[-1].split(".")[0].replace("Process_","")

		# test if vsi file are not aleardy process
		if not os.path.isfile(output_path+number+"_resultat.csv"):
		pathFile = os.path.dirname(fileVSI)
		try:
		print("Process Image number", number)
		AI_OnePile(number,pathFile, output_path,start="0",End="3",OptionS="9")
		except:
		print("erreur :",pathFile,number)
		pass

		print("End process")

		```

		%% Output

		Process Image number 195

		1/1 [==============================] - 1s 500ms/step
		1/1 [==============================] - 0s 381ms/step
		1/1 [==============================] - 0s 407ms/step
		End process