Coverage for local_installation_linux/mumott/pipelines/fbp_utilities.py: 88%

1import logging 

2 

3import numpy as np 

4 

5from mumott import Geometry 

6from mumott.data_handling.utilities import get_absorbances 

7from mumott.core.projection_stack import ProjectionStack 

8 

9logger = logging.getLogger(__name__) 

10 

11 

12def _get_unique_indices(angles: np.ndarray[float], 

13 indices: np.ndarray[int], 

14 tolerance: float = 1e-4) -> np.ndarray[float]: 

15 """Filters indices based on angles.""" 

16 angles = np.mod(angles, np.pi) 

17 out_indices = [indices[0]] 

18 out_angles = [angles[0]] 

19 

20 for angle, index in zip(angles, indices): 

21 angle_delta = np.abs(angle - np.array(out_angles)) 

22 if np.all((angle_delta > tolerance) & (angle_delta < np.pi - tolerance)): 

23 out_indices.append(index) 

24 out_angles.append(angle) 

25 return np.array(out_indices) 

26 

27 

28def _get_orthogonal_axis(geometry: Geometry, 

29 projection_index: int = 0, 

30 axis_string='inner'): 

31 """Retrieves index of axis orthogonal to inner or outer axis in geometry.""" 

32 if axis_string == 'inner': 

33 axis = geometry.inner_axes[projection_index] 

34 elif axis_string == 'outer': 34 ↛ 35line 34 didn't jump to line 35, because the condition on line 34 was never true

35 axis = geometry.outer_axes[projection_index] 

36 else: 

37 raise ValueError('axis_string must be either "inner" or "outer", ' 

38 f'but a value of "{axis_string}" was specified.') 

39 j_projection = np.dot(axis, geometry.j_direction_0) 

40 k_projection = np.dot(axis, geometry.k_direction_0) 

41 if not np.isclose(abs(j_projection + k_projection), 1): 

42 raise ValueError('Rotation axis must be orthogonal to the j or k-axis.') 

43 

44 if np.isclose(j_projection, 0): 

45 return 1 

46 else: 

47 return 2 

48 

49 

50def _get_filter(length: int, 

51 filter_type: str) -> np.ndarray[float]: 

52 """Retrieves a high-pass filter for FBP based on the string.""" 

53 u = np.fft.fftfreq(length) / (4 * np.pi) 

54 filter = abs(u) 

55 filter[abs(u) > 1 / (16 * np.pi)] = 0.0 

56 if filter_type.lower() == 'ram-lak': 

57 return filter 

58 elif filter_type.lower() == 'hann': 

59 filter *= 0.5 + 0.5 * np.cos(2 * np.pi * u) 

60 elif filter_type.lower() == 'hamming': 

61 filter *= 0.54 + 0.46 * np.cos(2 * np.pi * u) 

62 elif filter_type.lower() == 'shepp-logan': 

63 filter *= np.sinc(u) 

64 elif filter_type.lower() == 'cosine': 

65 filter *= np.cos(u * np.pi) 

66 else: 

67 raise ValueError(f'Unknown filter type: "{filter_type}"!' 

68 ' Permitted values are: "Ram-Lak", "Hamming", "Hann",' 

69 ' "cosine", and "Shepp-Logan".') 

70 return filter 

71 

72 

73def get_filtered_projections(projections: ProjectionStack, 

74 axis_string: str = 'inner', 

75 filter_type: str = 'Ram-Lak', 

76 normalization_percentile: float = 99.9, 

77 transmittivity_cutoff: tuple[float, float] = (None, None)) \ 

78 -> tuple[np.ndarray, np.ndarray]: 

79 """ 

80 Applies a high-pass filter to a selected subset of the 

81 absorbances for filtered back projection. 

82 

83 Parameters 

84 ---------- 

85 projections 

86 The :class:`ProjectionStack <mumott.data_handling.ProjectionStack>` to calculate the 

87 filtered projections from. 

88 axis_string 

89 Default is ``'inner'``, the value depends on how the sample is mounted to the holder. Typically, 

90 the inner axis is the rotation axis while the ``'outer'`` axis refers to the tilt axis. 

91 filter_string 

92 Default is ``'ram-lak'``, a high-pass filter. Other options are ``'Hamming'`` and ``'Hanning'``. 

93 normalization_percentile 

94 The normalization percentile to use for the transmittivity calculation. See 

95 :func:`get_transmittivities <mumott.data_handling.utilities.get_transmittivities>` for details. 

96 transmittivity_cutoff 

97 The cutoffs to use for the transmittivity calculation. See 

98 :func:`get_transmittivities <mumott.data_handling.utilities.get_transmittivities>` for details. 

99 

100 Returns 

101 ------- 

102 A tuple containing the filtered subset of the absorbances 

103 and the index of the axis orthogonal to the inner or outer axis. 

104 """ 

105 geometry = projections.geometry 

106 if axis_string == 'inner': 106 ↛ 109line 106 didn't jump to line 109, because the condition on line 106 was never false

107 tilt_angles = geometry.outer_angles_as_array 

108 rotation_angles = geometry.inner_angles_as_array 

109 elif axis_string == 'outer': 

110 tilt_angles = geometry.inner_angles_as_array 

111 rotation_angles = geometry.outer_angles_as_array 

112 else: 

113 raise ValueError(f'Unknown axis: {axis_string}, ' 

114 'please specify "inner" or "outer".') 

115 

116 # Check if we have any transmittivity cutoff values to consider 

117 cutoff_values = (1e-4, 1) # default 

118 for k in range(2): 

119 if transmittivity_cutoff[k] is not None: 119 ↛ 120line 119 didn't jump to line 120, because the condition on line 119 was never true

120 cutoff_values[k] = transmittivity_cutoff[k] 

121 

122 abs_dict = get_absorbances( 

123 projections.diode, 

124 normalize_per_projection=True, 

125 normalization_percentile=normalization_percentile, 

126 cutoff_values=cutoff_values) 

127 absorbances = abs_dict['absorbances'] 

128 # Find projections without any tilt and the rotation of each such projection 

129 no_tilt_indices = np.where(np.isclose(tilt_angles, 0))[0] 

130 no_tilt_angles = rotation_angles[no_tilt_indices] 

131 

132 # Remove 'duplicate' projections with equivalent or very similar rotation angles 

133 no_tilt_indices = _get_unique_indices(no_tilt_angles, no_tilt_indices) 

134 absorbances = absorbances[no_tilt_indices] 

135 filter_axis = _get_orthogonal_axis(geometry, 0, axis_string) 

136 filter = _get_filter(length=absorbances.shape[filter_axis], 

137 filter_type=filter_type) 

138 if filter_axis == 1: 138 ↛ 141line 138 didn't jump to line 141, because the condition on line 138 was never false

139 filter = filter.reshape(1, -1, 1, 1) 

140 else: 

141 filter = filter.reshape(1, 1, -1, 1) 

142 

143 # reduce weights over last index 

144 abs_mask = np.all(abs_dict['cutoff_mask'][no_tilt_indices] * 

145 projections.weights[no_tilt_indices] > 0., axis=-1).astype(float) 

146 absorbances *= abs_mask[..., None] # redundant axis needed for consistency reasons 

147 abs_fft = np.fft.fft(absorbances, axis=filter_axis) 

148 abs_fft *= filter 

149 filtered_absorbances = np.fft.ifft(abs_fft, axis=filter_axis).real 

150 return np.ascontiguousarray(filtered_absorbances), no_tilt_indices