Mercurial > repos > perssond > unmicst
comparison toolbox/imtools.py @ 0:6bec4fef6b2e draft
"planemo upload for repository https://github.com/ohsu-comp-bio/unmicst commit 73e4cae15f2d7cdc86719e77470eb00af4b6ebb7-dirty"
author | perssond |
---|---|
date | Fri, 12 Mar 2021 00:17:29 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:6bec4fef6b2e |
---|---|
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 |