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

1"""

2Library of SPAM functions for dealing with labelled images

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

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

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

8any later version.

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

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

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

13more details.

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

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

17"""

19import multiprocessing

21import matplotlib

22import matplotlib.pyplot as plt

23import numpy

24import progressbar

25import scipy.ndimage

26import scipy.spatial

27import spam.DIC

28import spam.filters

29from spam.label.labelToolkit import boundingBoxes as boundingBoxesCPP

30from spam.label.labelToolkit import centresOfMass as centresOfMassCPP

31from spam.label.labelToolkit import labelToFloat as labelToFloatCPP

32from spam.label.labelToolkit import momentOfInertia as momentOfInertiaCPP

33from spam.label.labelToolkit import relabel as relabelCPP

34from spam.label.labelToolkit import setVoronoi as setVoronoiCPP

35from spam.label.labelToolkit import tetPixelLabel as tetPixelLabelCPP

36from spam.label.labelToolkit import volumes as volumesCPP

38# Define a random colourmap for showing labels

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

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

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

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

43randomCmap = matplotlib.colors.ListedColormap(randomCmapVals)

44del randomCmapVals

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

48labelType = "<u4"

50# Global number of processes

51nProcessesDefault = multiprocessing.cpu_count()

54def boundingBoxes(lab):

55 """

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

58 Parameters

59 ----------

60 lab : 3D array of integers

61 Labelled volume, with lab.max() labels

63 Returns

64 -------

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

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

68 - Zmin, Zmax

69 - Ymin, Ymax

70 - Xmin, Xmax

72 Note

73 ----

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

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

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

77 1,1,1,1,1,1

79 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

81 """

82 # Catch 2D image, and pad

83 if lab.ndim == 2:

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

86 lab = lab.astype(labelType)

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

90 boundingBoxesCPP(lab, boundingBoxes)

92 return boundingBoxes

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

96 """

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

99 Parameters

100 ----------

101 lab : 3D array of integers

102 Labelled volume, with lab.max() labels

103

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

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

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

107

108 minVol : int, optional

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

110

111 Returns

112 -------

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

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

115 """

116 if boundingBoxes is None:

117 boundingBoxes = spam.label.boundingBoxes(lab)

118 if minVol is None:

119 minVol = 0

120 # Catch 2D image, and pad

121 if lab.ndim == 2:

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

123

124 lab = lab.astype(labelType)

125

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

127

128 centresOfMassCPP(lab, boundingBoxes, centresOfMass, minVol)

129

130 return centresOfMass

131

132

133def volumes(lab, boundingBoxes=None):

134 """

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

136

137 Parameters

138 ----------

139 lab : 3D array of integers

140 Labelled volume, with lab.max() labels

141

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

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

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

145

146 Returns

147 -------

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

149 This array contains the volume in voxels of each label

150 """

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

152

153 lab = lab.astype(labelType)

154

155 if boundingBoxes is None:

156 boundingBoxes = spam.label.boundingBoxes(lab)

157

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

159

160 volumesCPP(lab, boundingBoxes, volumes)

161

162 return volumes

163

164

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

166 """

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

168

169 Parameters

170 ----------

171 lab : 3D array of integers

172 Labelled volume, with lab.max() labels

173

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

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

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

177

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

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

180

181 Returns

182 -------

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

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

185 """

186

187 def vol2rad(volumes):

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

189

190 # If we have volumes, just go for it

191 if volumes is not None:

192 return vol2rad(volumes)

193

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

195 if boundingBoxes is None:

196 boundingBoxes = spam.label.boundingBoxes(lab)

197

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

199

200

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

202 """

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

204

205 Parameters

206 ----------

207 lab : 3D array of integers

208 Labelled volume, with lab.max() labels

209

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

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

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

213

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

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

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

217

218 minVol : int, optional

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

220 Default = default for spam.label.centresOfMass

221

222 Returns

223 -------

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

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

226

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

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

229 """

230 if boundingBoxes is None:

231 boundingBoxes = spam.label.boundingBoxes(lab)

232 if centresOfMass is None:

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

234

235 lab = lab.astype(labelType)

236

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

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

239

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

241

242 return [eigenValues, eigenVectors]

243

244

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

246 """

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

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

249

250 See appendix of for inital work:

251 "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).

252

253 Parameters

254 ----------

255 lab : 3D array of integers

256 Labelled volume, with lab.max() labels

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

258

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

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

261

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

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

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

265

266 enforceVolume = bool (default = True)

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

268 the same as the particle?

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

270

271 twoD : bool (default = False)

272 Are these in fact 2D ellipses?

273 Not implemented!!

274

275 Returns

276 -------

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

278 a, b, c lengths of particle in pixels

279

280 Note

281 -----

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

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

284

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

286

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

288

289 Function contributed by Takashi Matsushima (University of Tsukuba)

290 """

291 # Full ref:

292 # @misc{ikeda2000three,

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

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

295 # year={2000},

296 # publisher={De Gruyter}

297 # }

298

299 if volumes is None:

300 volumes = spam.label.volumes(lab)

301 if MOIeigenValues is None:

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

303

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

305

306 Ia = MOIeigenValues[:, 0]

307 Ib = MOIeigenValues[:, 1]

308 Ic = MOIeigenValues[:, 2]

309

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

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

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

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

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

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

316

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

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

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

320 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)

321 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)

322 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)

323

324 if enforceVolume:

325 # Compute volume of ellipse:

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

327 # filter zeros and infs

328 # print volumes.shape

329 # print ellipseVol.shape

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

331 # print volRatio

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

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

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

335

336 return ABCaxes

337

338

339def convertLabelToFloat(lab, vector):

