Statistiques
| Révision :

root / bin / normalise.py @ 21

Historique | Voir | Annoter | Télécharger (2,88 ko)

1
#!/usr/bin/env python
2

    
3
# before lounching the executable activate the timagetk conda environement
4
# conda activate florivar
5

    
6
# general imports :
7
from sys import path, argv
8
import os
9
import time
10
import numpy as np
11
from skimage.filters import threshold_otsu
12

    
13
from titk_tools.io import imread, imsave, SpatialImage
14
#from timagetk.plugins import linear_filtering
15
from scipy.ndimage import gaussian_filter
16

    
17
# multiprocessing
18
from multiprocessing import Pool
19
nproc = 20 # number of processors
20

    
21
print('=======================================')
22

    
23
def gauss_titk(img, radius):
24
        '''
25
        gaussian filter with radius in real units
26
        - if radius=a number : isotropic
27
        - if radius= a list of numbers : anisotropic
28
        '''
29
        img_filtered = linear_filtering(img, std_dev=radius, method='gaussian_smoothing')#, real=True)
30
        return img_filtered
31

    
32

    
33
def gauss_nd(img, radius):
34
        '''
35
        gaussian filter with radius in real units
36
        - if radius=a number : isotropic
37
        - if radius= a list of numbers : anisotropic
38
        '''
39
        img_filtered = gaussian_filter(img, radius)
40
        return img_filtered
41

    
42
def iterated_gauss_2d(img):
43
        vs=img.voxelsize
44
        rad=2*vs[2]/vs[0]
45
        #rad=10
46
        radius = [rad, rad, 0]
47
        img_filtered=gaussian_filter(img, radius)
48
        rad=rad/2
49
        while rad>1:
50
                print(rad)
51
                radius=[rad,rad,0]
52
                img_filtered = gaussian_filter(img_filtered, radius)
53
                rad=rad/2
54
        return SpatialImage(img_filtered, dtype='uint8', voxelsize=vs, origin=[0, 0, 0])
55

    
56

    
57
indir = argv[1]
58
T = int(argv[2])
59

    
60
newdir=indir+'_normalised-OtsuT'+str(T)
61
os.system("mkdir -p "+newdir)
62

    
63
print('Normalising images in the folder:',indir)
64
print('Common Otsu threshold for normalisation:',T)
65
print('Results will be written in the folder:',newdir)
66

    
67
pathdir=indir
68
listfiles=os.listdir(pathdir)
69

    
70
def treat_file(i):
71
        filename = listfiles[i]
72
        extension = filename.split('.')[-1]
73
        if extension=='tif':
74
                print('--------------------------')
75
                print(filename)
76
                im = imread(pathdir+'/'+filename)
77
                # filter for computing Otsu
78
                #imf = iterated_gauss_2d(im)
79
                imf=im
80
                if len(np.unique(imf))>0:#20:
81
                        iotsu = threshold_otsu(imf)
82
                        # normalise the original image
83
                        imin = np.min(im)+0.0
84
                        imax= np.max(im)+0.0
85
                        a = T/(iotsu-imin)
86
                        b = - imin*a
87
                        transformed_im = SpatialImage(a*im+b, dtype='uint16', voxelsize=im.voxelsize, origin=[0, 0, 0])
88
                        transformed_im[transformed_im>255]=255                 
89
                        transformed_im = SpatialImage(transformed_im, dtype='uint8', voxelsize=im.voxelsize, origin=[0, 0, 0])
90
                        tiotsu = threshold_otsu(transformed_im)
91
                        imf = iterated_gauss_2d(transformed_im)
92
                        tiotsuf = threshold_otsu(imf)
93
                        print('   otsu=',iotsu,'  otsu normalised=',tiotsu,'  otsu normalised filtered=',tiotsuf)
94
                        newname = 'norm'+str(T)+'-'+filename
95
                        newnamef = 'filtered-'+newname
96
                        #imsave(newdir+'/'+newname,transformed_im)
97
                        imsave(newdir+'/'+newnamef,imf)
98
        return 0
99

    
100

    
101
inputs = range(len(listfiles))
102

    
103

    
104
for i in range(len(inputs)):
105
        treat_file(i)
106

    
107

    
108
'''
109
pool = Pool(processes=nproc)
110
pool.map(treat_file, inputs)
111
pool.close()
112
pool.join()
113
'''