Coverage for /usr/local/lib/python3.10/site-packages/spam/label/label.py: 89%

1# Library of SPAM functions for dealing with labelled images

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.

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/>.

17import multiprocessing

19try:

20 multiprocessing.set_start_method("fork")

21except RuntimeError:

22 pass

24import random

26import matplotlib

27import matplotlib.pyplot as plt

28import numba

29import numpy

30import progressbar

31import scipy.ndimage

32import scipy.spatial

33import spam.DIC

34import spam.filters

35from spambind.label.labelToolkit import boundingBoxes as boundingBoxesCPP

36from spambind.label.labelToolkit import centresOfMass as centresOfMassCPP

37from spambind.label.labelToolkit import labelToFloat as labelToFloatCPP

38from spambind.label.labelToolkit import momentOfInertia as momentOfInertiaCPP

39from spambind.label.labelToolkit import relabel as relabelCPP

40from spambind.label.labelToolkit import setVoronoi as setVoronoiCPP

41from spambind.label.labelToolkit import tetPixelLabel as tetPixelLabelCPP

42from spambind.label.labelToolkit import volumes as volumesCPP

44# Define a random colourmap for showing labels

45# This is taken from https://gist.github.com/jgomezdans/402500

46randomCmapVals = numpy.random.rand(256, 3)

47randomCmapVals[0, :] = numpy.array([1.0, 1.0, 1.0])

48randomCmapVals[-1, :] = numpy.array([0.0, 0.0, 0.0])

49randomCmap = matplotlib.colors.ListedColormap(randomCmapVals)

50del randomCmapVals

53# If you change this, remember to change the typedef in tools/labelToolkit/labelToolkitC.hpp

54labelType = "<u4"

56# Global number of processes

57nProcessesDefault = multiprocessing.cpu_count()

60def boundingBoxes(lab):

61 """

62 Returns bounding boxes for labelled objects using fast C-code which runs a single time through lab

64 Parameters

65 ----------

66 lab : 3D array of integers

67 Labelled volume, with lab.max() labels

69 Returns

70 -------

71 boundingBoxes : lab.max()x6 array of ints

72 This array contains, for each label, 6 integers:

74 - Zmin, Zmax

75 - Ymin, Ymax

76 - Xmin, Xmax

78 Note

79 ----

80 Bounding boxes `are not slices` and so to extract the correct bounding box from a numpy array you should use:

81 lab[ Zmin:Zmax+1, Ymin:Ymax+1, Xmin:Xmax+1 ]

82 Otherwise said, the bounding box of a single-voxel object at 1,1,1 will be:

83 1,1,1,1,1,1

85 Also note: for labelled images where some labels are missing, the bounding box returned for this case will be obviously wrong: `e.g.`, Zmin = (z dimension-1) and Zmax = 0

87 """

88 # Catch 2D image, and pad

89 if lab.ndim == 2:

90 lab = lab[numpy.newaxis, ...]

92 lab = lab.astype(labelType)

94 boundingBoxes = numpy.zeros((lab.max() + 1, 6), dtype="<u2")

96 boundingBoxesCPP(lab, boundingBoxes)

98 return boundingBoxes

100

101def centresOfMass(lab, boundingBoxes=None, minVol=None):

102 """

103 Calculates (binary) centres of mass of each label in labelled image

104

105 Parameters

106 ----------

107 lab : 3D array of integers

108 Labelled volume, with lab.max() labels

109

110 boundingBoxes : lab.max()x6 array of ints, optional

111 Bounding boxes in format returned by ``boundingBoxes``.

112 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

113

114 minVol : int, optional

115 The minimum volume in vx to be treated, any object below this threshold is returned as 0

116

117 Returns

118 -------

119 centresOfMass : lab.max()x3 array of floats

120 This array contains, for each label, 3 floats, describing the centre of mass of each label in Z, Y, X order

121 """

122 if boundingBoxes is None:

123 boundingBoxes = spam.label.boundingBoxes(lab)

124 if minVol is None:

125 minVol = 0

126 # Catch 2D image, and pad

127 if lab.ndim == 2:

128 lab = lab[numpy.newaxis, ...]

129

130 lab = lab.astype(labelType)

131

132 centresOfMass = numpy.zeros((lab.max() + 1, 3), dtype="<f4")

133

134 centresOfMassCPP(lab, boundingBoxes, centresOfMass, minVol)

135

136 return centresOfMass

137

138

139def volumes(lab, boundingBoxes=None):

140 """

141 Calculates (binary) volumes each label in labelled image, using potentially slow numpy.where

142

143 Parameters

144 ----------

145 lab : 3D array of integers

146 Labelled volume, with lab.max() labels

147

148 boundingBoxes : lab.max()x6 array of ints, optional

149 Bounding boxes in format returned by ``boundingBoxes``.

150 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

151

152 Returns

153 -------

154 volumes : lab.max()x1 array of ints

155 This array contains the volume in voxels of each label

156 """

157 # print "label.toolkit.volumes(): Warning this is a crappy python implementation"

158

159 lab = lab.astype(labelType)

160

161 if boundingBoxes is None:

162 boundingBoxes = spam.label.boundingBoxes(lab)

163

164 volumes = numpy.zeros((lab.max() + 1), dtype="<u4")

165

166 volumesCPP(lab, boundingBoxes, volumes)

167

168 return volumes

169

170

171def equivalentRadii(lab, boundingBoxes=None, volumes=None):

172 """

173 Calculates (binary) equivalent sphere radii of each label in labelled image

174

175 Parameters

176 ----------

177 lab : 3D array of integers

178 Labelled volume, with lab.max() labels

179

180 boundingBoxes : lab.max()x6 array of ints, optional

181 Bounding boxes in format returned by ``boundingBoxes``.

182 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

183

184 volumes : lab.max()x1 array of ints

185 Vector contining volumes, if this is passed, the others are ignored

186

187 Returns

188 -------

189 equivRadii : lab.max()x1 array of floats

190 This array contains the equivalent sphere radius in pixels of each label

191 """

192

193 def vol2rad(volumes):

194 return ((3.0 * volumes) / (4.0 * numpy.pi)) ** (1.0 / 3.0)

195

196 # If we have volumes, just go for it

197 if volumes is not None:

198 return vol2rad(volumes)

199

200 # If we don't have bounding boxes, recalculate them

201 if boundingBoxes is None:

202 boundingBoxes = spam.label.boundingBoxes(lab)

203

204 return vol2rad(spam.label.volumes(lab, boundingBoxes=boundingBoxes))

205

206

207def momentOfInertia(lab, boundingBoxes=None, minVol=None, centresOfMass=None):

208 """

209 Calculates (binary) moments of inertia of each label in labelled image

210

211 Parameters

212 ----------

213 lab : 3D array of integers

214 Labelled volume, with lab.max() labels

215

216 boundingBoxes : lab.max()x6 array of ints, optional

217 Bounding boxes in format returned by ``boundingBoxes``.

218 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

219

220 centresOfMass : lab.max()x3 array of floats, optional

221 Centres of mass in format returned by ``centresOfMass``.

222 If not defined (Default = None), it is recomputed by running ``centresOfMass``

223

224 minVol : int, optional

225 The minimum volume in vx to be treated, any object below this threshold is returned as 0

226 Default = default for spam.label.centresOfMass

227

228 Returns

229 -------

230 eigenValues : lab.max()x3 array of floats

231 The values of the three eigenValues of the moment of inertia of each labelled shape

232

233 eigenVectors : lab.max()x9 array of floats

234 3 x Z,Y,X components of the three eigenValues in the order of the eigenValues

235 """

236 if boundingBoxes is None:

237 boundingBoxes = spam.label.boundingBoxes(lab)

238 if centresOfMass is None:

239 centresOfMass = spam.label.centresOfMass(lab, boundingBoxes=boundingBoxes, minVol=minVol)

240

241 lab = lab.astype(labelType)

242

243 eigenValues = numpy.zeros((lab.max() + 1, 3), dtype="<f4")

244 eigenVectors = numpy.zeros((lab.max() + 1, 9), dtype="<f4")

245

246 momentOfInertiaCPP(lab, boundingBoxes, centresOfMass, eigenValues, eigenVectors)

247

248 return [eigenValues, eigenVectors]

249

250

251def ellipseAxes(lab, volumes=None, MOIeigenValues=None, enforceVolume=True, twoD=False):

252 """

253 Calculates length of half-axes a,b,c of the ellipitic fit of the particle.

254 These are half-axes and so are comparable to the radius -- and not the diameter -- of the particle.

255

256 See appendix of for inital work:

257 "Three-dimensional study on the interconnection and shape of crystals in a graphic granite by X-ray CT and image analysis.", Ikeda, S., Nakano, T., & Nakashima, Y. (2000).

258

259 Parameters

260 ----------

261 lab : 3D array of integers

262 Labelled volume, with lab.max() labels

263 Note: This is not strictly necessary if volumes and MOI is given

264

265 volumes : 1D array of particle volumes (optional, default = None)

266 Volumes of particles (length of array = lab.max())

267

268 MOIeigenValues : lab.max()x3 array of floats, (optional, default = None)

269 Bounding boxes in format returned by ``boundingBoxes``.

270 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

271

272 enforceVolume = bool (default = True)

273 Should a, b and c be scaled to enforce the fitted ellipse volume to be

274 the same as the particle?

275 This causes eigenValues are no longer completely consistent with fitted ellipse

276

277 twoD : bool (default = False)

278 Are these in fact 2D ellipses?

279 Not implemented!!

280

281 Returns

282 -------

283 ABCaxes : lab.max()x3 array of floats

284 a, b, c lengths of particle in pixels

285

286 Note

287 -----

288 Our elliptic fit is not necessarily of the same volume as the original particle,

289 although by default we scale all axes linearly with `enforceVolumes` to enforce this condition.

290

291 Reminder: volume of an ellipse is (4/3)*pi*a*b*c

292

293 Useful check from TM: Ia = (4/15)*pi*a*b*c*(b**2+c**2)

294

295 Function contributed by Takashi Matsushima (University of Tsukuba)

296 """

297 # Full ref:

298 # @misc{ikeda2000three,

299 # title={Three-dimensional study on the interconnection and shape of crystals in a graphic granite by X-ray CT and image analysis},

300 # author={Ikeda, S and Nakano, T and Nakashima, Y},

301 # year={2000},

302 # publisher={De Gruyter}

303 # }

304

305 if volumes is None:

306 volumes = spam.label.volumes(lab)

307 if MOIeigenValues is None:

308 MOIeigenValues = spam.label.momentOfInertia(lab)[0]

309

310 ABCaxes = numpy.zeros((volumes.shape[0], 3))

311

312 Ia = MOIeigenValues[:, 0]

313 Ib = MOIeigenValues[:, 1]

314 Ic = MOIeigenValues[:, 2]

315

316 # Initial derivation -- has quite a different volume from the original particle

317 # Use the particle's V. This is a source of inconsistency,

318 # since the condition V = (4/3) * pi * a * b * c is not necessarily respected

319 # ABCaxes[:,2] = numpy.sqrt( numpy.multiply((5.0/(2.0*volumes.ravel())),( Ib + Ia - Ic ) ) )

320 # ABCaxes[:,1] = numpy.sqrt( numpy.multiply((5.0/(2.0*volumes.ravel())),( Ia + Ic - Ib ) ) )

