Coverage for /usr/local/lib/python3.8/site-packages/spam/filters/morphologicalOperations.py: 97%

1""" 

2Library of SPAM morphological functions. 

3Copyright (C) 2020 SPAM Contributors 

4 

5This program is free software: you can redistribute it and/or modify it 

6under the terms of the GNU General Public License as published by the Free 

7Software Foundation, either version 3 of the License, or (at your option) 

8any later version. 

9 

10This program is distributed in the hope that it will be useful, but WITHOUT 

11ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 

13more details. 

14 

15You should have received a copy of the GNU General Public License along with 

16this program. If not, see <http://www.gnu.org/licenses/>. 

17""" 

18 

19import numpy 

20import spam.helpers 

21import multiprocessing 

22import progressbar 

23import spam.label.label 

24import scipy.ndimage 

25import skimage.morphology 

26# operations on greyscale images 

27 

28# Global number of processes 

29nProcessesDefault = multiprocessing.cpu_count() 

30 

31 

32def greyDilation(im, nBytes=1): 

33 """ 

34 This function applies a dilation on a grey scale image 

35 

36 Parameters 

37 ----------- 

38 im: numpy array 

39 The input image (greyscale) 

40 nBytes: int, default=1 

41 Number of bytes used to substitute the values on the border. 

42 

43 Returns 

44 -------- 

45 numpy array 

46 The dilated image 

47 """ 

48 # Step 1: check type and dimension 

49 dim = len(im.shape) 

50 # Step 2: Determine substitution value 

51 sub = 2**(8 * nBytes) - 1 

52 # Step 3: apply dilation # x y z 

53 outputIm = im # +0 0 0 

54 outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, 1, 0, sub=sub)) # +1 0 0 

55 outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, -1, 0, sub=sub)) # -1 0 0 

56 outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, 1, 1, sub=sub)) # +0 1 0 

57 outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, -1, 1, sub=sub)) # +0 -1 0 

58 if dim == 3: 

59 outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, 1, 2, sub=sub)) # 0 0 1 

60 outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, -1, 2, sub=sub)) # 0 0 -1 

61 return outputIm 

62 

63 

64def greyErosion(im, nBytes=1): 

65 """ 

66 This function applies a erosion on a grey scale image 

67 

68 Parameters 

69 ----------- 

70 im: numpy array 

71 The input image (greyscale) 

72 nBytes: int, default=1 

73 Number of bytes used to substitute the values on the border. 

74 

75 Returns 

76 -------- 

77 numpy array 

78 The eroded image 

79 """ 

80 # Step 1: check type and dimension 

81 dim = len(im.shape) 

82 # Step 2: Determine substitution value 

83 sub = 2**(8 * nBytes) - 1 

84 # Step 1: apply erosion # x y z 

85 outputIm = im # 0 0 0 

86 outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, 1, 0, sub=sub)) # 1 0 0 

87 outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, -1, 0, sub=sub)) # -1 0 0 

88 outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, 1, 1, sub=sub)) # 0 1 0 

89 outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, -1, 1, sub=sub)) # 0 -1 0 

90 if dim == 3: 

91 outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, 1, 2, sub=sub)) # 0 0 1 

92 outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, -1, 2, sub=sub)) # 0 0 -1 

93 return outputIm 

94 

95 

96def greyMorphologicalGradient(im, nBytes=1): 

97 """ 

98 This function applies a morphological gradient on a grey scale image 

99 

100 Parameters 

101 ----------- 

102 im: numpy array 

103 The input image (greyscale) 

104 nBytes: int, default=1 

105 Number of bytes used to substitute the values on the border. 

106 

107 Returns 

108 -------- 

109 numpy array 

110 The morphologycal gradient of the image 

111 """ 

112 return (greyDilation(im, nBytes=nBytes) - im) 

113 

114 

115# operation on binary images 

116 

117def binaryDilation(im, sub=False): 

118 """ 

119 This function applies a dilation on a binary scale image 

120 

121 Parameters 

122 ----------- 

123 im: numpy array 

124 The input image (greyscale) 

125 sub: bool, default=False 

126 Subtitute value. 

127 

128 Returns 

129 -------- 

130 numpy array 

131 The dilated image 

132 """ 

133 # Step 0: import as bool 

134 im = im.astype(bool) 

135 # Step 1: check type and dimension 

136 dim = len(im.shape) 

137 # Step 1: apply dilation # x y z 

138 outputIm = im # 0 0 0 

139 outputIm = outputIm | spam.helpers.singleShift(im, 1, 0, sub=sub) # 1 0 0 

140 outputIm = outputIm | spam.helpers.singleShift(im, -1, 0, sub=sub) # -1 0 0 

