First commit

pattern_reco.py

+__author__ = "Benjamin GALLOIS"
+__license__ = "GPL"
+__version__ = "1.0.0"
+__pythonversion__ = "3.6.0"
+__maintainer__ = "Benjamin GALLOIS"
+__email__ = "benjamin.gallois@curie.fr"
+__status__ = "Production"
+
+
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy as sc
+import skimage
+from skimage import io
+from skimage.io.collection import ImageCollection,MultiImage
+import glob
+from skimage import filters as filters
+import os
+from skimage.morphology import disk
+from skimage.filters.rank import threshold
+from skimage.filters.rank import mean
+from skimage.filters.rank import enhance_contrast
+import pandas as pd
+
+
+
+def stack_importation(chemin):
+	'''This function imports a stack of images .TIF in a 3D matrix, it takes in arguments the path of the folder.'''
+	img = MultiImage(chemin)
+	size = len(img)
+	shape = img[0].shape
+
+	stack = np.zeros((shape[0],shape[1],size))
+	for i in range(size):
+		stack[:,:,i] = img[i]
+	return stack, size, shape
+
+
+def number_MIP(matrix,shape,number):
+	'''This function makes a N-maximums intensity projection of a 3D matrix, it takes in arguments the 3D matrix, the shape of the 3D matri in a tuple and the number of maximas to project.'''
+	projection = np.zeros(shape)
+	for x in range(shape[0]):
+		for y in range(shape[1]):
+			array_sort = np.sort(matrix[x, y, :])
+			projection[x, y] = np.sum(array_sort[-number::])
+	return projection
+
+
+
+def main(chemin,num_max):
+	'''It's the main function of the program, it takes the path to a folder where the stacks .TIFF are placed, it projects the N-maximums, sets the background to 0 and computes the mean intensity normalized by the numbers of pixels non sets to 0.
+It takes in argument the path of a folder and return the mean intensity.'''
+	(img, size, shape) = stack_importation(chemin)
+	projection_corrected = number_MIP(img,shape,num_max)
+	return projection_corrected
+
+
+
+
+def otsufilter(img):
+    '''This function applies an Otsu filter, it takes in arguments a matrix and return a mask where the background is set to 0 and the foreground to 1. Change nbins with the type of images (8bits = 256, 16bits = 65536'''
+    val = filters.threshold_otsu(img, nbins = 65536)
+    mask = img < val
+    mask = np.invert(mask)
+    mask = sc.ndimage.binary_fill_holes(mask)
+    return mask
+
+
+
+
+def analysis(image, maskOtsu, meanRadius, contrastRadius, thresholdRadius, withImage, chemin):
+    '''This function selects a pattern of maximum intensity and computes the mean intensity of the pattern and the mean intensity of the inverted pattern
+        Args :
+            image = matrix uint16 of the image
+            maskOtsu = mask of the embryo without the background
+            meanRadius = parameter which indicate the radius of the disk where is computed the local mean.
+            constrastRadius = parameter which indicate the radius of the disk where is computed the local constrast enhancement.
+            thresholdRadius = parameter which indicate the radius of the disk where is computed the local.
+            withImage = if 0 it doesn't save image, if 1 it saves the invert of the pattern and the fill image, if 2 it saves the pattern and the full image. 
+            chemin = pathway to the file where the image stack is located.
+            
+        return :
+            The mask of the pattern selected.
+            The mean intensity of the pattern.
+            The mean intensity of the inverted pattern.
+    '''
+    
+    maxima = mean(image, disk(meanRadius), mask=maskOtsu)
+    maxima = enhance_contrast(maxima, disk(contrastRadius), mask=maskOtsu)
+    mask3 = threshold(maxima,disk(thresholdRadius), mask=maskOtsu)
+    mask3 = mask3.astype(np.bool)
+    patternImage = image*maskOtsu*mask3
+    invertPatternImage = image*maskOtsu*np.invert(mask3)  
+    
+    
+    if withImage == 0:   
+        sum1 = np.sum(patternImage) / np.count_nonzero(patternImage)
+        sum2 = np.sum(invertPatternImage) / np.count_nonzero(invertPatternImage)
+    
+    if withImage == 1:
+        sum1 = np.sum(patternImage) / np.count_nonzero(patternImage)
+        sum2 = np.sum(invertPatternImage) / np.count_nonzero(invertPatternImage)
+        name2 =str("projected_"+str(chemin)+ "full.tif")
+        name3 =str("projected_"+str(chemin)+ "dark.tif")
+        io.imsave(os.path.join('Projection', name3), invertPatternImage)
+        io.imsave(os.path.join('Projection', name2), face)
+        
+    if withImage == 2:
+        sum1 = np.sum(patternImage) / np.count_nonzero(patternImage)
+        sum2 = np.sum(invertPatternImage) / np.count_nonzero(invertPatternImage)
+        name =str("projected_"+str(chemin)+".tif")
+        io.imsave(os.path.join('Projection', name), patternImage)
+
+
+    return mask3*maskOtsu, sum1, sum2
+
+
+
+
+def relativeIntensity(image,chemin):
+    '''This function selects a pattern with two successives application of the analysis function, and computes the pattern of the second analysis divided by the invert of the pattern of the first analysis.
+        Argz :
+            image = matrix of the image.
+            chemin = pathway of the stack image.
+        return :
+            ratio of mean intensities.
+    '''
+    
+    
+    image = image.astype(np.uint16)
+    mask = otsufilter(image)
+    c,d,e = analysis(image, mask,200,150,200, 0, chemin)
+    g,h,i = analysis(image, c , 60, 50, 500, 0, chemin)
+    return h/e
+
+
+
+pattern = ['Const', 'Inv-', ' ']
+std = []
+data = []
+means = []
+DATA = []
+NAME= []
+
+for i in range(3):
+    cheminFace = (glob.glob(pattern[i] + '/' + '*491.TIF*'))
+    relativeInt = []
+    for j in cheminFace:
+        face = main(j, 1)
+        NAME.append(j)
+        relativeInt.append(relativeIntensity(face,j))
+        DATA.append(j)
+        
+
+    data.append(relativeInt)
+    means.append(np.mean(relativeInt))
+    std.append(np.std(relativeInt))
+
+
+save = pd.Series(DATA, index=NAME)
+save.to_csv('out.csv', sep=',')
+
+
+
+'''ind = np.arange(3)  # the x locations for the groups
+width = 0.2     # the width of the bars
+fig, ax = plt.subplots()
+rects1 = ax.bar(ind, means, width, color='y', yerr=std)
+ax.set_xticks(ind + width / 2)
+ax.set_xticklabels(pattern)
+
+def label_diff(i,j,text,X,Y):
+	#This function draws a line between two bars of the bar chart and displays a text on top of it.It takes in arguments the i-th and the j-th bar of the bar chart, the text to display, the center of the bars as an array and the y value of the bars as an array
+	x = (X[i]+X[j])/2
+	y = 1.1*max(Y[i], Y[j])
+	dx = abs(X[i]-X[j])
+	props = {'arrowstyle':'|-|',\
+                 'lw':2}
+	ax.annotate('', xy=(X[i],y+0.2*(i+j)), xytext=(X[j],y+0.2*(i+j)), arrowprops=props)
+	ax.annotate(text, xy=(X[i]+dx/2,y+0.2*(i+j)+0.05), zorder=10)
+
+
+for i in range(3):
+    for j in range(3):
+        if j > 0 and i < j :
+            st = np.round(sc.stats.mannwhitneyu(data[i],data[j]),6)
+            stat = ('p = ' + str(st[1]))
+            print(stat)
+            label_diff(i,j,stat,ind+width*.5,means)
+plt.show()'''

