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

1# Library of SPAM morphological functions. 

2# Copyright (C) 2020 SPAM Contributors 

3# 

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

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

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

7# any later version. 

8# 

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

10# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

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

12# more details. 

13# 

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

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

16 

17import multiprocessing 

18 

19import numpy 

20import spam.helpers 

21 

22try: 

23 multiprocessing.set_start_method("fork") 

24except RuntimeError: 

25 pass 

26 

27import progressbar 

28import scipy.ndimage 

29import skimage.morphology 

30import spam.label 

31 

32# operations on greyscale images 

33 

34# Global number of processes 

35nProcessesDefault = multiprocessing.cpu_count() 

36 

37 

38def greyDilation(im, nBytes=1): 

39 """ 

40 This function applies a dilation on a grey scale image 

41 

42 Parameters 

43 ----------- 

44 im: numpy array 

45 The input image (greyscale) 

46 nBytes: int, default=1 

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

48 

49 Returns 

50 -------- 

51 numpy array 

52 The dilated image 

53 """ 

54 # Step 1: check type and dimension 

55 dim = len(im.shape) 

56 # Step 2: Determine substitution value 

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

58 # Step 3: apply dilation # x y z 

59 outputIm = im # +0 0 0 

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

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

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

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

64 if dim == 3: 

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

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

67 return outputIm 

68 

69 

70def greyErosion(im, nBytes=1): 

71 """ 

72 This function applies a erosion on a grey scale image 

73 

74 Parameters 

75 ----------- 

76 im: numpy array 

77 The input image (greyscale) 

78 nBytes: int, default=1 

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

80 

81 Returns 

82 -------- 

83 numpy array 

84 The eroded image 

85 """ 

86 # Step 1: check type and dimension 

87 dim = len(im.shape) 

88 # Step 2: Determine substitution value 

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

90 # Step 1: apply erosion # x y z 

91 outputIm = im # 0 0 0 

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

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

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

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

96 if dim == 3: 

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

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

99 return outputIm 

100 

101 

102def greyMorphologicalGradient(im, nBytes=1): 

103 """ 

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

105 

106 Parameters 

107 ----------- 

108 im: numpy array 

109 The input image (greyscale) 

110 nBytes: int, default=1 

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

112 

113 Returns 

114 -------- 

115 numpy array 

116 The morphologycal gradient of the image 

117 """ 

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

119 

120 

121# operation on binary images 

122 

123 

124def binaryDilation(im, sub=False): 

125 """ 

126 This function applies a dilation on a binary scale image 

127 

128 Parameters 

129 ----------- 

130 im: numpy array 

131 The input image (greyscale) 

132 sub: bool, default=False 

133 Subtitute value. 

134 

135 Returns 

136 -------- 

137 numpy array 

138 The dilated image 

139 """ 

140 # Step 0: import as bool 

141 im = im.astype(bool) 

142 # Step 1: check type and dimension 

143 dim = len(im.shape) 

144 # Step 1: apply dilation # x y z 

145 outputIm = im # 0 0 0 

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

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

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

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

150 if dim == 3: 

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

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

153 return numpy.array(outputIm).astype("<u1") 

154 

155 

156def binaryErosion(im, sub=False): 

157 """ 

158 This function applies a erosion on a binary scale image 

159 

160 Parameters 

161 ----------- 

162 im: numpy array 

163 The input image (greyscale) 

164 sub: bool, default=False 

165 Substitute value. 

166 

167 Returns 

168 -------- 

169 numpy array 

170 The eroded image 

171 """ 

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

173 return numpy.logical_not(binaryDilation(numpy.logical_not(im), sub=sub)).astype("<u1") 

174 

175 

176def binaryMorphologicalGradient(im, sub=False): 

177 """ 

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

179 

180 Parameters 

181 ----------- 

182 im: numpy array 

183 The input image (greyscale) 

184 nBytes: int, default=False 

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

186 

187 Returns 

188 -------- 

189 numpy array 

190 The morphologycal gradient of the image 

191 """ 

192 return (numpy.logical_xor(binaryDilation(im, sub=sub), im)).astype("<u1") 

193 

194 

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

196 """ 

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

198 

199 Parameters 

200 ----------- 

201 im: numpy.array 

202 The input binary image 

203 

204 marker: numpy.array or list 

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

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

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

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

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

210 

211 dmax: int, default=None 

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

213 

214 verbose: bool, default=False 

215 Verbose mode 

216 

217 Returns 

218 -------- 

219 numpy.array 

220 The reconstructed image 

221 """ 