340 """

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

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

343

344 Parameters

345 ----------

346 lab : 3D array of integers

347 Labelled volume, with lab.max() labels

348

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

350

351 Returns

352 -------

353 relabelled : 3D array of converted floats

354 """

355 lab = lab.astype(labelType)

356

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

358

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

360

361 labelToFloatCPP(lab, vector, relabelled)

362

363 return relabelled

364

365

366def makeLabelsSequential(lab):

367 """

368 This function fills gaps in labelled images,

369 by relabelling them to be sequential integers.

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

371

372 Parameters

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

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

375 An array of labels with 0 as the background

376

377 Returns

378 --------

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

380 An array of labels with 0 as the background

381 """

382 maxLabel = int(lab.max())

383 lab = lab.astype(labelType)

384

385 uniqueLabels = numpy.unique(lab)

386 # print uniqueLabels

387

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

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

390

391 relabelCPP(lab, relabelMap)

392

393 return lab

394

395

396def getLabel(

397 labelledVolume,

398 label,

399 boundingBoxes=None,

400 centresOfMass=None,

401 margin=0,

402 extractCube=False,

403 extractCubeSize=None,

404 maskOtherLabels=True,

405 labelDilate=0,

406 labelDilateMaskOtherLabels=False,

407):

408 """

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

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

411 Passing boundingBoxes and centresOfMass is highly recommended.

412

413 Parameters

414 ----------

415 labelVolume : 3D array of ints

416 3D Labelled volume

417

418 label : int

419 Label that we want information about

420

421 boundingBoxes : nLabels*2 array of ints, optional

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

423 Optional but highly recommended.

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

425

426 centresOfMass : nLabels*3 array of floats, optional

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

428 Optional but highly recommended.

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

430

431 extractCube : bool, optional

432 Return label subvolume in the middle of a cube?

433 Default = False

434

435 extractCubeSize : int, optional

436 half-size of cube to extract.

437 Default = calculate minimum cube

438

439 margin : int, optional

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

441 Default = 0

442

443 maskOtherLabels : bool, optional

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

445 If true, the mask is directly returned.

446 Default = True

447

448 labelDilate : int, optional

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

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

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

452 Requires ``maskOtherLabels``.

453 Default = 0

454

455 labelDilateMaskOtherLabels : bool, optional

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

457 Only pertinent for positive labelDilate values.

458 Default = False

459

460

461 Returns

462 -------

463 Dictionary containing:

464

465 Keys:

466 subvol : 3D array of bools or ints

467 subvolume from labelled image

468

469 slice : tuple of 3*slices

470 Slice used to extract subvol for the bounding box mode

471

472 sliceCube : tuple of 3*slices

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

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

475 and so it will be smaller than the returned cube

476

477 boundingBox : 1D numpy array of 6 ints

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

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

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

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

482

483 boundingBoxCube : 1D numpy array of 6 ints

484 Bounding box, including margin, in cube 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 centreOfMassABS : 3*float

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

491

492 centreOfMassREL : 3*float

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

494

495 volumeInitial : int

496 Volume of label (before dilating)

497

498 volumeDilated : int

499 Volume of label (after dilating, if requested)

500

501 """

502 import spam.mesh

503

504 if boundingBoxes is None:

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

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

507 boundingBoxes = spam.label.boundingBoxes(labelledVolume)

508

509 if centresOfMass is None:

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

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

512 centresOfMass = spam.label.centresOfMass(labelledVolume)

513

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

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

516 return

517 raise "No bounding boxes for this grain"

518

519 bbo = boundingBoxes[label]

520 com = centresOfMass[label]

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

522

523 # 1. Check if boundingBoxes are correct:

524 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):

525 pass

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

527

528 else:

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

530 output = {}

531 output["centreOfMassABS"] = com

532

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

534 if extractCube:

535 # Calculate offsets between centre of mass and bounding box

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

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

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

539

540 # If is none, assume closest fitting cube.

541 if extractCubeSize is not None:

542 if extractCubeSize < offset:

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

544 offset = int(extractCubeSize)

545

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

547 # if margin is not None:

548 offset += margin

549

550 offset = int(offset)

551

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

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

554

555 topOfSlice = numpy.array(

556 [

557 int(comRound[0] - offset),

558 int(comRound[1] - offset),

559 int(comRound[2] - offset),

560 ]

561 )

562 botOfSlice = numpy.array(

563 [

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

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

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

567 ]

568 )

569

570 labSubVol = spam.helpers.slicePadded(

571 labelledVolume,

572 [

573 topOfSlice[0],

574 botOfSlice[0],

575 topOfSlice[1],

576 botOfSlice[1],

577 topOfSlice[2],

578 botOfSlice[2],

579 ],

580 )

581

582 output["sliceCube"] = (

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

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

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

586 )

587

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

589 [

590 topOfSlice[0],

591 botOfSlice[0] - 1,

592 topOfSlice[1],

593 botOfSlice[1] - 1,

594 topOfSlice[2],

595 botOfSlice[2] - 1,

596 ]

597 )

598

599 output["centreOfMassREL"] = com - topOfSlice

600

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

602 else:

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

604 botOfSlice = numpy.array(

605 [

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

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

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

609 ]

610 )

611

612 labSubVol = spam.helpers.slicePadded(

613 labelledVolume,

614 [

615 topOfSlice[0],

616 botOfSlice[0],

617 topOfSlice[1],

618 botOfSlice[1],

619 topOfSlice[2],

620 botOfSlice[2],

621 ],

622 )

623

624 output["slice"] = (

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

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

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

628 )

629

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

631 [

632 topOfSlice[0],

633 botOfSlice[0] - 1,

634 topOfSlice[1],

635 botOfSlice[1] - 1,

636 topOfSlice[2],

637 botOfSlice[2] - 1,

638 ]