321 # ABCaxes[:,0] = numpy.sqrt( numpy.multiply((5.0/(2.0*volumes.ravel())),( Ic + Ib - Ia ) ) )

322

323 mask = numpy.logical_and(Ia != 0, numpy.isfinite(Ia))

324 # Calculate a, b and c: TM calculation 2018-03-30

325 # 2018-04-30 EA and MW: swap A and C so that A is the biggest

326 ABCaxes[mask, 2] = ((15.0 / (8.0 * numpy.pi)) * numpy.square(Ib[mask] + Ic[mask] - Ia[mask]) / numpy.sqrt((Ia[mask] - Ib[mask] + Ic[mask]) * (Ia[mask] + Ib[mask] - Ic[mask]))) ** (1.0 / 5.0)

327 ABCaxes[mask, 1] = ((15.0 / (8.0 * numpy.pi)) * numpy.square(Ic[mask] + Ia[mask] - Ib[mask]) / numpy.sqrt((Ib[mask] - Ic[mask] + Ia[mask]) * (Ib[mask] + Ic[mask] - Ia[mask]))) ** (1.0 / 5.0)

328 ABCaxes[mask, 0] = ((15.0 / (8.0 * numpy.pi)) * numpy.square(Ia[mask] + Ib[mask] - Ic[mask]) / numpy.sqrt((Ic[mask] - Ia[mask] + Ib[mask]) * (Ic[mask] + Ia[mask] - Ib[mask]))) ** (1.0 / 5.0)

329

330 if enforceVolume:

331 # Compute volume of ellipse:

332 ellipseVol = (4.0 / 3.0) * numpy.pi * ABCaxes[:, 0] * ABCaxes[:, 1] * ABCaxes[:, 2]

333 # filter zeros and infs

334 # print volumes.shape

335 # print ellipseVol.shape

336 volRatio = (volumes[mask] / ellipseVol[mask]) ** (1.0 / 3.0)

337 # print volRatio

338 ABCaxes[mask, 0] = ABCaxes[mask, 0] * volRatio

339 ABCaxes[mask, 1] = ABCaxes[mask, 1] * volRatio

340 ABCaxes[mask, 2] = ABCaxes[mask, 2] * volRatio

341

342 return ABCaxes

343

344

345def convertLabelToFloat(lab, vector):

346 """

347 Replaces all values of a labelled array with a given value.

348 Useful for visualising properties attached to labels, `e.g.`, sand grain displacements.

349

350 Parameters

351 ----------

352 lab : 3D array of integers

353 Labelled volume, with lab.max() labels

354

355 vector : a lab.max()x1 vector with values to replace each label with

356

357 Returns

358 -------

359 relabelled : 3D array of converted floats

360 """

361 lab = lab.astype(labelType)

362

363 relabelled = numpy.zeros_like(lab, dtype="<f4")

364

365 vector = vector.ravel().astype("<f4")

366

367 labelToFloatCPP(lab, vector, relabelled)

368

369 return relabelled

370

371

372def makeLabelsSequential(lab):

373 """

374 This function fills gaps in labelled images,

375 by relabelling them to be sequential integers.

376 Don't forget to recompute all your grain properties since the label numbers will change

377

378 Parameters

379 -----------

380 lab : 3D numpy array of ints ( of type spam.label.toolkit.labelType)

381 An array of labels with 0 as the background

382

383 Returns

384 --------

385 lab : 3D numpy array of ints ( of type spam.label.toolkit.labelType)

386 An array of labels with 0 as the background

387 """

388 maxLabel = int(lab.max())

389 lab = lab.astype(labelType)

390

391 uniqueLabels = numpy.unique(lab)

392 # print uniqueLabels

393

394 relabelMap = numpy.zeros((maxLabel + 1), dtype=labelType)

395 relabelMap[uniqueLabels] = range(len(uniqueLabels))

396

397 relabelCPP(lab, relabelMap)

398

399 return lab

400

401

402def getLabel(

403 labelledVolume,

404 label,

405 boundingBoxes=None,

406 centresOfMass=None,

407 margin=0,

408 extractCube=False,

409 extractCubeSize=None,

410 maskOtherLabels=True,

411 labelDilate=0,

412 labelDilateMaskOtherLabels=False,

413):

414 """

415 Helper function to extract labels from a labelled image/volume.

416 A dictionary is returned with the a subvolume around the particle.

417 Passing boundingBoxes and centresOfMass is highly recommended.

418

419 Parameters

420 ----------

421 labelVolume : 3D array of ints

422 3D Labelled volume

423

424 label : int

425 Label that we want information about

426

427 boundingBoxes : nLabels*2 array of ints, optional

428 Bounding boxes as returned by ``boundingBoxes``.

429 Optional but highly recommended.

430 If unset, bounding boxes are recalculated for every call.

431

432 centresOfMass : nLabels*3 array of floats, optional

433 Centres of mass as returned by ``centresOfMass``.

434 Optional but highly recommended.

435 If unset, centres of mass are recalculated for every call.

436

437 extractCube : bool, optional

438 Return label subvolume in the middle of a cube?

439 Default = False

440

441 extractCubeSize : int, optional

442 half-size of cube to extract.

443 Default = calculate minimum cube

444

445 margin : int, optional

446 Extract a ``margin`` pixel margin around bounding box or cube.

447 Default = 0

448

449 maskOtherLabels : bool, optional

450 In the returned subvolume, should other labels be masked?

451 If true, the mask is directly returned.

452 Default = True

453

454 labelDilate : int, optional

455 Number of times label should be dilated before returning it?

456 This can be useful for catching the outside/edge of an image.

457 ``margin`` should at least be equal to this value.

458 Requires ``maskOtherLabels``.

459 Default = 0

460

461 labelDilateMaskOtherLabels : bool, optional

462 Strictly cut the other labels out of the dilated image of the requested label?

463 Only pertinent for positive labelDilate values.

464 Default = False

465

466

467 Returns

468 -------

469 Dictionary containing:

470

471 Keys:

472 subvol : 3D array of bools or ints

473 subvolume from labelled image

474

475 slice : tuple of 3*slices

476 Slice used to extract subvol for the bounding box mode

477

478 sliceCube : tuple of 3*slices

479 Slice used to extract subvol for the cube mode, warning,

480 if the label is near the edge, this is the slice up to the edge,

481 and so it will be smaller than the returned cube

482

483 boundingBox : 1D numpy array of 6 ints

484 Bounding box, including margin, in bounding box mode. Contains:

485 [Zmin, Zmax, Ymin, Ymax, Xmin, Xmax]

486 Note: uses the same convention as spam.label.boundingBoxes, so

487 if you want to use this to extract your subvolume, add +1 to max

488

489 boundingBoxCube : 1D numpy array of 6 ints

490 Bounding box, including margin, in cube mode. Contains:

491 [Zmin, Zmax, Ymin, Ymax, Xmin, Xmax]

492 Note: uses the same convention as spam.label.boundingBoxes, so

493 if you want to use this to extract your subvolume, add +1 to max

494

495 centreOfMassABS : 3*float

496 Centre of mass with respect to ``labelVolume``

497

498 centreOfMassREL : 3*float

499 Centre of mass with respect to ``subvol``

500

501 volumeInitial : int

502 Volume of label (before dilating)

503

504 volumeDilated : int

505 Volume of label (after dilating, if requested)

506

507 """

508 import spam.mesh

509

510 if boundingBoxes is None:

511 print("\tlabel.toolkit.getLabel(): Bounding boxes not passed.")

512 print("\tThey will be recalculated for each label, highly recommend calculating outside this function")

513 boundingBoxes = spam.label.boundingBoxes(labelledVolume)

514

515 if centresOfMass is None:

516 print("\tlabel.toolkit.getLabel(): Centres of mass not passed.")

517 print("\tThey will be recalculated for each label, highly recommend calculating outside this function")

518 centresOfMass = spam.label.centresOfMass(labelledVolume)

519

520 # Check if there is a bounding box for this label:

521 if label >= boundingBoxes.shape[0]:

522 return

523 raise "No bounding boxes for this grain"

524

525 bbo = boundingBoxes[label]

526 com = centresOfMass[label]

527 comRound = numpy.floor(centresOfMass[label])

528

529 # 1. Check if boundingBoxes are correct:

530 if (bbo[0] == labelledVolume.shape[0] - 1) and (bbo[1] == 0) and (bbo[2] == labelledVolume.shape[1] - 1) and (bbo[3] == 0) and (bbo[4] == labelledVolume.shape[2] - 1) and (bbo[5] == 0):

531 pass

532 # print("\tlabel.toolkit.getLabel(): Label {} does not exist".format(label))

533

534 else:

535 # Define output dictionary since we'll add different things to it

536 output = {}

537 output["centreOfMassABS"] = com

538

539 # We have a bounding box, let's extract it.

540 if extractCube:

541 # Calculate offsets between centre of mass and bounding box

542 offsetTop = numpy.ceil(com - bbo[0::2])

543 offsetBot = numpy.ceil(com - bbo[0::2])

544 offset = numpy.max(numpy.hstack([offsetTop, offsetBot]))

545

546 # If is none, assume closest fitting cube.

547 if extractCubeSize is not None:

548 if extractCubeSize < offset:

549 print("\tlabel.toolkit.getLabel(): size of desired cube is smaller than minimum to contain label. Continuing anyway.")

550 offset = int(extractCubeSize)

551

552 # if a margin is set, add it to offset

553 # if margin is not None:

554 offset += margin

555

556 offset = int(offset)

557

558 # we may go outside the volume. Let's check this

559 labSubVol = numpy.zeros(3 * [2 * offset + 1])

560

561 topOfSlice = numpy.array(

562 [

563 int(comRound[0] - offset),

564 int(comRound[1] - offset),

565 int(comRound[2] - offset),

566 ]

567 )

568 botOfSlice = numpy.array(

569 [

570 int(comRound[0] + offset + 1),

571 int(comRound[1] + offset + 1),

572 int(comRound[2] + offset + 1),

573 ]

574 )

575

576 labSubVol = spam.helpers.slicePadded(

577 labelledVolume,

578 [

579 topOfSlice[0],

580 botOfSlice[0],

581 topOfSlice[1],

582 botOfSlice[1],

583 topOfSlice[2],

584 botOfSlice[2],

585 ],

586 )

587

588 output["sliceCube"] = (

589 slice(topOfSlice[0], botOfSlice[0]),

590 slice(topOfSlice[1], botOfSlice[1]),

591 slice(topOfSlice[2], botOfSlice[2]),

592 )

593

594 output["boundingBoxCube"] = numpy.array(

595 [

596 topOfSlice[0],

597 botOfSlice[0] - 1,

598 topOfSlice[1],

599 botOfSlice[1] - 1,

600 topOfSlice[2],

601 botOfSlice[2] - 1,

602 ]

603 )

604

605 output["centreOfMassREL"] = com - topOfSlice

606

607 # We have a bounding box, let's extract it.

608 else:

609 topOfSlice = numpy.array([int(bbo[0]) - margin, int(bbo[2]) - margin, int(bbo[4]) - margin])

610 botOfSlice = numpy.array(

611 [

612 int(bbo[1] + margin + 1),

613 int(bbo[3] + margin + 1),

614 int(bbo[5] + margin + 1),

615 ]

616 )

617

