Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Helping functions for analysis of lif microscopy images
"""
import os
import skimage.io as io
import skimage.transform
import numpy as np
import local_features as lf
import matplotlib.pyplot as plt
#%% General functions
# sort list of directories from selected string position
def sort_strings_customStart(strings, string_start_pos = 0):
"""Sorts list of strings
Input(s)
strings: list of strings
string_start_pos (default 0): defines the section of the string on which the ordering is based.
Output(s)
strings_sortered: list of sorted strings
Author: Monica Jane Emerson (monj@dtu.dk)"""
strings_roi = [i[string_start_pos:] for i in strings]
#key indicates that the indices (k) are sorted according to subjects_list[k]
ind = np.array(sorted(range(len(strings_roi)), key=lambda k: strings_roi[k]))
strings_sorted = [strings[i] for i in ind]
return strings_sorted
#Provide the contents of a directory sorted
def listdir_custom(directory, string_start_pos = 0, dir_flag = False, base_name = False):
'monj@dtu.dk'
if dir_flag:
if base_name:
list_dirs = [dI for dI in os.listdir(directory) if (os.path.isdir(os.path.join(directory,dI)) & dI[0:len(base_name)]==base_name)]
else:
list_dirs = [dI for dI in os.listdir(directory) if os.path.isdir(os.path.join(directory,dI))]
else:
if base_name:
list_dirs = [dI for dI in os.listdir(directory)if dI[0:len(base_name)]==base_name]
else:
list_dirs = [dI for dI in os.listdir(directory)]
listdirs_sorted = sort_strings_customStart(list_dirs,string_start_pos)
return listdirs_sorted
def flatten_list(list):
'monj@dtu.dk'
list_flat = [item for sublist in list for item in sublist]
return list_flat
#%% IO functions
#make directory, and subdirectories within, if they don't exist
def make_output_dirs(directory,subdirectories = False):
'monj@dtu.dk'
os.makedirs(directory, exist_ok = True)
if subdirectories:
for subdir in subdirectories:
os.makedirs(directory + subdir + '/', exist_ok = True)
# def make_output_dirs(directory,subdirectories):
# 'monj@dtu.dk'
# os.makedirs(directory, exist_ok = True)
# os.makedirs(directory + 'control/', exist_ok = True)
# for disease in diseases:
# os.makedirs(directory + disease + '/', exist_ok = True)
#Reads images starting with base_name from the subdirectories of the input directory
def read_max_imgs(dir_condition, base_name, sc_fac = 1, colour_mode = 'colour'):
'monj@dtu.dk'
sample_list = listdir_custom(dir_condition, string_start_pos = -4, dir_flag = True)
#print(sample_list)
max_img_list = []
for sample in sample_list:
sample_dir = dir_condition + '/' + sample + '/'
frame_list = listdir_custom(sample_dir, base_name = base_name)
#print(frame_list)
frame_img_list = []
for frame in frame_list:
frame_path = sample_dir + frame
#Option to load in bnw or colour
if colour_mode == 'bnw':
img = io.imread(frame_path, as_gray = True).astype('uint8')
frame_img_list += [skimage.transform.rescale(img, sc_fac, preserve_range = True).astype('uint8')]
else:
img = io.imread(frame_path).astype('uint8')
#print(img.dtype)
frame_img_list += [skimage.transform.rescale(img, sc_fac, preserve_range = True, multichannel=True).astype('uint8')]
max_img_list += [frame_img_list]
#print(frame_img_list[0].dtype)
return max_img_list
#%% Functions for intensity inspection and image preprocessing
# computes the maximum projection image from a directory of images
def get_max_img(in_dir, ext = '.png', n_img = 0):
file_names = [f for f in os.listdir(in_dir) if f.endswith(ext)]
file_names.sort()
if ( n_img < 1 ):
n_img = len(file_names)
img_in = io.imread(in_dir + file_names[0])
for i in range(1,n_img):
img_in = np.maximum(img_in, io.imread(in_dir + file_names[i]))
return img_in
# computes a list of maximum projection images from a list of directories
def compute_max_img_list(in_dir, ext, base_name, dir_out = ''):
"""by abda
Modified by monj"""
dir_list = [dI for dI in os.listdir(in_dir) if (os.path.isdir(os.path.join(in_dir,dI)) and dI[0:len(base_name)]==base_name)]
dir_list.sort()
max_img_list = []
for d in dir_list:
image_dir_in = in_dir + d + '/'
max_img = get_max_img(image_dir_in, ext)
if dir_out!='':
os.makedirs(dir_out, exist_ok = True)
io.imsave(dir_out + d + '.png', max_img.astype('uint8'))
max_img_list += [max_img]
return max_img_list
# One more level up from compute_max_img_list. Computes a list of lists of maximum
#projection images from a list of directories, each containing a list of directories
#with the set of images that should be combined into a maximum projection.
def comp_max_imgs(dir_condition, base_name, dir_out = ''):
'monj@dtu.dk'
dir_list_condition = listdir_custom(dir_condition, string_start_pos = -4, dir_flag = True)
max_img_list_condition = []
for directory in dir_list_condition:
dir_in = dir_condition + '/' + directory + '/'
if dir_out!= '':
subdir_out = dir_in.replace(dir_condition,dir_out)
else:
subdir_out = ''
max_img_condition = compute_max_img_list(dir_in, '.png', base_name, subdir_out)
max_img_list_condition+= [max_img_condition]
return max_img_list_condition
#TO DO: Eliminate histogram part from this function
def comp_std_perChannel(max_im_list, th_int, dir_condition):
'monj@dtu.dk'
sample_list = listdir_custom(dir_condition, string_start_pos = -4, dir_flag = True)
std_list = [[],[],[]]
for sample, frame_img_list in zip(sample_list, max_im_list):
for ind,img in enumerate(frame_img_list):
h, w, channels = img.shape
for channel in range(0,channels):
intensities = img[:,:,channel].ravel()
std_list[channel] += [(intensities[intensities>th_int]).std()]
return std_list
def intensity_spread_normalisation(img_list, th_int, mean_std, dir_condition, base_name, dir_results):
'monj@dtu.dk'
sample_list = listdir_custom(dir_condition, string_start_pos = -4, dir_flag = True)
#print(sample_list)
img_corr_list = []
for sample, sample_imgs in zip(sample_list, img_list):
sample_dir = dir_condition + '/' + sample + '/'
#print(sample_dir)
frame_list = listdir_custom(sample_dir, base_name = base_name)
#print(frame_list)
frame_img_corr_list = []
for frame,img in zip(frame_list,sample_imgs):
h, w, channels = img.shape
#img_corr = np.empty(img.shape,dtype = 'uint8')
img_corr = img
for channel in range(0,channels):
img_channel = img_corr[:,:,channel]
img_channel[img_channel>th_int] = img_channel[img_channel>th_int]*(mean_std/img_channel.std())
img_corr[:,:,channel] = img_channel
frame_img_corr_list += [img_corr]
os.makedirs(dir_results + '/' + sample, exist_ok = True)
io.imsave(dir_results + '/' + sample + '/' + frame +'.png', img_corr)
img_corr_list += [frame_img_corr_list]
return img_corr_list
#def plotHist_perChannel_imgset
def plotHist_perChannel_imgset_list(max_img_list, dir_condition, dir_results = '', name_tag = 'original'):
'monj@dtu.dk'
sample_list = listdir_custom(dir_condition, string_start_pos = -4, dir_flag = True)
for sample, frame_img_list in zip(sample_list, max_img_list):
fig, axs = plt.subplots(4,len(frame_img_list), figsize = (len(frame_img_list)*2,4*2))
plt.suptitle('Sample '+sample[-4:] + ', acq. date: ' + sample[:6])
for ind,img in enumerate(frame_img_list):
h, w, channels = img.shape
axs[0][ind].imshow(img)
for channel in range(0,channels):
intensities = img[:,:,channel].ravel()
axs[channel+1][ind].hist(intensities, bins = 50)
axs[channel+1][ind].set_aspect(1.0/axs[channel+1][ind].get_data_ratio())
if dir_results!= '':
plt.savefig(dir_results + '/' + sample[-4:] + '_' + name_tag + '_perChannelHistograms.png', dpi = 300)
plt.close(fig)
#def compare_imgpairs
def compare_imgpairs_list(list1_imgsets, list2_imgsets, dir_condition, dir_results = ''):
'monj@dtu.dk'
sample_list = listdir_custom(dir_condition, string_start_pos = -4, dir_flag = True)
for imgset1, imgset2, sample in zip(list1_imgsets, list2_imgsets, sample_list):
fig, axs = plt.subplots(2,len(imgset1), figsize = (2*len(imgset1),2*2),sharex=True,sharey=True)
plt.suptitle('Sample '+sample[-4:] + ', acq. date: ' + sample[:6])
for ind, img1, img2 in zip(range(0,len(imgset1)), imgset1, imgset2):
axs[0][ind].imshow(img1)
axs[1][ind].imshow(img2)
if dir_results!= '':
plt.savefig(dir_results + '/' + sample[-4:] + '_originalVScorrected.png', dpi=300)
plt.close(fig)
#%%Functions for Feature analysis
# computes the max projection image and the features into a list
def compute_max_img_feat(in_dir, ext, base_name, sigma, sc_fac, save = False, abs_intensity = True):
dir_list = [dI for dI in os.listdir(in_dir) if (os.path.isdir(os.path.join(in_dir,dI)) and dI[0:len(base_name)]==base_name)]
dir_list.sort()
max_img_list = []
for d in dir_list:
image_dir = in_dir + d + '/'
max_img = get_max_img(image_dir, ext)
if save:
dir_name_out_img = in_dir.replace('data','maxProjImages')
os.makedirs(dir_name_out_img, exist_ok = True)
io.imsave( dir_name_out_img + d + '.png', max_img.astype('uint8'))
max_img_list += [skimage.transform.rescale(max_img, sc_fac, multichannel=True)]
feat_list = []
for max_img in max_img_list:
r,c = max_img.shape[:2]
feat = np.zeros((r*c,45))
for i in range(0,3):
feat_tmp = lf.get_gauss_feat_im(max_img[:,:,i], sigma)
feat[:,i*15:(i+1)*15] = feat_tmp.reshape(-1,feat_tmp.shape[2])
if not(abs_intensity):
feat_list += [np.concatenate((feat[:,1:15],feat[:,16:30],feat[:,31:45]),axis = 1)]
else:
feat_list += [feat]
return max_img_list, feat_list
# computes a histogram of features from a kmeans object
def compute_assignment_hist(feat_list, kmeans, background_feat = None):
assignment_list = []
for feat in feat_list:
assignment_list += [kmeans.predict(feat)]
edges = np.arange(kmeans.n_clusters+1)-0.5 # histogram edges halfway between integers
hist = np.zeros(kmeans.n_clusters)
for a in assignment_list:
hist += np.histogram(a,bins=edges)[0]
sum_hist = np.sum(hist)
hist/= sum_hist
if background_feat is not None:
assignment_back = kmeans.predict(background_feat)
hist_back = np.histogram(assignment_back,bins=edges)[0]
return hist, assignment_list, hist_back, assignment_back
else:
return hist, assignment_list
# image to array of patches
def im2col(A, BSZ, stepsize=1, norm=False):
# Parameters
m,n = A.shape
s0, s1 = A.strides
nrows = m-BSZ[0]+1
ncols = n-BSZ[1]+1
shp = BSZ[0],BSZ[1],nrows,ncols
strd = s0,s1,s0,s1
out_view = np.lib.stride_tricks.as_strided(A, shape=shp, strides=strd)
out_view_shaped = out_view.reshape(BSZ[0]*BSZ[1],-1)[:,::stepsize]
if norm:
one_norm = np.sum(out_view_shaped,axis=0)
ind_zero_norm = np.where(one_norm !=0)
out_view_shaped[:,ind_zero_norm] = 255*out_view_shaped[:,ind_zero_norm]/one_norm[ind_zero_norm]
return out_view_shaped
# nd image to array of patches
def ndim2col(A, BSZ, stepsize=1, norm=False):
if(A.ndim == 2):
return im2col(A, BSZ, stepsize, norm)
else:
r,c,l = A.shape
patches = np.zeros((l*BSZ[0]*BSZ[1],(r-BSZ[0]+1)*(c-BSZ[1]+1)))
for i in range(l):
patches[i*BSZ[0]*BSZ[1]:(i+1)*BSZ[0]*BSZ[1],:] = im2col(A[:,:,i],BSZ,stepsize,norm)
return patches
# nd image to array of patches with mirror padding along boundaries
def ndim2col_pad(A, BSZ, stepsize=1, norm=False):
r,c = A.shape[:2]
if (A.ndim == 2):
l = 1
else:
l = A.shape[2]
tmp = np.zeros((r+BSZ[0]-1,c+BSZ[1]-1,l))
fhr = int(np.floor(BSZ[0]/2))
fhc = int(np.floor(BSZ[1]/2))
thr = int(BSZ[0]-fhr-1)
thc = int(BSZ[1]-fhc-1)
tmp[fhr:fhr+r,fhc:fhc+c,:] = A.reshape((r,c,l))
tmp[:fhr,:] = np.flip(tmp[fhr:fhr*2,:], axis=0)
tmp[fhr+r:,:] = np.flip(tmp[r:r+thr,:], axis=0)
tmp[:,:fhc] = np.flip(tmp[:,fhc:fhc*2], axis=1)
tmp[:,fhc+c:] = np.flip(tmp[:,c:c+thc], axis=1)
tmp = np.squeeze(tmp)
return ndim2col(tmp,BSZ,stepsize,norm)
# def plot_mapsAndimages(dir_condition, directory_list, map_img, max_img_list, base_name = 'frame', r = 1024, c = 1024):
# nr_list = 0
# for directory in directory_list:
# in_dir = dir_condition + directory + '/'
# dir_list = [dI for dI in os.listdir(in_dir) if (os.path.isdir(os.path.join(in_dir ,dI)) and dI[0:len(base_name)]==base_name)]
# print(dir_list)
# img_all_control = []
# img_all_control += map_img[nr_list:nr_list+len(dir_list)]
# img_all_control += max_img_list[nr_list:nr_list+len(dir_list)]
# fig, axs = plt.subplots(2,len(dir_list), sharex=True, sharey=True)
# nr = 0
# prob_sample = np.zeros((len(dir_list),1))
# for ax, img in zip(axs.ravel(), img_all_control):
# print(nr)
# if nr<len(dir_list):
# prob_img = sum(sum(img))/(r*c)
# ax.set_title('Pc '+str(round(prob_img,2)))
# prob_sample[nr] = prob_img
# #print(prob_sample)
# nr += 1
# if ( img.ndimg == 2 ):
# ax.imshow(skimage.transform.resize(img, (1024,1024)), cmap=plt.cm.bwr, vmin = 0, vmax = 1)
# else:
# ax.imshow(skimage.transform.resize(img.astype(np.uint8), (1024,1024)))
# plt.suptitle('Probability control of sample '+str(directory)+', avg:' + str(round(prob_sample.mean(),2))+ ' std:' + str(round(prob_sample.std(),2)))
# plt.savefig(dir_probs + '/control/' + '/probImControl_'+str(directory)+'_%dclusters_%dsigma%ddownscale_absInt%s.png'%(nr_clusters,sigma,1/sc_fac,abs_intensity), dpi=1000)
# plt.show()
# nr_list += len(dir_list)