639 )

640

641 output["centreOfMassREL"] = com - topOfSlice

642

643 # Get mask for this label

644 maskLab = labSubVol == label

645 volume = numpy.sum(maskLab)

646 output["volumeInitial"] = volume

647

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

649 if maskOtherLabels:

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

651 # iterations of the standard.

652 if labelDilate > 0:

653 if labelDilate >= margin:

654 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!")

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

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

657 if labelDilateMaskOtherLabels:

658 # remove voxels that are neither our label nor pore

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

660 if labelDilate < 0:

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

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

663

664 # Just overwrite "labSubVol"

665 labSubVol = maskLab

666 # Update volume output

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

668

669 output["subvol"] = labSubVol

670

671 return output

672

673

674def getImagettesLabelled(

675 lab1,

676 label,

677 Phi,

678 im1,

679 im2,

680 searchRange,

681 boundingBoxes,

682 centresOfMass,

683 margin=0,

684 labelDilate=0,

685 maskOtherLabels=True,

686 applyF="all",

687 volumeThreshold=100,

688):

689 """

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

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

692

693 Parameters

694 ----------

695 lab1 : 3D numpy array of ints

696 Labelled image containing nLabels

697

698 label : int

699 Label of interest

700

701 Phi : 4x4 numpy array of floats

702 Phi matrix representing the movement of imagette1,

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

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

705

706 im1 : 3D numpy array

707 This is the large input reference image of greyvalues

708

709 im2 : 3D numpy array

710 This is the large input deformed image of greyvalues

711

712 searchRange : 6-component numpy array of ints

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

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

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

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

717

718 boundingBoxes : nLabels*2 array of ints

719 Bounding boxes as returned by ``boundingBoxes``

720

721 centresOfMass : nLabels*3 array of floats

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

723

724 margin : int, optional

725 Margin around the grain to extract in pixels

726 Default = 0

727

728 labelDilate : int, optional

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

730 Default = 0

731

732 maskOtherLabels : bool, optional

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

734 If true, the mask is directly returned.

735 Default = True

736

737 applyF : string, optional

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

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

740 Default = 'all'

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

742

743 volumeThreshold : int, optional

744 Pixel volume of labels that are discarded

745 Default = 100

746

747 Returns

748 -------

749 Dictionary :

750

751 'imagette1' : 3D numpy array,

752

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

754

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

756

757 'returnStatus': int,

758 Describes success in extracting imagette1 and imagette2.

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

760

761 'pixelSearchOffset': 3-component list of ints

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

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

764 """

765 returnStatus = 1

766 imagette1 = None

767 imagette1mask = None

768 imagette2 = None

769

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

771 PhiNoDisp = Phi.copy()

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

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

774 if applyF == "rigid":

775 PhiNoDisp = spam.deformation.computeRigidPhi(PhiNoDisp)

776

777 gottenLabel = spam.label.getLabel(

778 lab1,

779 label,

780 extractCube=False,

781 boundingBoxes=boundingBoxes,

782 centresOfMass=centresOfMass,

783 margin=labelDilate + margin,

784 maskOtherLabels=True,

785 labelDilate=labelDilate,

786 labelDilateMaskOtherLabels=maskOtherLabels,

787 )

788

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

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

791 returnStatus = -7

792

793 else:

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

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

796

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

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

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

800

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

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

803 imagette1mask = (

804 spam.DIC.applyPhi(

805 gottenLabel["subvol"] > 0,

806 PhiNoDisp,

807 PhiCentre=gottenLabel["centreOfMassREL"],

808 interpolationOrder=0,

809 )

810 > 0

811 )

812 elif applyF == "no":

813 imagette1mask = gottenLabel["subvol"]

814 else:

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

816

817 maskette1vol = numpy.sum(imagette1mask)

818

819 if maskette1vol > volumeThreshold:

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

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

822 startStopIm2 = [

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

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

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

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

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

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

829 ]

830

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

832

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

834

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

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

837 returnStatus = -5

838 else:

839 # Failed volume condition

840 returnStatus = -5

841

842 return {

843 "imagette1": imagette1,

844 "imagette1mask": imagette1mask,

845 "imagette2": imagette2,

846 "returnStatus": returnStatus,

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

848 }

849

850

851def labelsOnEdges(lab):

852 """

853 Return labels on edges of volume

854

855 Parameters

856 ----------

857 lab : 3D numpy array of ints

858 Labelled volume

859

860 Returns

861 -------

862 uniqueLabels : list of ints

863 List of labels on edges

864 """

865

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

867

868 uniqueLabels = []

869

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

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

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

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

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

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

876

877 # Flatten list of lists:

878 # 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

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

880

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

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

883

884 return uniqueLabels.astype(labelType)

885

886

887def removeLabels(lab, listOfLabelsToRemove):

888 """

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

890

891 Parameters

892 ----------

893 lab : 3D numpy array of ints

894 Labelled volume

895

896 listOfLabelsToRemove : list-like of ints

897 Labels to remove

898

899 Returns

900 -------

901 lab : 3D numpy array of ints

902 Labelled volume with desired labels blanked

903

904 Note

905 ----

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

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

908 """

909 lab = lab.astype(labelType)

910

911 # define a vector with sequential ints

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

913

914 # Remove the ones that have been asked for

915 for label in listOfLabelsToRemove:

916 arrayOfLabels[label] = 0

917

918 relabelCPP(lab, arrayOfLabels)

919

920 return lab

921

922

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

