Coverage for local_installation_linux/mumott/data_handling/utilities/data_processing.py: 78%
30 statements
« prev ^ index » next coverage.py v7.3.2, created at 2024-08-11 23:08 +0000
1import logging
3import numpy as np
4from numpy.typing import NDArray
5from typing import Any, Dict, Tuple
8logger = logging.getLogger(__name__)
11def get_transmittivities(diode: NDArray,
12 normalize_per_projection: bool = False,
13 normalization_percentile: float = 99.9,
14 cutoff_values: Tuple[float, float] = (1e-4, 1.0)) -> Dict[str, Any]:
15 r""" Calculates the transmittivity from the diode, i.e., the fraction of transmitted
16 intensity relative to a high percentile.
18 Notes
19 -----
21 Diode readouts may be given in various formats such as a count or a current.
22 When doing absorption tomography, one is generally interested in the fraction of
23 transmitted intensity. Since we do not generally have access to the incoming flux,
24 or its theoretical readout at complete transmission, we can instead normalize
25 the diode readout based on the largest values, where the beam has only passed through
26 some air. We thus want to compute
28 .. math::
29 T(i, j, k) = \frac{I_T(i, j, k)}{I_0}
31 where :math:`I_T(i, j, k)` is the diode readout value at projection :math:`i`, and pixel :math:`(j, k)`
32 with the approximation :math:`I_0 \approx \text{max}(I_T(i, j, k))`.
34 To avoid routine normalization based on individual spurious readouts
35 (from, e.g., hot pixels), by default the normalization is done based on the
36 99.9th percentile rather than the strict maximum.
37 The normalized values are then clipped to the interval specified
38 by :attr:`cutoff_values`, by default ``(1e-4, 1.0)``. A mask is returned which masks out
39 any values outside this range, which can be useful to mask out spurious readouts.
41 If the transmittivities are to be used to normalize :term:`SAXS` data, one should
42 leave the :attr:`normalize_per_projection` option at ``False``, because the :term:`SAXS` data
43 also scales with the incoming flux, and thus we want any variations in flux between
44 projections to be preserved in the transmittivities.
46 However, if the transmittivities are to be used for transmission (or absorption) tomography,
47 then the :attr:`normalize_per_projection`
48 option should be set to ``True``. Since we are interested in the transmittivity
49 of the sample irrespective of the incoming flux, we are therefore better off assuming
50 the flux is constant over each projection. This corresponds to the slightly modified computation
52 .. math::
53 T(i, j, k) = \frac{I_T(i, j, k)}{I_0(i)}
55 with the approximation :math:`I_0(i) \approx \text{max}_{j, k}(I_T(i, j, k))`,
56 with the understanding that we take the maximum value for each projection :math:`i`.
58 Parameters
59 ----------
60 diode
61 An array of diode readouts.
62 normalize_per_projection
63 If ``True``, the diode will be normalized projection-wise.
64 This is the appropriate choice for absorption tomography.
65 For :term:`SAXS`, it is preferable to normalize the diode across the entire set
66 of diode measurements in order to account for possible variations in flux.
67 Default value is ``False``.
68 normalization_percentile
69 The percentile of values in either the entire set of diode measurements or each projection
70 (depending on :attr:`normalize_per_projection`) to use for normalization.
71 The default value is ``99.9``. Values above this range will be clipped. If you are
72 certain that you do not have spuriously large diode readout values, you can specify
73 ``100.`` as the percentile instead.
74 cutoff_values
75 The cutoffs to use for the transmittivity calculation. Default value is
76 ``(1e-4, 1.0)``, i.e., one part in ten thousand for the lower bound, and
77 a hundred percent for the upper bound. For values outside of this range, it may
78 be desirable to mask them out during *any* calculation. For this purpose,
79 a ``cutoff_mask`` is included in the return dictionary with the same shape as
80 the weights in :attr:`projections`.
81 In some cases, you may wish to specify other bounds. For example, if you know that
82 your sample is embedded in a substrate which reduces the maximum possible transmittivity,
83 you may wish to lower the upper bound. If you know that your sample has extremely low
84 transmittivity (perhaps compensated with a very long exposure time), then you can set
85 the lower cutoff even lower.
86 The cutoffs must lie within the open interval ``(0, 1]``. A lower bound of ``0``
87 is not permitted since this would lead to an invalid absorbance.
89 Returns
90 -------
91 A dictionary with three entries, ``transmittivity``, ``cutoff_mask_lower``,
92 and ``cutoff_mask_upper``.
93 """
94 if cutoff_values[0] <= 0 or cutoff_values[1] > 1: 94 ↛ 95line 94 didn't jump to line 95, because the condition on line 94 was never true
95 raise ValueError('cutoff_values must lie in the open interval (0, 1], but'
96 f' cutoff values of {cutoff_values} were specified!')
97 # These should already be copies, but for future-proofing.
98 transmittivity = diode.copy()
99 cutoff_mask_lower = np.ones_like(diode)[..., None]
100 cutoff_mask_upper = np.ones_like(diode)[..., None]
101 if normalize_per_projection: 101 ↛ 111line 101 didn't jump to line 111, because the condition on line 101 was never false
102 for i, t in enumerate(transmittivity):
103 # Normalize each projection's maximum value to 1, cutting off values outside the range.
104 normalization_value = np.percentile(t, normalization_percentile)
105 transmittivity[i] = t * np.reciprocal(normalization_value)
106 cutoff_mask_lower[i] = (transmittivity[i, ..., None] >= cutoff_values[0])
107 cutoff_mask_upper[i] = (transmittivity[i, ..., None] <= cutoff_values[1])
108 np.clip(transmittivity[i], *cutoff_values, out=transmittivity[i])
109 else:
110 # Normalize the maximum value of all projections to 1, cutting off values outside the range.
111 normalization_value = np.percentile(transmittivity, normalization_percentile)
112 transmittivity *= np.reciprocal(normalization_value)
113 cutoff_mask_lower = (transmittivity[..., None] >= cutoff_values[0])
114 cutoff_mask_upper = (transmittivity[..., None] <= cutoff_values[1])
115 np.clip(transmittivity, *cutoff_values, out=transmittivity)
117 return dict(transmittivity=transmittivity, cutoff_mask_lower=cutoff_mask_lower,
118 cutoff_mask_upper=cutoff_mask_upper)
121def get_absorbances(diode: NDArray, **kwargs) -> Dict[str, Any]:
122 r""" Calculates the absorbance based on the transmittivity of the diode data.
124 Notes
125 -----
126 The absorbance is defined as the negative base-10 logarithm of the transmittivity.
127 Specifically,
129 .. math::
130 A(i, j, k) = -\log_{10}(T(i, j, k))
132 where :math:`T` is the transmittivity, normalized to the open interval :math:`(0, 1]`.
133 It can be inferred from this formula why :math:`T(i, j, k)` must not have values which
134 are equal to or smaller than :math:`0`, as that would give a non-finite absorbance.
135 Similarly, values greater than :math:`1` would result in physically impossible negative
136 absorbances.
138 The transmittivity is calculated directly from diode readouts, which may or may not already
139 be normalized and clipped. When using already normalized diode values, it is best
140 to set the keyword argument :attr:`normalization_percentile` to ``100``, and the
141 argument :attr:`cutoff_values` to whatever cutoff your normalization used.
143 See :func:`get_transmittivities` for more details on the transmittivity calculation.
145 Parameters
146 ----------
147 diode
148 An array of diode readouts.
149 kwargs
150 Keyword arguments which are passed on to :func:`get_transmittivities`.
152 Returns
153 -------
154 A dictionary with the absorbance, and the union of the cutoff masks from
155 :func:`get_transmittivities`.
156 """
158 transmittivity_dictionary = get_transmittivities(diode, **kwargs)
159 transmittivity = transmittivity_dictionary['transmittivity']
160 absorbances = -np.log10(transmittivity)
161 assert np.all(np.isfinite(absorbances))
162 return dict(absorbances=absorbances[..., None],
163 cutoff_mask=(transmittivity_dictionary['cutoff_mask_lower'] *
164 transmittivity_dictionary['cutoff_mask_upper']))