618 labSubVol = spam.helpers.slicePadded(

619 labelledVolume,

620 [

621 topOfSlice[0],

622 botOfSlice[0],

623 topOfSlice[1],

624 botOfSlice[1],

625 topOfSlice[2],

626 botOfSlice[2],

627 ],

628 )

629

630 output["slice"] = (

631 slice(topOfSlice[0], botOfSlice[0]),

632 slice(topOfSlice[1], botOfSlice[1]),

633 slice(topOfSlice[2], botOfSlice[2]),

634 )

635

636 output["boundingBox"] = numpy.array(

637 [

638 topOfSlice[0],

639 botOfSlice[0] - 1,

640 topOfSlice[1],

641 botOfSlice[1] - 1,

642 topOfSlice[2],

643 botOfSlice[2] - 1,

644 ]

645 )

646

647 output["centreOfMassREL"] = com - topOfSlice

648

649 # Get mask for this label

650 maskLab = labSubVol == label

651 volume = numpy.sum(maskLab)

652 output["volumeInitial"] = volume

653

654 # if we should mask, just return the mask.

655 if maskOtherLabels:

656 # 2019-09-07 EA: changing dilation/erosion into a single pass by a spherical element, rather than repeated

657 # iterations of the standard.

658 if labelDilate > 0:

659 if labelDilate >= margin:

660 print("\tlabel.toolkit.getLabel(): labelDilate requested with a margin smaller than or equal to the number of times to dilate. I hope you know what you're doing!")

661 strucuringElement = spam.mesh.structuringElement(radius=labelDilate, order=2, dim=3)

662 maskLab = scipy.ndimage.binary_dilation(maskLab, structure=strucuringElement, iterations=1)

663 if labelDilateMaskOtherLabels:

664 # remove voxels that are neither our label nor pore

665 maskLab[numpy.logical_and(labSubVol != label, labSubVol != 0)] = 0

666 if labelDilate < 0:

667 strucuringElement = spam.mesh.structuringElement(radius=-1 * labelDilate, order=2, dim=3)

668 maskLab = scipy.ndimage.binary_erosion(maskLab, structure=strucuringElement, iterations=1)

669

670 # Just overwrite "labSubVol"

671 labSubVol = maskLab

672 # Update volume output

673 output["volumeDilated"] = labSubVol.sum()

674

675 output["subvol"] = labSubVol

676

677 return output

678

679

680def getImagettesLabelled(

681 lab1,

682 label,

683 Phi,

684 im1,

685 im2,

686 searchRange,

687 boundingBoxes,

688 centresOfMass,

689 margin=0,

690 labelDilate=0,

691 maskOtherLabels=True,

692 applyF="all",

693 volumeThreshold=100,

694):

695 """

696 This function is responsible for extracting correlation windows ("imagettes") from two larger images (im1 and im2) with the help of a labelled im1.

697 This is generally to do image correlation, this function will be used for spam-ddic and pixelSearch modes.

698

699 Parameters

700 ----------

701 lab1 : 3D numpy array of ints

702 Labelled image containing nLabels

703

704 label : int

705 Label of interest

706

707 Phi : 4x4 numpy array of floats

708 Phi matrix representing the movement of imagette1,

709 if not equal to `I`, imagette1 is deformed by the non-translation parts of Phi (F)

710 and the displacement is added to the search range (see below)

711

712 im1 : 3D numpy array

713 This is the large input reference image of greyvalues

714

715 im2 : 3D numpy array

716 This is the large input deformed image of greyvalues

717

718 searchRange : 6-component numpy array of ints

719 This defines where imagette2 should be extracted with respect to imagette1's position in im1.

720 The 6 components correspond to [ Zbot Ztop Ybot Ytop Xbot Xtop ].

721 If Z, Y and X values are the same, then imagette2 will be displaced and the same size as imagette1.

722 If 'bot' is lower than 'top', imagette2 will be larger in that dimension

723

724 boundingBoxes : nLabels*2 array of ints

725 Bounding boxes as returned by ``boundingBoxes``

726

727 centresOfMass : nLabels*3 array of floats

728 Centres of mass as returned by ``centresOfMass``

729

730 margin : int, optional

731 Margin around the grain to extract in pixels

732 Default = 0

733

734 labelDilate : int, optional

735 How much to dilate the label before computing the mask?

736 Default = 0

737

738 maskOtherLabels : bool, optional

739 In the returned subvolume, should other labels be masked?

740 If true, the mask is directly returned.

741 Default = True

742

743 applyF : string, optional

744 If a non-identity Phi is passed, should the F be applied to the returned imagette1?

745 Options are: 'all', 'rigid', 'no'

746 Default = 'all'

747 Note: as of January 2021, it seems to make more sense to have this as 'all' for pixelSearch, and 'no' for local DIC

748

749 volumeThreshold : int, optional

750 Pixel volume of labels that are discarded

751 Default = 100

752

753 Returns

754 -------

755 Dictionary :

756

757 'imagette1' : 3D numpy array,

758

759 'imagette1mask': 3D numpy array of same size as imagette1 or None,

760

761 'imagette2': 3D numpy array, bigger or equal size to imagette1

762

763 'returnStatus': int,

764 Describes success in extracting imagette1 and imagette2.

765 If == 1 success, otherwise negative means failure.

766

767 'pixelSearchOffset': 3-component list of ints

768 Coordinates of the top of the pixelSearch range in im1, i.e., the displacement that needs to be

769 added to the raw pixelSearch output to make it a im1 -> im2 displacement

770 """

771 returnStatus = 1

772 imagette1 = None

773 imagette1mask = None

774 imagette2 = None

775

776 intDisplacement = numpy.round(Phi[0:3, 3]).astype(int)

777 PhiNoDisp = Phi.copy()

778 # PhiNoDisp[0:3,-1] -= intDisplacement

779 PhiNoDisp[0:3, -1] = numpy.zeros(3)

780 if applyF == "rigid":

781 PhiNoDisp = spam.deformation.computeRigidPhi(PhiNoDisp)

782

783 gottenLabel = spam.label.getLabel(

784 lab1,

785 label,

786 extractCube=False,

787 boundingBoxes=boundingBoxes,

788 centresOfMass=centresOfMass,

789 margin=labelDilate + margin,

790 maskOtherLabels=True,

791 labelDilate=labelDilate,

792 labelDilateMaskOtherLabels=maskOtherLabels,

793 )

794

795 # In case the label is missing or the Phi is duff

796 if gottenLabel is None or not numpy.all(numpy.isfinite(Phi)):

797 returnStatus = -7

798

799 else:

800 # Maskette 1 is either a boolean array if args.MASK

801 # otherwise it contains ints i.e., labels

802

803 # Use new padded slicer, to remain aligned with getLabel['subvol']

804 # + add 1 on the "max" side for bounding box -> slice

805 imagette1 = spam.helpers.slicePadded(im1, gottenLabel["boundingBox"] + numpy.array([0, 1, 0, 1, 0, 1]))

806

807 if applyF == "all" or applyF == "rigid":

808 imagette1 = spam.DIC.applyPhi(imagette1, PhiNoDisp, PhiCentre=gottenLabel["centreOfMassREL"])

809 imagette1mask = (

810 spam.DIC.applyPhi(

811 gottenLabel["subvol"] > 0,

812 PhiNoDisp,

813 PhiCentre=gottenLabel["centreOfMassREL"],

814 interpolationOrder=0,

815 )

816 > 0

817 )

818 elif applyF == "no":

819 imagette1mask = gottenLabel["subvol"]

820 else:

821 print("spam.label.getImagettesLabelled(): unknown option for applyF options are: ['all', 'rigid', 'no']")

822

823 maskette1vol = numpy.sum(imagette1mask)

824

825 if maskette1vol > volumeThreshold:

826 # 2020-09-25 OS and EA: Prepare startStop array for imagette 2 to be extracted with new

827 # slicePadded, this should solved "Boss: failed imDiff" and RS=-5 forever

828 startStopIm2 = [

829 int(gottenLabel["boundingBox"][0] - margin - max(labelDilate, 0) + searchRange[0] + intDisplacement[0]),

830 int(gottenLabel["boundingBox"][1] + margin + max(labelDilate, 0) + searchRange[1] + intDisplacement[0] + 1),

831 int(gottenLabel["boundingBox"][2] - margin - max(labelDilate, 0) + searchRange[2] + intDisplacement[1]),

832 int(gottenLabel["boundingBox"][3] + margin + max(labelDilate, 0) + searchRange[3] + intDisplacement[1] + 1),

833 int(gottenLabel["boundingBox"][4] - margin - max(labelDilate, 0) + searchRange[4] + intDisplacement[2]),

834 int(gottenLabel["boundingBox"][5] + margin + max(labelDilate, 0) + searchRange[5] + intDisplacement[2] + 1),

835 ]

836

837 imagette2 = spam.helpers.slicePadded(im2, startStopIm2)

838

839 # imagette2imagette1sizeDiff = numpy.array(imagette2.shape) - numpy.array(imagette1.shape)

840

841 # If all of imagette2 is nans it fell outside im2 (or in any case it's going to be difficult to correlate)

842 if numpy.all(numpy.isnan(imagette2)):

843 returnStatus = -5

844 else:

845 # Failed volume condition

846 returnStatus = -5

847

848 return {

849 "imagette1": imagette1,

850 "imagette1mask": imagette1mask,

851 "imagette2": imagette2,

852 "returnStatus": returnStatus,

853 "pixelSearchOffset": searchRange[0::2] - numpy.array([max(labelDilate, 0)] * 3) - margin + intDisplacement,

854 }

855

856

857def labelsOnEdges(lab):

858 """

859 Return labels on edges of volume

860

861 Parameters

862 ----------

863 lab : 3D numpy array of ints

864 Labelled volume

865

866 Returns

867 -------

868 uniqueLabels : list of ints

869 List of labels on edges

870 """

871

872 numpy.arange(lab.max() + 1)

873

874 uniqueLabels = []

875

876 uniqueLabels.append(numpy.unique(lab[:, :, 0]))

877 uniqueLabels.append(numpy.unique(lab[:, :, -1]))

878 uniqueLabels.append(numpy.unique(lab[:, 0, :]))

879 uniqueLabels.append(numpy.unique(lab[:, -1, :]))

880 uniqueLabels.append(numpy.unique(lab[0, :, :]))

881 uniqueLabels.append(numpy.unique(lab[-1, :, :]))

882

883 # Flatten list of lists:

884 # https://stackoverflow.com/questions/952914/making-a-flat-list-out-of-list-of-lists-in-python?utm_medium=organic&utm_source=google_rich_qa&utm_campaign=google_rich_qa

885 uniqueLabels = [item for sublist in uniqueLabels for item in sublist]

886

887 # There might well be labels that appears on multiple faces of the cube, remove them

888 uniqueLabels = numpy.unique(numpy.array(uniqueLabels))

889

890 return uniqueLabels.astype(labelType)

891

892

893def removeLabels(lab, listOfLabelsToRemove):

894 """

895 Resets a list of labels to zero in a labelled volume.

896

897 Parameters

898 ----------

899 lab : 3D numpy array of ints

900 Labelled volume

901

902 listOfLabelsToRemove : list-like of ints

903 Labels to remove

904

905 Returns

906 -------

907 lab : 3D numpy array of ints

908 Labelled volume with desired labels blanked

909

910 Note

911 ----

912 You might want to use `makeLabelsSequential` after using this function,

913 but don't forget to recompute all your grain properties since the label numbers will change

914 """