924 """

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

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

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

928

929 See:

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

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

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

933

934 and

935

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

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

938

939 Parameters

940 -----------

941 lab: 3D numpy array of labelTypes

942 Labelled image

943

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

945 Euclidean distance map of the pores.

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

947

948 maxPoreRadius: int (optional, default = 10)

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

950

951 Returns

952 --------

953 lab: 3D numpy array of labelTypes

954 Image labelled with set voronoi labels

955 """

956 if poreEDT is None:

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

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

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

960

961 lab = lab.astype(labelType)

962 labOut = numpy.zeros_like(lab)

963 maxPoreRadius = int(maxPoreRadius)

964

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

966 # This precomutation saves a lot of time

967 # Local grid of values, centred at zero

968 gridD = numpy.mgrid[

969 -maxPoreRadius : maxPoreRadius + 1,

970 -maxPoreRadius : maxPoreRadius + 1,

971 -maxPoreRadius : maxPoreRadius + 1,

972 ]

973

974 # Compute distances from centre

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

976 sortedIndices = numpy.argsort(Rarray)

977

978 # Array to hold sorted points

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

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

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

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

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

984 del gridD

985

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

987 # entry point into coords at the nearest integer values

988 sortedDistances = Rarray[sortedIndices]

989 indices = []

990 n = 0

991 i = 0

992 while i <= maxPoreRadius + 1:

993 if sortedDistances[n] >= i:

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

995 indices.append(n)

996 i += 1

997 n += 1

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

999

1000 # Call C++ code

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

1002

1003 return labOut

1004

1005

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

1007 """

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

1009

1010 Parameters

1011 ----------

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

1013

1014 points: number of points x 3 array of floats

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

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

1017

1018 connectivity: number of tetrahedra x 4 array of integers

1019 Connectivity matrix between points that define tetrahedra.

1020 Each line defines a tetrahedron whose number is the line number + 1.

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

1022 nThreads: int (optional, default=1)

1023 The number of threads used for the cpp parallelisation.

1024

1025 Returns

1026 -------

1027 3D array of ints, shape = dims

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

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

1030 """

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

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

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

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

1035

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

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

1038 # bold choice but that's ok

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

1040

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

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

1043

1044 tetPixelLabelCPP(lab, connectivity, points, nThreads)

1045

1046 return lab

1047

1048

1049def labelTetrahedraForScipyDelaunay(dims, delaunay):

1050 """

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

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

1053 when using `labelledTetrahedra`

1054

1055 Parameters

1056 ----------

1057 dims: tuple

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

1059

1060 delaunay: "delaunay" object

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

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

1063 the position of zero.

1064

1065 Returns

1066 -------

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

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

1069 """

1070

1071 # Big matrix of points poisitions

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

1073

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

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

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

1077

1078 del mgrid

1079

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

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

1082

1083 return lab

1084

1085

1086def filterIsolatedCells(array, struct, size):

1087 """

1088 Return array with completely isolated single cells removed

1089

1090 Parameters

1091 ----------

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

1093 Array with completely isolated single cells

1094

1095 struct: 3-D binary array

1096 Structure array for generating unique regions

1097

1098 size: integer

1099 Size of the isolated cells to exclude

1100 (Number of Voxels)

1101

1102 Returns

1103 -------

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

1105 Array with minimum region size > size

1106

1107 Notes

1108 -----

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

1110 """

1111

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

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

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

1115 areaMask = idSizes <= size

1116 filteredArray[areaMask[idRegions]] = 0

1117

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

1119 array = filteredArray * array

1120

1121 return array

1122

1123

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

1125 """

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

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

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

1129

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

1131

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

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

1134

1135 Parameters

1136 ----------

1137 lab : 3D array of integers

1138 Labelled volume, with lab.max() labels

1139

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

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

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

1143

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

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

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

1147

1148 gaussianFilterSigma : float, optional

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

1150 Default = 0.75

1151

1152 minVol : int, optional

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

1154 Default = 256 voxels

1155

1156 Returns

1157 -------

1158 trueSphericity : lab.max() array of floats

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

1160

1161 Notes

1162 -----

1163 Function contributed by Riccardo Rorato (UPC Barcelona)

1164

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

1166 """

1167 import skimage.measure

1168

1169 lab = lab.astype(labelType)

1170

1171 if boundingBoxes is None:

1172 boundingBoxes = spam.label.boundingBoxes(lab)

1173 if centresOfMass is None:

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

1175

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

1177

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

1179

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

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

1182 # Extract grain

1183 GL = spam.label.getLabel(

1184 lab,

1185 label,

1186 boundingBoxes=boundingBoxes,

1187 centresOfMass=centresOfMass,

1188 extractCube=True,

1189 margin=2,

1190 maskOtherLabels=True,

1191 )

1192 # Gaussian smooth

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

1194 # mesh edge

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

1196 # compute surface

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

1198 # compute psi

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

1200 return trueSphericity

1201

1202

1203def convexVolume(

1204 lab,

1205 boundingBoxes=None,

1206 centresOfMass=None,

1207 volumes=None,

1208 nProcesses=nProcessesDefault,

1209 verbose=True,

1210):

1211 """

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

1213 list with the convex volume of each particle.

1214

1215 Parameters

1216 ----------

1217 lab : 3D array of integers

1218 Labelled volume, with lab.max() labels

1219

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

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

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

1223

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

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

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

1227

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

1229 Volumes in format returned by ``volumes``

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

1231

1232 nProcesses : integer (optional, default = nProcessesDefault)

1233 Number of processes for multiprocessing

1234 Default = number of CPUs in the system

1235

1236 verbose : boolean (optional, default = False)

1237 True for printing the evolution of the process

1238 False for not printing the evolution of process

1239

1240 Returns

1241 --------

1242

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

1244

1245 Note

1246 ----

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

1248

1249 """

1250 lab = lab.astype(labelType)

