Coverage for local_installation_linux/mumott/core/geometry.py: 95%
635 statements
« prev ^ index » next coverage.py v7.3.2, created at 2024-08-11 23:08 +0000
1from __future__ import annotations
2import logging
3import tarfile
4import tempfile
5import json
6import os
7import codecs
8import numpy as np
9from typing import NamedTuple, Union
10from scipy.spatial.transform import Rotation
11from mumott.core.hashing import list_to_hash
12from mumott.core.probed_coordinates import ProbedCoordinates
13from numpy.typing import NDArray
15logger = logging.getLogger(__name__)
18class GeometryTuple(NamedTuple):
19 """Tuple for passing and returning projection-wise geometry information.
20 This is a helper class used by :class:`Geometry`.
22 Attributes
23 ----------
24 rotation
25 Rotation matrix. If the :attr:`angle` and :attr:`axis` arguments
26 are given, this matrix should correspond to the
27 matrix given by R_outer @ R_inner, where R_inner
28 is defined by a rotation by :attr:`inner_angle` about
29 :attr:`inner_axis`, and similarly for R_outer.
30 j_offset
31 Offset to align projection in the direction j.
32 k_offset
33 Offset to align projection in the direction k.
34 inner_angle
35 Angle of rotation about :attr:`inner_axis` in radians.
36 outer_angle
37 Angle of rotation about :attr:`outer_axis` in radians.
38 inner_axis
39 Inner rotation axis.
40 outer_axis
41 Outer rotation axis.
42 """
43 rotation: np.ndarray[float] = np.eye(3, dtype=float)
44 j_offset: float = float(0)
45 k_offset: float = float(0)
46 inner_angle: float = None
47 outer_angle: float = None
48 inner_axis: np.ndarray[float] = None
49 outer_axis: np.ndarray[float] = None
51 def __hash__(self) -> int:
52 to_hash = [self.rotation.ravel(), self.j_offset, self.k_offset,
53 self.inner_angle, self.outer_angle, self.inner_axis,
54 self.outer_axis]
55 return int(list_to_hash(to_hash), 16)
57 def __str__(self) -> str:
58 wdt = 74
59 s = []
60 s += ['-' * wdt]
61 s += ['GeometryTuple'.center(wdt)]
62 s += ['-' * wdt]
63 with np.printoptions(threshold=4, precision=5, linewidth=60, edgeitems=1):
64 ss = ', '.join([f'{r}' for r in self.rotation])
65 s += ['{:18} : {}'.format('rotation', ss)]
66 s += ['{:18} : {}'.format('j_offset', self.j_offset)]
67 s += ['{:18} : {}'.format('k_offset', self.k_offset)]
68 s += ['{:18} : {}'.format('inner_angle', self.inner_angle)]
69 s += ['{:18} : {}'.format('outer_angle', self.outer_angle)]
70 ss = ', '.join([f'{r}' for r in self.inner_axis])
71 s += ['{:18} : {}'.format('inner_axis', ss)]
72 ss = ', '.join([f'{r}' for r in self.outer_axis])
73 s += ['{:18} : {}'.format('outer_axis', ss)]
74 s += ['{:18} : {}'.format('Hash', hex(hash(self))[2:8])]
75 s += ['-' * wdt]
76 return '\n'.join(s)
78 def _repr_html_(self) -> str:
79 s = []
80 s += ['<h3>GeometryTuple</h3>']
81 s += ['<table border="1" class="dataframe">']
82 s += ['<thead><tr><th style="text-align: left;">Field</th><th>Size</th><th>Data</th></tr></thead>']
83 s += ['<tbody>']
84 with np.printoptions(threshold=4, edgeitems=2, precision=2, linewidth=40):
85 s += ['<tr><td style="text-align: left;">Rotation</td>']
86 s += [f'<td>{self.rotation.shape}</td><td>{self.rotation}</td></tr>']
87 s += ['<tr><td style="text-align: left;">j_offset</td>']
88 s += [f'<td>{1}</td><td>{self.j_offset:.4f}</td></tr>']
89 s += ['<tr><td style="text-align: left;">k_offset</td>']
90 s += [f'<td>{1}</td><td>{self.k_offset:.4f}</td></tr>']
91 s += ['<tr><td style="text-align: left;">inner_angle</td>']
92 s += [f'<td>{1}</td>'
93 f'<td>{self.inner_angle}']
94 s += ['<tr><td style="text-align: left;">outer_angle</td>']
95 s += [f'<td>{1}</td>'
96 f'<td>{self.outer_angle}</td>']
97 s += ['<tr><td style="text-align: left;">inner_axis</td>']
98 s += [f'<td>{self.inner_axis.shape}</td><td>{self.inner_axis}</td></tr>']
99 s += ['<tr><td style="text-align: left;">outer_axis</td>']
100 s += [f'<td>{self.outer_axis.shape}</td><td>{self.outer_axis}</td></tr>']
101 s += [f'<td>{6}</td>'
102 f'<td>{hex(hash(self))[2:8]} (hash)</td></tr>']
103 s += ['</tbody>']
104 s += ['</table>']
105 return '\n'.join(s)
108class Geometry:
109 """ Stores information about the system geometry.
110 Instances of this class are used by
111 :class:`DataContainer <mumott.data_handling.DataContainer>`
112 and :class:`ProjectionStack <mumott.core.projection_stack.ProjectionStack>` to
113 maintain geometry information.
114 They can be stored as a file using the :meth:`write` method.
115 This allows one to (re)create a :class:`Geometry` instance from an earlier
116 and overwrite the geometry information read by a
117 :class:`DataContainer <mumott.data_handling.DataContainer>` instance.
118 This is useful, for example, in the context of alignment.
120 Parameters
121 ----------
122 filename
123 Name of file from which to read geometry information.
124 Defaults to ``None``, in which case the instance is created with
125 default parameters.
126 """
127 def __init__(self, filename: str = None):
128 self._rotations = []
129 self._j_offsets = []
130 self._k_offsets = []
131 self._p_direction_0 = np.array([0, 1, 0]).astype(float)
132 self._j_direction_0 = np.array([1, 0, 0]).astype(float)
133 self._k_direction_0 = np.array([0, 0, 1]).astype(float)
134 self._detector_direction_origin = np.array([1, 0, 0]).astype(float)
135 self._detector_direction_positive_90 = np.array([0, 0, 1]).astype(float)
136 self._projection_shape = np.array([0, 0]).astype(np.int32)
137 self._volume_shape = np.array([0, 0, 0]).astype(np.int32)
138 self._detector_angles = np.array([]).astype(float)
139 self._two_theta = np.array([0.0])
140 self._reconstruction_rotations = []
141 self._system_rotations = []
142 self._inner_angles = []
143 self._outer_angles = []
144 self._inner_axes = []
145 self._outer_axes = []
146 self._full_circle_covered = False
147 if filename is not None:
148 self.read(filename)
150 def write(self, filename: str) -> None:
151 """Method for writing the current state of a :class:`Geometry` instance to file.
153 Notes
154 -----
155 Any rotations in :attr:`reconstruction_rotations` and :attr:`system_rotations`
156 will be applied to the :attr:`rotations` and system vectors respectively prior to writing.
158 Parameters
159 ----------
160 filename
161 Name of output file.
162 """
163 to_write = dict(_rotations=self.rotations_as_array.tolist(),
164 _j_offsets=self._j_offsets,
165 _k_offsets=self._k_offsets,
166 p_direction_0=self.p_direction_0.tolist(),
167 j_direction_0=self.j_direction_0.tolist(),
168 k_direction_0=self.k_direction_0.tolist(),
169 detector_direction_origin=self.detector_direction_origin.tolist(),
170 detector_direction_positive_90=self.detector_direction_positive_90.tolist(),
171 two_theta=self.two_theta.tolist(),
172 projection_shape=self.projection_shape.tolist(),
173 volume_shape=self.volume_shape.tolist(),
174 detector_angles=self.detector_angles.tolist(),
175 full_circle_covered=[bool(self.full_circle_covered)],
176 _inner_angles=self._inner_angles,
177 _outer_angles=self._outer_angles,
178 _inner_axes=[j.tolist() if j is not None else None for j in self._inner_axes],
179 _outer_axes=[j.tolist() if j is not None else None for j in self._outer_axes],
180 checksum=[hash(self)])
181 with tarfile.open(name=filename, mode='w') as tar_file:
182 for key, item in to_write.items():
183 temp_file = tempfile.NamedTemporaryFile(delete=False)
184 temp_file.close()
185 with codecs.open(temp_file.name, 'w', encoding='utf-8') as tf:
186 json.dump(item, tf)
187 tf.flush()
188 with open(temp_file.name, 'rb') as tt:
189 tar_info = tar_file.gettarinfo(arcname=key, fileobj=tt)
190 tar_file.addfile(tar_info, tt)
191 os.remove(temp_file.name)
193 def read(self, filename: str) -> None:
194 """Method for reading the current state of a :class:`Geometry` instance from file.
196 Parameters
197 ----------
198 filename
199 Name of input file.
200 """
201 to_read = ['_rotations',
202 '_j_offsets',
203 '_k_offsets',
204 'p_direction_0',
205 'j_direction_0',
206 'k_direction_0',
207 'detector_direction_origin',
208 'detector_direction_positive_90',
209 'two_theta',
210 'projection_shape',
211 'volume_shape',
212 'detector_angles',
213 'full_circle_covered',
214 '_inner_angles',
215 '_outer_angles',
216 '_inner_axes',
217 '_outer_axes',
218 'checksum']
219 with tarfile.open(name=filename, mode='r') as tar_file:
220 for key in to_read:
221 temp_file = tempfile.NamedTemporaryFile(delete=False)
222 try:
223 temp_file.write(tar_file.extractfile(key).read())
224 except KeyError:
225 logger.warning(f'Key {key} not found!')
226 temp_file.close()
227 continue
228 temp_file.close()
229 with codecs.open(temp_file.name, 'r', encoding='utf-8') as f:
230 text = f.read()
231 data_as_list = json.loads(text)
232 for i, entry in enumerate(data_as_list):
233 if entry == 'null': 233 ↛ 234line 233 didn't jump to line 234, because the condition on line 233 was never true
234 data_as_list[i] = None
235 if key == 'checksum':
236 checksum = data_as_list[0]
237 elif key in ('_rotations', '_inner_axes', '_outer_axes'):
238 setattr(self, key, [])
239 for entry in data_as_list:
240 # necessary check to avoid [array(None), ...] structure
241 if entry is not None:
242 entry = np.array(entry)
243 getattr(self, key).append(entry)
244 elif key in ('_j_offsets', '_k_offsets', '_inner_angles', '_outer_angles'):
245 setattr(self, key, data_as_list)
246 elif key in ('full_circle_covered'):
247 setattr(self, key, data_as_list[0])
248 else:
249 setattr(self, key, np.array(data_as_list))
250 if checksum != hash(self): 250 ↛ 251line 250 didn't jump to line 251, because the condition on line 250 was never true
251 logger.warning(f'Checksum does not match! Checksum is {checksum},'
252 f' but hash(self) is {hash(self)}. This may be due to'
253 ' version differences, but please proceed with caution!')
255 def rotate_reconstruction(self,
256 A: np.ndarray[float] = None,
257 axis: np.ndarray[float] = None,
258 angle: np.ndarray[float] = None):
259 r""" Rotates the reconstruction geometry. The given rotation matrix will modify the rotation
260 matrix of each projection by multiplication from the right, such that
262 .. math ::
263 R_i' = R_i A
265 where :math:`R_i` is the rotation matrix of projection :math:`i` and :math:`A` is the rotation matrix.
266 For each projection, the system vectors are then rotated by
268 .. math ::
269 v_i = (R_i')^T v = A^T R_i^T v
271 where :math:`v` corresponds to e.g., :attr:`p_direction_0`.
273 Notes
274 -----
275 It is not possible to directly modify :attr:`rotations` after adding a reconstruction rotation.
277 Parameters
278 ----------
279 A
280 A 3-by-3 rotation matrix. If not given, then :attr:`axis` and :attr:`angle` must be provided.
281 axis
282 An axis, given as a unit length 3-vector, about which the rotation is defined. Not used
283 if :attr:`A` is provided.
284 angle
285 The angle in radians of the rotation about :attr:`axis`. Not used if :attr:`A` is provided.
286 """
287 if A is None:
288 A = Rotation.from_rotvec(axis * angle / np.linalg.norm(axis)).as_matrix()
289 elif axis is not None or angle is not None:
290 logger.warning('A provided along with axis and/or angle; axis/angle will be ignored!')
292 self._reconstruction_rotations.append(A)
294 def rotate_system_vectors(self,
295 A: np.ndarray[float] = None,
296 axis: np.ndarray[float] = None,
297 angle: np.ndarray[float] = None):
298 r""" Rotates the system vectors. The given rotation matrix will modify the system vectors by
300 .. math ::
301 v' = A v
303 where :math:`v` is a system vector, e.g., :attr:`p_direction_0`, and :math:`A` is the rotation matrix.
304 For each projection, the system vectors are then rotated by
306 .. math ::
307 v_i = R_i^T A v
309 where :math:`R_i` corresponds to :attr:`rotations` for projection :math:`i`.
311 Notes
312 -----
313 It is not possible to directly modify the system vectors after adding a system rotation.
315 Parameters
316 ----------
317 A
318 A 3-by-3 rotation matrix. If not given, then :attr:`axis` and :attr:`angle` must be provided.
319 axis
320 An axis, given as a 3-vector, about which a rotation can be defined. Not used
321 if :attr:`A` is provided.
322 angle
323 The angle in radians of the rotation about :attr:`axis`. Not used if :attr:`A` is provided.
324 """
325 if A is None:
326 A = Rotation.from_rotvec(axis * angle / np.linalg.norm(axis)).as_matrix()
327 elif axis is not None or angle is not None:
328 logger.warning('A provided along with axis and/or angle; axis/angle will be ignored!')
330 self._system_rotations.append(A)
332 def append(self, value: GeometryTuple) -> None:
333 """ Appends projection-wise geometry data provided as a
334 :class:`GeometryTuple <mumott.core.geometry.GeometryTuple>`. """
335 self._rotations.append(value.rotation)
336 self._j_offsets.append(value.j_offset)
337 self._k_offsets.append(value.k_offset)
338 self._inner_angles.append(value.inner_angle)
339 self._outer_angles.append(value.outer_angle)
340 self._inner_axes.append(value.inner_axis)
341 self._outer_axes.append(value.outer_axis)
343 def insert(self, key: int, value: GeometryTuple) -> None:
344 """ Inserts projection-wise data handed via a
345 :class:`GeometryTuple <mumott.core.geometry.GeometryTuple>`. """
346 self._rotations.insert(key, value.rotation)
347 self._j_offsets.insert(key, value.j_offset)
348 self._k_offsets.insert(key, value.k_offset)
349 self._inner_angles.insert(key, value.inner_angle)
350 self._outer_angles.insert(key, value.outer_angle)
351 self._inner_axes.insert(key, value.inner_axis)
352 self._outer_axes.insert(key, value.outer_axis)
354 def __setitem__(self, key: int, value: GeometryTuple) -> None:
355 """ Sets projection-wise data handed via a :class:`GeometryTuple`."""
356 self._rotations[key] = value.rotation
357 self._j_offsets[key] = value.j_offset
358 self._k_offsets[key] = value.k_offset
359 self._inner_angles[key] = value.inner_angle
360 self._outer_angles[key] = value.outer_angle
361 self._inner_axes[key] = value.inner_axis
362 self._outer_axes[key] = value.outer_axis
364 def __getitem__(self, key: int) -> GeometryTuple:
365 """ Returns projection-wise data as a :class:`GeometryTuple`."""
366 return GeometryTuple(rotation=self.rotations[key],
367 j_offset=self._j_offsets[key],
368 k_offset=self._k_offsets[key],
369 inner_angle=self._inner_angles[key],
370 outer_angle=self._outer_angles[key],
371 inner_axis=self._inner_axes[key],
372 outer_axis=self._outer_axes[key])
374 def __delitem__(self, key: int) -> None:
375 del self._rotations[key]
376 del self._j_offsets[key]
377 del self._k_offsets[key]
378 del self._inner_angles[key]
379 del self._outer_angles[key]
380 del self._inner_axes[key]
381 del self._outer_axes[key]
383 def _get_probed_coordinates(self) -> NDArray[np.float_]:
384 """
385 Calculates and returns the probed polar and azimuthal coordinates on the unit sphere at
386 each angle of projection and for each detector segment in the geometry of the system.
387 """
388 n_proj = len(self)
389 n_seg = len(self.detector_angles)
390 probed_directions_zero_rot = np.zeros((n_seg, 3, 3))
391 # Impose symmetry if needed.
392 if not self.full_circle_covered: 392 ↛ 395line 392 didn't jump to line 395, because the condition on line 392 was never false
393 shift = np.pi
394 else:
395 shift = 0
396 det_bin_middles_extended = np.copy(self.detector_angles)
397 det_bin_middles_extended = np.insert(det_bin_middles_extended, 0,
398 det_bin_middles_extended[-1] + shift)
399 det_bin_middles_extended = np.append(det_bin_middles_extended, det_bin_middles_extended[1] + shift)
401 for ii in range(n_seg):
403 # Check if the interval from the previous to the next bin goes over the -pi +pi discontinuity
404 before = det_bin_middles_extended[ii]
405 now = det_bin_middles_extended[ii + 1]
406 after = det_bin_middles_extended[ii + 2]
408 if abs(before - now + 2 * np.pi) < abs(before - now): 408 ↛ 409line 408 didn't jump to line 409, because the condition on line 408 was never true
409 before = before + 2 * np.pi
410 elif abs(before - now - 2 * np.pi) < abs(before - now):
411 before = before - 2 * np.pi
413 if abs(now - after + 2 * np.pi) < abs(now - after): 413 ↛ 414line 413 didn't jump to line 414, because the condition on line 413 was never true
414 after = after - 2 * np.pi
415 elif abs(now - after - 2 * np.pi) < abs(now - after): 415 ↛ 416line 415 didn't jump to line 416, because the condition on line 415 was never true
416 after = after + 2 * np.pi
418 # Generate a linearly spaced set of angles covering the detector segment
419 start = 0.5 * (before + now)
420 end = 0.5 * (now + after)
421 angles = np.linspace(start, end, 3)
423 # Make the zero-rotation-projection vectors corresponding to the given angles
424 probed_directions_zero_rot[ii, :, :] = np.cos(angles[:, np.newaxis]) * \
425 self.detector_direction_origin[np.newaxis, :]
426 probed_directions_zero_rot[ii, :, :] += np.sin(angles[:, np.newaxis]) * \
427 self.detector_direction_positive_90[np.newaxis, :]
429 # Do wide-angle rotation
430 n_twotheta = len(self.two_theta)
431 twotheta = np.repeat(self.two_theta, n_seg)[:, np.newaxis, np.newaxis]
432 probed_directions_zero_rot = np.tile(probed_directions_zero_rot, (n_twotheta, 1, 1))
434 probed_directions_zero_rot = probed_directions_zero_rot * np.cos(twotheta/2)\
435 - np.sin(twotheta/2) * self.p_direction_0
437 # Initialize array for vectors
438 probed_direction_vectors = np.zeros((n_proj,
439 n_seg * n_twotheta,
440 3,
441 3), dtype=np.float64)
442 # Calculate all the rotations
443 probed_direction_vectors[...] = \
444 np.einsum('kij,mli->kmlj',
445 self.rotations,
446 probed_directions_zero_rot)
447 great_circle_offsets = np.einsum('kij,mli->kmlj',
448 self.rotations,
449 -np.sin(twotheta / 2) * self.p_direction_0)
450 return ProbedCoordinates(probed_direction_vectors, great_circle_offsets)
452 def delete_projections(self) -> None:
453 """ Delete all projections."""
454 self._rotations = []
455 self._j_offsets = []
456 self._k_offsets = []
457 self._inner_angles = []
458 self._outer_angles = []
459 self._inner_axes = []
460 self._outer_axes = []
462 def _get_reconstruction_rotation(self) -> np.ndarray[float]:
463 """ Internal method for composing reconstruction rotations. """
464 reconstruction_rotation = np.eye(3, dtype=float)
465 for r in self.reconstruction_rotations:
466 reconstruction_rotation = reconstruction_rotation @ r
467 return reconstruction_rotation
469 def _get_system_rotation(self) -> np.ndarray[float]:
470 """ Internal method for composing system rotations. """
471 system_rotation = np.eye(3, dtype=float)
472 for r in self.system_rotations:
473 system_rotation = r @ system_rotation
474 return system_rotation
476 @property
477 def system_rotations(self) -> list[np.ndarray[float]]:
478 """ list of rotation matrices sequentially applied to the basis vectors of the system. """
479 return self._system_rotations
481 @system_rotations.setter
482 def system_rotations(self, value: list[np.ndarray[float]]) -> list[np.ndarray[float]]:
483 self._system_rotations = list(value)
485 @property
486 def reconstruction_rotations(self) -> list[np.ndarray[float]]:
487 """ list of rotation matrices sequentially applied to the reconstruction geometry of the system. """
488 return self._reconstruction_rotations
490 @reconstruction_rotations.setter
491 def reconstruction_rotations(self, value: list[np.ndarray[float]]) -> list[np.ndarray[float]]:
492 self._reconstruction_rotations = list(value)
494 @property
495 def rotations(self) -> list[np.ndarray[float]]:
496 """ Rotation matrices for the experimental rotation corresponding to each projection of data."""
497 if len(self.reconstruction_rotations) > 0:
498 reconstruction_rotation = self._get_reconstruction_rotation()
499 return [r @ reconstruction_rotation for r in self._rotations]
501 return self._rotations
503 @property
504 def rotations_as_array(self) -> np.ndarray[float]:
505 """ Rotation matrices corresponding to each projection of data as an array."""
506 if len(self) == 0:
507 return np.array([])
508 return np.stack(list(self.rotations), axis=0)
510 @rotations.setter
511 def rotations(self, value: Union[list, np.ndarray[float]]) -> None:
512 if len(self._reconstruction_rotations) > 0:
513 raise ValueError('Cannot modify rotations when reconstruction '
514 'rotations are in use.')
515 self._rotations = list(value)
517 @property
518 def p_direction_0(self) -> np.ndarray[float]:
519 """ The projection direction when no experimental rotation is applied."""
520 if len(self._system_rotations) > 0:
521 system_rotation = self._get_system_rotation()
522 return system_rotation @ self._p_direction_0
524 return self._p_direction_0
526 @p_direction_0.setter
527 def p_direction_0(self, value: np.ndarray[float]) -> None:
528 if len(self.system_rotations) > 0:
529 raise ValueError('Cannot modify system vectors when system '
530 'rotations are in use.')
531 if np.size(value) != 3:
532 raise ValueError('The size of the new value must be 3, but '
533 f'the provided value has size {np.size(value)}')
534 self._p_direction_0[...] = value
536 @property
537 def j_direction_0(self) -> np.ndarray[float]:
538 """ The direction corresponding to the first index in each projection
539 when no experimental rotation is applied."""
540 if len(self._system_rotations) > 0:
541 system_rotation = np.eye(3)
542 for r in self.system_rotations:
543 system_rotation = system_rotation @ r
544 return system_rotation @ self._j_direction_0
546 return self._j_direction_0
548 @j_direction_0.setter
549 def j_direction_0(self, value: np.ndarray[float]) -> None:
550 if len(self.system_rotations) > 0:
551 raise ValueError('Cannot modify system vectors when system '
552 'rotations are in use.')
553 if np.size(value) != 3:
554 raise ValueError('The size of the new value must be 3, but '
555 f'the provided value has size {np.size(value)}')
556 self._j_direction_0[...] = value
558 @property
559 def k_direction_0(self) -> np.ndarray[float]:
560 """ The direction corresponding to the second index in each projection
561 when no experimental rotation is applied."""
562 if len(self._system_rotations) > 0:
563 system_rotation = self._get_system_rotation()
564 return system_rotation @ self._k_direction_0
566 return self._k_direction_0
568 @k_direction_0.setter
569 def k_direction_0(self, value: np.ndarray[float]) -> None:
570 if len(self.system_rotations) > 0:
571 raise ValueError('Cannot modify system vectors when system '
572 'rotations are in use.')
573 if np.size(value) != 3:
574 raise ValueError('The size of the new value must be 3, but '
575 f'the provided value has size {np.size(value)}')
576 self._k_direction_0[...] = value
578 @property
579 def detector_direction_origin(self) -> np.ndarray[float]:
580 """ The direction at which the angle on the detector is zero,
581 when no experimental rotation is applied."""
582 if len(self._system_rotations) > 0:
583 system_rotation = self._get_system_rotation()
584 return system_rotation @ self._detector_direction_origin
586 return self._detector_direction_origin
588 @detector_direction_origin.setter
589 def detector_direction_origin(self, value: np.ndarray[float]) -> None:
590 if len(self.system_rotations) > 0:
591 raise ValueError('Cannot modify system vectors when system '
592 'rotations are in use.')
593 if np.size(value) != 3:
594 raise ValueError('The size of the new value must be 3, but '
595 f'the provided value has size {np.size(value)}')
596 self._detector_direction_origin[...] = value
598 @property
599 def detector_direction_positive_90(self) -> np.ndarray[float]:
600 """ Rotation matrices corresponding to each projection of data."""
601 if len(self._system_rotations) > 0:
602 system_rotation = self._get_system_rotation()
603 return system_rotation @ self._detector_direction_positive_90
605 return self._detector_direction_positive_90
607 @detector_direction_positive_90.setter
608 def detector_direction_positive_90(self, value: np.ndarray[float]) -> None:
609 if len(self.system_rotations) > 0:
610 raise ValueError('Cannot modify system vectors when system '
611 'rotations are in use.')
612 if np.size(value) != 3:
613 raise ValueError('The size of the new value must be 3, but '
614 f'the provided value has size {np.size(value)}')
615 self._detector_direction_positive_90[...] = value
617 @property
618 def j_offsets(self) -> list[float]:
619 """Offsets to align projection in the direction j."""
620 return self._j_offsets
622 @property
623 def j_offsets_as_array(self) -> np.ndarray[float]:
624 """Offsets to align projection in the direction j as an array."""
625 if len(self._j_offsets) == 0:
626 return np.array([])
627 return np.stack(self.j_offsets, axis=0)
629 @j_offsets.setter
630 def j_offsets(self, value: Union[list[float], np.ndarray[float]]) -> None:
631 self._j_offsets = list(value)
633 @property
634 def outer_angles(self) -> list[float]:
635 """Rotation angles for inner rotation, in radians."""
636 return list(self._outer_angles)
638 @property
639 def outer_angles_as_array(self) -> np.ndarray[float]:
640 """Rotation angles for inner rotations, in radians, as an array."""
641 if len(self._outer_angles) == 0:
642 return np.array([])
643 return np.stack(self.outer_angles, axis=0)
645 @outer_angles.setter
646 def outer_angles(self, value: Union[list[float], np.ndarray[float]]) -> None:
647 self._outer_angles = list(value)
648 self._update_rotations()
650 @property
651 def inner_angles(self) -> list[float]:
652 """Rotation angles for inner rotation, in radians."""
653 return self._inner_angles
655 @property
656 def probed_coordinates(self) -> ProbedCoordinates:
657 """ An array of 3-vectors with the (x, y, z)-coordinates
658 on the reciprocal space map probed by the method.
659 Structured as ``(N, K, 3, 3)``, where ``N``
660 is the number of projections, ``K`` is the number of
661 detector segments, the second-to-last axis contains
662 start-, mid-, and endpoints, and the last axis contains the
663 (x, y, z)-coordinates.
665 Notes
666 -----
667 The number of detector segments is
668 `len(geometry.detecor_angles)*len(geometry.two_theta)`
669 i.e. the product of the number of two_theta bins times the number of
670 azimuthal bins. As a default, only on two theta bin is used.
671 When several two_theta bins are used, the second index corresponds
672 to a raveled array, where the azimuthal is the fast index and
673 two theta is the slow index.
674 """
675 return self._get_probed_coordinates()
677 def _get_hashable_axes_and_angles(self) -> tuple[np.ndarray[float], ...]:
678 """ Internal method for getting hashable ste of axes and angles,
679 as well as checking if set is valid or if it contains any ``None``
680 entries."""
681 attributes = [self._inner_angles,
682 self._outer_angles,
683 self._inner_axes,
684 self._outer_axes]
685 array_props = [self.inner_angles_as_array,
686 self.outer_angles_as_array,
687 self.inner_axes_as_array,
688 self.outer_axes_as_array]
689 return_values = []
690 for elements, arrays in zip(attributes, array_props):
691 for a in elements:
692 if a is None:
693 return_values.append(None)
694 break
695 else:
696 return_values.append(arrays)
697 return tuple(return_values)
699 @property
700 def inner_angles_as_array(self) -> np.ndarray[float]:
701 """Rotation angles for inner rotations, in radians, as an array."""
702 if len(self._inner_angles) == 0:
703 return np.array([])
704 return np.stack(self.inner_angles, axis=0)
706 @inner_angles.setter
707 def inner_angles(self, value: Union[list[float], np.ndarray]) -> None:
708 self._inner_angles = [j for j in value]
709 self._update_rotations()
711 @property
712 def inner_axes(self) -> list[np.ndarray[float]]:
713 """Inner rotation axes. All axes can be set
714 at once using a single array with three entries."""
715 return self._inner_axes
717 @property
718 def inner_axes_as_array(self) -> np.ndarray[float]:
719 """Inner rotation axes as an array."""
720 if len(self._inner_axes) == 0:
721 return np.array([])
722 return np.stack([j if j is not None else [None, None, None] for j in self._inner_axes], axis=0)
724 @inner_axes.setter
725 def inner_axes(self, value: Union[list[float], np.ndarray]) -> None:
726 value = np.array(value)
727 if value.ndim == 1:
728 if value.shape != (3,):
729 raise ValueError('inner_axes may be set using either '
730 'a list/array of size-3 arrays or an array '
731 'of shape (3,), but the provided array has shape '
732 f'{value.shape}.')
733 self._inner_axes = [value for _ in self._inner_axes]
734 elif value.ndim == 2:
735 if len(value) != len(self):
736 raise ValueError('If inner_axes is set using a list/array of '
737 'size-3 arrays, then it must be of the same length '
738 f'as the Geometry instance ({len(self)}), '
739 f'but it is of length {len(value)}.')
740 if value[0].size != 3:
741 raise ValueError('inner_axes may be set using either '
742 'a list/array of size-3 arrays or an array '
743 'of shape 3, but the provided array has shape '
744 f'{value.shape}.')
745 self._inner_axes = [j for j in value]
746 else:
747 raise ValueError('inner_axes must be set either with a list/array of '
748 f'shape ({len(self)}, 3) or with an array of shape '
749 f'(3,), but the provided array has shape {value.shape}!')
750 self._update_rotations()
752 @property
753 def outer_axes(self) -> list[np.ndarray[float]]:
754 """Inner rotation axes. All axes can be set
755 at once using a single array with three entries."""
756 return self._outer_axes
758 @property
759 def outer_axes_as_array(self) -> np.ndarray[float]:
760 """Outer rotation axes as an array."""
761 if len(self._outer_axes) == 0:
762 return np.array([])
763 return np.stack([j if j is not None else [None, None, None] for j in self._outer_axes], axis=0)
765 @outer_axes.setter
766 def outer_axes(self, value: Union[list[np.ndarray[float]], np.ndarray[float]]) -> None:
767 value = np.array(value)
768 if value.ndim == 1:
769 if value.shape != (3,):
770 raise ValueError('outer_axes may be set using either '
771 'a list/array of size-3 arrays or an array '
772 'of shape (3,), but the provided array has shape '
773 f'{value.shape}.')
774 self._outer_axes = [value for j in self._outer_axes]
775 elif value.ndim == 2:
776 if len(value) != len(self):
777 raise ValueError('If outer_axes is set using a list/array of '
778 'size-3 arrays, then it must be of the same length '
779 f'as the Geometry instance ({len(self)}), but it is of length {len(value)}.')
780 if value[0].size != 3:
781 raise ValueError('outer_axes may be set using either '
782 'a list/array of size-3 arrays or an array '
783 'of shape 3, but the provided array has shape '
784 f'{value.shape}.')
785 self._outer_axes = [j for j in value]
786 else:
787 raise ValueError('outer_axes must be set either with a list/array of '
788 f'shape ({len(self)}, 3) or with an array of shape '
789 f'(3,), but the provided array has shape {value.shape}!')
790 self._update_rotations()
792 def _update_rotations(self) -> None:
793 """Internal method for updating rotations based on changes
794 to inner and outer axes. """
795 can_update_rotations = np.all([h is not None for h in self._get_hashable_axes_and_angles()])
796 if not can_update_rotations:
797 logger.info('None values found in some axis or angle entries,'
798 ' rotations not updated.')
799 return
801 for i in range(len(self)):
802 R_inner = Rotation.from_rotvec(self.inner_angles[i] * self.inner_axes[i]).as_matrix()
803 R_outer = Rotation.from_rotvec(self.outer_angles[i] * self.outer_axes[i]).as_matrix()
804 self.rotations[i] = R_outer @ R_inner
806 @property
807 def k_offsets(self) -> list[float]:
808 """Offsets to align projection in the direction k."""
809 return self._k_offsets
811 @property
812 def k_offsets_as_array(self) -> np.ndarray[float]:
813 """Offsets to align projection in the direction k as an array."""
814 if len(self._k_offsets) == 0:
815 return np.array([])
816 return np.stack(self.k_offsets, axis=0)
818 @property
819 def hash_rotations(self) -> str:
820 """ A blake2b hash of :attr:`rotations_as_array`. """
821 return list_to_hash([self.rotations_as_array])
823 @property
824 def hash_j_offsets(self) -> str:
825 """ A blake2b hash of :attr:`j_offsets_as_array`. """
826 return list_to_hash([self.j_offsets_as_array])
828 @property
829 def hash_k_offsets(self) -> str:
830 """ A blake2b hash of :attr:`k_offsets_as_array`. """
831 return list_to_hash([self.k_offsets_as_array])
833 @property
834 def hash_inner_angles(self) -> str:
835 """ A blake2b hash of :attr:`inner_angle`. """
836 return list_to_hash(self.inner_angles)
838 @property
839 def hash_outer_angles(self) -> str:
840 """ A blake2b hash of :attr:`outer_anglesy`. """
841 return list_to_hash(self.outer_angles)
843 @property
844 def hash_inner_axes(self) -> str:
845 """ A blake2b hash of :attr:`inner_axes`. """
846 return list_to_hash(self.inner_axes)
848 @property
849 def hash_outer_axes(self) -> str:
850 """ A blake2b hash of :attr:`outer_axes`. """
851 return list_to_hash(self.outer_axes)
853 @property
854 def projection_shape(self) -> NDArray[int]:
855 """ 2D shape of the raster-scan. 1st element is the number of steps in the
856 j-direction and the second is the number of steps in the k-direction."""
857 return self._projection_shape
859 @projection_shape.setter
860 def projection_shape(self, value: NDArray[int]) -> None:
861 if type(value) is not tuple:
862 value = tuple(value)
864 if len(value) != 2:
865 raise ValueError('Length of projection_shape must be exactly 2.')
867 if not all(isinstance(item, (int, np.integer)) for item in value):
868 first_wrong = next((item for item in value if type(item) is not int)) 868 ↛ exitline 868 didn't finish the generator expression on line 868
869 raise TypeError(f'{type(first_wrong)} cannot be interpreted as an integer.')
871 self._projection_shape = np.array(value).astype(np.int32)
873 @property
874 def volume_shape(self) -> NDArray[int]:
875 """ 3D shape of the reconstruction voxel array. 1st element is the number of points
876 along the x-direction. 2nd is y and 3rd is z."""
877 return self._volume_shape
879 @volume_shape.setter
880 def volume_shape(self, value: NDArray[int]) -> None:
881 if type(value) is not tuple:
882 value = tuple(value)
884 if len(value) != 3:
885 raise ValueError('Length of volume_shape must be exactly 2.')
887 if not all(isinstance(item, (int, np.integer)) for item in value):
888 first_wrong = next((item for item in value if type(item) is not int)) 888 ↛ exitline 888 didn't finish the generator expression on line 888
889 raise TypeError(f'{type(first_wrong)} cannot be interpreted as an integer.')
891 self._volume_shape = np.array(value).astype(np.int32)
893 @property
894 def detector_angles(self) -> NDArray(float):
895 """ Azimuthal angles of detector segments in radians.
896 One-dimensional sequence of center positions"""
897 return self._detector_angles
899 @detector_angles.setter
900 def detector_angles(self, value: NDArray(float)) -> None:
901 value = np.array(value).astype(float)
902 if np.ndim(value) != 1:
903 raise ValueError('Detector angles must be a one-dimensional sequence.')
904 self._detector_angles = value
906 @property
907 def two_theta(self) -> NDArray(float):
908 """ Scattering angle in radians. Can be list of angles, if multiple
909 radial bins are used."""
910 return self._two_theta
912 @two_theta.setter
913 def two_theta(self, value: NDArray(float)) -> None:
915 value = np.array(value)
917 if np.ndim(value) == 0:
918 value = value.flatten()
919 elif np.ndim(value) > 1:
920 raise ValueError('Only scalars or one-dimensional sequences are valid values of two_theta.')
922 self._two_theta = value.astype(float)
924 @property
925 def full_circle_covered(self) -> bool:
926 """ Whether the azimuthal bins cover a half-circle of the detector (False)
927 or the full circle (True). """
928 return self._full_circle_covered
930 @full_circle_covered.setter
931 def full_circle_covered(self, value: bool) -> None:
932 self._full_circle_covered = bool(value)
934 @k_offsets.setter
935 def k_offsets(self, value: Union[list[float], np.ndarray[float]]) -> None:
936 self._k_offsets = list(value)
938 def __hash__(self) -> int:
939 to_hash = [self.rotations_as_array,
940 self.j_offsets_as_array,
941 self.k_offsets_as_array,
942 self.p_direction_0, self.j_direction_0, self.k_direction_0,
943 self.detector_direction_origin, self.detector_direction_positive_90,
944 self.two_theta,
945 self.projection_shape, self.volume_shape,
946 self.detector_angles, self.full_circle_covered,
947 *self._get_hashable_axes_and_angles()]
948 return int(list_to_hash(to_hash), 16)
950 def __len__(self) -> int:
951 return len(self._rotations)
953 def _get_str_representation(self, max_lines: int = 25) -> str:
954 """ Retrieves a string representation of the object with specified
955 maximum number of lines.
957 Parameters
958 ----------
959 max_lines
960 The maximum number of lines to return.
961 """
962 wdt = 74
963 s = []
964 s += ['-' * wdt]
965 s += ['Geometry'.center(wdt)]
966 s += ['-' * wdt]
967 with np.printoptions(threshold=3, edgeitems=1, precision=3, linewidth=60):
968 s += ['{:18} : {}'.format('hash_rotations',
969 self.hash_rotations[:6])]
970 s += ['{:18} : {}'.format('hash_j_offsets',
971 self.hash_j_offsets[:6])]
972 s += ['{:18} : {}'.format('hash_k_offsets',
973 self.hash_k_offsets[:6])]
974 s += ['{:18} : {}'.format('p_direction_0', self.p_direction_0)]
975 s += ['{:18} : {}'.format('j_direction_0', self.j_direction_0)]
976 s += ['{:18} : {}'.format('k_direction_0', self.k_direction_0)]
977 s += ['{:18} : {}'.format('hash_inner_angles', self.hash_inner_angles[:6])]
978 s += ['{:18} : {}'.format('hash_outer_angles', self.hash_outer_angles[:6])]
979 s += ['{:18} : {}'.format('hash_inner_axes', self.hash_inner_axes[:6])]
980 s += ['{:18} : {}'.format('hash_outer_axes', self.hash_outer_axes[:6])]
981 s += ['{:18} : {}'.format('detector_direction_origin', self.detector_direction_origin)]
982 s += ['{:18} : {}'.format('detector_direction_positive_90', self.detector_direction_positive_90)]
983 s += ['{:18} : {}°'.format('two_theta', np.rad2deg(self.two_theta))]
984 s += ['{:18} : {}'.format('projection_shape', self.projection_shape)]
985 s += ['{:18} : {}'.format('volume_shape', self.volume_shape)]
986 s += ['{:18} : {}'.format('detector_angles', self.detector_angles)]
987 s += ['{:18} : {}'.format('full_circle_covered', self.full_circle_covered)]
988 s += ['-' * wdt]
989 truncated_s = []
990 leave_loop = False
991 while not leave_loop:
992 line = s.pop(0).split('\n')
993 for split_line in line:
994 if split_line != '': 994 ↛ 996line 994 didn't jump to line 996, because the condition on line 994 was never false
995 truncated_s += [split_line]
996 if len(truncated_s) > max_lines - 2: 996 ↛ 997line 996 didn't jump to line 997, because the condition on line 996 was never true
997 if split_line != '...':
998 truncated_s += ['...']
999 if split_line != ('=' * wdt):
1000 truncated_s += ['=' * wdt]
1001 leave_loop = True
1002 break
1003 if len(s) == 0:
1004 leave_loop = True
1005 return '\n'.join(truncated_s)
1007 def __str__(self) -> str:
1008 return self._get_str_representation()
1010 def _get_html_representation(self, max_lines: int = 25) -> str:
1011 """ Retrieves an html representation of the object with specified
1012 maximum number of lines.
1014 Parameters
1015 ----------
1016 max_lines
1017 The maximum number of lines to return.
1018 """
1019 s = []
1020 s += ['<h3>Geometry</h3>']
1021 s += ['<table border="1" class="dataframe">']
1022 s += ['<thead><tr><th style="text-align: left;">Field</th><th>Size</th><th>Data</th></tr></thead>']
1023 s += ['<tbody>']
1024 with np.printoptions(threshold=3, edgeitems=1, precision=2, linewidth=40):
1025 s += ['<tr><td style="text-align: left;">rotations</td>']
1026 s += [f'<td>{len(self.rotations)}</td>'
1027 f'<td>{self.hash_rotations[:6]} (hash)</td></tr>']
1028 s += ['<tr><td style="text-align: left;">j_offsets</td>']
1029 s += [f'<td>{len(self.j_offsets)}</td>'
1030 f'<td>{self.hash_j_offsets[:6]} (hash)</td></tr>']
1031 s += ['<tr><td style="text-align: left;">k_offsets</td>']
1032 s += [f'<td>{len(self.k_offsets)}</td>'
1033 f'<td>{self.hash_k_offsets[:6]} (hash)</td></tr>']
1034 s += ['<tr><td style="text-align: left;">p_direction_0</td>']
1035 s += [f'<td>{len(self.p_direction_0)}</td><td>{self.p_direction_0}</td></tr>']
1036 s += ['<tr><td style="text-align: left;">j_direction_0</td>']
1037 s += [f'<td>{len(self.j_direction_0)}</td><td>{self.j_direction_0}</td></tr>']
1038 s += ['<tr><td style="text-align: left;">k_direction_0</td>']
1039 s += [f'<td>{len(self.k_direction_0)}</td><td>{self.k_direction_0}</td></tr>']
1040 s += ['<tr><td style="text-align: left;">inner_angles</td>']
1041 s += [f'<td>{len(self.inner_angles)}</td>'
1042 f'<td>{self.hash_inner_angles[:6]} (hash)</td></tr>']
1043 s += ['<tr><td style="text-align: left;">outer_angles</td>']
1044 s += [f'<td>{len(self.outer_angles)}</td>'
1045 f'<td>{self.hash_outer_angles[:6]} (hash)</td></tr>']
1046 s += ['<tr><td style="text-align: left;">inner_axes</td>']
1047 s += [f'<td>{len(self.inner_axes)}</td>'
1048 f'<td>{self.hash_inner_axes[:6]} (hash)</td></tr>']
1049 s += ['<tr><td style="text-align: left;">outer_axes</td>']
1050 s += [f'<td>{len(self.outer_axes)}</td>'
1051 f'<td>{self.hash_outer_axes[:6]} (hash)</td></tr>']
1052 s += ['<tr><td style="text-align: left;">detector_direction_origin</td>']
1053 s += [f'<td>{len(self.detector_direction_origin)}</td>'
1054 f'<td>{self.detector_direction_origin}</td></tr>']
1055 s += ['<tr><td style="text-align: left;">detector_direction_positive_90</td>']
1056 s += [f'<td>{len(self.detector_direction_positive_90)}</td>'
1057 f'<td>{self.detector_direction_positive_90}</td></tr>']
1058 s += ['<tr><td style="text-align: left;">two_theta</td>']
1059 s += [f'<td>{len(self.two_theta)}</td>'
1060 '<td>' + f'{self.two_theta * 180 / np.pi}' + r'${}^{\circ}$</td>']
1061 s += ['<tr><td style="text-align: left;">projection_shape</td>']
1062 s += [f'<td>{len(self.projection_shape)}</td><td>{self.projection_shape}</td></tr>']
1063 s += ['<tr><td style="text-align: left;">volume_shape</td>']
1064 s += [f'<td>{len(self.volume_shape)}</td><td>{self.volume_shape}</td></tr>']
1065 s += ['<tr><td style="text-align: left;">detector_angles</td>']
1066 s += [f'<td>{len(self.detector_angles)}</td><td>{self.detector_angles}</td></tr>']
1067 s += ['<tr><td style="text-align: left;">full_circle_covered</td>']
1068 s += [f'<td>{1}</td><td>{self.full_circle_covered}</td></tr>']
1069 s += ['</tbody>']
1070 s += ['</table>']
1071 truncated_s = []
1072 line_count = 0
1073 leave_loop = False
1074 while not leave_loop:
1075 line = s.pop(0).split('\n')
1076 for split_line in line:
1077 truncated_s += [split_line]
1078 if '</tr>' in split_line:
1079 line_count += 1
1080 # Catch if last line had ellipses
1081 last_tr = split_line
1082 if line_count > max_lines - 1: 1082 ↛ 1083line 1082 didn't jump to line 1083, because the condition on line 1082 was never true
1083 if last_tr != '<tr><td style="text-align: left;">...</td></tr>':
1084 truncated_s += ['<tr><td style="text-align: left;">...</td></tr>']
1085 truncated_s += ['</tbody>']
1086 truncated_s += ['</table>']
1087 leave_loop = True
1088 break
1089 if len(s) == 0:
1090 leave_loop = True
1091 return '\n'.join(truncated_s)
1093 def _repr_html_(self) -> str:
1094 return self._get_html_representation()