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

49 statements  

« prev     ^ index     » next       coverage.py v7.2.3, created at 2023-11-22 13:26 +0000

1# Library of SPAM functions for measuring a porosity field 

2# Copyright (C) 2020 SPAM Contributors 

3# 

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

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

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

7# any later version. 

8# 

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

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

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

12# more details. 

13# 

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

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

16 

17 

18import matplotlib.pyplot as plt 

19import numpy 

20import spam.DIC.grid as grid 

21from spam.measurements.measurementsToolkit import ( 

22 porosityFieldBinary as porosityFieldBinaryCPP, 

23) 

24 

25 

26def porosityField( 

27 image, 

28 poreValue=None, 

29 solidValue=None, 

30 maskValue=None, 

31 nodeSpacing="auto", 

32 hws=10, 

33 showGraph=False, 

34): 

35 """ 

36 This function calculates a porosity field on a greyscale image, outputting a field in % porosity. 

37 A list of half-window sizes can be given, allowing for porosity measurement stability to to evaluate easily. 

38 

39 Parameters 

40 ---------- 

41 image : 3D numpy array 

42 

43 poreValue : float 

44 The grey value corresponding to the peak of the pores. 

45 Grey values of this value as considered 100% pore. 

46 In the case of 2 peaks in the pores (having water and air), should put the principal one. 

47 

48 solidValue : float 

49 The grey value corresponding to the peak of the solids. 

50 Grey values of this value as considered 100% solid. 

51 

52 maskValue : float, optional 

53 Indicate whether there is a special value in the input image 

54 that corrseponds to a mask which indicates voxels that are not taken 

55 into account in the porosity calculation. 

56 Default = None 

57 

58 hws : int or list, optional 

59 half-window size for cubic measurement window 

60 Default = 10 

61 

62 nodeSpacing : int, optional 

63 spacing of nodes in pixels 

64 Default = 10 pixels 

65 

66 showGraph : bool, optional 

67 If a list of hws is input, show a graph of the porosity evolution 

68 Default = False 

69 

70 Returns 

71 ------- 

72 3D (or 4D) numpy array for constant half-window size (list of half-window sizes) 

73 

74 """ 

75 

76 # The C function for porosity needs an 8-bit image with greyvalues 

77 if poreValue is None or solidValue is None: 

78 print( 

79 "spam.measurements.porosityField(): Need to give me the greyvalue of the pores and solids" 

80 ) 

81 return 

82 

83 if poreValue == solidValue: 

84 print( 

85 "spam.measurements.porosityField(): Need to give me the correct greyvalue of the pores and solids, where greyvalue of pores always smaller than that of solids" 

86 ) 

87 return 

88 

89 imPorosity = numpy.zeros_like(image, dtype="<u1") 

90 # Rescale image to imPorosity between poreValue and solidValue, ! 

91 if poreValue < solidValue: 

92 imPorosity[image > poreValue] = ( 

93 numpy.true_divide(100, solidValue - poreValue) 

94 ) * (solidValue - image[image > poreValue]) 

95 imPorosity[image <= poreValue] = 100 

96 imPorosity[image >= solidValue] = 0 

97 

98 if poreValue > solidValue: 

99 imPorosity[image > solidValue] = ( 

100 numpy.true_divide(100, poreValue - solidValue) 

101 ) * (poreValue - image[image > solidValue]) 

102 imPorosity[image <= solidValue] = 100 

103 imPorosity[image >= poreValue] = 0 

104 

105 # If a mask value is given, apply it 

106 if maskValue is not None: 

107 imPorosity[image == maskValue] = 255 

108 

109 if nodeSpacing == "auto": 

110 nodeSpacing = min(image.shape) / 2 

111 

112 # if hws == "auto": 10 

113 if type(hws) == int: 

114 hws = [hws] 

115 

116 points, layout = grid.makeGrid(imPorosity.shape, nodeSpacing) 

117 

118 pointsArray = [] 

119 hwsArray = [] 

120 

121 for p in points: 

122 for i in hws: 

123 pointsArray.append(p) 

124 hwsArray.append(i) 

125 

126 pointsArray = numpy.array(pointsArray).astype("<i4") 

127 hwsArray = numpy.array(hwsArray).astype("<i4") 

128 porosities = numpy.zeros_like(hwsArray, dtype="<f4") 

129 

130 # call C function on scaled 8 bit image 

131 porosityFieldBinaryCPP(imPorosity, pointsArray, hwsArray, porosities) 

132 

133 if len(hws) == 1: 

134 porosities = porosities.reshape((layout[0], layout[1], layout[2])) 

135 

136 else: 

137 porosities = porosities.reshape((layout[0], layout[1], layout[2], len(hws))) 

138 

139 if showGraph is True: 

140 for z in range(porosities.shape[0]): 

141 for y in range(porosities.shape[1]): 

142 for x in range(porosities.shape[2]): 

143 plt.plot(hws, porosities[z, y, x, :]) 

144 plt.xlabel("Half-window size") 

145 plt.ylabel("Porosity (%)") 

146 plt.show() 

147 

148 return porosities