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
« 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/>.
18import matplotlib.pyplot as plt
19import numpy
20import spam.DIC.grid as grid
21from spam.measurements.measurementsToolkit import (
22 porosityFieldBinary as porosityFieldBinaryCPP,
23)
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.
39 Parameters
40 ----------
41 image : 3D numpy array
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.
48 solidValue : float
49 The grey value corresponding to the peak of the solids.
50 Grey values of this value as considered 100% solid.
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
58 hws : int or list, optional
59 half-window size for cubic measurement window
60 Default = 10
62 nodeSpacing : int, optional
63 spacing of nodes in pixels
64 Default = 10 pixels
66 showGraph : bool, optional
67 If a list of hws is input, show a graph of the porosity evolution
68 Default = False
70 Returns
71 -------
72 3D (or 4D) numpy array for constant half-window size (list of half-window sizes)
74 """
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
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
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
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
105 # If a mask value is given, apply it
106 if maskValue is not None:
107 imPorosity[image == maskValue] = 255
109 if nodeSpacing == "auto":
110 nodeSpacing = min(image.shape) / 2
112 # if hws == "auto": 10
113 if type(hws) == int:
114 hws = [hws]
116 points, layout = grid.makeGrid(imPorosity.shape, nodeSpacing)
118 pointsArray = []
119 hwsArray = []
121 for p in points:
122 for i in hws:
123 pointsArray.append(p)
124 hwsArray.append(i)
126 pointsArray = numpy.array(pointsArray).astype("<i4")
127 hwsArray = numpy.array(hwsArray).astype("<i4")
128 porosities = numpy.zeros_like(hwsArray, dtype="<f4")
130 # call C function on scaled 8 bit image
131 porosityFieldBinaryCPP(imPorosity, pointsArray, hwsArray, porosities)
133 if len(hws) == 1:
134 porosities = porosities.reshape((layout[0], layout[1], layout[2]))
136 else:
137 porosities = porosities.reshape((layout[0], layout[1], layout[2], len(hws)))
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()
148 return porosities