suppress cache data

%% 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 numpy

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

```

%% Output

    __init__.py (155): A NumPy version >=1.18.5 and <1.25.0 is required for this version of SciPy (detected version 1.16.6

%% 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: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,numberstack=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

    numberstack : 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(numberstack)+"_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 numberstack is None:

        numberstack=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(numberstack)+'.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(numberstack)):

        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)

        df = pd.DataFrame()

        df['class']=details['class_id']

        df['coord']=[Polygon(numpy.squeeze(x)) for x in numpy.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'])

        """

        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 *****************************************************************************************

            threshold = details['class_id']==1

            #nbClasse1=numpy.sum(threshold)

            arr = numpy.insert(threshold, 0, False, axis=None)

            Class1 = arr[label]

            threshold = details['class_id']==2

            #nbClasse2=numpy.sum(threshold)

            arr = numpy.insert(threshold, 0, False, axis=None)

            Class2 = arr[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 = numpy.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(numberstack)+"_dark_binary.tif",compression="tiff_deflate", save_all=True, append_images=multilist1[1:])

        multilist2[0].save(output_path + number+'_'+start+'-'+str(numberstack)+"_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(numberstack)+"_dark_label.tif",compression="tiff_deflate", save_all=True, append_images=multilist1[1:])

        multilist4[0].save(output_path + number+'_'+start+'-'+str(numberstack)+"_light_label.tif",compression="tiff_deflate", save_all=True, append_images=multilist2[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,"0","3")

        except:

            print("erreur :",pathFile,number)

            pass

print("End process")

```

%% Output

    ['Examples/Process_13.vsi']

    Process Image number 13

    End process

%% Cell type:code id: tags:

``` python

```

%% Cell type:code id: tags:

``` python

from __future__ import print_function

from collections import OrderedDict

from io import BytesIO

import random

import struct

import gzip

import geojson

from shapely.geometry import shape

from shapely.strtree import STRtree

from shapely.geometry import Point

from shapely.geometry import Polygon

from skimage.draw import polygon

import numpy

import multiprocessing as mp

import subprocess, os, skimage,mahotas, cv2

from skimage import io, filters, morphology, util, transform, segmentation, measure,io,feature,exposure

import math

#from skimage.segmentation import clear_border

import pandas as pd

import glob, scipy

import itertools

from skimage.filters.rank import entropy

import javabridge

import bioformats

from xml import etree as et

```

%% Cell type:code id: tags:

``` python

def drawnucleiFinal(SizeImage,allshape):

    """

    SizeImage : size image

    allshape : coordinates of objects

    """

    TotalUnique=numpy.zeros((SizeImage[3],SizeImage[2]),numpy.uint16())

    compteurU=1

    Dictionnaire={}

    for objet in allshape:

        ROI=numpy.zeros((SizeImage[3],SizeImage[2]),numpy.uint8())

        if objet ["geometry"]["type"]=="Polygon":

            # Create a Polygon from the coordinates

            poly = Polygon(objet["geometry"]["coordinates"][0])

            "conversion du polygon en numarray arrondi"

            contour=numpy.uint32(numpy.around(numpy.asarray(poly.exterior.coords)))

            "creation de la matrice de soustraction pour revenir à un depart à 0,0"

            try:

                bc_arr, bc_row_means = numpy.broadcast_arrays(contour, numpy.array([SizeImage[0],SizeImage[1]]))

            except:

                print('erreur')

            polygone=bc_arr-bc_row_means

            "dessin du polygon sur l'image pour faire un masque"

            c,r=numpy.split(polygone,2,axis=1)

            rr, cc = polygon(r.reshape(r.shape[0]), c.reshape(c.shape[0]))

            tab=numpy.where(rr>=SizeImage[3])

            rr[tab]=SizeImage[3]-1

            tab=numpy.where(cc>=SizeImage[2])

            cc[tab]=SizeImage[2]-1

            try:

                ROI[rr,cc]=255

                Tab=numpy.where(ROI==255)

                if objet["properties"]["classification"]["name"]=="sombre":

                    Dictionnaire[compteurU]=1

                    TotalUnique[Tab]=compteurU

                    compteurU+=1

                if objet["properties"]["classification"]["name"]=="clair":

                    Dictionnaire[compteurU]=2

                    TotalUnique[Tab]=compteurU

                    compteurU+=1

            except :

                print("erreur 2")

                pass

    return TotalUnique,Dictionnaire

```

%% Cell type:code id: tags:

``` python

"dossier contenant les contours dessinées sur qupath et exportés en json"

Liste=sorted(glob.glob("Examples/qupathProject/ExportJSON/*.json"))

"dossier des images d'origines"

dossierimage="Examples/"

"parcours des fichiers de contours"

for json_detections in Liste:

        namefile=json_detections[:-5]

        #namefiletif=json_detections[:-5].split("/")[-1]+".tif"

        namefiletif=json_detections[:-19].split("/")[-1]

        with open(json_detections) as g:

                alldetections = geojson.load(g)

        imagecouleur=io.imread(dossierimage+namefiletif)

        "Contour design, labeling and dictionary creation for classification"

        TotalUnique,my_dict=drawnucleiFinal([0,0,1544,1038],alldetections)

        "lecture de l'image couleur d'origine"

        imagecouleur=io.imread(dossierimage+namefiletif)

        nombre=0

        "decoupage de l'image en carré de 256x256"

        "on fait 24 images avec la taille de l'image "

        for x in range(0,6):

                for y in range(0,4):

                        Dictionnaire={}

                        imagcropC=imagecouleur[(7+y*256):(7+(y+1)*256),(4+x*256):(4+(x+1)*256)]

                        imagcropL=TotalUnique[(7+y*256):(7+(y+1)*256),(4+x*256):(4+(x+1)*256)]

                        labelunique=numpy.unique(imagcropL)

                        "sauvegarde du crop couleur et label associés"

                        io.imsave(namefile+"_"+str(nombre)+"_TU_C.tif",imagcropC,check_contrast=False)

                        io.imsave(namefile+"_"+str(nombre)+"_TU_L.tif",imagcropL,check_contrast=False)

                        for l in labelunique[1:]:

                                Dictionnaire[int(l)]=my_dict.get(l)

                        "sauvegarde du dico ave les labels indivs por les images 256*256"

                        tf1 = open(namefile+"_"+str(nombre)+"_TU.json", "w")

                        geojson.dump(Dictionnaire,tf1)

                        tf1.close()

                        nombre+=1

```

%% Cell type:code id: tags:

``` python

```