1251

1252 # Compute boundingBoxes if needed

1253 if boundingBoxes is None:

1254 boundingBoxes = spam.label.boundingBoxes(lab)

1255 # Compute centresOfMass if needed

1256 if centresOfMass is None:

1257 centresOfMass = spam.label.centresOfMass(lab)

1258 # Compute volumes if needed

1259 if volumes is None:

1260 volumes = spam.label.volumes(lab)

1261 # Compute number of labels

1262 nLabels = lab.max()

1263

1264 # Result array

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

1266

1267 if verbose:

1268 widgets = [

1269 progressbar.FormatLabel(""),

1270 " ",

1271 progressbar.Bar(),

1272 " ",

1273 progressbar.AdaptiveETA(),

1274 ]

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

1276 pbar.start()

1277 finishedNodes = 0

1278

1279 # Function for convex volume

1280 global computeConvexVolume

1281

1282 def computeConvexVolume(label):

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

1284 subvol = labelI["subvol"]

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

1286 try:

1287 hull = scipy.spatial.ConvexHull(points)

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

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

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

1291 hullIm = numpy.zeros(subvol.shape)

1292 hullIm[out_idx] = 1

1293 hullVol = spam.label.volumes(hullIm)

1294 return label, hullVol[-1]

1295 except Exception:

1296 return label, 0

1297

1298 # Run multiprocessing

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

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

1301 if verbose:

1302 finishedNodes += 1

1303 pbar.update(finishedNodes)

1304 convexVolume[returns[0]] = returns[1]

1305 pool.close()

1306 pool.join()

1307

1308 if verbose:

1309 pbar.finish()

1310

1311 return convexVolume

1312

1313

1314def moveLabels(

1315 lab,

1316 PhiField,

1317 returnStatus=None,

1318 boundingBoxes=None,

1319 centresOfMass=None,

1320 margin=3,

1321 PhiCOM=True,

1322 threshold=0.5,

1323 labelDilate=0,

1324 nProcesses=nProcessesDefault,

1325):

1326 """

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

1328

1329 Parameters

1330 -----------

1331 lab : 3D numpy array

1332 Labelled image

1333

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

1335 Spatial field of Phis

1336

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

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

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

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

1341

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

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

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

1345

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

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

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

1349

1350 margin : int, optional

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

1352 Default = 3

1353

1354 PhiCOM : bool, optional

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

1356 Default = True

1357

1358 threshold : float, optional

1359 Threshold to keep interpolated voxels in the binary image.

1360 Default = 0.5

1361

1362 labelDilate : int, optional

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

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

1365 Default = 0

1366

1367 nProcesses : integer (optional, default = nProcessesDefault)

1368 Number of processes for multiprocessing

1369 Default = number of CPUs in the system

1370

1371 Returns

1372 --------

1373 labOut : 3D numpy array

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

1375

1376 """

1377

1378 # Check for boundingBoxes

1379 if boundingBoxes is None:

1380 boundingBoxes = spam.label.boundingBoxes(lab)

1381 # Check for centresOfMass

1382 if centresOfMass is None:

1383 centresOfMass = spam.label.centresOfMass(lab)

1384 # Create output label image

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

1386 # Get number of labels

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

1388 numberOfLabelsToMove = 0

1389 labelsToMove = []

1390 # Add the labels to move

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

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

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

1394 pass

1395 else: # Add the particles

1396 labelsToMove.append(label)

1397 numberOfLabelsToMove += 1

1398

1399 # Function for moving labels

1400 global funMoveLabels

1401

1402 def funMoveLabels(label):

1403 getLabelReturn = spam.label.getLabel(

1404 lab,

1405 label,

1406 labelDilate=labelDilate,

1407 margin=margin,

1408 boundingBoxes=boundingBoxes,

1409 centresOfMass=centresOfMass,

1410 extractCube=True,

1411 )

1412 # Check that the label exist

1413 if getLabelReturn is not None:

1414 # Get Phi field

1415 Phi = PhiField[label].copy()

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

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

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

1419 # Check that the displacement exist

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

1421 # Just move binary label

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

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

1424 if PhiCOM:

1425 labSubvolDefInterp = spam.DIC.applyPhi(

1426 getLabelReturn["subvol"],

1427 Phi=Phi,

1428 interpolationOrder=1,

1429 PhiCentre=getLabelReturn["centreOfMassREL"],

1430 )

1431 else:

1432 labSubvolDefInterp = spam.DIC.applyPhi(

1433 getLabelReturn["subvol"],

1434 Phi=Phi,

1435 interpolationOrder=1,

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

1437 )

1438

1439 # "death mask"

1440 labSubvolDefMask = labSubvolDefInterp >= threshold

1441

1442 del labSubvolDefInterp

1443 # Get the boundary of the cube

1444 topOfSlice = numpy.array(

1445 [

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

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

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

1449 ]

1450 )

1451

1452 botOfSlice = numpy.array(

1453 [

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

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

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

1457 ]

1458 )

1459

1460 topOfSliceCrop = numpy.array(

1461 [

1462 max(topOfSlice[0], 0),

1463 max(topOfSlice[1], 0),

1464 max(topOfSlice[2], 0),

1465 ]

1466 )

1467 botOfSliceCrop = numpy.array(

1468 [

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

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

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

1472 ]

1473 )

1474 # Update grainSlice with disp

1475 grainSlice = (

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

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

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

1479 )

1480

1481 # Update labSubvolDefMask

1482 labSubvolDefMaskCrop = labSubvolDefMask[

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

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

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

1486 ]

1487 return label, grainSlice, labSubvolDefMaskCrop, 1

1488

1489 # Nan displacement, run away

1490 else:

