Coverage for /usr/local/lib/python3.10/site-packages/spam/label/contacts.py: 89%
371 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-24 17:26 +0000
1"""
2Library of SPAM functions for dealing with contacts between particles
3Copyright (C) 2020 SPAM Contributors
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
21try:
22 multiprocessing.set_start_method("fork")
23except RuntimeError:
24 pass
26import numpy
27import progressbar
28import scipy.ndimage # for the contact detection
29import skimage.feature # for the maxima of the EDM
30import spam.label
31from spambind.label.labelToolkit import labelContacts
33contactType = "<u4"
34# Global number of processes
35nProcessesDefault = multiprocessing.cpu_count()
38def contactingLabels(lab, labelsList=None, areas=False, boundingBoxes=None, centresOfMass=None):
39 """
40 This function returns contacting labels for a given label or list of labels,
41 and optionally the number of voxels involved in the contact.
42 This is designed for an interpixel watershed where there are no watershed lines,
43 and so contact detection is done by dilating the particle in question once and
44 seeing what other labels in comes in contact with.
46 Parameters
47 -----------
48 lab : 3D numpy array of ints ( of type spam.label.toolkit.labelType)
49 An array of labels with 0 as the background
51 labels : int or a list of labels
52 Labels for which we should compute neighbours
54 areas : bool, optional (default = False)
55 Compute contact "areas"? *i.e.*, number of voxels
56 Careful, this is not a physically correct quantity, and
57 is measured only as the overlap from ``label`` to its
58 neightbours. See ``c ontactPoint`` for something
59 more meaningful
61 boundingBoxes : lab.max()x6 array of ints, optional
62 Bounding boxes in format returned by ``boundingBoxes``.
63 If not defined (Default = None), it is recomputed by running ``boundingBoxes``
65 centresOfMass : lab.max()x3 array of floats, optional
66 Centres of mass in format returned by ``centresOfMass``.
67 If not defined (Default = None), it is recomputed by running ``centresOfMass``
69 Returns
70 --------
71 For a single "labels", a list of contacting labels.
72 if areas == True, then a list of "areas" is also returned.
74 For multiple labels, as above but a lsit of lists in the order of the labels
76 Note
77 -----
78 Because of dilation, this function is not very safe at the edges,
79 and will throw an exception
80 """
81 # Default state for single
82 single = False
84 if boundingBoxes is None:
85 boundingBoxes = spam.label.boundingBoxes(lab)
86 if centresOfMass is None:
87 centresOfMass = spam.label.centresOfMass(lab, boundingBoxes=boundingBoxes)
89 if labelsList is None:
90 labelsList = list(range(0, lab.max() + 1))
92 # I guess there's only one label
93 if type(labelsList) != list:
94 labelsList = [labelsList]
95 single = True
97 contactingLabels = []
98 contactingAreas = []
100 for label in labelsList:
101 # catch zero and return it in order to keep arrays aligned
102 if label == 0:
103 contactingLabels.append([0])
104 if areas:
105 contactingAreas.append([0])
106 else:
107 p1 = spam.label.getLabel(
108 lab,
109 label,
110 boundingBoxes=boundingBoxes,
111 centresOfMass=centresOfMass,
112 margin=2,
113 )
114 p2 = spam.label.getLabel(
115 lab,
116 label,
117 boundingBoxes=boundingBoxes,
118 centresOfMass=centresOfMass,
119 margin=2,
120 labelDilate=1,
121 )
123 dilOnly = numpy.logical_xor(p2["subvol"], p1["subvol"])
125 labSlice = lab[
126 p1["slice"][0].start : p1["slice"][0].stop,
127 p1["slice"][1].start : p1["slice"][1].stop,
128 p1["slice"][2].start : p1["slice"][2].stop,
129 ]
131 if dilOnly.shape == labSlice.shape:
132 intersection = dilOnly * labSlice
134 counts = numpy.unique(intersection, return_counts=True)
136 # Print counts starting from the 1th column since we'll always have a ton of zeros
137 # print "\tLabels:\n\t",counts[0][1:]
138 # print "\tCounts:\n\t",counts[1][1:]
140 contactingLabels.append(counts[0][1:])
141 if areas:
142 contactingAreas.append(counts[1][1:])
144 else:
145 # The particle is near the edge, pad it
147 padArray = numpy.zeros((3, 2)).astype(int)
148 sliceArray = numpy.zeros((3, 2)).astype(int)
149 for i, sl in enumerate(p1["slice"]):
150 if sl.start < 0:
151 padArray[i, 0] = -1 * sl.start
152 sliceArray[i, 0] = 0
153 else:
154 sliceArray[i, 0] = sl.start
155 if sl.stop > lab.shape[i]:
156 padArray[i, 1] = sl.stop - lab.shape[i]
157 sliceArray[i, 1] = lab.shape[i]
158 else:
159 sliceArray[i, 1] = sl.stop
160 labSlice = lab[
161 sliceArray[0, 0] : sliceArray[0, 1],
162 sliceArray[1, 0] : sliceArray[1, 1],
163 sliceArray[2, 0] : sliceArray[2, 1],
164 ]
165 labSlicePad = numpy.pad(labSlice, padArray)
166 intersection = dilOnly * labSlicePad
167 counts = numpy.unique(intersection, return_counts=True)
168 contactingLabels.append(counts[0][1:])
169 if areas:
170 contactingAreas.append(counts[1][1:])
172 # Flatten if it's a list with only one object
173 if single:
174 contactingLabels = contactingLabels[0]
175 if areas:
176 contactingAreas = contactingAreas[0]
177 # Now return things
178 if areas:
179 return contactingLabels, contactingAreas
180 else:
181 return contactingLabels
184def contactPoints(
185 lab,
186 contactPairs,
187 returnContactZones=False,
188 boundingBoxes=None,
189 centresOfMass=None,
190 verbose=False,
191):
192 """
193 Get the point, or area of contact between contacting particles.
194 This is done by looking at the overlap of a 1-dilation of particle1 onto particle 2
195 and vice versa.
197 Parameters
198 -----------
199 lab : 3D numpy array of ints ( of type spam.label.toolkit.labelType)
200 An array of labels with 0 as the background
202 contactPairs : list (or list of list) of labels
203 Pairs of labels for which we should compute neighbours
205 returnContactZones : bool, optional (default = False)
206 Output a labelled volume where contact zones are labelled according to the order
207 of contact pairs?
208 If False, the centres of mass of the contacts are returned
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``
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``
218 Returns
219 --------
220 if returnContactZones:
221 lab : 3D numpy array of ints ( of type spam.label.toolkit.labelType)
222 else:
223 Z Y X Centres of mass of contacts
224 """
225 # detect list of lists with code from: https://stackoverflow.com/questions/5251663/determine-if-a-list-contains-other-lists
226 # if not any(isinstance(el, list) for el in contactPairs ):
227 # contactPairs = [ contactPairs ]
228 # different approach:
229 contactPairs = numpy.array(contactPairs)
230 # Pad with extra dim if only one contact pair
231 if len(contactPairs.shape) == 1:
232 contactPairs = contactPairs[numpy.newaxis, ...]
234 if boundingBoxes is None:
235 boundingBoxes = spam.label.boundingBoxes(lab)
236 if centresOfMass is None:
237 centresOfMass = spam.label.centresOfMass(lab, boundingBoxes=boundingBoxes)
239 # To this volume will be added each contact "area", labelled
240 # with the contact pair number (+1)
241 analysisVolume = numpy.zeros_like(lab, dtype=numpy.uint32)
242 if verbose:
243 widgets = [
244 progressbar.FormatLabel(""),
245 " ",
246 progressbar.Bar(),
247 " ",
248 progressbar.AdaptiveETA(),
249 ]
250 pbar = progressbar.ProgressBar(widgets=widgets, maxval=len(contactPairs) + 1)
251 pbar.start()
252 finishedNodes = 0
253 for n, contactPair in enumerate(contactPairs):
254 # Do both orders in contacts
255 if len(contactPair) != 2:
256 print("spam.label.contacts.contactPoints(): contact pair does not have 2 elements")
257 for label1, label2 in [contactPair, contactPair[::-1]]:
258 # print( "Label1: {} Label2: {}".format( label1, label2 ) )
259 p1 = spam.label.getLabel(
260 lab,
261 label1,
262 boundingBoxes=boundingBoxes,
263 centresOfMass=centresOfMass,
264 margin=2,
265 )
266 p2 = spam.label.getLabel(
267 lab,
268 label1,
269 boundingBoxes=boundingBoxes,
270 centresOfMass=centresOfMass,
271 margin=2,
272 labelDilate=1,
273 )
275 dilOnly = numpy.logical_xor(p2["subvol"], p1["subvol"])
277 labSlice = (
278 slice(p1["slice"][0].start, p1["slice"][0].stop),
279 slice(p1["slice"][1].start, p1["slice"][1].stop),
280 slice(p1["slice"][2].start, p1["slice"][2].stop),
281 )
283 labSubvol = lab[labSlice]
285 if dilOnly.shape == labSubvol.shape:
286 intersection = dilOnly * labSubvol
288 analysisVolume[labSlice][intersection == label2] = n + 1
290 else:
291 raise Exception("\tdilOnly and labSlice are not the same size... getLabel probably did some safe cutting around the edges")
292 if verbose:
293 finishedNodes += 1
294 widgets[0] = progressbar.FormatLabel("{}/{} ".format(finishedNodes, len(contactPairs)))
295 pbar.update(finishedNodes)
297 if returnContactZones:
298 return analysisVolume
299 else:
300 return spam.label.centresOfMass(analysisVolume)
303def labelledContacts(lab, maximumCoordinationNumber=20):
304 """
305 Uniquely names contacts based on grains.
307 Parameters
308 -----------
309 lab : 3D numpy array of ints ( of type spam.label.toolkit.labelType)
310 An array of labels with 0 as the background
312 maximumCoordinationNumber : int (optional, default = 20)
313 Maximum number of neighbours to consider per grain
315 Returns
316 --------
317 An array, containing:
318 contactVolume : array of ints
319 3D array where contact zones are uniquely labelled
321 Z : array of ints
322 Vector of length lab.max()+1 contaning the coordination number
323 (number of touching labels for each label)
325 contactTable : array of ints
326 2D array of lab.max()+1 by 2*maximumCoordinationNumber
327 containing, per grain: touching grain label, contact number
329 contactingLabels : array of ints
330 2D array of numberOfContacts by 2
331 containing, per contact: touching grain label 1, touching grain label 2
332 """
333 contacts = numpy.zeros_like(lab, dtype=contactType)
334 Z = numpy.zeros((lab.max() + 1), dtype="<u1")
335 contactTable = numpy.zeros((lab.max() + 1, maximumCoordinationNumber * 2), dtype=contactType)
336 contactingLabels = numpy.zeros(((lab.max() + 1) * maximumCoordinationNumber, 2), dtype=spam.label.labelType)
338 labelContacts(lab, contacts, Z, contactTable, contactingLabels)
340 # removed the first row of contactingLabels, as it is 0, 0
341 return [contacts, Z, contactTable, contactingLabels[1 : contacts.max() + 1, :]]
344def fetchTwoGrains(volLab, labels, volGrey=None, boundingBoxes=None, padding=0, size_exclude=5):
345 """
346 Fetches the sub-volume of two grains from a labelled image
348 Parameters
349 ----------
350 volLab : 3D array of integers
351 Labelled volume, with lab.max() labels
353 labels : 1x2 array of integers
354 the two labels that should be contained in the subvolume
356 volGrey : 3D array
357 Grey-scale volume
359 boundingBoxes : lab.max()x6 array of ints, optional
360 Bounding boxes in format returned by ``boundingBoxes``.
361 If not defined (Default = None), it is recomputed by running ``boundingBoxes``
363 padding : integer
364 padding of the subvolume
365 for some purpose it might be benefitial to have a subvolume
366 with a boundary of zeros
368 Returns
369 -------
370 subVolLab : 3D array of integers
371 labelled sub-volume containing the two input labels
373 subVolBin : 3D array of integers
374 binary subvolume
376 subVolGrey : 3D array
377 grey-scale subvolume
378 """
380 # check if bounding boxes are given
381 if boundingBoxes is None:
382 # print "bounding boxes are not given. calculating ...."
383 boundingBoxes = spam.label.boundingBoxes(volLab)
384 # else:
385 # print "bounding boxes are given"
386 lab1, lab2 = labels
387 # Define output dictionary since we'll add different things to it
388 output = {}
389 # get coordinates of the big bounding box
390 startZ = min(boundingBoxes[lab1, 0], boundingBoxes[lab2, 0]) - padding
391 stopZ = max(boundingBoxes[lab1, 1], boundingBoxes[lab2, 1]) + padding
392 startY = min(boundingBoxes[lab1, 2], boundingBoxes[lab2, 2]) - padding
393 stopY = max(boundingBoxes[lab1, 3], boundingBoxes[lab2, 3]) + padding
394 startX = min(boundingBoxes[lab1, 4], boundingBoxes[lab2, 4]) - padding
395 stopX = max(boundingBoxes[lab1, 5], boundingBoxes[lab2, 5]) + padding
396 output["slice"] = (
397 slice(startZ, stopZ + 1),
398 slice(startY, stopY + 1),
399 slice(startX, stopX + 1),
400 )
401 subVolLab = volLab[
402 output["slice"][0].start : output["slice"][0].stop,
403 output["slice"][1].start : output["slice"][1].stop,
404 output["slice"][2].start : output["slice"][2].stop,
405 ]
407 subVolLab_A = numpy.where(subVolLab == lab1, lab1, 0)
408 subVolLab_B = numpy.where(subVolLab == lab2, lab2, 0)
409 subVolLab = subVolLab_A + subVolLab_B
411 struc = numpy.zeros((3, 3, 3))
412 struc[1, 1:2, :] = 1
413 struc[1, :, 1:2] = 1
414 struc[0, 1, 1] = 1
415 struc[2, 1, 1] = 1
416 subVolLab = spam.label.filterIsolatedCells(subVolLab, struc, size_exclude)
417 output["subVolLab"] = subVolLab
418 subVolBin = numpy.where(subVolLab != 0, 1, 0)
419 output["subVolBin"] = subVolBin
420 if volGrey is not None:
421 subVolGrey = volGrey[
422 output["slice"][0].start : output["slice"][0].stop,
423 output["slice"][1].start : output["slice"][1].stop,
424 output["slice"][2].start : output["slice"][2].stop,
425 ]
426 subVolGrey = subVolGrey * subVolBin
427 output["subVolGrey"] = subVolGrey
429 return output
432def localDetection(subVolGrey, localThreshold, radiusThresh=None):
433 """
434 local contact refinement
435 checks whether two particles are in contact with a local threshold,
436 that is higher than the global one used for binarisation
438 Parameters
439 ----------
440 subVolGrey : 3D array
441 Grey-scale volume
443 localThreshold : integer or float, same type as the 3D array
444 threshold for binarisation of the subvolume
446 radiusThresh : integer, optional
447 radius for excluding patches that might not be relevant,
448 e.g. noise can lead to patches that are not connected with the grains in contact
449 the patches are calculated with ``equivalentRadii()``
450 Default is None and such patches are not excluded from the subvolume
452 Returns
453 -------
454 CONTACT : boolean
455 if True, that particles appear in contact for the local threshold
457 Note
458 ----
459 see https://doi.org/10.1088/1361-6501/aa8dbf for further information
460 """
462 CONTACT = False
463 subVolBin = ((subVolGrey > localThreshold) * 1).astype("uint8")
464 if radiusThresh is not None:
465 # clean the image of isolated voxels or pairs of voxels
466 # due to higher threshold they could exist
467 subVolLab, numObjects = scipy.ndimage.label(subVolBin)
468 if numObjects > 1:
469 radii = spam.label.equivalentRadii(subVolLab)
470 labelsToRemove = numpy.where(radii < radiusThresh)
471 if len(labelsToRemove[0]) > 1:
472 subVolLab = spam.label.removeLabels(subVolLab, labelsToRemove)
473 subVolBin = ((subVolLab > 0) * 1).astype("uint8")
475 # fill holes
476 subVolBin = scipy.ndimage.binary_fill_holes(subVolBin).astype("uint8")
477 labeledArray, numObjects = scipy.ndimage.label(subVolBin, structure=None, output=None)
478 if numObjects == 1:
479 CONTACT = True
481 return CONTACT
484def contactOrientations(volBin, volLab, watershed="ITK", peakDistance=1, markerShrink=True, verbose=False):
485 """
486 determines contact normal orientation between two particles
487 uses the random walker implementation from skimage
489 Parameters
490 ----------
491 volBin : 3D array
492 binary volume containing two particles in contact
494 volLab : 3D array of integers
495 labelled volume containing two particles in contact
497 watershed : string
498 sets the basis for the determination of the orientation
499 options are "ITK" for the labelled image from the input,
500 "RW" for a further segmentation by the random walker
501 Default = "ITK"
503 peakDistance : int, optional
504 Distance in pixels used in skimage.feature.peak_local_max
505 Default value is 1
507 markerShrink : boolean
508 shrinks the markers to contain only one voxel.
509 For large markers, the segmentation might yield strange results depending on the shape.
510 Default = True
512 verbose : boolean (optional, default = False)
513 True for printing the evolution of the process
514 False for not printing the evolution of process
517 Returns
518 -------
519 contactNormal : 1x3 array
520 contact normal orientation in z,y,x coordinates
521 the vector is flipped such that the z coordinate is positive
523 len(coordIntervox) : integer
524 the number of points for the principal component analysis
525 indicates the quality of the fit
527 notTreatedContact : boolean
528 if False that contact is not determined
529 if the contact consisted of too few voxels a fit is not possible or feasible
530 """
531 notTreatedContact = False
532 from skimage.segmentation import random_walker
534 if watershed == "ITK":
535 # ------------------------------
536 # ITK watershed for orientations
537 # ------------------------------
538 # need to relabel, because the _contactPairs functions looks for 1 and 2
539 labels = numpy.unique(volLab)
540 volLab[volLab == labels[1]] = 1
541 volLab[volLab == labels[2]] = 2
543 contactITK = _contactPairs(volLab)
545 if len(contactITK) <= 2:
546 # print('WARNING: not enough contacting voxels (beucher)... aborting this calculation')
547 notTreatedContact = True
548 return numpy.zeros(3), 0, notTreatedContact
550 coordIntervox = contactITK[:, 0:3]
551 # Determining the contact orientation! using PCA
552 contactNormal = _contactNormals(coordIntervox)
554 if watershed == "RW":
555 # Get Labels
556 label1, label2 = numpy.unique(volLab)[1:]
557 # Create sub-domains
558 subVolBinA = numpy.where(volLab == label1, 1, 0)
559 subVolBinB = numpy.where(volLab == label2, 1, 0)
560 volBin = subVolBinA + subVolBinB
561 # Calculate distance map
562 distanceMapA = scipy.ndimage.distance_transform_edt(subVolBinA)
563 distanceMapB = scipy.ndimage.distance_transform_edt(subVolBinB)
564 # Calculate local Max
565 localMaxiPointsA = skimage.feature.peak_local_max(
566 distanceMapA,
567 min_distance=peakDistance,
568 exclude_border=False,
569 labels=subVolBinA,
570 )
571 localMaxiPointsB = skimage.feature.peak_local_max(
572 distanceMapB,
573 min_distance=peakDistance,
574 exclude_border=False,
575 labels=subVolBinB,
576 )
577 # 2022-09-26: Dropping indices, and so rebuilding image of peaks
578 localMaxiA = numpy.zeros_like(distanceMapA, dtype=bool)
579 localMaxiB = numpy.zeros_like(distanceMapB, dtype=bool)
580 localMaxiA[tuple(localMaxiPointsA.T)] = True
581 localMaxiB[tuple(localMaxiPointsB.T)] = True
582 # Label the Peaks
583 struc = numpy.ones((3, 3, 3))
584 markersA, numMarkersA = scipy.ndimage.label(localMaxiA, structure=struc)
585 markersB, numMarkersB = scipy.ndimage.label(localMaxiB, structure=struc)
586 if numMarkersA == 0 or numMarkersB == 0:
587 notTreatedContact = True
588 if verbose:
589 print("spam.contacts.contactOrientations: No grain markers found for this contact. Setting the contact as non-treated")
590 return numpy.zeros(3), 0, notTreatedContact
591 # Shrink the markers
592 # The index property should be a list with all the labels encountered, exlucing the label for the backgroun (0).
593 centerOfMassA = scipy.ndimage.center_of_mass(markersA, labels=markersA, index=list(numpy.unique(markersA)[1:]))
594 centerOfMassB = scipy.ndimage.center_of_mass(markersB, labels=markersB, index=list(numpy.unique(markersB)[1:]))
595 # Calculate the value of the distance map for the markers
596 marksValA = []
597 marksValB = []
598 for x in range(len(centerOfMassA)):
599 marksValA.append(
600 distanceMapA[
601 int(centerOfMassA[x][0]),
602 int(centerOfMassA[x][1]),
603 int(centerOfMassA[x][2]),
604 ]
605 )
606 for x in range(len(centerOfMassB)):
607 marksValB.append(
608 distanceMapB[
609 int(centerOfMassB[x][0]),
610 int(centerOfMassB[x][1]),
611 int(centerOfMassB[x][2]),
612 ]
613 )
614 # Sort the markers coordinates acoording to their distance map value
615 # First element correspond to the coordinates of the marker with the higest value in the distance map
616 centerOfMassA = numpy.flipud([x for _, x in sorted(zip(marksValA, centerOfMassA))])
617 centerOfMassB = numpy.flipud([x for _, x in sorted(zip(marksValB, centerOfMassB))])
618 marksValA = numpy.flipud(sorted(marksValA))
619 marksValB = numpy.flipud(sorted(marksValB))
620 # Select markers
621 markers = numpy.zeros(volLab.shape)
622 markers[int(centerOfMassA[0][0]), int(centerOfMassA[0][1]), int(centerOfMassA[0][2])] = 1.0
623 markers[int(centerOfMassB[0][0]), int(centerOfMassB[0][1]), int(centerOfMassB[0][2])] = 2.0
624 # set the void phase of the binary image to -1, excludes this phase from calculation of the random walker (saves time)
625 volRWbin = volBin.astype(float)
626 markers[~volRWbin.astype(bool)] = -1
627 if numpy.max(numpy.unique(markers)) != 2:
628 notTreatedContact = True
629 if verbose:
630 print("spam.contacts.contactOrientations: The grain markers where wrongly computed. Setting the contact as non-treated")
631 return numpy.zeros(3), 0, notTreatedContact
632 probMaps = random_walker(volRWbin, markers, beta=80, mode="cg_mg", return_full_prob=True)
633 # reset probability of voxels in the void region to 0! (-1 right now, as they were not taken into account!)
634 probMaps[:, ~volRWbin.astype(bool)] = 0
635 # create a map of the probabilities of belonging to either label 1 oder label 2 (which is just the inverse map)
636 labRW1 = probMaps[0].astype(numpy.float32)
637 # label the image depending on the probabilty of each voxel belonging to marker = 1, save one power watershed
638 labRW = numpy.array(labRW1)
639 labRW[labRW > 0.5] = 1
640 labRW[(labRW < 0.5) & (labRW > 0)] = 2
641 # seed of the second marker has to be labelled by 2 as well
642 labRW[markers == 2] = 2
644 contactVox = _contactPairs(labRW)
645 if len(contactVox) <= 2:
646 # print 'WARNING: not enough contacting voxels (rw).... aborting this calculation'
647 notTreatedContact = True
648 if verbose:
649 print("spam.contacts.contactOrientations: Not enough contacting voxels, aborting this calculation. Setting the contact as non-treated")
650 return numpy.zeros(3), 0, notTreatedContact
652 # get subvoxel precision - using the probability values at the contacting voxels!
653 coordIntervox = _contactPositions(contactVox, labRW1)
654 # Determining the contact orientation! using PCA
655 contactNormal = _contactNormals(coordIntervox)
657 return contactNormal[0], len(coordIntervox), notTreatedContact
660def _contactPairs(lab):
661 """
662 finds the voxels involved in the contact
663 i.e. the voxels of one label in direct contact with another label
665 Parameters
666 ----------
667 lab : 3D array of integers
668 labelled volume containing two particles in contact
669 labels of the considered particles are 1 and 2
671 Returns
672 -------
673 contact : (len(contactVoxels))x4 array
674 z,y,x positions and the label of the voxel
675 """
676 dimensions = numpy.shape(lab)
677 dimZ, dimY, dimX = dimensions[0], dimensions[1], dimensions[2]
679 # position of the voxels labelled by 1 ... row 1 = coord index 1, row 2 = coord 2, ...
680 positionLabel = numpy.array(numpy.where(lab == 1))
681 # initialize the array for the contact positions and labels!
682 contact = numpy.array([[], [], [], []]).transpose()
684 # loop over all voxels labeled with 1
685 for x in range(0, len(positionLabel.transpose())):
686 pix = numpy.array([positionLabel[0][x], positionLabel[1][x], positionLabel[2][x], 1])
687 # pix ... positions (axis 0,1,2) of the first voxel containing 1
688 # x axis neighbor
689 if positionLabel[0][x] < dimZ - 1: # if the voxel is still in the image!
690 # if neighbor in pos. x-direction is 2 then save the actual voxel and the neighbor!!
691 if lab[positionLabel[0][x] + 1, positionLabel[1][x], positionLabel[2][x]] == 2:
692 vx1 = numpy.array(
693 [
694 positionLabel[0][x] + 1,
695 positionLabel[1][x],
696 positionLabel[2][x],
697 2,
698 ]
699 )
700 # stacks the data, adds a row to the data, so you get (contact, pix, vx2)
701 contact = numpy.vstack((contact, pix))
702 contact = numpy.vstack((contact, vx1))
703 # this condition prevents from getting stuck at the border of the image, where the x-value cannot be reduced!
704 if positionLabel[0][x] != 0:
705 # if neighbor in neg. x-direction is 2 then save the actual voxel and the neighbor!!
706 if lab[positionLabel[0][x] - 1, positionLabel[1][x], positionLabel[2][x]] == 2:
707 vx2 = numpy.array(
708 [
709 positionLabel[0][x] - 1,
710 positionLabel[1][x],
711 positionLabel[2][x],
712 2,
713 ]
714 )
715 contact = numpy.vstack((contact, pix))
716 contact = numpy.vstack((contact, vx2))
717 # y axis neighbor
718 if positionLabel[1][x] < dimY - 1:
719 if lab[positionLabel[0][x], positionLabel[1][x] + 1, positionLabel[2][x]] == 2:
720 vy1 = numpy.array(
721 [
722 positionLabel[0][x],
723 positionLabel[1][x] + 1,
724 positionLabel[2][x],
725 2,
726 ]
727 )
728 contact = numpy.vstack((contact, pix))
729 contact = numpy.vstack((contact, vy1))
730 if positionLabel[1][x] != 0:
731 if lab[positionLabel[0][x], positionLabel[1][x] - 1, positionLabel[2][x]] == 2:
732 vy2 = numpy.array(
733 [
734 positionLabel[0][x],
735 positionLabel[1][x] - 1,
736 positionLabel[2][x],
737 2,
738 ]
739 )
740 contact = numpy.vstack((contact, pix))
741 contact = numpy.vstack((contact, vy2))
742 # z axis neighbor
743 if positionLabel[2][x] < dimX - 1:
744 if lab[positionLabel[0][x], positionLabel[1][x], positionLabel[2][x] + 1] == 2:
745 vz1 = numpy.array(
746 [
747 positionLabel[0][x],
748 positionLabel[1][x],
749 positionLabel[2][x] + 1,
750 2,
751 ]
752 )
753 contact = numpy.vstack((contact, pix))
754 contact = numpy.vstack((contact, vz1))
755 if positionLabel[2][x] != 0:
756 if lab[positionLabel[0][x], positionLabel[1][x], positionLabel[2][x] - 1] == 2:
757 vz2 = numpy.array(
758 [
759 positionLabel[0][x],
760 positionLabel[1][x],
761 positionLabel[2][x] - 1,
762 2,
763 ]
764 )
765 contact = numpy.vstack((contact, pix))
766 contact = numpy.vstack((contact, vz2))
768 return contact
771def _contactPositions(contact, probMap):
772 """
773 determines the position of the points that define the contact area
774 for the random walker probability map
775 the 50 percent probability plane is considered as the area of contact
776 linear interpolation is used to determine the 50 percent probability area
778 Parameters
779 ----------
780 contact : (len(contactVoxels))x4 array
781 z,y,x positions and the label of the voxel
782 as returned by the function contactPairs()
784 probMap : 3D array of floats
785 probability map of one of the labels
786 as determined by the random walker
788 Returns
789 -------
790 coordIntervox : (len(coordIntervox))x3 array
791 z,y,x positions of the 50 percent probability area
792 """
793 coordIntervox = numpy.array([[], [], []]).transpose()
794 for x in range(0, len(contact), 2): # call only contact pairs (increment 2)
795 prob1 = probMap[int(contact[x][0]), int(contact[x][1]), int(contact[x][2])]
796 # pick the probability value of the concerning contact voxel!
797 prob2 = probMap[int(contact[x + 1][0]), int(contact[x + 1][1]), int(contact[x + 1][2])]
798 if prob2 - prob1 != 0:
799 add = float((0.5 - prob1) / (prob2 - prob1))
800 else:
801 add = 0.5 # if both voxels have the same probability (0.5), the center is just in between them!
802 # check whether the next contact voxel is x-axis (0) neighbor or not
803 # x axis neighbor
804 if contact[x][0] != contact[x + 1][0]:
805 # 2 possibilities: a coordinates > or < b coordinates
806 if contact[x][0] < contact[x + 1][0]: # find the contact point in subvoxel resolution!
807 midX = numpy.array([contact[x][0] + add, contact[x][1], contact[x][2]])
808 coordIntervox = numpy.vstack((coordIntervox, midX))
809 else:
810 midX = numpy.array([contact[x][0] - add, contact[x][1], contact[x][2]])
811 coordIntervox = numpy.vstack((coordIntervox, midX))
812 # y axis neighbor
813 elif contact[x][1] != contact[x + 1][1]:
814 if contact[x][1] < contact[x + 1][1]:
815 midY = numpy.array([contact[x][0], contact[x][1] + add, contact[x][2]])
816 coordIntervox = numpy.vstack((coordIntervox, midY))
817 else:
818 midY = numpy.array([contact[x][0], contact[x][1] - add, contact[x][2]])
819 coordIntervox = numpy.vstack((coordIntervox, midY))
820 # z axis neighbor
821 else:
822 if contact[x][2] < contact[x + 1][2]:
823 midZ = numpy.array([contact[x][0], contact[x][1], contact[x][2] + add])
824 coordIntervox = numpy.vstack((coordIntervox, midZ))
825 else:
826 midZ = numpy.array([contact[x][0], contact[x][1], contact[x][2] - add])
827 coordIntervox = numpy.vstack((coordIntervox, midZ))
829 return coordIntervox
832def _contactNormals(dataSet):
833 """
834 determines the contact normal orientation
835 based on the intervoxel positions determined by the function contactPositions()
836 uses a principal component analysis
837 the smallest eigenvector of the covariance matrix is considered as the contact orientation
839 Parameters
840 ----------
841 dataSet : (len(dataSet))x3 array
842 z,y,x positions of the 50 percent probability area of the random walker
844 Returns
845 -------
846 contactNormal : 1x3 array
847 z,y,x directions of the normalised contact normal orientation
848 """
849 covariance = numpy.cov(dataSet, rowvar=0, bias=0)
851 # step 3: calculate the eigenvectors and eigenvalues
852 # eigenvalues are not necessarily ordered ...
853 # colum[:,i] corresponds to eig[i]
854 eigVal, eigVec = numpy.linalg.eig(covariance)
856 # look for the eigenvector corresponding to the minimal eigenvalue!
857 minEV = eigVec[:, numpy.argmin(eigVal)]
859 # vector orientation
860 contactNormal = numpy.zeros((1, 3))
861 if minEV[2] < 0:
862 contactNormal[0, 0], contactNormal[0, 1], contactNormal[0, 2] = (
863 -minEV[0],
864 -minEV[1],
865 -minEV[2],
866 )
867 else:
868 contactNormal[0, 0], contactNormal[0, 1], contactNormal[0, 2] = (
869 minEV[0],
870 minEV[1],
871 minEV[2],
872 )
874 return contactNormal
877def _markerCorrection(
878 markers,
879 numMarkers,
880 distanceMap,
881 volBin,
882 peakDistance=5,
883 struc=numpy.ones((3, 3, 3)),
884):
885 """
886 corrects the number of markers used for the random walker segmentation of two particles
887 if too many markers are found, the markers with the largest distance are considered
888 if too few markers are found, the minimum distance between parameters is reduced
890 Parameters
891 ----------
892 markers : 3D array of integers
893 volume containing the potential markers for the random walker
895 numMarkers : integer
896 number of markers found so far
898 distanceMap : 3D array of floats
899 euclidean distance map
901 volBin : 3D array of integers
902 binary volume
904 peakDistance : integer
905 peak distance as used in skimage.feature.peak_local_max
906 Default value is 15 px
908 struc=numpy.ones((3,3,3)) : 3x3x3 array of integers
909 structuring element for the labelling of the markers
910 Default element is a 3x3x3 array of ones
912 Returns
913 -------
914 markers : 3-D array of integers
915 volume containing the markers
916 """
917 counter = 0
918 # If there are too many markers, slowly increse peakDistance until it's the size of the image, if this overshoots the next bit will bring it back down
919 if numpy.amax(markers) > 2:
920 while numpy.amax(markers) > 2 and peakDistance < min(volBin.shape):
921 peakDistance = max(1, int(peakDistance * 1.3))
922 localMaxiPoints = skimage.feature.peak_local_max(
923 distanceMap,
924 min_distance=peakDistance,
925 exclude_border=False,
926 labels=volBin,
927 num_peaks=2,
928 )
929 localMaxi = numpy.zeros_like(distanceMap)
930 localMaxi[tuple(localMaxiPoints.T)] = True
931 markers, numMarkers = scipy.ndimage.label(localMaxi, structure=struc)
933 # If there are no markers, decrease peakDistance
934 while numpy.any(markers) is False or numpy.amax(markers) < 2:
935 if counter > 10:
936 markers = numpy.zeros(markers.shape)
937 return markers
938 peakDistance = max(1, int(peakDistance / 1.3))
939 localMaxiPoints = skimage.feature.peak_local_max(
940 distanceMap,
941 min_distance=peakDistance,
942 exclude_border=False,
943 labels=volBin,
944 num_peaks=2,
945 )
946 localMaxi = numpy.zeros_like(distanceMap, dtype=bool)
947 localMaxi[tuple(localMaxiPoints.T)] = True
948 markers, numMarkers = scipy.ndimage.label(localMaxi, structure=struc)
949 counter = counter + 1
951 if numpy.amax(markers) > 2:
952 centerOfMass = scipy.ndimage.center_of_mass(markers, labels=markers, index=range(1, numMarkers + 1))
953 distances = numpy.zeros((numMarkers, numMarkers))
954 for i in range(0, numMarkers):
955 for j in range(i, numMarkers):
956 distances[i, j] = numpy.linalg.norm(numpy.array(centerOfMass[i]) - numpy.array(centerOfMass[j]))
958 maxDistance = numpy.amax(distances)
959 posMaxDistance = numpy.where(distances == maxDistance)
961 # let's choose the markers with the maximum distance
962 markers1 = numpy.where(markers == posMaxDistance[0][0] + 1, 1, 0)
963 markers2 = numpy.where(markers == posMaxDistance[1][0] + 1, 2, 0)
965 localMaxi = markers1 + markers2
966 markers, numMarkers = scipy.ndimage.label(localMaxi, structure=struc)
968 return markers
971def localDetectionAssembly(
972 volLab,
973 volGrey,
974 contactList,
975 localThreshold,
976 boundingBoxes=None,
977 nProcesses=nProcessesDefault,
978 radiusThresh=4,
979):
980 """
981 Local contact refinement of a set of contacts
982 checks whether two particles are in contact with a local threshold using ``localDetection()``
984 Parameters
985 ----------
986 volLab : 3D array
987 Labelled volume
989 volGrey : 3D array
990 Grey-scale volume
992 contactList : (ContactNumber)x2 array of integers
993 contact list with grain labels of each contact in a seperate row,
994 as given by ``spam.label.contacts.labelledContacts`` in the list contactingLabels
996 localThreshold : integer or float, same type as the 3D array
997 threshold for binarisation of the subvolume
999 boundingBoxes : lab.max()x6 array of ints, optional
1000 Bounding boxes in format returned by ``boundingBoxes``.
1001 If not defined (Default = None), it is recomputed by running ``boundingBoxes``
1003 nProcesses : integer (optional, default = nProcessesDefault)
1004 Number of processes for multiprocessing
1005 Default = number of CPUs in the system
1007 radiusThresh : integer, optional
1008 radius in px for excluding patches that might not be relevant,
1009 e.g. noise can lead to patches that are not connected with the grains in contact
1010 the patches are calculated with ``equivalentRadii()``
1011 Default is 4
1013 Returns
1014 -------
1015 contactListRefined : (ContactNumber)x2 array of integers
1016 refined contact list, based on the chosen local threshold
1018 Note
1019 ----
1020 see https://doi.org/10.1088/1361-6501/aa8dbf for further information
1021 """
1022 import progressbar
1024 # check if bounding boxes are given
1025 if boundingBoxes is None:
1026 # print "bounding boxes are not given. calculating ...."
1027 boundingBoxes = spam.label.boundingBoxes(volLab)
1029 contactListRefined = []
1030 numberOfJobs = len(contactList)
1031 # print ("\n\tLocal contact refinement")
1032 # print ("\tApplying a local threshold of ", localThreshold, " to each contact.")
1033 # print ("\tNumber of contacts for treatment ", numberOfJobs)
1034 # print ("")
1036 ##########################################
1038 # Function for localDetectionAssembly
1039 global _multiprocessingLocalDetectionAssembly
1041 def _multiprocessingLocalDetectionAssembly(job):
1042 grainA, grainB = contactList[job].astype("int")
1043 labels = [grainA, grainB]
1044 subset = fetchTwoGrains(volLab, labels, volGrey, boundingBoxes)
1045 contact = localDetection(subset["subVolGrey"], localThreshold, radiusThresh)
1046 if contact is True:
1047 # print ("we are in contact!")
1048 return job, grainA, grainB
1050 # Create progressbar
1051 widgets = [
1052 progressbar.FormatLabel(""),
1053 " ",
1054 progressbar.Bar(),
1055 " ",
1056 progressbar.AdaptiveETA(),
1057 ]
1058 pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfJobs)
1059 pbar.start()
1060 finishedNodes = 0
1062 # Run multiprocessing
1063 with multiprocessing.Pool(processes=nProcesses) as pool:
1064 for returns in pool.imap_unordered(_multiprocessingLocalDetectionAssembly, range(0, numberOfJobs)):
1065 finishedNodes += 1
1066 widgets[0] = progressbar.FormatLabel("{}/{} ".format(finishedNodes, numberOfJobs))
1067 pbar.update(finishedNodes)
1068 if returns is not None:
1069 contactListRefined.append([returns[1], returns[2]])
1070 pool.close()
1071 pool.join()
1073 # End progressbar
1074 pbar.finish()
1076 return numpy.asarray(contactListRefined)
1079def contactOrientationsAssembly(
1080 volLab,
1081 volGrey,
1082 contactList,
1083 watershed="ITK",
1084 peakDistance=5,
1085 boundingBoxes=None,
1086 nProcesses=nProcessesDefault,
1087 verbose=False,
1088):
1089 """
1090 Determines contact normal orientation in an assembly of touching particles
1091 uses either directly the labelled image or the random walker implementation from skimage
1093 Parameters
1094 ----------
1095 volLab : 3D array
1096 Labelled volume
1098 volGrey : 3D array
1099 Grey-scale volume
1101 contactList : (ContactNumber)x2 array of integers
1102 contact list with grain labels of each contact in a seperate row,
1103 as given by ``spam.label.contacts.labelledContacts`` in the list contactingLabels
1104 or by ``spam.label.contacts.localDetectionAssembly()``
1106 watershed : string, optional
1107 sets the basis for the determination of the orientation
1108 options are "ITK" for the labelled image from the input,
1109 "RW" for a further segmentation by the random walker
1110 Default = "ITK"
1112 peakDistance : int, optional
1113 Distance in pixels used in skimage.feature.peak_local_max
1114 Default value is 5
1116 boundingBoxes : lab.max()x6 array of ints, optional
1117 Bounding boxes in format returned by ``boundingBoxes``.
1118 If not defined (Default = None), it is recomputed by running ``boundingBoxes``
1120 nProcesses : integer (optional, default = nProcessesDefault)
1121 Number of processes for multiprocessing
1122 Default = number of CPUs in the system
1124 verbose : boolean (optional, default = False)
1125 True for printing the evolution of the process
1126 False for not printing the evolution of process
1128 Returns
1129 -------
1130 contactOrientations : (numContacts)x6 array
1131 contact normal orientation for every contact
1132 [labelGrainA, labelGrainB, orientationZ, orientationY, orientationX, intervoxel positions for PCA]
1134 notTreatedContact : boolean
1135 if False that contact is not determined
1136 if the contact consisted of too few voxels a fit is not possible or feasible
1137 """
1139 # check if bounding boxes are given
1140 if boundingBoxes is None:
1141 # print ("bounding boxes are not given. calculating ....")
1142 boundingBoxes = spam.label.boundingBoxes(volLab)
1144 contactOrientations = []
1145 numberOfJobs = len(contactList)
1146 # print ("\n\tDetermining the contact orientations of an assembly of particles")
1147 # print ("\tNumber of contacts ", numberOfJobs)
1148 # print ("")
1150 ##########################################
1152 # Function for ContactOrientationsAssembly
1153 global _multiprocessingContactOrientationsAssembly
1155 def _multiprocessingContactOrientationsAssembly(job):
1156 grainA, grainB = contactList[job, 0:2].astype("int")
1157 labels = [grainA, grainB]
1158 subset = fetchTwoGrains(volLab, labels, volGrey, boundingBoxes)
1159 contactNormal, intervox, NotTreatedContact = spam.label.contactOrientations(
1160 subset["subVolBin"],
1161 subset["subVolLab"],
1162 watershed,
1163 peakDistance=peakDistance,
1164 verbose=verbose,
1165 )
1166 # TODO work on not treated contacts -- output them!
1167 return (
1168 grainA,
1169 grainB,
1170 contactNormal[0],
1171 contactNormal[1],
1172 contactNormal[2],
1173 intervox,
1174 )
1176 # Create progressbar
1177 widgets = [
1178 progressbar.FormatLabel(""),
1179 " ",
1180 progressbar.Bar(),
1181 " ",
1182 progressbar.AdaptiveETA(),
1183 ]
1184 pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfJobs)
1185 pbar.start()
1186 finishedNodes = 0
1188 # Run multiprocessing
1189 with multiprocessing.Pool(processes=nProcesses) as pool:
1190 for returns in pool.imap_unordered(_multiprocessingContactOrientationsAssembly, range(0, numberOfJobs)):
1191 finishedNodes += 1
1192 widgets[0] = progressbar.FormatLabel("{}/{} ".format(finishedNodes, numberOfJobs))
1193 pbar.update(finishedNodes)
1194 if returns is not None:
1195 contactOrientations.append(
1196 [
1197 returns[0],
1198 returns[1],
1199 returns[2],
1200 returns[3],
1201 returns[4],
1202 returns[5],
1203 ]
1204 )
1205 # result[2], result[3], result[4], result[5], result[6], result[7]
1206 pool.close()
1207 pool.join()
1209 # End progressbar
1210 pbar.finish()
1212 return numpy.asarray(contactOrientations)