pattern_reco_actine.py

+__author__ = "Benjamin GALLOIS"
+__license__ = "GPL"
+__version__ = "1.0.0"
+__pythonversion__ = "3.6.0"
+__maintainer__ = "Benjamin GALLOIS"
+__email__ = "benjamin.gallois@curie.fr"
+__status__ = "Production"
+
+
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy as sc
+import skimage
+from skimage import io
+from skimage.io.collection import ImageCollection,MultiImage
+import glob
+from skimage import filters as filters
+import os
+from skimage.morphology import disk
+from skimage.filters.rank import threshold
+from skimage.filters.rank import mean
+from skimage.filters.rank import enhance_contrast
+import pandas as pd
+
+
+def nameReco(cheminA, cheminBList):
+    for i, j in enumerate(cheminBList):
+        if cheminA[5:-8] == j[len('/run/media/benjamin/BENJAMIN/17-03-23/'):-8]:
+            return j
+            break
+
+
+    
+
+
+
+def stack_importation(chemin):
+	'''This function imports a stack of images .TIF in a 3D matrix, it takes in arguments the path of the folder.'''
+	img = MultiImage(chemin)
+	size = len(img)
+	shape = img[0].shape
+
+	stack = np.zeros((shape[0],shape[1],size))
+	for i in range(size):
+		stack[:,:,i] = img[i]
+	return stack, size, shape
+
+
+def number_MIP(matrix,shape,number):
+	'''This function makes a N-maximums intensity projection of a 3D matrix, it takes in arguments the 3D matrix, the shape of the 3D matri in a tuple and the number of maximas to project.'''
+	projection = np.zeros(shape)
+	for x in range(shape[0]):
+		for y in range(shape[1]):
+			array_sort = np.sort(matrix[x, y, :])
+			projection[x, y] = np.sum(array_sort[-number::])
+	return projection
+
+
+
+def main(chemin,num_max):
+	'''It's the main function of the program, it takes the path to a folder where the stacks .TIFF are placed, it projects the N-maximums, sets the background to 0 and computes the mean intensity normalized by the numbers of pixels non sets to 0.
+It takes in argument the path of a folder and return the mean intensity.'''
+	(img, size, shape) = stack_importation(chemin)
+	projection_corrected = number_MIP(img,shape,num_max)
+	return projection_corrected
+
+
+
+
+def otsufilter(img):
+    '''This function applies an Otsu filter, it takes in arguments a matrix and return a mask where the background is set to 0 and the foreground to 1. Change nbins with the type of images (8bits = 256, 16bits = 65536'''
+    val = filters.threshold_otsu(img, nbins = 65536)
+    mask = img < val
+    mask = np.invert(mask)
+    mask = sc.ndimage.binary_fill_holes(mask)
+    return mask
+
+
+
+
+def analysis(image, maskOtsu, meanRadius, contrastRadius, thresholdRadius, withImage, chemin):
+    '''This function selects a pattern of maximum intensity and computes the mean intensity of the pattern and the mean intensity of the inverted pattern
+        Args :
+            image = matrix uint16 of the image
+            maskOtsu = mask of the embryo without the background
+            meanRadius = parameter which indicate the radius of the disk where is computed the local mean.
+            constrastRadius = parameter which indicate the radius of the disk where is computed the local constrast enhancement.
+            thresholdRadius = parameter which indicate the radius of the disk where is computed the local.
+            withImage = if 0 it doesn't save image, if 1 it saves the invert of the pattern and the fill image, if 2 it saves the pattern and the full image. 
+            chemin = pathway to the file where the image stack is located.
+            
+        return :
+            The mask of the pattern selected.
+            The mean intensity of the pattern.
+            The mean intensity of the inverted pattern.
+    '''
+    
+    maxima = mean(image, disk(meanRadius), mask=maskOtsu)
+    maxima = enhance_contrast(maxima, disk(contrastRadius), mask=maskOtsu)
+    mask3 = threshold(maxima,disk(thresholdRadius), mask=maskOtsu)
+    mask3 = mask3.astype(np.bool)
+    patternImage = image*maskOtsu*mask3
+    invertPatternImage = image*maskOtsu*np.invert(mask3)  
+    
+    
+    if withImage == 0:   
+        sum1 = np.sum(patternImage) / np.count_nonzero(patternImage)
+        sum2 = np.sum(invertPatternImage) / np.count_nonzero(invertPatternImage)
+    
+    if withImage == 1:
+        sum1 = np.sum(patternImage) / np.count_nonzero(patternImage)
+        sum2 = np.sum(invertPatternImage) / np.count_nonzero(invertPatternImage)
+        name2 =str("projected_"+str(chemin)+ "full.tif")
+        name3 =str("projected_"+str(chemin)+ "dark.tif")
+        io.imsave(os.path.join('Projection', name3), invertPatternImage)
+        io.imsave(os.path.join('Projection', name2), face)
+        
+    if withImage == 2:
+        sum1 = np.sum(patternImage) / np.count_nonzero(patternImage)
+        sum2 = np.sum(invertPatternImage) / np.count_nonzero(invertPatternImage)
+        name =str("projected_"+str(chemin)+".tif")
+        io.imsave(os.path.join('Projection', name), patternImage)
+
+
+    return mask3*maskOtsu, sum1, sum2, mask3
+
+
+
+
+def relativeIntensity(image, chemin, cheminAct):
+    '''This function selects a pattern with two successives application of the analysis function, and computes the pattern of the second analysis divided by the invert of the pattern of the first analysis.
+        Argz :
+            image = matrix of the image.
+            chemin = pathway of the stack image.
+            cheminActine = pathway of the stack image in the actine channel.
+        return :
+            ratio of mean intensities.
+    '''
+    
+    
+    image = image.astype(np.uint16)
+    mask = otsufilter(image)
+    c,d,e, f = analysis(image, mask,200,150,200, 0, chemin)
+    g,h,i, k = analysis(image, c , 60, 50, 500, 0, chemin)
+    print(cheminAct)
+    imageActine = main(cheminAct, 1)
+    
+    if imageActine.shape == image.shape:
+        invertActinePattern = np.invert(f)*mask*imageActine
+        actinePattern = imageActine*g
+        return (np.sum(actinePattern) / np.count_nonzero(actinePattern))/(np.sum(invertActinePattern) / np.count_nonzero(invertActinePattern))
+    else:
+        return 0
+
+
+pattern = ['Cons', 'Inv-', ' ']
+std = []
+data = []
+means = []
+DATA = []
+NAME= []
+
+for i in pattern:
+    cheminFace = (glob.glob(i + '/' + '*491.TIF*'))
+    cheminActine = (glob.glob('/run/media/benjamin/BENJAMIN/17-03-23/*561.TIF*'))
+    relativeInt = []
+    for k,j in enumerate(cheminFace):
+        actine = nameReco(j, cheminActine)
+        face = main(j, 1)
+        NAME.append(i+j)
+        relativeInt.append(relativeIntensity(face,j,actine))
+        DATA.append(relativeInt[k])
+        
+
+    data.append(relativeInt)
+    means.append(np.mean(relativeInt))
+    std.append(np.std(relativeInt))
+
+
+save = pd.Series(DATA, index=NAME)
+save.to_csv('out.csv', sep=',')
\ No newline at end of file