915 lab = lab.astype(labelType)

916

917 # define a vector with sequential ints

918 arrayOfLabels = numpy.arange(lab.max() + 1, dtype=labelType)

919

920 # Remove the ones that have been asked for

921 for label in listOfLabelsToRemove:

922 arrayOfLabels[label] = 0

923

924 relabelCPP(lab, arrayOfLabels)

925

926 return lab

927

928

929def setVoronoi(lab, poreEDT=None, maxPoreRadius=10):

930 """

931 This function computes an approximate set Voronoi for a given labelled image.

932 This is a voronoi which does not have straight edges, and which necessarily

933 passes through each contact point, so it is respectful of non-spherical grains.

934

935 See:

936 Schaller, F. M., Kapfer, S. C., Evans, M. E., Hoffmann, M. J., Aste, T., Saadatfar, M., ... & Schroder-Turk, G. E. (2013).

937 Set Voronoi diagrams of 3D assemblies of aspherical particles. Philosophical Magazine, 93(31-33), 3993-4017.

938 https://doi.org/10.1080/14786435.2013.834389

939

940 and

941

942 Weis, S., Schonhofer, P. W., Schaller, F. M., Schroter, M., & Schroder-Turk, G. E. (2017).

943 Pomelo, a tool for computing Generic Set Voronoi Diagrams of Aspherical Particles of Arbitrary Shape. In EPJ Web of Conferences (Vol. 140, p. 06007). EDP Sciences.

944

945 Parameters

946 -----------

947 lab: 3D numpy array of labelTypes

948 Labelled image

949

950 poreEDT: 3D numpy array of floats (optional, default = None)

951 Euclidean distance map of the pores.

952 If not given, it is computed by scipy.ndimage.distance_transform_edt

953

954 maxPoreRadius: int (optional, default = 10)

955 Maximum pore radius to be considered (this threshold is for speed optimisation)

956

957 Returns

958 --------

959 lab: 3D numpy array of labelTypes

960 Image labelled with set voronoi labels

961 """

962 if poreEDT is None:

963 # print( "\tlabel.toolkit.setVoronoi(): Calculating the Euclidean Distance Transform of the pore with" )

964 # print "\t\tscipy.ndimage.distance_transform_edt, this takes a lot of memory"

965 poreEDT = scipy.ndimage.distance_transform_edt(lab == 0).astype("<f4")

966

967 lab = lab.astype(labelType)

968 labOut = numpy.zeros_like(lab)

969 maxPoreRadius = int(maxPoreRadius)

970

971 # Prepare sorted distances in a cube to fit a maxPoreRadius.

972 # This precomutation saves a lot of time

973 # Local grid of values, centred at zero

974 gridD = numpy.mgrid[

975 -maxPoreRadius : maxPoreRadius + 1,

976 -maxPoreRadius : maxPoreRadius + 1,

977 -maxPoreRadius : maxPoreRadius + 1,

978 ]

979

980 # Compute distances from centre

981 Rarray = numpy.sqrt(numpy.square(gridD[0]) + numpy.square(gridD[1]) + numpy.square(gridD[2])).ravel()

982 sortedIndices = numpy.argsort(Rarray)

983

984 # Array to hold sorted points

985 coords = numpy.zeros((len(Rarray), 3), dtype="<i4")

986 # Fill in with Z, Y, X points in order of distance to centre

987 coords[:, 0] = gridD[0].ravel()[sortedIndices]

988 coords[:, 1] = gridD[1].ravel()[sortedIndices]

989 coords[:, 2] = gridD[2].ravel()[sortedIndices]

990 del gridD

991

992 # Now define a simple array (by building a list) that gives the linear

993 # entry point into coords at the nearest integer values

994 sortedDistances = Rarray[sortedIndices]

995 indices = []

996 n = 0

997 i = 0

998 while i <= maxPoreRadius + 1:

999 if sortedDistances[n] >= i:

1000 # indices.append( [ i, n ] )

1001 indices.append(n)

1002 i += 1

1003 n += 1

1004 indices = numpy.array(indices).astype("<i4")

1005

1006 # Call C++ code

1007 setVoronoiCPP(lab, poreEDT.astype("<f4"), labOut, coords, indices)

1008

1009 return labOut

1010

1011

1012def labelTetrahedra(dims, points, connectivity, nThreads=1):

1013 """

1014 Labels voxels corresponding to tetrahedra according to a connectivity matrix and node points

1015

1016 Parameters

1017 ----------

1018 dims: tuple representing z,y,x dimensions of the desired labelled output

1019

1020 points: number of points x 3 array of floats

1021 List of points that define the vertices of the tetrahedra in Z,Y,X format.

1022 These points are referred to by line number in the connectivity array

1023

1024 connectivity: number of tetrahedra x 4 array of integers

1025 Connectivity matrix between points that define tetrahedra.

1026 Each line defines a tetrahedron whose number is the line number.

1027 Each line contains 4 integers that indicate the 4 points in the nodePos array.

1028

1029 nThreads: int (optional, default=1)

1030 The number of threads used for the cpp parallelisation.

1031

1032 Returns

1033 -------

1034 3D array of ints, shape = dims

1035 Labelled 3D volume where voxels are numbered according to the tetrahedron number they fall inside of

1036 # WARNING: Voxels outside of the mesh get a value of #tetrahedra + 1

1037 """

1038 assert len(dims) == 3, "spam.label.labelTetrahedra(): dim is not length 3"

1039 assert points.shape[1] == 3, "spam.label.labelTetrahedra(): points doesn't have 3 colums"

1040 assert connectivity.shape[1] == 4, "spam.label.labelTetrahedra(): connectivity doesn't have 4 colums"

1041 assert points.shape[0] >= connectivity.max(), "spam.label.labelTetrahedra(): connectivity should not refer to points numbers biggest than the number of rows in points"

1042

1043 dims = numpy.array(dims).astype("<u2")

1044 # WARNING: here we set the background to be number of tetra + 1

1045 # bold choice but that's ok

1046 lab = numpy.ones(tuple(dims), dtype=labelType) * connectivity.shape[0] + 1

1047

1048 connectivity = connectivity.astype("<u4")

1049 points = points.astype("<f4")

1050

1051 tetPixelLabelCPP(lab, connectivity, points, nThreads)

1052

1053 return lab

1054

1055

1056def labelTetrahedraForScipyDelaunay(dims, delaunay):

1057 """

1058 Labels voxels corresponding to tetrahedra coming from scipy.spatial.Delaunay

1059 Apparently the cells are not well-numbered, which causes a number of zeros

1060 when using `labelledTetrahedra`

1061

1062 Parameters

1063 ----------

1064 dims: tuple

1065 represents z,y,x dimensions of the desired labelled output

1066

1067 delaunay: "delaunay" object

1068 Object returned by scipy.spatial.Delaunay( centres )

1069 Hint: If using label.toolkit.centresOfMass( ), do centres[1:] to remove

1070 the position of zero.

1071

1072 Returns

1073 -------

1074 lab: 3D array of ints, shape = dims

1075 Labelled 3D volume where voxels are numbered according to the tetrahedron number they fall inside of

1076 """

1077

1078 # Big matrix of points poisitions

1079 points = numpy.zeros((dims[0] * dims[1] * dims[2], 3))

1080

1081 mgrid = numpy.mgrid[0 : dims[0], 0 : dims[1], 0 : dims[2]]

1082 for i in [0, 1, 2]:

1083 points[:, i] = mgrid[i].ravel()

1084

1085 del mgrid

1086

1087 lab = numpy.ones(tuple(dims), dtype=labelType) * delaunay.nsimplex + 1

1088 lab = delaunay.find_simplex(points).reshape(dims)

1089

1090 return lab

1091

1092

1093@numba.njit(cache=True, parallel=True)

1094def labelTriangles(dims, points, connectivity):

1095 """

1096 Labels pixels corresponding to triangles according to a connectivity matrix and node points

1097

1098 Parameters

1099 ----------

1100 dims: tuple representing y,x dimensions of the desired labelled output

1101

1102 points: number of points x 2 array of floats

1103 List of points that define the vertices of the triangles in Y,X format.

1104 These points are referred to by line number in the connectivity array

1105

1106 connectivity: number of triangles x 3 array of integers

1107 Connectivity matrix between points that define triangles.

1108 Each line defines a triangle whose number is the line number.

1109 Each line contains 3 integers that indicate the 3 points in the nodePos array.

1110

1111 nThreads: int (optional, default=1)

1112 The number of threads used for the cpp parallelisation.

1113

1114 Returns

1115 -------

1116 2D array of ints, shape = dims

1117 Labelled 2D image where pixels are numbered according to the triangle number they fall inside of

1118 # WARNING: Voxels outside of the mesh get a value of #tetrahedra + 1

1119 """

1120 assert len(dims) == 2

1121 assert len(points.shape) == 2

1122 assert points.shape[1] == 2

1123 assert len(connectivity.shape) == 2

1124 assert connectivity.shape[1] == 3

1125

1126 lab = numpy.zeros(dims, dtype=numpy.uint16)

1127 # lab = numpy.zeros(dims)

1128 lab[:, :] = connectivity.shape[0]

1129

1130 for ntri in numba.prange(connectivity.shape[0]):

1131 # for ntri in range(connectivity.shape[0]):

1132 tri = connectivity[ntri]

1133 # print(f"{ntri=} {tri=}")

1134 # step 1: bounding box

1135 yxLocal = numpy.zeros((3, 2))

1136 yxLocalMin = numpy.zeros(2)

1137 yxLocalMin[:] = 65535

1138 yxLocalMax = numpy.zeros(2)

1139 yxLocalMax[:] = 0

1140 for nnode, node in enumerate(tri):

1141 yxLocal[nnode] = points[node]

1142 if yxLocalMin[0] > points[node][0]:

1143 yxLocalMin[0] = points[node][0]

1144 if yxLocalMin[1] > points[node][1]:

1145 yxLocalMin[1] = points[node][1]

1146 if yxLocalMax[0] < points[node][0]:

1147 yxLocalMax[0] = points[node][0]

1148 if yxLocalMax[1] < points[node][1]:

1149 yxLocalMax[1] = points[node][1]

1150 yxLocalMin = numpy.floor(yxLocalMin)

1151 yxLocalMax = numpy.ceil(yxLocalMax)

1152

1153 for y in range(int(yxLocalMin[0]), int(yxLocalMax[0])):

1154 for x in range(int(yxLocalMin[1]), int(yxLocalMax[1])):

1155 if _pointInTriangle(yxLocal[0], yxLocal[1], yxLocal[2], numpy.array([y, x])):

1156 lab[y, x] = ntri

1157

1158 return lab

1159

1160

1161@numba.njit(cache=True)

1162def _triangleSide(p1, p2, p):

1163 return (p2[0] - p1[0]) * (p[1] - p1[1]) + (-p2[1] + p1[1]) * (p[0] - p1[0])

1164

1165

1166@numba.njit(cache=True)

1167def _pointInTriangle(p1, p2, p3, p):

1168 b1 = _triangleSide(p1, p2, p) >= 0

1169 b2 = _triangleSide(p2, p3, p) >= 0

1170 b3 = _triangleSide(p3, p1, p) >= 0

1171 return b1 and b2 and b3

1172

1173

1174def filterIsolatedCells(array, struct, size):