141 outputIm = outputIm | spam.helpers.singleShift(im, 1, 1, sub=sub) # 0 1 0 

142 outputIm = outputIm | spam.helpers.singleShift(im, -1, 1, sub=sub) # 0 -1 0 

143 if dim == 3: 

144 outputIm = outputIm | spam.helpers.singleShift(im, 1, 2, sub=sub) # 0 0 1 

145 outputIm = outputIm | spam.helpers.singleShift(im, -1, 2, sub=sub) # 0 0 -1 

146 return numpy.array(outputIm).astype('<u1') 

147 

148 

149def binaryErosion(im, sub=False): 

150 """ 

151 This function applies a erosion on a binary scale image 

152 

153 Parameters 

154 ----------- 

155 im: numpy array 

156 The input image (greyscale) 

157 sub: bool, default=False 

158 Substitute value. 

159 

160 Returns 

161 -------- 

162 numpy array 

163 The eroded image 

164 """ 

165 # Step 1: apply erosion with dilation --> erosion = ! dilation( ! image ) 

166 return numpy.logical_not(binaryDilation(numpy.logical_not(im), sub=sub)).astype('<u1') 

167 

168 

169def binaryMorphologicalGradient(im, sub=False): 

170 """ 

171 This function applies a morphological gradient on a binary scale image 

172 

173 Parameters 

174 ----------- 

175 im: numpy array 

176 The input image (greyscale) 

177 nBytes: int, default=False 

178 Number of bytes used to substitute the values on the border. 

179 

180 Returns 

181 -------- 

182 numpy array 

183 The morphologycal gradient of the image 

184 """ 

185 return (numpy.logical_xor(binaryDilation(im, sub=sub), im)).astype('<u1') 

186 

187 

188def binaryGeodesicReconstruction(im, marker, dmax=None, verbose=False): 

189 """ 

190 Calculate the geodesic reconstruction of a binary image with a given marker 

191 

192 Parameters 

193 ----------- 

194 im: numpy.array 

195 The input binary image 

196 

197 marker: numpy.array or list 

198 If numpy array: direct input of the marker (must be the size of im) 

199 If list: description of the plans of the image considered as the marker 

200 | ``[1, 0]`` plan defined by all voxels at ``x1=0`` 

201 | ``[0, -1]`` plan defined by all voxels at ``x0=x0_max`` 

202 | ``[0, 0, 2, 5]`` plans defined by all voxels at ``x0=0`` and ``x2=5`` 

203 

204 dmax: int, default=None 

205 The maximum geodesic distance. If None, the reconstruction is complete. 

206 

207 verbose: bool, default=False 

208 Verbose mode 

209 

210 Returns 

211 -------- 

212 numpy.array 

213 The reconstructed image 

214 """ 

215 from spam.errors import InputError 

216 

217 # Step 1: Define marker 

218 if isinstance(marker, list): 

219 # marker based on list of plans 

220 if len(marker) % 2: 

221 raise InputError("marker", explanation="len(marker) must be a multiple of 2") 

222 

223 plans = marker[:] 

224 marker = numpy.zeros(im.shape, dtype=bool) 

