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