1175 """

1176 Return array with completely isolated single cells removed

1177

1178 Parameters

1179 ----------

1180 array: 3-D (labelled or binary) array

1181 Array with completely isolated single cells

1182

1183 struct: 3-D binary array

1184 Structure array for generating unique regions

1185

1186 size: integer

1187 Size of the isolated cells to exclude

1188 (Number of Voxels)

1189

1190 Returns

1191 -------

1192 filteredArray: 3-D (labelled or binary) array

1193 Array with minimum region size > size

1194

1195 Notes

1196 -----

1197 function from: http://stackoverflow.com/questions/28274091/removing-completely-isolated-cells-from-python-array

1198 """

1199

1200 filteredArray = ((array > 0) * 1).astype("uint8")

1201 idRegions, numIDs = scipy.ndimage.label(filteredArray, structure=struct)

1202 idSizes = numpy.array(scipy.ndimage.sum(filteredArray, idRegions, range(numIDs + 1)))

1203 areaMask = idSizes <= size

1204 filteredArray[areaMask[idRegions]] = 0

1205

1206 filteredArray = ((filteredArray > 0) * 1).astype("uint8")

1207 array = filteredArray * array

1208

1209 return array

1210

1211

1212def trueSphericity(lab, boundingBoxes=None, centresOfMass=None, gaussianFilterSigma=0.75, minVol=256):

1213 """

1214 Calculates the degree of True Sphericity (psi) for all labels, as per:

1215 "Sphericity measures of sand grains" Rorato et al., Engineering Geology, 2019

1216 and originlly proposed in: "Volume, shape, and roundness of rock particles", Waddell, The Journal of Geology, 1932.

1217

1218 True Sphericity (psi) = Surface area of equivalent sphere / Actual surface area

1219

1220 The actual surface area is computed by extracting each particle with getLabel, a Gaussian smooth of 0.75 is applied

1221 and the marching cubes algorithm from skimage is used to mesh the surface and compute the surface area.

1222

1223 Parameters

1224 ----------

1225 lab : 3D array of integers

1226 Labelled volume, with lab.max() labels

1227

1228 boundingBoxes : lab.max()x6 array of ints, optional

1229 Bounding boxes in format returned by ``boundingBoxes``.

1230 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

1231

1232 centresOfMass : lab.max()x3 array of floats, optional

1233 Centres of mass in format returned by ``centresOfMass``.

1234 If not defined (Default = None), it is recomputed by running ``centresOfMass``

1235

1236 gaussianFilterSigma : float, optional

1237 Sigma of the Gaussian filter used to smooth the binarised shape

1238 Default = 0.75

1239

1240 minVol : int, optional

1241 The minimum volume in vx to be treated, any object below this threshold is returned as 0

1242 Default = 256 voxels

1243

1244 Returns

1245 -------

1246 trueSphericity : lab.max() array of floats

1247 The values of the degree of true sphericity for each particle

1248

1249 Notes

1250 -----

1251 Function contributed by Riccardo Rorato (UPC Barcelona)

1252

1253 Due to numerical errors, this value can be >1, it should be clipped at 1.0

1254 """

1255 import skimage.measure

1256

1257 lab = lab.astype(labelType)

1258

1259 if boundingBoxes is None:

1260 boundingBoxes = spam.label.boundingBoxes(lab)

1261 if centresOfMass is None:

1262 centresOfMass = spam.label.centresOfMass(lab, boundingBoxes=boundingBoxes, minVol=minVol)

1263

1264 trueSphericity = numpy.zeros((lab.max() + 1), dtype="<f4")

1265

1266 sphereSurfaceArea = 4.0 * numpy.pi * (equivalentRadii(lab, boundingBoxes=boundingBoxes) ** 2)

1267

1268 for label in range(1, lab.max() + 1):

1269 if not (centresOfMass[label] == numpy.array([0.0, 0.0, 0.0])).all():

1270 # Extract grain

1271 GL = spam.label.getLabel(

1272 lab,

1273 label,

1274 boundingBoxes=boundingBoxes,

1275 centresOfMass=centresOfMass,

1276 extractCube=True,

1277 margin=2,

1278 maskOtherLabels=True,

1279 )

1280 # Gaussian smooth

1281 grainCubeFiltered = scipy.ndimage.gaussian_filter(GL["subvol"].astype("<f4"), sigma=gaussianFilterSigma)

1282 # mesh edge

1283 verts, faces, _, _ = skimage.measure.marching_cubes(grainCubeFiltered, level=0.5)

1284 # compute surface

1285 surfaceArea = skimage.measure.mesh_surface_area(verts, faces)

1286 # compute psi

1287 trueSphericity[label] = sphereSurfaceArea[label] / surfaceArea

1288 return trueSphericity

1289

1290

1291def convexVolume(

1292 lab,

1293 boundingBoxes=None,

1294 centresOfMass=None,

1295 volumes=None,

1296 nProcesses=nProcessesDefault,

1297 verbose=True,

1298):

1299 """

1300 This function compute the convex hull of each label of the labelled image and return a

1301 list with the convex volume of each particle.

1302

1303 Parameters

1304 ----------

1305 lab : 3D array of integers

1306 Labelled volume, with lab.max() labels

1307

1308 boundingBoxes : lab.max()x6 array of ints, optional

1309 Bounding boxes in format returned by ``boundingBoxes``.

1310 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

1311

1312 centresOfMass : lab.max()x3 array of floats, optional

1313 Centres of mass in format returned by ``centresOfMass``.

1314 If not defined (Default = None), it is recomputed by running ``centresOfMass``

1315

1316 volumes : lab.max()x1 array of ints

1317 Volumes in format returned by ``volumes``

1318 If not defined (Default = None), it is recomputed by running ``volumes``

1319

1320 nProcesses : integer (optional, default = nProcessesDefault)

1321 Number of processes for multiprocessing

1322 Default = number of CPUs in the system

1323

1324 verbose : boolean (optional, default = False)

1325 True for printing the evolution of the process

1326 False for not printing the evolution of process

1327

1328 Returns

1329 --------

1330

1331 convexVolume : lab.max()x1 array of floats with the convex volume.

1332

1333 Note

1334 ----

1335 convexVolume can only be computed for particles with volume greater than 3 voxels. If it is not the case, it will return 0.

1336

1337 """

1338 lab = lab.astype(labelType)

1339

1340 # Compute boundingBoxes if needed

1341 if boundingBoxes is None:

1342 boundingBoxes = spam.label.boundingBoxes(lab)

1343 # Compute centresOfMass if needed

1344 if centresOfMass is None:

1345 centresOfMass = spam.label.centresOfMass(lab)

1346 # Compute volumes if needed

1347 if volumes is None:

1348 volumes = spam.label.volumes(lab)

1349 # Compute number of labels

1350 nLabels = lab.max()

1351

1352 # Result array

1353 convexVolume = numpy.zeros(nLabels + 1, dtype="float")

1354

1355 if verbose:

1356 widgets = [

1357 progressbar.FormatLabel(""),

1358 " ",

1359 progressbar.Bar(),

1360 " ",

1361 progressbar.AdaptiveETA(),

1362 ]

1363 pbar = progressbar.ProgressBar(widgets=widgets, maxval=nLabels)

1364 pbar.start()

1365 finishedNodes = 0

1366

1367 # Function for convex volume

1368 global computeConvexVolume

1369

1370 def computeConvexVolume(label):

1371 labelI = spam.label.getLabel(lab, label, boundingBoxes=boundingBoxes, centresOfMass=centresOfMass)

1372 subvol = labelI["subvol"]

1373 points = numpy.transpose(numpy.where(subvol))

1374 try:

1375 hull = scipy.spatial.ConvexHull(points)

1376 deln = scipy.spatial.Delaunay(points[hull.vertices])

1377 idx = numpy.stack(numpy.indices(subvol.shape), axis=-1)

1378 out_idx = numpy.nonzero(deln.find_simplex(idx) + 1)

1379 hullIm = numpy.zeros(subvol.shape)

1380 hullIm[out_idx] = 1

1381 hullVol = spam.label.volumes(hullIm)

1382 return label, hullVol[-1]

1383 except Exception:

1384 return label, 0

1385

1386 # Run multiprocessing

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

1388 for returns in pool.imap_unordered(computeConvexVolume, range(1, nLabels + 1)):

1389 if verbose:

1390 finishedNodes += 1

1391 pbar.update(finishedNodes)

1392 convexVolume[returns[0]] = returns[1]

1393 pool.close()

1394 pool.join()

1395

1396 if verbose:

1397 pbar.finish()

1398

1399 return convexVolume

1400

1401

1402def moveLabels(

1403 lab,

1404 PhiField,

1405 returnStatus=None,

1406 boundingBoxes=None,

1407 centresOfMass=None,

1408 margin=3,

1409 PhiCOM=True,

1410 threshold=0.5,

1411 labelDilate=0,

1412 nProcesses=nProcessesDefault,

1413):

1414 """

1415 This function applies a discrete Phi field (from DDIC?) over a labelled image.

1416

1417 Parameters

1418 -----------

1419 lab : 3D numpy array

1420 Labelled image

1421

1422 PhiField : (multidimensional x 4 x 4 numpy array of floats)

1423 Spatial field of Phis

1424

1425 returnStatus : lab.max()x1 array of ints, optional

1426 Array with the return status for each label (usually returned by ``spam-ddic``)

1427 If not defined (Default = None), all the labels will be moved

1428 If returnStatus[i] == 2, the label will be moved, otherwise is omitted and erased from the final image

1429

1430 boundingBoxes : lab.max()x6 array of ints, optional

1431 Bounding boxes in format returned by ``boundingBoxes``.

1432 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

1433

1434 centresOfMass : lab.max()x3 array of floats, optional

1435 Centres of mass in format returned by ``centresOfMass``.

1436 If not defined (Default = None), it is recomputed by running ``centresOfMass``

1437

1438 margin : int, optional

1439 Margin, in pixels, to take in each label.

1440 Default = 3

1441

1442 PhiCOM : bool, optional

1443 Apply Phi to centre of mass of particle?, otherwise it will be applied in the middle of the particle\'s bounding box.

1444 Default = True

1445

1446 threshold : float, optional

1447 Threshold to keep interpolated voxels in the binary image.

1448 Default = 0.5

1449

1450 labelDilate : int, optional

1451 Number of times label should be dilated/eroded before returning it.

1452 If ``labelDilate > 0`` a dilated label is returned, while ``labelDilate < 0`` returns an eroded label.

1453 Default = 0

1454

1455 nProcesses : integer (optional, default = nProcessesDefault)

1456 Number of processes for multiprocessing

1457 Default = number of CPUs in the system

1458

1459 Returns

1460 --------

1461 labOut : 3D numpy array

1462 New labelled image with the labels moved by the deformations established by the PhiField.

1463

1464 Note

1465 ----

1466 When using more than one process (nProcesses > 1), the order of label updates is not guaranteed due to the use of imap_unordered. As a result, small differences in the final output may occur, especially in cases where labels overlap or touch.

1467

1468 """

1469

1470 # Check for boundingBoxes

1471 if boundingBoxes is None:

1472 boundingBoxes = spam.label.boundingBoxes(lab)

1473 # Check for centresOfMass

1474 if centresOfMass is None:

1475 centresOfMass = spam.label.centresOfMass(lab)

1476 # Create output label image

1477 labOut = numpy.zeros_like(lab, dtype=spam.label.labelType)