225 for i in range(len(plans) // 2): 

226 direction = plans[2 * i] 

227 distance = plans[2 * i + 1] 

228 if len(im.shape) == 2: 

229 if direction == 0: 

230 marker[distance, :] = im[distance, :] 

231 elif direction == 1: 

232 marker[:, distance] = im[:, distance] 

233 else: 

234 raise InputError("marker", explanation=f"Wrong marker plan direction {direction}") 

235 

236 elif len(im.shape) == 3: 

237 if direction == 0: 

238 marker[distance, :, :] = im[distance, :, :] 

239 elif direction == 1: 

240 marker[:, distance, :] = im[:, distance, :] 

241 elif direction == 2: 

242 marker[:, :, distance] = im[:, :, distance] 

243 else: 

244 raise InputError("marker", explanation=f"Wrong marker plan direction {direction}") 

245 

246 else: 

247 raise InputError("marker", explanation=f"Image dimension should be 2 or 3, not {len(im.shape)}") 

248 

249 if verbose: 

250 print(f"binaryGeodesicReconstruction: marker -> set plan in direction {direction} at distance {distance}") 

251 

252 elif isinstance(marker, numpy.ndarray): 

253 # direct input of the marker 

254 if im.shape != marker.shape: 

255 raise InputError("marker", explanation="im and marker must have same shape") 

256 else: 

257 raise InputError("marker", explanation="must be a numpy array or a list") 

258 

259 # Step 2: first dilation and intersection 

260 r1 = binaryDilation(marker) & im 

261 r2 = r1 

262 r1 = binaryDilation(r2) & im 

263 d = 1 

264 dmax = numpy.inf if dmax is None else dmax 

265 if verbose: 

266 print(f'binaryGeodesicReconstruction: geodesic distance = {d} (sum = {r1.sum()} & {r2.sum()})') 

267 

268 # binary dilation until: 

269 # geodesic distance reach dmax 

270 # reconstuction complete 

271 while not numpy.array_equal(r1, r2) and d < dmax: 

272 r2 = r1 

273 r1 = binaryDilation(r2) & im 

274 d += 1 

275 if verbose: 

276 print(f'binaryGeodesicReconstruction: geodesic distance = {d} (sum = {r1.sum()} & {r2.sum()})') 

277 

278 return r1 # send the reconstructed image 

279 

280 

281def directionalErosion(bwIm, vect, a, c, nProcesses=nProcessesDefault, verbose = False): 

282 """ 

283 This functions performs direction erosion over the binarized image using 

284 an ellipsoidal structuring element over a range of directions. It is highly 

285 recommended that the structuring element is slightly smaller than the  

286 expected particle (50% smaller in each axis is a fair guess) 

287 

288 Parameters 

289 ----------- 

290 bwIm : 3D numpy array 

291 Binarized image to perform the erosion 

292 

293 vect : list of n elements, each element correspond to a 1X3 array of floats 

294 List of directional vectors for the structuring element 

295 

296 a : int or float 

297 Length of the secondary semi-axis of the structuring element in px 

298 

299 c : int or float 

300 Lenght of the principal semi-axis of the structuring element in px 

301 

302 nProcesses : integer (optional, default = nProcessesDefault) 

303 Number of processes for multiprocessing 

304 Default = number of CPUs in the system 

305 

306 verbose : boolean, optional (Default = False) 

307 True for printing the evolution of the process 

308 False for not printing the evolution of process 

309 

310 Returns 

311 -------- 

312 imEroded : 3D boolean array 

313 Booean array with the result of the erosion 

314 

315 Note 

316 ----- 

317 Taken from https://sbrisard.github.io/posts/20150930-orientation_correlations_among_rice_grains-06.html 

318 

319 """ 

320 

321 

322 #Check if the directional vector input is a list 

323 if isinstance(vect,list) == False: 

324 print("spam.contacts.directionalErosion: The directional vector must be a list") 

325 return 

326 

327 numberOfJobs = len( vect ) 

328 imEroded = numpy.zeros(bwIm.shape) 

329 

330 # Function for directionalErosion 

331 global funDirectionalErosion 

332 def funDirectionalErosion(job): 

333 maxDim = numpy.max([a,c]) 

334 spheroid = spam.kalisphera.makeBlurryNoisySpheroid([maxDim,maxDim,maxDim], 

335 [numpy.floor(maxDim/2), numpy.floor(maxDim/2), numpy.floor(maxDim/2)], 

336 [a,c], 

337 vect[job], 

338 background=0, 

339 foreground=1) 

340 imEroded_i = scipy.ndimage.binary_erosion(bwIm, structure = spheroid) 

341 return imEroded_i 

342 

343 if verbose: 

344 widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar.AdaptiveETA()] 

345 pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfJobs) 

346 pbar.start() 

347 finishedNodes = 0 

348 

349 # Run multiprocessing 

350 with multiprocessing.Pool(processes=nProcesses) as pool: 

351 for returns in pool.imap_unordered(funDirectionalErosion, range(0, numberOfJobs)): 

352 if verbose: 

353 finishedNodes += 1 

354 widgets[0] = progressbar.FormatLabel("{}/{} ".format(finishedNodes, numberOfJobs)) 

355 pbar.update(finishedNodes) 

356 imEroded = imEroded + returns 

357 pool.close() 

358 pool.join() 

359 

360 if verbose: 

361 pbar.finish() 

362 

363 return imEroded 

364 

365 

366def morphologicalReconstruction(im, selem = skimage.morphology.ball(1)): 

367 """ 

368 This functions performs a morphological reconstruction (greyscale opening followed by greyscale closing). The ouput image presents less variability in the greyvalues inside each phase, without modifying the original 

369 shape of the objects of the image. 

370 - 

371 

372 Parameters 

373 ----------- 

374 im : 3D numpy array 

375 Greyscale image to perform the reconstuction 

376 

377 selem : structuring element, optional 

378 Structuring element 

379 Default = None 

380 

381 Returns 

382 -------- 

383 imReconstructed : 3D boolean array 

384 Greyscale image after the reconstuction 

385 

386 """ 

387 

388 # Perform the opening 

389 imOpen = scipy.ndimage.grey_opening(im, footprint=selem) 

390 # Perform the closing 

391 imReconstructed = (scipy.ndimage.grey_closing(imOpen, footprint=selem)).astype(numpy.float32) 

392 

393 return imReconstructed