1491 return label, 0, 0, -1

1492 # Got None from getLabel()

1493 else:

1494 return label, 0, 0, -1

1495

1496 # Create progressbar

1497 widgets = [

1498 progressbar.FormatLabel(""),

1499 " ",

1500 progressbar.Bar(),

1501 " ",

1502 progressbar.AdaptiveETA(),

1503 ]

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

1505 pbar.start()

1506 finishedNodes = 0

1507 # Run multiprocessing

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

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

1510 finishedNodes += 1

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

1512 pbar.update(finishedNodes)

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

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

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

1516 pool.close()

1517 pool.join()

1518

1519 # End progressbar

1520 pbar.finish()

1521

1522 return labOut

1523

1524

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

1526 """

1527 This function erodes a labelled image.

1528

1529 Parameters

1530 -----------

1531 lab : 3D numpy array

1532 Labelled image

1533

1534 erosion : int, optional

1535 Erosion level

1536

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

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

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

1540

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

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

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

1544

1545 nProcesses : integer (optional, default = nProcessesDefault)

1546 Number of processes for multiprocessing

1547 Default = number of CPUs in the system

1548

1549 Returns

1550 --------

1551 erodeImage : 3D numpy array

1552 New labelled image with the eroded labels.

1553

1554 Note

1555 ----

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

1557

1558 """

1559 # Get number of labels

1560 numberOfLabels = lab.max()

1561 # Create the Empty Phi field

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

1563 # Setup Phi as the identity

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

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

1566 # Use moveLabels

1567 erodeImage = spam.label.moveLabels(

1568 lab,

1569 PhiField,

1570 boundingBoxes=boundingBoxes,

1571 centresOfMass=centresOfMass,

1572 margin=1,

1573 PhiCOM=True,

1574 threshold=0.5,

1575 labelDilate=-erosion,

1576 nProcesses=nProcesses,

1577 )

1578 return erodeImage

1579

1580

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

1582 """

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

1584

1585 Parameters

1586 -----------

1587 lab : 3D numpy array

1588 Labelled image

1589

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

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

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

1593

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

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

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

1597

1598 Returns

1599 --------

1600 labOut : 3D numpy array

1601 New labelled image.

1602

1603 Note

1604 ----

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

1606

1607 """

1608

1609 # Check for boundingBoxes

1610 if boundingBoxes is None:

1611 boundingBoxes = spam.label.boundingBoxes(lab)

1612 # Check for centresOfMass

1613 if centresOfMass is None:

1614 centresOfMass = spam.label.centresOfMass(lab)

1615 # Create output label image

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

1617 # Get number of labels

1618 numberOfLabels = lab.max()

1619 # Create progressbar

1620 widgets = [

1621 progressbar.FormatLabel(""),

1622 " ",

1623 progressbar.Bar(),

1624 " ",

1625 progressbar.AdaptiveETA(),

1626 ]

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

1628 pbar.start()

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

1630 # Get label

1631 getLabelReturn = spam.label.getLabel(

1632 lab,

1633 i,

1634 labelDilate=0,

1635 margin=3,

1636 boundingBoxes=boundingBoxes,

1637 centresOfMass=centresOfMass,

1638 maskOtherLabels=False,

1639 )

1640 # Get subvolume

1641 subVol = getLabelReturn["subvol"]

1642 # Transform to binary

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

1644 # Mask out all the other labels

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

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

1647 subVolMaskOtherLabel = subVolBinMask - subVolBinMaskLabel

1648 # Fill holes with convex volume

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

1650 hull = scipy.spatial.ConvexHull(points)

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

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

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

1654 hullIm = numpy.zeros(subVol.shape)

1655 hullIm[out_idx] = 1

1656 hullIm = hullIm > 0

1657 # Identify added voxels

1658 subVolAdded = hullIm - subVolBinMaskLabel

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

1660 subVolWrongAdded = subVolAdded * subVolMaskOtherLabel

1661 # Remove wrong filling areas

1662 subVolCorrect = (hullIm - subVolWrongAdded) > 0

1663 # Get slice

1664 grainSlice = (

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

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

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

1668 )

1669 # Add it to the output file

1670 labOut[grainSlice][subVolCorrect] = i

1671 # Update the progressbar

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

1673 pbar.update(i)

1674

1675 return labOut

1676

1677

1678def getNeighbours(

1679 lab,

1680 listOfLabels,

1681 method="getLabel",

1682 neighboursRange=None,

1683 centresOfMass=None,

1684 boundingBoxes=None,

1685):

1686 """

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

1688

1689 Parameters

1690 -----------

1691 lab : 3D numpy array

1692 Labelled image

1693

1694 listOfLabels : list of ints

1695 List of labels to which the neighbours will be computed

1696

1697 method : string

1698 Method to compute the neighbours.

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

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

1701 Default = 'getLabel'

1702

1703 neighboursRange : int

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

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

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

1707

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

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

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

1711

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

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

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

1715

1716 Returns

1717 --------

1718 neighbours : list

1719 List with the neighbours for each label in listOfLabels.

1720

1721 """

1722 # Create result list

1723 neighbours = []

1724 # Compute centreOfMass if needed

1725 if centresOfMass is None:

1726 centresOfMass = spam.label.centresOfMass(lab)

1727 # Compute boundingBoxes if needed

1728 if boundingBoxes is None:

1729 boundingBoxes = spam.label.boundingBoxes(lab)

1730 # Compute Radii

1731 radii = spam.label.equivalentRadii(lab)

1732 if method == "getLabel":

1733 # Compute neighboursRange if needed

1734 if neighboursRange is None:

1735 neighboursRange = numpy.mean(radii)

1736 # Compute for each label in the list of labels