222 from spam.errors import InputError 

223 

224 # Step 1: Define marker 

225 if isinstance(marker, list): 

226 # marker based on list of plans 

227 if len(marker) % 2: 

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

229 

230 plans = marker[:] 

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

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

233 direction = plans[2 * i] 

234 distance = plans[2 * i + 1] 

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

236 if direction == 0: 

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

238 elif direction == 1: 

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

240 else: 

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

242 

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

244 if direction == 0: 

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

246 elif direction == 1: 

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

248 elif direction == 2: 

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

250 else: 

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

252 

253 else: 

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

255 

256 if verbose: 

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

258 

259 elif isinstance(marker, numpy.ndarray): 

260 # direct input of the marker 

261 if im.shape != marker.shape: 

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

263 else: 

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

265 

266 # Step 2: first dilation and intersection 

267 r1 = binaryDilation(marker) & im 

268 r2 = r1 

269 r1 = binaryDilation(r2) & im 

270 d = 1 

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

272 if verbose: 

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

274 

275 # binary dilation until: 

276 # geodesic distance reach dmax 

277 # reconstuction complete 

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

279 r2 = r1 

280 r1 = binaryDilation(r2) & im 

281 d += 1 

282 if verbose: 

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

284 

285 return r1 # send the reconstructed image 

286 

287 

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

289 """ 

290 This functions performs direction erosion over the binarized image using 

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

292 recommended that the structuring element is slightly smaller than the 

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

294 

295 Parameters 

296 ----------- 

297 bwIm : 3D numpy array 

298 Binarized image to perform the erosion 

299 

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

301 List of directional vectors for the structuring element 

302 

303 a : int or float 

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

305 

306 c : int or float 

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

308 

309 nProcesses : integer (optional, default = nProcessesDefault) 

310 Number of processes for multiprocessing 

311 Default = number of CPUs in the system 

312 

313 verbose : boolean, optional (Default = False) 

314 True for printing the evolution of the process 

315 False for not printing the evolution of process 

316 

317 Returns 

318 -------- 

319 imEroded : 3D boolean array 

320 Booean array with the result of the erosion 

321 

322 Note 

323 ----- 

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

325 

326 """ 

327 

328 # Check if the directional vector input is a list 

329 if isinstance(vect, list) is False: 

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

331 return 

332 

333 numberOfJobs = len(vect) 

334 imEroded = numpy.zeros(bwIm.shape) 

335 

336 # Function for directionalErosion 

337 global _multiprocessingDirectionalErosion 

338 

339 def _multiprocessingDirectionalErosion(job): 

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

341 spheroid = spam.kalisphera.makeBlurryNoisySpheroid( 

342 [maxDim, maxDim, maxDim], [numpy.floor(maxDim / 2), numpy.floor(maxDim / 2), numpy.floor(maxDim / 2)], [a, c], vect[job], background=0, foreground=1 

343 ) 

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

345 return imEroded_i 

346 

347 if verbose: 

348 widgets = [progressbar.FormatLabel(""), " ", progressbar.Bar(), " ", progressbar.AdaptiveETA()] 

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

350 pbar.start() 

351 finishedNodes = 0 

352 

353 # Run multiprocessing 

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

355 for returns in pool.imap_unordered(_multiprocessingDirectionalErosion, range(0, numberOfJobs)): 

356 if verbose: 

357 finishedNodes += 1 

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

359 pbar.update(finishedNodes) 

360 imEroded = imEroded + returns 

361 pool.close() 

362 pool.join() 

363 

364 if verbose: 

365 pbar.finish() 

366 

367 return imEroded 

368 

369 

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

371 """ 

372 This functions performs a morphological reconstruction (greyscale opening followed by greyscale closing). 

373 The ouput image presents less variability in the greyvalues inside each phase, without modifying the original 

374 shape of the objects of the image. 

375 - 

376 

377 Parameters 

378 ----------- 

379 im : 3D numpy array 

380 Greyscale image to perform the reconstuction 

381 

382 selem : structuring element, optional 

383 Structuring element 

384 Default = None 

385 

386 Returns 

387 -------- 

388 imReconstructed : 3D boolean array 

389 Greyscale image after the reconstuction 

390 

391 """ 

392 

393 # Perform the opening 

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

395 # Perform the closing 

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

397 

398 return imReconstructed