1478 # Get number of labels

1479 numberOfLabels = min(lab.max(), PhiField.shape[0] - 1)

1480 numberOfLabelsToMove = 0

1481 labelsToMove = []

1482 # Add the labels to move

1483 for label in range(1, numberOfLabels + 1):

1484 # Skip the particles if the returnStatus == 2 and returnStatus != None

1485 if type(returnStatus) == numpy.ndarray and returnStatus[label] != 2:

1486 pass

1487 else: # Add the particles

1488 labelsToMove.append(label)

1489 numberOfLabelsToMove += 1

1490

1491 # Function for moving labels

1492 global funMoveLabels

1493

1494 def funMoveLabels(label):

1495 getLabelReturn = spam.label.getLabel(

1496 lab,

1497 label,

1498 labelDilate=labelDilate,

1499 margin=margin,

1500 boundingBoxes=boundingBoxes,

1501 centresOfMass=centresOfMass,

1502 extractCube=True,

1503 )

1504 # Check that the label exist

1505 if getLabelReturn is not None:

1506 # Get Phi field

1507 Phi = PhiField[label].copy()

1508 # Phi will be split into a local part and a part of floored displacements

1509 disp = numpy.floor(Phi[0:3, -1]).astype(int)

1510 Phi[0:3, -1] -= disp

1511 # Check that the displacement exist

1512 if numpy.isfinite(disp).sum() == 3:

1513 # Just move binary label

1514 # Need to do backtracking here to avoid holes in the NN interpolation

1515 # Here we will cheat and do order 1 and re-threshold full pixels

1516 if PhiCOM:

1517 labSubvolDefInterp = spam.DIC.applyPhi(

1518 getLabelReturn["subvol"],

1519 Phi=Phi,

1520 interpolationOrder=1,

1521 PhiCentre=getLabelReturn["centreOfMassREL"],

1522 )

1523 else:

1524 labSubvolDefInterp = spam.DIC.applyPhi(

1525 getLabelReturn["subvol"],

1526 Phi=Phi,

1527 interpolationOrder=1,

1528 PhiCentre=(numpy.array(getLabelReturn["subvol"].shape) - 1) / 2.0,

1529 )

1530

1531 # "death mask"

1532 labSubvolDefMask = labSubvolDefInterp >= threshold

1533

1534 del labSubvolDefInterp

1535 # Get the boundary of the cube

1536 topOfSlice = numpy.array(

1537 [

1538 getLabelReturn["boundingBoxCube"][0] + disp[0],

1539 getLabelReturn["boundingBoxCube"][2] + disp[1],

1540 getLabelReturn["boundingBoxCube"][4] + disp[2],

1541 ]

1542 )

1543

1544 botOfSlice = numpy.array(

1545 [

1546 getLabelReturn["boundingBoxCube"][1] + disp[0],

1547 getLabelReturn["boundingBoxCube"][3] + disp[1],

1548 getLabelReturn["boundingBoxCube"][5] + disp[2],

1549 ]

1550 )

1551

1552 topOfSliceCrop = numpy.array(

1553 [

1554 max(topOfSlice[0], 0),

1555 max(topOfSlice[1], 0),

1556 max(topOfSlice[2], 0),

1557 ]

1558 )

1559 botOfSliceCrop = numpy.array(

1560 [

1561 min(botOfSlice[0], lab.shape[0]),

1562 min(botOfSlice[1], lab.shape[1]),

1563 min(botOfSlice[2], lab.shape[2]),

1564 ]

1565 )

1566 # Update grainSlice with disp

1567 grainSlice = (

1568 slice(topOfSliceCrop[0], botOfSliceCrop[0]),

1569 slice(topOfSliceCrop[1], botOfSliceCrop[1]),

1570 slice(topOfSliceCrop[2], botOfSliceCrop[2]),

1571 )

1572

1573 # Update labSubvolDefMask

1574 labSubvolDefMaskCrop = labSubvolDefMask[

1575 topOfSliceCrop[0] - topOfSlice[0] : labSubvolDefMask.shape[0] - 1 + botOfSliceCrop[0] - botOfSlice[0],

1576 topOfSliceCrop[1] - topOfSlice[1] : labSubvolDefMask.shape[1] - 1 + botOfSliceCrop[1] - botOfSlice[1],

1577 topOfSliceCrop[2] - topOfSlice[2] : labSubvolDefMask.shape[2] - 1 + botOfSliceCrop[2] - botOfSlice[2],

1578 ]

1579 return label, grainSlice, labSubvolDefMaskCrop, 1

1580

1581 # Nan displacement, run away

1582 else:

1583 return label, 0, 0, -1

1584 # Got None from getLabel()

1585 else:

1586 return label, 0, 0, -1

1587

1588 # Create progressbar

1589 widgets = [

1590 progressbar.FormatLabel(""),

1591 " ",

1592 progressbar.Bar(),

1593 " ",

1594 progressbar.AdaptiveETA(),

1595 ]

1596 pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfLabels)

1597 pbar.start()

1598 finishedNodes = 0

1599 # Run multiprocessing

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

1601 for returns in pool.imap_unordered(funMoveLabels, labelsToMove):

1602 finishedNodes += 1

1603 # widgets[0] = progressbar.FormatLabel("{}/{} ".format(finishedNodes, numberOfLabels))

1604 pbar.update(finishedNodes)

1605 # Set voxels in slice to the value of the label if not in greyscale mode

1606 if returns[0] > 0 and returns[3] == 1:

1607 labOut[returns[1]][returns[2]] = returns[0]

1608 pool.close()

1609 pool.join()

1610

1611 # End progressbar

1612 pbar.finish()

1613

1614 return labOut

1615

1616

1617def erodeLabels(lab, erosion=1, boundingBoxes=None, centresOfMass=None, nProcesses=nProcessesDefault):

1618 """

1619 This function erodes a labelled image.

1620

1621 Parameters

1622 -----------

1623 lab : 3D numpy array

1624 Labelled image

1625

1626 erosion : int, optional

1627 Erosion level

1628

1629 boundingBoxes : lab.max()x6 array of ints, optional

1630 Bounding boxes in format returned by ``boundingBoxes``.

1631 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

1632

1633 centresOfMass : lab.max()x3 array of floats, optional

1634 Centres of mass in format returned by ``centresOfMass``.

1635 If not defined (Default = None), it is recomputed by running ``centresOfMass``

1636

1637 nProcesses : integer (optional, default = nProcessesDefault)

1638 Number of processes for multiprocessing

1639 Default = number of CPUs in the system

1640

1641 Returns

1642 --------

1643 erodeImage : 3D numpy array

1644 New labelled image with the eroded labels.

1645

1646 Note

1647 ----

1648 The function makes use of spam.label.moveLabels() to generate the eroded image.

1649

1650 """

1651 # Get number of labels

1652 numberOfLabels = lab.max()

1653 # Create the Empty Phi field

1654 PhiField = numpy.zeros((numberOfLabels + 1, 4, 4))

1655 # Setup Phi as the identity

1656 for i in range(0, numberOfLabels + 1, 1):

1657 PhiField[i] = numpy.eye(4)

1658 # Use moveLabels

1659 erodeImage = spam.label.moveLabels(

1660 lab,

1661 PhiField,

1662 boundingBoxes=boundingBoxes,

1663 centresOfMass=centresOfMass,

1664 margin=1,

1665 PhiCOM=True,

1666 threshold=0.5,

1667 labelDilate=-erosion,

1668 nProcesses=nProcesses,

1669 )

1670 return erodeImage

1671

1672

1673def convexFillHoles(lab, boundingBoxes=None, centresOfMass=None):

1674 """

1675 This function fills the holes computing the convex volume around each label.

1676

1677 Parameters

1678 -----------

1679 lab : 3D numpy array

1680 Labelled image

1681

1682 boundingBoxes : lab.max()x6 array of ints, optional

1683 Bounding boxes in format returned by ``boundingBoxes``.

1684 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

1685

1686 centresOfMass : lab.max()x3 array of floats, optional

1687 Centres of mass in format returned by ``centresOfMass``.

1688 If not defined (Default = None), it is recomputed by running ``centresOfMass``

1689

1690 Returns

1691 --------

1692 labOut : 3D numpy array

1693 New labelled image.

1694

1695 Note

1696 ----

1697 The function works nicely for convex particles. For non-convex particles, it will alter the shape.

1698

1699 """

1700

1701 # Check for boundingBoxes

1702 if boundingBoxes is None:

1703 boundingBoxes = spam.label.boundingBoxes(lab)

1704 # Check for centresOfMass

1705 if centresOfMass is None:

1706 centresOfMass = spam.label.centresOfMass(lab)

1707 # Create output label image

1708 labOut = numpy.zeros_like(lab, dtype=spam.label.labelType)

1709 # Get number of labels

1710 numberOfLabels = lab.max()

1711 # Create progressbar

1712 widgets = [

1713 progressbar.FormatLabel(""),

1714 " ",

1715 progressbar.Bar(),

1716 " ",

1717 progressbar.AdaptiveETA(),

1718 ]

1719 pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfLabels)

1720 pbar.start()

1721 for i in range(1, numberOfLabels + 1, 1):

1722 # Get label

1723 getLabelReturn = spam.label.getLabel(

1724 lab,

1725 i,

1726 labelDilate=0,

1727 margin=3,

1728 boundingBoxes=boundingBoxes,

1729 centresOfMass=centresOfMass,

1730 maskOtherLabels=False,

1731 )

1732 # Get subvolume

1733 subVol = getLabelReturn["subvol"]

1734 # Transform to binary

1735 subVolBinMask = (subVol > 0).astype(int)

1736 # Mask out all the other labels

1737 subVolBinMaskLabel = numpy.where(subVol == i, 1, 0).astype(int)

1738 # Mask only the current label - save all the other labels

1739 subVolMaskOtherLabel = subVolBinMask - subVolBinMaskLabel

1740 # Fill holes with convex volume

1741 points = numpy.transpose(numpy.where(subVolBinMaskLabel))

1742 hull = scipy.spatial.ConvexHull(points)

1743 deln = scipy.spatial.Delaunay(points[hull.vertices])

1744 idx = numpy.stack(numpy.indices(subVol.shape), axis=-1)

1745 out_idx = numpy.nonzero(deln.find_simplex(idx) + 1)

1746 hullIm = numpy.zeros(subVol.shape)

1747 hullIm[out_idx] = 1

1748 hullIm = hullIm > 0

1749 # Identify added voxels

1750 subVolAdded = hullIm - subVolBinMaskLabel

1751 # Identify the wrong voxels - they are inside other labels

1752 subVolWrongAdded = subVolAdded * subVolMaskOtherLabel

1753 # Remove wrong filling areas

1754 subVolCorrect = (hullIm - subVolWrongAdded) > 0

1755 # Get slice

1756 grainSlice = (

1757 slice(getLabelReturn["slice"][0].start, getLabelReturn["slice"][0].stop),

1758 slice(getLabelReturn["slice"][1].start, getLabelReturn["slice"][1].stop),

1759 slice(getLabelReturn["slice"][2].start, getLabelReturn["slice"][2].stop),

1760 )

1761 # Add it to the output file

1762 labOut[grainSlice][subVolCorrect] = i

1763 # Update the progressbar

1764 widgets[0] = progressbar.FormatLabel("{}/{} ".format(i, numberOfLabels))

1765 pbar.update(i)

1766