1737 for label in listOfLabels:

1738 getLabelReturn = spam.label.getLabel(

1739 lab,

1740 label,

1741 labelDilate=neighboursRange,

1742 margin=neighboursRange,

1743 boundingBoxes=boundingBoxes,

1744 centresOfMass=centresOfMass,

1745 maskOtherLabels=False,

1746 )

1747 # Get subvolume

1748 subVol = getLabelReturn["subvol"]

1749 # Get neighbours

1750 neighboursLabel = numpy.unique(subVol)

1751 # Remove label and 0 from the list of neighbours

1752 neighboursLabel = neighboursLabel[~numpy.in1d(neighboursLabel, label)]

1753 neighboursLabel = neighboursLabel[~numpy.in1d(neighboursLabel, 0)]

1754 # Add the neighbours to the list

1755 neighbours.append(neighboursLabel)

1756

1757 elif method == "mesh":

1758 # Compute neighboursRange if needed

1759 if neighboursRange is None:

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

1761 # Get connectivity matrix

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

1763 # Compute for each label in the list of labels

1764 for label in listOfLabels:

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

1766 # Remove label from the list of neighbours

1767 neighboursLabel = neighboursLabel[~numpy.in1d(neighboursLabel, label)]

1768 # Add the neighbours to the list

1769 neighbours.append(neighboursLabel)

1770 else:

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

1772

1773 return neighbours

1774

1775

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

1777 """

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

1779

1780 Parameters

1781 -----------

1782 lab : 3D numpy array

1783 Labelled image

1784

1785 nProcesses : integer (optional, default = nProcessesDefault)

1786 Number of processes for multiprocessing

1787 Default = number of CPUs in the system

1788

1789 verbose : boolean (optional, default = False)

1790 True for printing the evolution of the process

1791

1792 Returns

1793 --------

1794 underSegCoeff : lab.max() array of floats

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

1796

1797 Note

1798 ----

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

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

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

1802

1803 """

1804 # Compute the volume

1805 vol = spam.label.volumes(lab)

1806 # Compute the convex volume

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

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

1809 vol[0] = 1

1810 # Compute the underSegmentation Coefficient

1811 underSegCoeff = convexVol / vol

1812 # Set the coefficient of the void to 0

1813 underSegCoeff[0] = 0

1814 return underSegCoeff

1815

1816

1817def detectOverSegmentation(lab):

1818 """

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

1820 and a characteristic length of the particle.

1821

1822 Parameters

1823 -----------

1824 lab : 3D numpy array

1825 Labelled image

1826

1827 Returns

1828 --------

1829 overSegCoeff : lab.max() array of floats

1830 An array of float values with the oversegmentation coefficient

1831

1832 sharedLabel : lab.max() array of floats

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

1834

1835 Note

1836 ----

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

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

1839

1840 """

1841 # Get the labels

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

1843 # Compute the volumes

1844 vol = spam.label.volumes(lab)

1845 # Compute the eq diameter

1846 eqDiam = spam.label.equivalentRadii(lab)

1847 # Compute the areas

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

1849 # Create result list

1850 overSegCoeff = []

1851 sharedLabel = []

1852 for label in labels:

1853 if label == 0:

1854 overSegCoeff.append(0)

1855 sharedLabel.append(0)

1856 else:

1857 # Check if there are contacting areas and volumes

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

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

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

1861 # Get the label for the max contacting area

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

1863 # Compute the coefficient

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

1865 # Add the label

1866 sharedLabel.append(maxLabel)

1867 else:

1868 overSegCoeff.append(0)

1869 sharedLabel.append(0)

1870 overSegCoeff = numpy.array(overSegCoeff)

1871 sharedLabel = numpy.array(sharedLabel)

1872 return overSegCoeff, sharedLabel

1873

1874

1875def fixUndersegmentation(

1876 lab,

1877 imGrey,

1878 targetLabels,

1879 underSegCoeff,

1880 boundingBoxes=None,

1881 centresOfMass=None,

1882 imShowProgress=False,

1883 verbose=True,

1884):

1885 """

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

1887

1888 Parameters

1889 -----------

1890 lab : 3D numpy array

1891 Labelled image

1892

1893 imGrey : 3D numpy array

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

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

1896

1897 targetLabels : int or a list of labels

1898 List of target labels to solve undersegmentation

1899

1900 underSegCoeff : lab.max() array of floats

1901 Undersegmentation coefficient as returned by ``detectUnderSegmentation``

1902

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

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

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

1906

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

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

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

1910 Labelled image

1911

1912 imShowProgress : bool, optional

1913 Graphical interface to observe the process for each label.

1914 Default = False

1915

1916 verbose : boolean (optional, default = False)

1917 True for printing the evolution of the process

1918

1919 Returns

1920 --------

1921 lab : 3D numpy array

1922 Labelled image after running ``makeLabelsSequential``

1923 """

1924

1925 # Usual checks

1926 if boundingBoxes is None:

1927 boundingBoxes = spam.label.boundingBoxes(lab)

1928 if centresOfMass is None:

1929 centresOfMass = spam.label.centresOfMass(lab)

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

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

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

1933 return

1934 # Start counters

1935 labelCounter = numpy.max(lab)

1936 labelDummy = numpy.zeros(lab.shape)

1937 successCounter = 0

1938 finishedLabels = 0

1939 if verbose:

1940 widgets = [

1941 progressbar.FormatLabel(""),

1942 " ",

1943 progressbar.Bar(),

1944 " ",

1945 progressbar.AdaptiveETA(),

1946 ]

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

1948 pbar.start()

1949 # Main loop

1950 for label in targetLabels:

1951 # Get the subset

