Coverage for local_installation_linux / mumott / output_handling / reconstruction_derived_quantities.py: 99%
82 statements
« prev ^ index » next coverage.py v7.13.0, created at 2025-12-26 11:25 +0000
1from dataclasses import dataclass
2import h5py
3import importlib.resources
4import os
5import shutil
6import numpy as np
7from matplotlib.colors import hsv_to_rgb
10def get_sorted_eigenvectors(tensors_array: np.array):
11 """ Calculate eigenvectors and eigenvalues of an array of
12 3 by 3 matrices and sort according to increasing eigenvalues.
14 Parameters
15 ----------
16 tensors_array
17 numpy array containing the 3 by 3 tensors with the tensor
18 indicies as the two last.
20 Returns
21 -------
22 eigenvalues
23 numpy array containing the eigenvalues of the tensors where
24 the last dimension indexes the eigenvalues. Smallest eigenvalues
25 come first.
26 eigenvectors
27 numpy array containing the eigenvalues of the tensors where
28 the last dimension indexes the eigenvalues and the second
29 to last dimension is the vector index.
30 """
32 volume_shape = tensors_array.shape[:-2]
34 # Compute and sort eigenvectors
35 w, v = np.linalg.eigh(tensors_array.reshape(-1, 3, 3))
36 sorting = np.argsort(w, axis=1).reshape(-1, 3, 1)
37 v = v.transpose(0, 2, 1)
38 v = np.take_along_axis(v, sorting, axis=1)
39 v = v.transpose(0, 2, 1)
40 v = v / np.sqrt(np.sum(v ** 2, axis=1).reshape(-1, 1, 3))
41 eigenvectors = v.reshape(volume_shape + (3, 3,))
42 eigenvalues = np.sort(w, axis=-1).reshape(volume_shape + (3,))
44 # Flip eigenvectors to have a positive z-component
45 for ii in range(3):
46 whereflip = eigenvectors[..., 2, ii] < 0.0
47 eigenvectors[whereflip, :, ii] = -eigenvectors[whereflip, :, ii]
49 return eigenvalues, eigenvectors
52@dataclass
53class ReconstructionDerivedQuantities:
54 """ A number of useful quantities that have been computed from the coefficients of a
55 reconstruction.
57 Attributes
58 ----------
59 volume_shape : tuple
60 The shape if the reconstructed volume.
61 mean_intensity : np.array
62 The mean intensity of the reconstructed functions over the unit sphere for each voxel.
63 fractional_anisostropy : np.array
64 A measure of the relative amount of anisotropic scattering in each voxel.
65 eigenvector_1
66 eigenvector_2
67 eigenvector_3 : np.array
68 The eingenvectors of the second order tensor component of the function.
69 Sorted by ascending eigenvalue.
70 eigenvalue_1
71 eigenvalue_2
72 eigenvalue_3 : np.array
73 The eingenvalues of the second order tensor component of the function.
74 Sorted by ascending eigenvalue.
75 second_moment_tensor : np.array
76 The second-moment tensor, which includes both the zero-th and second-order
77 parts of the reconstructed functions.
78 """
80 volume_shape: tuple
81 mean_intensity: np.array
82 fractional_anisotropy: np.array
83 eigenvector_1: np.array
84 eigenvector_2: np.array
85 eigenvector_3: np.array
86 eigenvalue_1: np.array
87 eigenvalue_2: np.array
88 eigenvalue_3: np.array
89 second_moment_tensor: np.array
91 def write(self, filename: str) -> None:
92 """ Save the derived reconstruction quantities in both an HDF5 file and a Paraview readable
93 XDMF file.
95 Parameters
96 ----------
97 filename : str
98 The filename to save the data in. The extension `.h5` and `.xdmf` will be added
99 to the filename.
100 """
102 if filename.endswith('.h5') or filename.lower().endswith('.xdmf'): 102 ↛ 108line 102 didn't jump to line 108 because the condition on line 102 was always true
104 filename, _ = os.path.splitext(filename)
105 basename = os.path.basename(filename)
107 # Make .h5 data file
108 with h5py.File(filename + '.h5', 'w') as file:
110 file.create_dataset('mean_intensity', data=self.mean_intensity.transpose((2, 1, 0)))
111 file.create_dataset('fractional_anisotropy', data=self.fractional_anisotropy.transpose((2, 1, 0)))
113 file.create_dataset('eigenvector_1', data=self.eigenvector_1.transpose((2, 1, 0, 3)))
114 file.create_dataset('eigenvector_2', data=self.eigenvector_2.transpose((2, 1, 0, 3)))
115 file.create_dataset('eigenvector_3', data=self.eigenvector_3.transpose((2, 1, 0, 3)))
117 file.create_dataset('eigenvector_1_rgb',
118 data=get_colors_on_halfsphere(self.eigenvector_1).transpose((2, 1, 0, 3)))
119 file.create_dataset('eigenvector_2_rgb',
120 data=get_colors_on_halfsphere(self.eigenvector_2).transpose((2, 1, 0, 3)))
121 file.create_dataset('eigenvector_3_rgb',
122 data=get_colors_on_halfsphere(self.eigenvector_3).transpose((2, 1, 0, 3)))
124 file.create_dataset('eigenvalue_1', data=self.eigenvalue_1.transpose((2, 1, 0)))
125 file.create_dataset('eigenvalue_2', data=self.eigenvalue_2.transpose((2, 1, 0)))
126 file.create_dataset('eigenvalue_3', data=self.eigenvalue_3.transpose((2, 1, 0)))
128 file.create_dataset('second_moment_tensor',
129 data=self.second_moment_tensor.transpose((2, 1, 0, 3, 4)))
131 # Make XDMF file for paraview
132 with open(filename + '.xdmf', 'w') as file:
134 # Header
135 file.write("""<?xml version="1.0" ?>
136<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
137<Xdmf Version="2.0">
138<Domain>
139<Grid Name="Structured Grid" GridType="Uniform">
140""")
141 file.write('<Topology TopologyType="3DCoRectMesh" NumberOfElements="' +
142 f'{self.volume_shape[2]} {self.volume_shape[1]} {self.volume_shape[0]}"/>\n')
143 file.write(""" <Geometry GeometryType="Origin_DxDyDz">
144 <DataItem Name="Origin" Dimensions="3" NumberType="Float" Precision="8" Format="XML">
145 0.5 0.5 0.5
146 </DataItem>
147 <DataItem Name="Spacing" Dimensions="3" NumberType="Float" Precision="8" Format="XML">
148 1.0 1.0 1.0
149 </DataItem>
150</Geometry>
152""")
154 # Scalars
155 for data_string in [
156 'mean_intensity',
157 'fractional_anisotropy',
158 'eigenvalue_1',
159 'eigenvalue_2',
160 'eigenvalue_3',]:
161 file.write(f'<Attribute Name="{data_string}" AttributeType="Scalar" Center="Node">\n')
162 file.write(' <DataItem Dimensions="' +
163 f'{self.volume_shape[2]} {self.volume_shape[1]} {self.volume_shape[0]}"' +
164 ' NumberType="Float" Precision="8" Format="HDF">\n')
165 file.write(f' {basename}.h5:/{data_string}')
166 file.write("""
167 </DataItem>
168</Attribute>
170""")
172 # Vectors
173 for data_string in [
174 'eigenvector_1',
175 'eigenvector_2',
176 'eigenvector_3',
177 'eigenvector_1_rgb',
178 'eigenvector_2_rgb',
179 'eigenvector_3_rgb',]:
180 file.write(f'<Attribute Name="{data_string}" AttributeType="Vector" Center="Node">\n')
181 file.write(' <DataItem Dimensions="' +
182 f'{self.volume_shape[2]} {self.volume_shape[1]} {self.volume_shape[0]} 3"' +
183 ' NumberType="Float" Precision="8" Format="HDF">\n')
184 file.write(f' {basename}.h5:/{data_string}')
185 file.write("""
186 </DataItem>
187</Attribute>
189""")
191 file.write("""</Grid>
192</Domain>
193</Xdmf>
194""")
196 # Color legend
197 with importlib.resources.path(__package__, 'color_map.png') as p:
198 path = p
199 shutil.copyfile(path, 'direction_colormap.png')
202def get_colors_on_halfsphere(vectors: np.array) -> np.array:
203 """ A colourmap on the unit halfsphere with inversion symmetry.
204 Vectors along the `y` direction are grey and the 'x-z' equator goes through the
205 hsv hue cycle. The +y+z and -y-z sections are made brighter than the
206 +y-z and -y+z sections such that the pairs (x, y, 0) and (x, -y, 0) are the
207 only non-equivalent points given the same colour.
209 Parameters
210 ----------
211 unit-vectors
212 Array where the last dimension is of length 3, corresponding to the vector index.
215 Returns
216 -------
217 Array of the same shape as the input where the third dimension is an RGB value in
218 floating point format.
220 """
222 vectors = np.copy(vectors)
223 whereflip = vectors[..., 1] < 0
224 vectors[whereflip, :] = -vectors[whereflip, :]
226 theta = np.arccos(vectors[..., 1])
227 phi = np.arctan2(vectors[..., 2], vectors[..., 0])
228 hue = ((phi) % (np.pi))/np.pi
229 saturation = (np.arctan(theta/2) / np.arctan(np.pi/4))**2
230 modifier = -np.sin(phi)*np.sin(2*theta)**2
231 modifier = np.sign(modifier) * np.abs(modifier)
232 value = np.ones(theta.shape)*0.7 + 0.2*modifier
234 hsv = np.stack([hue, saturation, value], axis=-1)
235 rgb_colours = hsv_to_rgb(hsv)
236 return rgb_colours