1767 return labOut

1768

1769

1770def getNeighbours(

1771 lab,

1772 listOfLabels,

1773 method="getLabel",

1774 neighboursRange=None,

1775 centresOfMass=None,

1776 boundingBoxes=None,

1777):

1778 """

1779 This function computes the neighbours for a list of labels.

1780

1781 Parameters

1782 -----------

1783 lab : 3D numpy array

1784 Labelled image

1785

1786 listOfLabels : list of ints

1787 List of labels to which the neighbours will be computed

1788

1789 method : string

1790 Method to compute the neighbours.

1791 'getLabel' : The neighbours are the labels inside the subset obtained through spam.getLabel()

1792 'mesh' : The neighbours are computed using a tetrahedral connectivity matrix

1793 Default = 'getLabel'

1794

1795 neighboursRange : int

1796 Parameter controlling the search range to detect neighbours for each method.

1797 For 'getLabel', it correspond to the size of the subset. Default = meanRadii

1798 For 'mesh', it correspond to the size of the alpha shape used for carving the mesh. Default = 5*meanDiameter.

1799

1800 boundingBoxes : lab.max()x6 array of ints, optional

1801 Bounding boxes in format returned by ``boundingBoxes``.

1802 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

1803

1804 centresOfMass : lab.max()x3 array of floats, optional

1805 Centres of mass in format returned by ``centresOfMass``.

1806 If not defined (Default = None), it is recomputed by running ``centresOfMass``

1807

1808 Returns

1809 --------

1810 neighbours : list

1811 List with the neighbours for each label in listOfLabels.

1812

1813 """

1814 # Create result list

1815 neighbours = []

1816 # Compute centreOfMass if needed

1817 if centresOfMass is None:

1818 centresOfMass = spam.label.centresOfMass(lab)

1819 # Compute boundingBoxes if needed

1820 if boundingBoxes is None:

1821 boundingBoxes = spam.label.boundingBoxes(lab)

1822 # Compute Radii

1823 radii = spam.label.equivalentRadii(lab)

1824 if method == "getLabel":

1825 # Compute neighboursRange if needed

1826 if neighboursRange is None:

1827 neighboursRange = numpy.mean(radii)

1828 # Compute for each label in the list of labels

1829 for label in listOfLabels:

1830 getLabelReturn = spam.label.getLabel(

1831 lab,

1832 label,

1833 labelDilate=neighboursRange,

1834 margin=neighboursRange,

1835 boundingBoxes=boundingBoxes,

1836 centresOfMass=centresOfMass,

1837 maskOtherLabels=False,

1838 )

1839 # Get subvolume

1840 subVol = getLabelReturn["subvol"]

1841 # Get neighbours

1842 neighboursLabel = numpy.unique(subVol)

1843 # Remove label and 0 from the list of neighbours

1844 neighboursLabel = neighboursLabel[~numpy.isin(neighboursLabel, label)]

1845 neighboursLabel = neighboursLabel[~numpy.isin(neighboursLabel, 0)]

1846 # Add the neighbours to the list

1847 neighbours.append(neighboursLabel)

1848

1849 elif method == "mesh":

1850 # Compute neighboursRange if needed

1851 if neighboursRange is None:

1852 neighboursRange = 5 * 2 * numpy.mean(radii)

1853 # Get connectivity matrix

1854 conn = spam.mesh.triangulate(centresOfMass, weights=radii**2, alpha=neighboursRange)

1855 # Compute for each label in the list of labels

1856 for label in listOfLabels:

1857 neighboursLabel = numpy.unique(conn[numpy.where(numpy.sum(conn == label, axis=1))])

1858 # Remove label from the list of neighbours

1859 neighboursLabel = neighboursLabel[~numpy.isin(neighboursLabel, label)]

1860 # Add the neighbours to the list

1861 neighbours.append(neighboursLabel)

1862 else:

1863 print("spam.label.getNeighbours(): Wrong method, aborting")

1864

1865 return neighbours

1866

1867

1868def detectUnderSegmentation(lab, nProcesses=nProcessesDefault, verbose=True):

1869 """

1870 This function computes the coefficient of undersegmentation for each particle, defined as the ratio of the convex volume and the actual volume.

1871

1872 Parameters

1873 -----------

1874 lab : 3D numpy array

1875 Labelled image

1876

1877 nProcesses : integer (optional, default = nProcessesDefault)

1878 Number of processes for multiprocessing

1879 Default = number of CPUs in the system

1880

1881 verbose : boolean (optional, default = False)

1882 True for printing the evolution of the process

1883

1884 Returns

1885 --------

1886 underSegCoeff : lab.max() array of floats

1887 An array of float values that suggests the respective labels are undersegmentated.

1888

1889 Note

1890 ----

1891 For perfect convex particles, any coefficient higher than 1 should be interpreted as a particle with undersegmentation problems.

1892 However, for natural materials the threshold to define undersegmentation varies.

1893 It is suggested to plot the histogram of the undersegmentation coefficient and select the threshold accordingly.

1894

1895 """

1896 # Compute the volume

1897 vol = spam.label.volumes(lab)

1898 # Compute the convex volume

1899 convexVol = spam.label.convexVolume(lab, verbose=verbose, nProcesses=nProcesses)

1900 # Set the volume of the void to 0 to avoid the division by zero error

1901 vol[0] = 1

1902 # Compute the underSegmentation Coefficient

1903 underSegCoeff = convexVol / vol

1904 # Set the coefficient of the void to 0

1905 underSegCoeff[0] = 0

1906 return underSegCoeff

1907

1908

1909def detectOverSegmentation(lab):

1910 """

1911 This function computes the coefficient of oversegmentation for each particle, defined as the ratio between a characteristic lenght of the maximum contact area

1912 and a characteristic length of the particle.

1913

1914 Parameters

1915 -----------

1916 lab : 3D numpy array

1917 Labelled image

1918

1919 Returns

1920 --------

1921 overSegCoeff : lab.max() array of floats

1922 An array of float values with the oversegmentation coefficient

1923

1924 sharedLabel : lab.max() array of floats

1925 Array of floats with the the oversegmentation coefficient neighbours - label that share the maximum contact area

1926

1927 Note

1928 ----

1929 The threshold to define oversegmentation is dependent on each material and conditions of the test.

1930 It is suggested to plot the histogram of the oversegmentation coefficient and select the threshold accordingly.

1931

1932 """

1933 # Get the labels

1934 labels = list(range(0, lab.max() + 1))

1935 # Compute the volumes

1936 vol = spam.label.volumes(lab)

1937 # Compute the eq diameter

1938 eqDiam = spam.label.equivalentRadii(lab)

1939 # Compute the areas

1940 contactLabels = spam.label.contactingLabels(lab, areas=True)

1941 # Create result list

1942 overSegCoeff = []

1943 sharedLabel = []

1944 for label in labels:

1945 if label == 0:

1946 overSegCoeff.append(0)

1947 sharedLabel.append(0)

1948 else:

1949 # Check if there are contacting areas and volumes

1950 if len(contactLabels[1][label]) > 0 and vol[label] > 0:

1951 # We have areas on the list, compute the area

1952 maxArea = numpy.max(contactLabels[1][label])

1953 # Get the label for the max contacting area

1954 maxLabel = contactLabels[0][label][numpy.argmax(contactLabels[1][label])]

1955 # Compute the coefficient

1956 overSegCoeff.append(maxArea * eqDiam[label] / vol[label])

1957 # Add the label

1958 sharedLabel.append(maxLabel)

1959 else:

1960 overSegCoeff.append(0)

1961 sharedLabel.append(0)

1962 overSegCoeff = numpy.array(overSegCoeff)

1963 sharedLabel = numpy.array(sharedLabel)

1964 return overSegCoeff, sharedLabel

1965

1966

1967def fixUndersegmentation(

1968 lab,

1969 imGrey,

1970 targetLabels,

1971 underSegCoeff,

1972 boundingBoxes=None,

1973 centresOfMass=None,

1974 imShowProgress=False,

1975 verbose=True,

1976):

1977 """

1978 This function fixes undersegmentation problems, by performing a watershed with a higher local threshold for the problematic labels.

1979

1980 Parameters

1981 -----------

1982 lab : 3D numpy array

1983 Labelled image

1984

1985 imGrey : 3D numpy array

1986 Normalised greyscale of the labelled image, with a greyscale range between 0 and 1 and with void/solid peaks at 0.25 and 0.75, respectively.

1987 You can use helpers.histogramTools.findHistogramPeaks and helpers.histogramTools.histogramNorm to obtain a normalized greyscale image.

1988

1989 targetLabels : int or a list of labels

1990 List of target labels to solve undersegmentation

1991

1992 underSegCoeff : lab.max() array of floats

1993 Undersegmentation coefficient as returned by ``detectUnderSegmentation``

1994

1995 boundingBoxes : lab.max()x6 array of ints, optional

1996 Bounding boxes in format returned by ``boundingBoxes``.

1997 If not defined (Default = None), it is recomputed by running ``boundingBoxes``

1998

1999 centresOfMass : lab.max()x3 array of floats, optional

2000 Centres of mass in format returned by ``centresOfMass``.

2001 If not defined (Default = None), it is recomputed by running ``centresOfMass``boundingBoxes : 3D numpy array

2002 Labelled image

2003

2004 imShowProgress : bool, optional

2005 Graphical interface to observe the process for each label.

2006 Default = False

2007

2008 verbose : boolean (optional, default = False)

2009 True for printing the evolution of the process

2010

2011 Returns

2012 --------

2013 lab : 3D numpy array

2014 Labelled image after running ``makeLabelsSequential``

2015 """

2016

2017 # Usual checks

2018 if boundingBoxes is None:

2019 boundingBoxes = spam.label.boundingBoxes(lab)

2020 if centresOfMass is None:

2021 centresOfMass = spam.label.centresOfMass(lab)

2022 # Check if imGrey is normalised (limits [0,1])

2023 if imGrey.max() > 1 or imGrey.min() < 0:

2024 print("\n spam.label.fixUndersegmentation(): imGrey is not normalised. Limits exceed [0,1]")

2025 return

2026 # Start counters

2027 labelCounter = numpy.max(lab)

2028 labelDummy = numpy.zeros(lab.shape)

2029 successCounter = 0

2030 finishedLabels = 0

2031 if verbose:

2032 widgets = [

2033 progressbar.FormatLabel(""),

2034 " ",

2035 progressbar.Bar(),

2036 " ",

2037 progressbar.AdaptiveETA(),

2038 ]

2039 pbar = progressbar.ProgressBar(widgets=widgets, maxval=len(targetLabels))

2040 pbar.start()

2041 # Main loop

2042 for label in targetLabels:

2043 # Get the subset

2044 labelData = spam.label.getLabel(

2045 lab,

2046 label,

2047 margin=5,

2048 boundingBoxes=boundingBoxes,

2049 centresOfMass=centresOfMass,

2050 extractCube=True,

2051 )

2052 # Get the slice on the greyscale

2053 imGreySlice = imGrey[

2054 labelData["sliceCube"][0].start : labelData["sliceCube"][0].stop,

2055 labelData["sliceCube"][1].start : labelData["sliceCube"][1].stop,

2056 labelData["sliceCube"][2].start : labelData["sliceCube"][2].stop,

2057 ]

2058 # Mask the imGreySubset

2059 greySubset = imGreySlice * labelData["subvol"]