1952 labelData = spam.label.getLabel(

1953 lab,

1954 label,

1955 margin=5,

1956 boundingBoxes=boundingBoxes,

1957 centresOfMass=centresOfMass,

1958 extractCube=True,

1959 )

1960 # Get the slice on the greyscale

1961 imGreySlice = imGrey[

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

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

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

1965 ]

1966 # Mask the imGreySubset

1967 greySubset = imGreySlice * labelData["subvol"]

1968 # Create seeds

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

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

1971 # Do we have seeds?

1972 if numpy.max(seeds) < 1:

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

1974 # We shouldn't change this label

1975 passBool = "Decline"

1976 else:

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

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

1979 # Run again the underSegCoeff for the subset

1980 res = detectUnderSegmentation(imLabSubset, verbose=False)

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

1982 if len(res) > 2:

1983 # We have at least one label

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

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

1986 # We can modify this label

1987 passBool = "Accept"

1988 successCounter += 1

1989 # Remove the label from the original label image

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

1991 # Assign the new labels to the grains

1992 # Create a subset to fill with the new labels

1993 imLabSubsetNew = numpy.zeros(imLabSubset.shape)

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

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

1996 labelCounter += 1

1997 # Create a disposable dummy sample to allocate the grains

1998 labelDummyUnit = numpy.zeros(lab.shape)

1999 # Alocate the grains

2000 labelDummyUnit[

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

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

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

2004 ] = imLabSubsetNew

2005 # Add the grains

2006 labelDummy = labelDummy + labelDummyUnit

2007 else:

2008 # We shouldn't change this label

2009 passBool = "Decline"

2010 else:

2011 # We shouldn't change this label

2012 passBool = "Decline"

2013 if imShowProgress:

2014 # Enter graphical mode

2015 # Change the labels to show different colourss

2016 fig = plt.figure()

2017 # Plot

2018 plt.subplot(3, 2, 1)

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

2020 plt.imshow(

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

2022 cmap="Greys_r",

2023 )

2024 plt.subplot(3, 2, 2)

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

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

2027 plt.subplot(3, 2, 3)

2028 plt.imshow(

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

2030 cmap="Greys_r",

2031 )

2032 plt.subplot(3, 2, 4)

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

2034 plt.subplot(3, 2, 5)

2035 plt.imshow(

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

2037 cmap="Greys_r",

2038 )

2039 plt.subplot(3, 2, 6)

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

2041 fig.suptitle(

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

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

2044 fontsize="xx-large",

2045 )

2046 plt.show()

2047 if verbose:

2048 finishedLabels += 1

2049 pbar.update(finishedLabels)

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

2051 lab = lab + labelDummy

2052 # Update the labels

2053 lab = spam.label.makeLabelsSequential(lab)

2054 if verbose:

2055 pbar.finish()

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

2057 return lab

2058

2059

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

2061 """

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

2063

2064 Parameters

2065 -----------

2066 lab : 3D numpy array

2067 Labelled image

2068

2069 targetLabels : int or a list of labels

2070 List of target labels to solve oversegmentation

2071

2072 sharedLabel : lab.max() array of floats

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

2074

2075 imShowProgress : bool, optional

2076 Graphical interface to observe the process for each label.

2077 Default = False

2078

2079 verbose : boolean (optional, default = False)

2080 True for printing the evolution of the process

2081

2082 Returns

2083 --------

2084 lab : 3D numpy array

2085 Labelled image after running ``makeLabelsSequential``

2086

2087 """

2088 # Start counters

2089 labelDummy = numpy.zeros(lab.shape)

2090 finishedLabelsCounter = 0

2091 finishedLabels = []

2092 if verbose:

2093 widgets = [

2094 progressbar.FormatLabel(""),

2095 " ",

2096 progressbar.Bar(),

2097 " ",

2098 progressbar.AdaptiveETA(),

2099 ]

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

2101 pbar.start()

2102 # Main loop

2103 for labelA in targetLabels:

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

2105 if labelA in finishedLabels:

2106 # It is already on the list, move on

2107 pass

2108 else:

2109 # Get the touching label

2110 labelB = sharedLabel[labelA]

2111 # Add then to the list

2112 finishedLabels.append(labelA)

2113 finishedLabels.append(labelB)

2114 # Get the subset of the two labels

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

2116 # Change the labelB by labelA in the subset

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

2118 # Create a disposable dummy sample to allocate the grains

2119 labelDummyUnit = numpy.zeros(lab.shape)

2120 # Alocate the grains

2121 labelDummyUnit[

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

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

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

2125 ] = subVolLabNew

2126 # Add the grains

2127 labelDummy = labelDummy + labelDummyUnit

2128 # Remove the label from the original label image

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

2130 # Enter graphical mode

2131 if imShowProgress:

2132 # Change the labels to show different colourss

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

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

2135 fig = plt.figure()

2136 plt.subplot(3, 2, 1)

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

2138 plt.imshow(

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

2140 cmap="cubehelix",

2141 )

2142 plt.subplot(3, 2, 2)

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

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

2145 plt.subplot(3, 2, 3)

2146 plt.imshow(

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

2148 cmap="cubehelix",

2149 )

2150 plt.subplot(3, 2, 4)

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

2152 plt.subplot(3, 2, 5)

2153 plt.imshow(

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

2155 cmap="cubehelix",

2156 )

2157 plt.subplot(3, 2, 6)

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

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

2160 plt.show()

2161 if verbose:

2162 finishedLabelsCounter += 1

2163 pbar.update(finishedLabelsCounter)

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

2165 lab = lab + labelDummy

2166 # Update the labels

2167 lab = spam.label.makeLabelsSequential(lab)

2168 if verbose:

2169 pbar.finish()

2170

2171 return lab