Mercurial > repos > perssond > coreograph

comparison toolbox/imtools.py @ 0:99308601eaa6 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/ohsu-comp-bio/UNetCoreograph commit fb90660a1805b3f68fcff80d525b5459c3f7dfd6-dirty"
author perssond
date Wed, 19 May 2021 21:34:38 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:99308601eaa6
1 import matplotlib.pyplot as plt
2 import tifffile
3 import os
4 import numpy as np
5 from skimage import io as skio
6 from scipy.ndimage import *
7 from skimage.morphology import *
8 from skimage.transform import resize
9
10 def tifread(path):
11 return tifffile.imread(path)
12
13 def tifwrite(I,path):
14 tifffile.imsave(path, I)
15
16 def imshow(I,**kwargs):
17 if not kwargs:
18 plt.imshow(I,cmap='gray')
19 else:
20 plt.imshow(I,**kwargs)
21
22 plt.axis('off')
23 plt.show()
24
25 def imshowlist(L,**kwargs):
26 n = len(L)
27 for i in range(n):
28 plt.subplot(1, n, i+1)
29 if not kwargs:
30 plt.imshow(L[i],cmap='gray')
31 else:
32 plt.imshow(L[i],**kwargs)
33 plt.axis('off')
34 plt.show()
35
36 def imread(path):
37 return skio.imread(path)
38
39 def imwrite(I,path):
40 skio.imsave(path,I)
41
42 def im2double(I):
43 if I.dtype == 'uint16':
44 return I.astype('float64')/65535
45 elif I.dtype == 'uint8':
46 return I.astype('float64')/255
47 elif I.dtype == 'float32':
48 return I.astype('float64')
49 elif I.dtype == 'float64':
50 return I
51 else:
52 print('returned original image type: ', I.dtype)
53 return I
54
55 def size(I):
56 return list(I.shape)
57
58 def imresizeDouble(I,sizeOut): # input and output are double
59 return resize(I,(sizeOut[0],sizeOut[1]),mode='reflect')
60
61 def imresize3Double(I,sizeOut): # input and output are double
62 return resize(I,(sizeOut[0],sizeOut[1],sizeOut[2]),mode='reflect')
63
64 def imresizeUInt8(I,sizeOut): # input and output are UInt8
65 return np.uint8(resize(I.astype(float),(sizeOut[0],sizeOut[1]),mode='reflect',order=0))
66
67 def imresize3UInt8(I,sizeOut): # input and output are UInt8
68 return np.uint8(resize(I.astype(float),(sizeOut[0],sizeOut[1],sizeOut[2]),mode='reflect',order=0))
69
70 def normalize(I):
71 m = np.min(I)
72 M = np.max(I)
73 if M > m:
74 return (I-m)/(M-m)
75 else:
76 return I
77
78 def snormalize(I):
79 m = np.mean(I)
80 s = np.std(I)
81 if s > 0:
82 return (I-m)/s
83 else:
84 return I
85
86 def cat(a,I,J):
87 return np.concatenate((I,J),axis=a)
88
89 def imerode(I,r):
90 return binary_erosion(I, disk(r))
91
92 def imdilate(I,r):
93 return binary_dilation(I, disk(r))
94
95 def imerode3(I,r):
96 return morphology.binary_erosion(I, ball(r))
97
98 def imdilate3(I,r):
99 return morphology.binary_dilation(I, ball(r))
100
101 def sphericalStructuralElement(imShape,fRadius):
102 if len(imShape) == 2:
103 return disk(fRadius,dtype=float)
104 if len(imShape) == 3:
105 return ball(fRadius,dtype=float)
106
107 def medfilt(I,filterRadius):
108 return median_filter(I,footprint=sphericalStructuralElement(I.shape,filterRadius))
109
110 def maxfilt(I,filterRadius):
111 return maximum_filter(I,footprint=sphericalStructuralElement(I.shape,filterRadius))
112
113 def minfilt(I,filterRadius):
114 return minimum_filter(I,footprint=sphericalStructuralElement(I.shape,filterRadius))
115
116 def ptlfilt(I,percentile,filterRadius):
117 return percentile_filter(I,percentile,footprint=sphericalStructuralElement(I.shape,filterRadius))
118
119 def imgaussfilt(I,sigma,**kwargs):
120 return gaussian_filter(I,sigma,**kwargs)
121
122 def imlogfilt(I,sigma,**kwargs):
123 return -gaussian_laplace(I,sigma,**kwargs)
124
125 def imgradmag(I,sigma):
126 if len(I.shape) == 2:
127 dx = imgaussfilt(I,sigma,order=[0,1])
128 dy = imgaussfilt(I,sigma,order=[1,0])
129 return np.sqrt(dx**2+dy**2)
130 if len(I.shape) == 3:
131 dx = imgaussfilt(I,sigma,order=[0,0,1])
132 dy = imgaussfilt(I,sigma,order=[0,1,0])
133 dz = imgaussfilt(I,sigma,order=[1,0,0])
134 return np.sqrt(dx**2+dy**2+dz**2)
135
136 def localstats(I,radius,justfeatnames=False):
137 ptls = [10,30,50,70,90]
138 featNames = []
139 for i in range(len(ptls)):
140 featNames.append('locPtl%d' % ptls[i])
141 if justfeatnames == True:
142 return featNames
143 sI = size(I)
144 nFeats = len(ptls)
145 F = np.zeros((sI[0],sI[1],nFeats))
146 for i in range(nFeats):
147 F[:,:,i] = ptlfilt(I,ptls[i],radius)
148 return F
149
150 def localstats3(I,radius,justfeatnames=False):
151 ptls = [10,30,50,70,90]
152 featNames = []
153 for i in range(len(ptls)):
154 featNames.append('locPtl%d' % ptls[i])
155 if justfeatnames == True:
156 return featNames
157 sI = size(I)
158 nFeats = len(ptls)
159 F = np.zeros((sI[0],sI[1],sI[2],nFeats))
160 for i in range(nFeats):
161 F[:,:,:,i] = ptlfilt(I,ptls[i],radius)
162 return F
163
164 def imderivatives(I,sigmas,justfeatnames=False):
165 if type(sigmas) is not list:
166 sigmas = [sigmas]
167 derivPerSigmaFeatNames = ['d0','dx','dy','dxx','dxy','dyy','normGrad','normHessDiag']
168 if justfeatnames == True:
169 featNames = [];
170 for i in range(len(sigmas)):
171 for j in range(len(derivPerSigmaFeatNames)):
172 featNames.append('derivSigma%d%s' % (sigmas[i],derivPerSigmaFeatNames[j]))
173 return featNames
174 nDerivativesPerSigma = len(derivPerSigmaFeatNames)
175 nDerivatives = len(sigmas)*nDerivativesPerSigma
176 sI = size(I)
177 D = np.zeros((sI[0],sI[1],nDerivatives))
178 for i in range(len(sigmas)):
179 sigma = sigmas[i]
180 dx = imgaussfilt(I,sigma,order=[0,1])
181 dy = imgaussfilt(I,sigma,order=[1,0])
182 dxx = imgaussfilt(I,sigma,order=[0,2])
183 dyy = imgaussfilt(I,sigma,order=[2,0])
184 D[:,:,nDerivativesPerSigma*i ] = imgaussfilt(I,sigma)
185 D[:,:,nDerivativesPerSigma*i+1] = dx
186 D[:,:,nDerivativesPerSigma*i+2] = dy
187 D[:,:,nDerivativesPerSigma*i+3] = dxx
188 D[:,:,nDerivativesPerSigma*i+4] = imgaussfilt(I,sigma,order=[1,1])
189 D[:,:,nDerivativesPerSigma*i+5] = dyy
190 D[:,:,nDerivativesPerSigma*i+6] = np.sqrt(dx**2+dy**2)
191 D[:,:,nDerivativesPerSigma*i+7] = np.sqrt(dxx**2+dyy**2)
192 return D
193 # derivatives are indexed by the last dimension, which is good for ML features but not for visualization,
194 # in which case the expected dimensions are [plane,channel,y(row),x(col)]; to obtain that ordering, do
195 # D = np.moveaxis(D,[0,3,1,2],[0,1,2,3])
196
197 def imderivatives3(I,sigmas,justfeatnames=False):
198 if type(sigmas) is not list:
199 sigmas = [sigmas]
200
201 derivPerSigmaFeatNames = ['d0','dx','dy','dz','dxx','dxy','dxz','dyy','dyz','dzz','normGrad','normHessDiag']
202
203 # derivPerSigmaFeatNames = ['d0','normGrad','normHessDiag']
204
205 if justfeatnames == True:
206 featNames = [];
207 for i in range(len(sigmas)):
208 for j in range(len(derivPerSigmaFeatNames)):
209 featNames.append('derivSigma%d%s' % (sigmas[i],derivPerSigmaFeatNames[j]))
210 return featNames
211 nDerivativesPerSigma = len(derivPerSigmaFeatNames)
212 nDerivatives = len(sigmas)*nDerivativesPerSigma
213 sI = size(I)
214 D = np.zeros((sI[0],sI[1],sI[2],nDerivatives)) # plane, channel, y, x
215 for i in range(len(sigmas)):
216 sigma = sigmas[i]
217 dx = imgaussfilt(I,sigma,order=[0,0,1]) # z, y, x
218 dy = imgaussfilt(I,sigma,order=[0,1,0])
219 dz = imgaussfilt(I,sigma,order=[1,0,0])
220 dxx = imgaussfilt(I,sigma,order=[0,0,2])
221 dyy = imgaussfilt(I,sigma,order=[0,2,0])
222 dzz = imgaussfilt(I,sigma,order=[2,0,0])
223
224 D[:,:,:,nDerivativesPerSigma*i ] = imgaussfilt(I,sigma)
225 D[:,:,:,nDerivativesPerSigma*i+1 ] = dx
226 D[:,:,:,nDerivativesPerSigma*i+2 ] = dy
227 D[:,:,:,nDerivativesPerSigma*i+3 ] = dz
228 D[:,:,:,nDerivativesPerSigma*i+4 ] = dxx
229 D[:,:,:,nDerivativesPerSigma*i+5 ] = imgaussfilt(I,sigma,order=[0,1,1])
230 D[:,:,:,nDerivativesPerSigma*i+6 ] = imgaussfilt(I,sigma,order=[1,0,1])
231 D[:,:,:,nDerivativesPerSigma*i+7 ] = dyy
232 D[:,:,:,nDerivativesPerSigma*i+8 ] = imgaussfilt(I,sigma,order=[1,1,0])
233 D[:,:,:,nDerivativesPerSigma*i+9 ] = dzz
234 D[:,:,:,nDerivativesPerSigma*i+10] = np.sqrt(dx**2+dy**2+dz**2)
235 D[:,:,:,nDerivativesPerSigma*i+11] = np.sqrt(dxx**2+dyy**2+dzz**2)
236
237 # D[:,:,:,nDerivativesPerSigma*i ] = imgaussfilt(I,sigma)
238 # D[:,:,:,nDerivativesPerSigma*i+1 ] = np.sqrt(dx**2+dy**2+dz**2)
239 # D[:,:,:,nDerivativesPerSigma*i+2 ] = np.sqrt(dxx**2+dyy**2+dzz**2)
240 return D
241 # derivatives are indexed by the last dimension, which is good for ML features but not for visualization,
242 # in which case the expected dimensions are [plane,y(row),x(col)]; to obtain that ordering, do
243 # D = np.moveaxis(D,[2,0,1],[0,1,2])
244
245 def imfeatures(I=[],sigmaDeriv=1,sigmaLoG=1,locStatsRad=0,justfeatnames=False):
246 if type(sigmaDeriv) is not list:
247 sigmaDeriv = [sigmaDeriv]
248 if type(sigmaLoG) is not list:
249 sigmaLoG = [sigmaLoG]
250 derivFeatNames = imderivatives([],sigmaDeriv,justfeatnames=True)
251 nLoGFeats = len(sigmaLoG)
252 locStatsFeatNames = []
253 if locStatsRad > 1:
254 locStatsFeatNames = localstats([],locStatsRad,justfeatnames=True)
255 nLocStatsFeats = len(locStatsFeatNames)
256 if justfeatnames == True:
257 featNames = derivFeatNames
258 for i in range(nLoGFeats):
259 featNames.append('logSigma%d' % sigmaLoG[i])
260 for i in range(nLocStatsFeats):
261 featNames.append(locStatsFeatNames[i])
262 return featNames
263 nDerivFeats = len(derivFeatNames)
264 nFeatures = nDerivFeats+nLoGFeats+nLocStatsFeats
265 sI = size(I)
266 F = np.zeros((sI[0],sI[1],nFeatures))
267 F[:,:,:nDerivFeats] = imderivatives(I,sigmaDeriv)
268 for i in range(nLoGFeats):
269 F[:,:,nDerivFeats+i] = imlogfilt(I,sigmaLoG[i])
270 if locStatsRad > 1:
271 F[:,:,nDerivFeats+nLoGFeats:] = localstats(I,locStatsRad)
272 return F
273
274 def imfeatures3(I=[],sigmaDeriv=2,sigmaLoG=2,locStatsRad=0,justfeatnames=False):
275 if type(sigmaDeriv) is not list:
276 sigmaDeriv = [sigmaDeriv]
277 if type(sigmaLoG) is not list:
278 sigmaLoG = [sigmaLoG]
279 derivFeatNames = imderivatives3([],sigmaDeriv,justfeatnames=True)
280 nLoGFeats = len(sigmaLoG)
281 locStatsFeatNames = []
282 if locStatsRad > 1:
283 locStatsFeatNames = localstats3([],locStatsRad,justfeatnames=True)
284 nLocStatsFeats = len(locStatsFeatNames)
285 if justfeatnames == True:
286 featNames = derivFeatNames
287 for i in range(nLoGFeats):
288 featNames.append('logSigma%d' % sigmaLoG[i])
289 for i in range(nLocStatsFeats):
290 featNames.append(locStatsFeatNames[i])
291 return featNames
292 nDerivFeats = len(derivFeatNames)
293 nFeatures = nDerivFeats+nLoGFeats+nLocStatsFeats
294 sI = size(I)
295 F = np.zeros((sI[0],sI[1],sI[2],nFeatures))
296 F[:,:,:,:nDerivFeats] = imderivatives3(I,sigmaDeriv)
297 for i in range(nLoGFeats):
298 F[:,:,:,nDerivFeats+i] = imlogfilt(I,sigmaLoG[i])
299 if locStatsRad > 1:
300 F[:,:,:,nDerivFeats+nLoGFeats:] = localstats3(I,locStatsRad)
301 return F
302
303 def stack2list(S):
304 L = []
305 for i in range(size(S)[2]):
306 L.append(S[:,:,i])
307 return L
308
309 def thrsegment(I,wsBlr,wsThr): # basic threshold segmentation
310 G = imgaussfilt(I,sigma=(1-wsBlr)+wsBlr*5) # min 1, max 5
311 M = G > wsThr
312 return M