2060 # Create seeds

2061 # 2021-08-02 GP: Maybe this can be changed by just a serie of binary erosion?

2062 seeds = spam.label.watershed(greySubset >= 0.75)

2063 # Do we have seeds?

2064 if numpy.max(seeds) < 1:

2065 # The threshold was too harsh on the greySubset and there are no seeds

2066 # We shouldn't change this label

2067 passBool = "Decline"

2068 else:

2069 # We have at least one seed, Run watershed again with markers

2070 imLabSubset = spam.label.watershed(labelData["subvol"], markers=seeds)

2071 # Run again the underSegCoeff for the subset

2072 res = detectUnderSegmentation(imLabSubset, verbose=False)

2073 # Safety check - do we have any labels at all?

2074 if len(res) > 2:

2075 # We have at least one label

2076 # Check if it should pass or not - is the new underSegCoeff of all the new labels less than the original coefficient?

2077 if all(map(lambda x: x < underSegCoeff[label], res[1:])):

2078 # We can modify this label

2079 passBool = "Accept"

2080 successCounter += 1

2081 # Remove the label from the original label image

2082 lab = spam.label.removeLabels(lab, [label])

2083 # Assign the new labels to the grains

2084 # Create a subset to fill with the new labels

2085 imLabSubsetNew = numpy.zeros(imLabSubset.shape)

2086 for newLab in numpy.unique(imLabSubset[imLabSubset != 0]):

2087 imLabSubsetNew = numpy.where(imLabSubset == newLab, labelCounter + 1, imLabSubsetNew)

2088 labelCounter += 1

2089 # Create a disposable dummy sample to allocate the grains

2090 labelDummyUnit = numpy.zeros(lab.shape)

2091 # Alocate the grains

2092 labelDummyUnit[

2093 labelData["sliceCube"][0].start : labelData["sliceCube"][0].stop,

2094 labelData["sliceCube"][1].start : labelData["sliceCube"][1].stop,

2095 labelData["sliceCube"][2].start : labelData["sliceCube"][2].stop,

2096 ] = imLabSubsetNew

2097 # Add the grains

2098 labelDummy = labelDummy + labelDummyUnit

2099 else:

2100 # We shouldn't change this label

2101 passBool = "Decline"

2102

2103 if imShowProgress:

2104 # Enter graphical mode

2105 # Change the labels to show different colourss

2106 fig = plt.figure()

2107 # Plot

2108 plt.subplot(3, 2, 1)

2109 plt.gca().set_title("Before")

2110 plt.imshow(

2111 labelData["subvol"][labelData["subvol"].shape[0] // 2, :, :],

2112 cmap="Greys_r",

2113 )

2114 plt.subplot(3, 2, 2)

2115 plt.gca().set_title("After")

2116 plt.imshow(imLabSubset[imLabSubset.shape[0] // 2, :, :], cmap="cubehelix")

2117 plt.subplot(3, 2, 3)

2118 plt.imshow(

2119 labelData["subvol"][:, labelData["subvol"].shape[1] // 2, :],

2120 cmap="Greys_r",

2121 )

2122 plt.subplot(3, 2, 4)

2123 plt.imshow(imLabSubset[:, imLabSubset.shape[1] // 2, :], cmap="cubehelix")

2124 plt.subplot(3, 2, 5)

2125 plt.imshow(

2126 labelData["subvol"][:, :, labelData["subvol"].shape[2] // 2],

2127 cmap="Greys_r",

2128 )

2129 plt.subplot(3, 2, 6)

2130 plt.imshow(imLabSubset[:, :, imLabSubset.shape[2] // 2], cmap="cubehelix")

2131 fig.suptitle(

2132 # r"Label {}. Status: $\bf{}$".format(label, passBool), # breaks for matplotlib 3.7.0

2133 f"Label {label}. Status: {passBool}",

2134 fontsize="xx-large",

2135 )

2136 plt.show()

2137 else:

2138 # We shouldn't change this label

2139 passBool = "Decline"

2140 if verbose:

2141 finishedLabels += 1

2142 pbar.update(finishedLabels)

2143 # We finish, lets add the new grains to the labelled image

2144 lab = lab + labelDummy

2145 # Update the labels

2146 lab = spam.label.makeLabelsSequential(lab)

2147 if verbose:

2148 pbar.finish()

2149 print(f"\n spam.label.fixUndersegmentation(): From {len(targetLabels)} target labels, {successCounter} were modified")

2150 return lab

2151

2152

2153def fixOversegmentation(lab, targetLabels, sharedLabel, verbose=True, imShowProgress=False):

2154 """

2155 This function fixes oversegmentation problems, by merging each target label with its oversegmentation coefficient neighbour.

2156

2157 Parameters

2158 -----------

2159 lab : 3D numpy array

2160 Labelled image

2161

2162 targetLabels : int or a list of labels

2163 List of target labels to solve oversegmentation

2164

2165 sharedLabel : lab.max() array of floats

2166 List ofoversegmentation coefficient neighbour as returned by ``detectOverSegmentation``

2167

2168 imShowProgress : bool, optional

2169 Graphical interface to observe the process for each label.

2170 Default = False

2171

2172 verbose : boolean (optional, default = False)

2173 True for printing the evolution of the process

2174

2175 Returns

2176 --------

2177 lab : 3D numpy array

2178 Labelled image after running ``makeLabelsSequential``

2179

2180 """

2181 # Start counters

2182 labelDummy = numpy.zeros(lab.shape)

2183 finishedLabelsCounter = 0

2184 finishedLabels = []

2185 if verbose:

2186 widgets = [

2187 progressbar.FormatLabel(""),

2188 " ",

2189 progressbar.Bar(),

2190 " ",

2191 progressbar.AdaptiveETA(),

2192 ]

2193 pbar = progressbar.ProgressBar(widgets=widgets, maxval=len(targetLabels))

2194 pbar.start()

2195 # Main loop

2196 for labelA in targetLabels:

2197 # Verify that the label is not on the finished list

2198 if labelA in finishedLabels:

2199 # It is already on the list, move on

2200 pass

2201 else:

2202 # Get the touching label

2203 labelB = sharedLabel[labelA]

2204 # Add then to the list

2205 finishedLabels.append(labelA)

2206 finishedLabels.append(labelB)

2207 # Get the subset of the two labels

2208 subset = spam.label.fetchTwoGrains(lab, [labelA, labelB])

2209 # Change the labelB by labelA in the subset

2210 subVolLabNew = numpy.where(subset["subVolLab"] == labelB, labelA, subset["subVolLab"])

2211 # Create a disposable dummy sample to allocate the grains

2212 labelDummyUnit = numpy.zeros(lab.shape)

2213 # Alocate the grains

2214 labelDummyUnit[

2215 subset["slice"][0].start : subset["slice"][0].stop,

2216 subset["slice"][1].start : subset["slice"][1].stop,

2217 subset["slice"][2].start : subset["slice"][2].stop,

2218 ] = subVolLabNew

2219 # Add the grains

2220 labelDummy = labelDummy + labelDummyUnit

2221 # Remove the label from the original label image

2222 lab = spam.label.removeLabels(lab, [labelA, labelB])

2223 # Enter graphical mode

2224 if imShowProgress:

2225 # Change the labels to show different colourss

2226 subVolLabNorm = numpy.where(subset["subVolLab"] == labelA, 1, subset["subVolLab"])

2227 subVolLabNorm = numpy.where(subset["subVolLab"] == labelB, 2, subVolLabNorm)

2228 fig = plt.figure()

2229 plt.subplot(3, 2, 1)

2230 plt.gca().set_title("Before")

2231 plt.imshow(

2232 subVolLabNorm[subset["subVolLab"].shape[0] // 2, :, :],

2233 cmap="cubehelix",

2234 )

2235 plt.subplot(3, 2, 2)

2236 plt.gca().set_title("After")

2237 plt.imshow(subVolLabNew[subVolLabNew.shape[0] // 2, :, :], cmap="cubehelix")

2238 plt.subplot(3, 2, 3)

2239 plt.imshow(

2240 subVolLabNorm[:, subset["subVolLab"].shape[1] // 2, :],

2241 cmap="cubehelix",

2242 )

2243 plt.subplot(3, 2, 4)

2244 plt.imshow(subVolLabNew[:, subVolLabNew.shape[1] // 2, :], cmap="cubehelix")

2245 plt.subplot(3, 2, 5)

2246 plt.imshow(

2247 subVolLabNorm[:, :, subset["subVolLab"].shape[2] // 2],

2248 cmap="cubehelix",

2249 )

2250 plt.subplot(3, 2, 6)

2251 plt.imshow(subVolLabNew[:, :, subVolLabNew.shape[2] // 2], cmap="cubehelix")

2252 fig.suptitle("Label {} and {}".format(labelA, labelB), fontsize="xx-large")

2253 plt.show()

2254 if verbose:

2255 finishedLabelsCounter += 1

2256 pbar.update(finishedLabelsCounter)

2257 # We finish, lets add the new grains to the labelled image

2258 lab = lab + labelDummy

2259 # Update the labels

2260 lab = spam.label.makeLabelsSequential(lab)

2261 if verbose:

2262 pbar.finish()

2263

2264 return lab

2265

2266

2267@numba.njit(cache=True)

2268def _updateLabels(lab, newLabels):

2269 """

2270 This function uses numba to go through all the voxels of a label image and assign a new label based on the newLabels list.

2271

2272 Parameters

2273 -----------

2274 lab : 3D array of integers

2275 Labelled volume, with lab.max() labels

2276

2277 newLabels : 1D array of integers

2278 Array with the order of the new labels

2279

2280 Returns

2281 --------

2282 lab : 3D array of integers

2283 Labelled volume, with lab.max() labels

2284

2285 """

2286

2287 # Loop over the image to change the values

2288 for z in range(lab.shape[0]):

2289 for y in range(lab.shape[1]):

2290 for x in range(lab.shape[2]):

2291 if lab[z, y, x] != 0:

2292 lab[z, y, x] = newLabels[lab[z, y, x]]

2293 return lab

2294

2295

2296def shuffleLabels(lab):

2297 """

2298 This function re-assigns randomly the labels of a label image, usually for visualisation purposes.

2299

2300 Parameters

2301 -----------

2302 lab : 3D array of integers

2303 Labelled volume, with lab.max() labels

2304

2305 Returns

2306 --------

2307 lab : 3D array of integers

2308 Labelled volume, with lab.max() labels

2309

2310 """

2311 # Ge the list of labels and a copy

2312 labels = numpy.unique(lab).tolist()

2313 labelsList = labels.copy()

2314 # Empty list ot fill

2315 newLabels = []

2316 # Loop over the labels

2317 for i in range(len(labelsList)):

2318 if i == 0:

2319 # Add the 0 to the list

2320 newLabels.append(0)

2321 # Remove 0 from the list

2322 labelsList.remove(0)

2323 else:

2324 # Generate a random number from a list length

2325 pos = random.randrange(len(labelsList))

2326 # Add the label to the list

2327 newLabels.append(labelsList[pos])

2328 # Remove 0 from the list

2329 labelsList.remove(labelsList[pos])

2330

2331 # Transform them into array for easier indexing

2332 labels = numpy.asarray(labels)

2333 newLabels = numpy.asarray(newLabels)

2334

2335 # Update the labels using the numba function

2336 lab = _updateLabels(lab, newLabels)

2337

2338 return lab