diff --git a/doc_sources/apidocs/corrct/corrct.alignment.cone_beam.md b/doc_sources/apidocs/corrct/corrct.alignment.cone_beam.md new file mode 100644 index 0000000..101264a --- /dev/null +++ b/doc_sources/apidocs/corrct/corrct.alignment.cone_beam.md @@ -0,0 +1,293 @@ +# {py:mod}`corrct.alignment.cone_beam` + +```{py:module} corrct.alignment.cone_beam +``` + +```{autodoc2-docstring} corrct.alignment.cone_beam +:allowtitles: +``` + +## Module Contents + +### Classes + +````{list-table} +:class: autosummary longtable +:align: left + +* - {py:obj}`ConeBeamGeometry ` + - ```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry + :summary: + ``` +* - {py:obj}`FitConeBeamGeometry ` + - ```{autodoc2-docstring} corrct.alignment.cone_beam.FitConeBeamGeometry + :summary: + ``` +```` + +### Functions + +````{list-table} +:class: autosummary longtable +:align: left + +* - {py:obj}`_class_to_json ` + - ```{autodoc2-docstring} corrct.alignment.cone_beam._class_to_json + :summary: + ``` +* - {py:obj}`_get_rot_axis_angle_deg ` + - ```{autodoc2-docstring} corrct.alignment.cone_beam._get_rot_axis_angle_deg + :summary: + ``` +* - {py:obj}`tune_acquisition_geometry ` + - ```{autodoc2-docstring} corrct.alignment.cone_beam.tune_acquisition_geometry + :summary: + ``` +```` + +### API + +````{py:function} _class_to_json(obj: object) -> str +:canonical: corrct.alignment.cone_beam._class_to_json + +```{autodoc2-docstring} corrct.alignment.cone_beam._class_to_json +``` +```` + +````{py:function} _get_rot_axis_angle_deg(center_1_vu: collections.abc.Sequence[float] | numpy.typing.NDArray, center_2_vu: collections.abc.Sequence[float] | numpy.typing.NDArray, decimals: int | None = 4, dtype: numpy.typing.DTypeLike = np.float32) -> float +:canonical: corrct.alignment.cone_beam._get_rot_axis_angle_deg + +```{autodoc2-docstring} corrct.alignment.cone_beam._get_rot_axis_angle_deg +``` +```` + +`````{py:class} ConeBeamGeometry +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry +``` + +````{py:attribute} theta_deg +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.theta_deg +:type: float +:value: > + 0.0 + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.theta_deg +``` + +```` + +````{py:attribute} phi_deg +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.phi_deg +:type: float +:value: > + 0.0 + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.phi_deg +``` + +```` + +````{py:attribute} eta_deg +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.eta_deg +:type: float +:value: > + 0.0 + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.eta_deg +``` + +```` + +````{py:attribute} D_pix +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.D_pix +:type: float +:value: > + 0.0 + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.D_pix +``` + +```` + +````{py:attribute} R_pix +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.R_pix +:type: float +:value: > + 0.0 + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.R_pix +``` + +```` + +````{py:attribute} v0_pix +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.v0_pix +:type: float +:value: > + 0.0 + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.v0_pix +``` + +```` + +````{py:attribute} u0_pix +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.u0_pix +:type: float +:value: > + 0.0 + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.u0_pix +``` + +```` + +````{py:attribute} det_size_v_pix +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.det_size_v_pix +:type: int +:value: > + 0 + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.det_size_v_pix +``` + +```` + +````{py:attribute} det_size_u_pix +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.det_size_u_pix +:type: int +:value: > + 0 + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.det_size_u_pix +``` + +```` + +````{py:attribute} pix_size_um +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.pix_size_um +:type: float +:value: > + 0.0 + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.pix_size_um +``` + +```` + +````{py:method} __str__() -> str +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.__str__ + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.__str__ +``` + +```` + +````{py:method} get_prj_geom(translate_z_to_center: bool = True) -> corrct.models.ProjectionGeometry +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.get_prj_geom + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.get_prj_geom +``` + +```` + +````{py:method} get_vol_geom(up_sampling: int = 1) -> corrct.models.VolumeGeometry +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.get_vol_geom + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.get_vol_geom +``` + +```` + +````{py:method} update(field: str, val: float, is_relative: bool = True, decimals: int | None = 3) -> corrct.alignment.cone_beam.ConeBeamGeometry +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.update + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.update +``` + +```` + +````{py:method} get_tuning_params(field: str, val_range: collections.abc.Sequence[float] | numpy.typing.NDArray, is_relative: bool = True) -> collections.abc.Sequence[corrct.alignment.cone_beam.ConeBeamGeometry] +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.get_tuning_params + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.get_tuning_params +``` + +```` + +````{py:method} to_json() -> str +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.to_json + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.to_json +``` + +```` + +````{py:method} from_json(data_json: str) -> None +:canonical: corrct.alignment.cone_beam.ConeBeamGeometry.from_json + +```{autodoc2-docstring} corrct.alignment.cone_beam.ConeBeamGeometry.from_json +``` + +```` + +````` + +`````{py:class} FitConeBeamGeometry(prj_size_vu: collections.abc.Sequence[int] | numpy.typing.NDArray, points_ell1: collections.abc.Sequence[collections.abc.Sequence[float]] | numpy.typing.NDArray, points_ell2: collections.abc.Sequence[collections.abc.Sequence[float]] | numpy.typing.NDArray, points_axis: collections.abc.Sequence[collections.abc.Sequence[float]] | numpy.typing.NDArray | None = None, pix_size_um: float | None = None, use_l1_norm: bool = False, verbose: bool = True, plot_result: bool = False) +:canonical: corrct.alignment.cone_beam.FitConeBeamGeometry + +```{autodoc2-docstring} corrct.alignment.cone_beam.FitConeBeamGeometry +``` + +```{rubric} Initialization +``` + +```{autodoc2-docstring} corrct.alignment.cone_beam.FitConeBeamGeometry.__init__ +``` + +````{py:attribute} acq_geom +:canonical: corrct.alignment.cone_beam.FitConeBeamGeometry.acq_geom +:type: corrct.alignment.cone_beam.ConeBeamGeometry +:value: > + None + +```{autodoc2-docstring} corrct.alignment.cone_beam.FitConeBeamGeometry.acq_geom +``` + +```` + +````{py:method} _initialize(use_l1_norm: bool) -> None +:canonical: corrct.alignment.cone_beam.FitConeBeamGeometry._initialize + +```{autodoc2-docstring} corrct.alignment.cone_beam.FitConeBeamGeometry._initialize +``` + +```` + +````{py:method} fit(r: float, e: float = 1, meas_D_pix: float | None = None) -> corrct.alignment.cone_beam.ConeBeamGeometry +:canonical: corrct.alignment.cone_beam.FitConeBeamGeometry.fit + +```{autodoc2-docstring} corrct.alignment.cone_beam.FitConeBeamGeometry.fit +``` + +```` + +````{py:method} _fit_distance_det2src(ellipse_1: corrct.alignment.fitting.Ellipse, ellipse_2: corrct.alignment.fitting.Ellipse, e: float = 1) -> float +:canonical: corrct.alignment.cone_beam.FitConeBeamGeometry._fit_distance_det2src +:staticmethod: + +```{autodoc2-docstring} corrct.alignment.cone_beam.FitConeBeamGeometry._fit_distance_det2src +``` + +```` + +````` + +````{py:function} tune_acquisition_geometry(acq_geom_init: corrct.alignment.cone_beam.ConeBeamGeometry, data: numpy.typing.NDArray, angles_rot_rad: collections.abc.Sequence[float] | numpy.typing.NDArray, params: dict[str, collections.abc.Sequence[float] | numpy.typing.NDArray], data_mask: numpy.typing.NDArray | None = None, verbose: bool = True) -> corrct.alignment.cone_beam.ConeBeamGeometry +:canonical: corrct.alignment.cone_beam.tune_acquisition_geometry + +```{autodoc2-docstring} corrct.alignment.cone_beam.tune_acquisition_geometry +``` +```` diff --git a/doc_sources/apidocs/corrct/corrct.alignment.fitting.md b/doc_sources/apidocs/corrct/corrct.alignment.fitting.md index 7471d49..de7ae21 100644 --- a/doc_sources/apidocs/corrct/corrct.alignment.fitting.md +++ b/doc_sources/apidocs/corrct/corrct.alignment.fitting.md @@ -9,6 +9,22 @@ ## Module Contents +### Classes + +````{list-table} +:class: autosummary longtable +:align: left + +* - {py:obj}`Trajectory ` + - ```{autodoc2-docstring} corrct.alignment.fitting.Trajectory + :summary: + ``` +* - {py:obj}`Ellipse ` + - ```{autodoc2-docstring} corrct.alignment.fitting.Ellipse + :summary: + ``` +```` + ### Functions ````{list-table} @@ -59,6 +75,22 @@ - ```{autodoc2-docstring} corrct.alignment.fitting.refine_max_position_2d :summary: ``` +* - {py:obj}`fit_parabola_min ` + - ```{autodoc2-docstring} corrct.alignment.fitting.fit_parabola_min + :summary: + ``` +* - {py:obj}`fit_ellipse_center ` + - ```{autodoc2-docstring} corrct.alignment.fitting.fit_ellipse_center + :summary: + ``` +* - {py:obj}`fit_ellipse_parameters ` + - ```{autodoc2-docstring} corrct.alignment.fitting.fit_ellipse_parameters + :summary: + ``` +* - {py:obj}`fit_ellipse ` + - ```{autodoc2-docstring} corrct.alignment.fitting.fit_ellipse + :summary: + ``` ```` ### Data @@ -175,3 +207,175 @@ ```{autodoc2-docstring} corrct.alignment.fitting.refine_max_position_2d ``` ```` + +````{py:function} fit_parabola_min(fun_x: numpy.typing.ArrayLike | numpy.typing.NDArray, fun_vals: numpy.typing.ArrayLike | numpy.typing.NDArray, scale: typing.Literal[linear, log] = 'linear', decimals: int = 2) -> tuple[float, float, tuple[numpy.typing.NDArray, numpy.typing.NDArray] | None] +:canonical: corrct.alignment.fitting.fit_parabola_min + +```{autodoc2-docstring} corrct.alignment.fitting.fit_parabola_min +``` +```` + +````{py:function} fit_ellipse_center(prj_points_vu: numpy.typing.NDArray, rescale: bool = True, use_l1_norm: bool = False, decimals: int | None = 2) -> numpy.typing.NDArray +:canonical: corrct.alignment.fitting.fit_ellipse_center + +```{autodoc2-docstring} corrct.alignment.fitting.fit_ellipse_center +``` +```` + +````{py:function} fit_ellipse_parameters(prj_points_vu: numpy.typing.NDArray, rescale: bool = True, use_l1_norm: bool = False) -> tuple[float, float, float, float, float] +:canonical: corrct.alignment.fitting.fit_ellipse_parameters + +```{autodoc2-docstring} corrct.alignment.fitting.fit_ellipse_parameters +``` +```` + +`````{py:class} Trajectory +:canonical: corrct.alignment.fitting.Trajectory + +Bases: {py:obj}`abc.ABC` + +```{autodoc2-docstring} corrct.alignment.fitting.Trajectory +``` + +````{py:method} __call__(uus: collections.abc.Sequence[float] | numpy.typing.NDArray) -> collections.abc.Sequence[numpy.typing.NDArray] +:canonical: corrct.alignment.fitting.Trajectory.__call__ +:abstractmethod: + +```{autodoc2-docstring} corrct.alignment.fitting.Trajectory.__call__ +``` + +```` + +````` + +`````{py:class} Ellipse(a: float, b: float, c: float, u: float, v: float, c_vu: corrct.alignment.fitting.NDArrayFloat) +:canonical: corrct.alignment.fitting.Ellipse + +Bases: {py:obj}`corrct.alignment.fitting.Trajectory` + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse +``` + +```{rubric} Initialization +``` + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.__init__ +``` + +````{py:attribute} a +:canonical: corrct.alignment.fitting.Ellipse.a +:type: float +:value: > + None + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.a +``` + +```` + +````{py:attribute} b +:canonical: corrct.alignment.fitting.Ellipse.b +:type: float +:value: > + None + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.b +``` + +```` + +````{py:attribute} c +:canonical: corrct.alignment.fitting.Ellipse.c +:type: float +:value: > + None + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.c +``` + +```` + +````{py:attribute} u +:canonical: corrct.alignment.fitting.Ellipse.u +:type: float +:value: > + None + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.u +``` + +```` + +````{py:attribute} v +:canonical: corrct.alignment.fitting.Ellipse.v +:type: float +:value: > + None + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.v +``` + +```` + +````{py:attribute} c_vu +:canonical: corrct.alignment.fitting.Ellipse.c_vu +:type: corrct.alignment.fitting.NDArrayFloat +:value: > + None + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.c_vu +``` + +```` + +````{py:method} __repr__() -> str +:canonical: corrct.alignment.fitting.Ellipse.__repr__ + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.__repr__ +``` + +```` + +````{py:property} extremes_u +:canonical: corrct.alignment.fitting.Ellipse.extremes_u +:type: tuple[float, float] + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.extremes_u +``` + +```` + +````{py:property} center_vu +:canonical: corrct.alignment.fitting.Ellipse.center_vu +:type: numpy.typing.NDArray + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.center_vu +``` + +```` + +````{py:property} parameters +:canonical: corrct.alignment.fitting.Ellipse.parameters +:type: tuple[float, float, float, float, float] + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.parameters +``` + +```` + +````{py:method} __call__(uus: numpy.typing.ArrayLike | numpy.typing.NDArray) -> collections.abc.Sequence[numpy.typing.NDArray] +:canonical: corrct.alignment.fitting.Ellipse.__call__ + +```{autodoc2-docstring} corrct.alignment.fitting.Ellipse.__call__ +``` + +```` + +````` + +````{py:function} fit_ellipse(prj_points_vu: numpy.typing.ArrayLike | numpy.typing.NDArray, rescale: bool = True, use_l1_norm: bool = False) -> corrct.alignment.fitting.Ellipse +:canonical: corrct.alignment.fitting.fit_ellipse + +```{autodoc2-docstring} corrct.alignment.fitting.fit_ellipse +``` +```` diff --git a/doc_sources/apidocs/corrct/corrct.alignment.markers.md b/doc_sources/apidocs/corrct/corrct.alignment.markers.md new file mode 100644 index 0000000..a388aa5 --- /dev/null +++ b/doc_sources/apidocs/corrct/corrct.alignment.markers.md @@ -0,0 +1,103 @@ +# {py:mod}`corrct.alignment.markers` + +```{py:module} corrct.alignment.markers +``` + +```{autodoc2-docstring} corrct.alignment.markers +:allowtitles: +``` + +## Module Contents + +### Classes + +````{list-table} +:class: autosummary longtable +:align: left + +* - {py:obj}`MarkerTrackingVisualizer ` + - ```{autodoc2-docstring} corrct.alignment.markers.MarkerTrackingVisualizer + :summary: + ``` +```` + +### Functions + +````{list-table} +:class: autosummary longtable +:align: left + +* - {py:obj}`cm2inch ` + - ```{autodoc2-docstring} corrct.alignment.markers.cm2inch + :summary: + ``` +* - {py:obj}`track_marker ` + - ```{autodoc2-docstring} corrct.alignment.markers.track_marker + :summary: + ``` +* - {py:obj}`create_marker_disk ` + - ```{autodoc2-docstring} corrct.alignment.markers.create_marker_disk + :summary: + ``` +```` + +### API + +````{py:function} cm2inch(dims: collections.abc.Sequence[float] | numpy.typing.NDArray) -> tuple[float] +:canonical: corrct.alignment.markers.cm2inch + +```{autodoc2-docstring} corrct.alignment.markers.cm2inch +``` +```` + +````{py:function} track_marker(prj_data: numpy.typing.NDArray, marker_vu: numpy.typing.NDArray, stack_axis: int = -2) -> numpy.typing.NDArray +:canonical: corrct.alignment.markers.track_marker + +```{autodoc2-docstring} corrct.alignment.markers.track_marker +``` +```` + +````{py:function} create_marker_disk(data_shape_vu: collections.abc.Sequence[int] | numpy.typing.NDArray, radius: float, super_sampling: int = 5, conv: bool = True) -> numpy.typing.NDArray +:canonical: corrct.alignment.markers.create_marker_disk + +```{autodoc2-docstring} corrct.alignment.markers.create_marker_disk +``` +```` + +`````{py:class} MarkerTrackingVisualizer(fitted_positions_vu: numpy.typing.NDArray, images: numpy.typing.NDArray, marker: numpy.typing.NDArray, trajectory: corrct.alignment.fitting.Trajectory | None = None) +:canonical: corrct.alignment.markers.MarkerTrackingVisualizer + +```{autodoc2-docstring} corrct.alignment.markers.MarkerTrackingVisualizer +``` + +```{rubric} Initialization +``` + +```{autodoc2-docstring} corrct.alignment.markers.MarkerTrackingVisualizer.__init__ +``` + +````{py:method} _update() -> None +:canonical: corrct.alignment.markers.MarkerTrackingVisualizer._update + +```{autodoc2-docstring} corrct.alignment.markers.MarkerTrackingVisualizer._update +``` + +```` + +````{py:method} _key_event(evnt) -> None +:canonical: corrct.alignment.markers.MarkerTrackingVisualizer._key_event + +```{autodoc2-docstring} corrct.alignment.markers.MarkerTrackingVisualizer._key_event +``` + +```` + +````{py:method} _scroll_event(evnt) -> None +:canonical: corrct.alignment.markers.MarkerTrackingVisualizer._scroll_event + +```{autodoc2-docstring} corrct.alignment.markers.MarkerTrackingVisualizer._scroll_event +``` + +```` + +````` diff --git a/doc_sources/apidocs/corrct/corrct.alignment.md b/doc_sources/apidocs/corrct/corrct.alignment.md index e6ac67b..9865b01 100644 --- a/doc_sources/apidocs/corrct/corrct.alignment.md +++ b/doc_sources/apidocs/corrct/corrct.alignment.md @@ -14,7 +14,9 @@ :maxdepth: 1 corrct.alignment.centering +corrct.alignment.cone_beam corrct.alignment.fitting +corrct.alignment.markers corrct.alignment.shifts ``` diff --git a/doc_sources/apidocs/corrct/corrct.models.md b/doc_sources/apidocs/corrct/corrct.models.md index d4e1a41..8e19a3c 100644 --- a/doc_sources/apidocs/corrct/corrct.models.md +++ b/doc_sources/apidocs/corrct/corrct.models.md @@ -63,6 +63,10 @@ - ```{autodoc2-docstring} corrct.models.get_vol_geom_from_volume :summary: ``` +* - {py:obj}`plot_projection_geometry ` + - ```{autodoc2-docstring} corrct.models.plot_projection_geometry + :summary: + ``` ```` ### Data @@ -485,3 +489,10 @@ Bases: {py:obj}`corrct.models.Geometry` ```{autodoc2-docstring} corrct.models.get_vol_geom_from_volume ``` ```` + +````{py:function} plot_projection_geometry(prj_geom: corrct.models.ProjectionGeometry, vol_geom: corrct.models.VolumeGeometry) +:canonical: corrct.models.plot_projection_geometry + +```{autodoc2-docstring} corrct.models.plot_projection_geometry +``` +```` diff --git a/doc_sources/cone_beam_calibration_tutorial.md b/doc_sources/cone_beam_calibration_tutorial.md new file mode 100644 index 0000000..924b785 --- /dev/null +++ b/doc_sources/cone_beam_calibration_tutorial.md @@ -0,0 +1,159 @@ +# Cone-Beam Geometry Calibration + +In this tutorial, we'll demonstrate how to use `corrct`'s cone-beam geometry calibration routines. +This process is crucial for accurate X-ray tomography reconstructions, especially when dealing with cone-beam geometry setups. +This calibration procedure is based on: + +- Noo, F., Clackdoyle, R., Mennessier, C., White, T. A. & Roney, T. J. (2000). Phys. Med. Biol. 45, 3489–3508. + doi: 10.1088/0031-9155/45/11/327 + +This procedure relies on certain assumptions of the geometry (e.g., non-degeneracy of certain parameters), and it has +limitations with respect to its accuracy or ability to determine certain parameters (i.e., one of the detector tilts). +More advanced procedures have been published in the literature more recently, which address its pitfalls. +However, it is quite straight-forward to be carried out, and it requires minimal data to work. +The remaining uncertainties can be fine tuned with the data-driven approach presented here below. + +In particular, this procedure assumes that you provide two different tomography scans of a rotating marker, +on circular orbits around the rotation axis of the sample coordinates. +Ideally, these two trajectories should be on opposite sides of the origin along the rotation axis (not necessarily at the exact same distance). +The projections of these two orbits will form two distinct ellipses over the detector. + +To follow this tutorial, you will need to download the demonstration dataset from [Zenodo](https://doi.org/10.5281/zenodo.20559974), and place it a sub-directory `./data/` of the current working directory. + +## Setup and Data Loading + +First, let's import the necessary libraries and load our example dataset: + +```python +from pathlib import Path +import h5py +import numpy as np +from numpy.typing import NDArray + +from corrct.alignment.cone_beam import FitConeBeamGeometry, tune_acquisition_geometry +from corrct.alignment.markers import create_marker_disk, track_marker +from corrct.models import plot_projection_geometry + +def _get_data(fid: h5py.File, data_path: str) -> NDArray: + dataset = fid[data_path] + if isinstance(dataset, h5py.Dataset): + return np.array(dataset[()]) + else: + raise ValueError(f"Path: {data_path}, is not a h5py.Dataset, but a {dataset} instead") + +def _load_data(fname: str | Path) -> dict[str, NDArray]: + with h5py.File(fname) as fid: + return {k: _get_data(fid, f"/{k}") for k in fid.keys()} + +try: + data = _load_data("./data/calibration_scans.h5") +except FileNotFoundError as exc: + raise ValueError("Please download the example dataset from https://doi.org/10.5281/zenodo.20559974") from exc +``` + +This dataset is composed of two different tomographic scans, with 60 angles each. In each scan, we record the motion of a +sphere in a circular orbit at radius of 3 millimeters around the rotation axis of the sample rotation stage. + +## Sphere Trajectory Tracking + +Before we can proceed with the geometry calibration, we need to identify the sphere position over the detector at each angle. +The sphere trajectory over the detector will be an ellipse, corresponding to the projection of each circular trajectory from the +sample space. + +We'll start by creating a marker disk: + +```python +prj_size_vu = (data["scan_1"].shape[0], data["scan_1"].shape[2]) + +probe = create_marker_disk(prj_size_vu, 3.5) +``` + +:::{note} + The `create_marker_disk` function is a convenience function for objects that look like projected solid spheres. + For other types of objects, you might want use one real image of the object (or an average of them). +::: + +We can then track its position in our projection images: +```python +pos_l, pos_u = (track_marker(imgs, probe) for imgs in (data["scan_1"], data["scan_2"])) +``` + +## Initial Geometry Estimation + +Next, we'll estimate the initial cone-beam geometry parameters: + +```python +pixel_size_um = float(data["pixel_size_um"]) +orbit_radius_pix = float(data["orbit_radius_um"]) / pixel_size_um + +fit_cb_geom = FitConeBeamGeometry( + prj_size_vu, points_ell1=pos_u, points_ell2=pos_l, pix_size_um=pixel_size_um, plot_result=True +) +acq_geom = fit_cb_geom.fit(r=orbit_radius_pix) + +print(acq_geom) +``` + +The radius of the circular trajectory should be known and provided in pixels. + +:::{note} + While the attribute `pix_size_um` is optional, we suggest passing it, as it will, as it will enable the fitting routines + to automatically present the fitted distances in both pixels and micrometers. +::: + +If the switch `plot_result` is `True`, this procedure will display the trajectory of the ellipses, both before `eta` fitting and after. + +![eta-det](images/alignment_cone-beam_ellipse-eta-fit.png) + +## Geometry Refinement + +For more accurate results, it is possible to refine our geometry estimation, by minimizing the residual of the reconstructions: + +```python +imgs_t = (data["scan_1"] + data["scan_2"]).astype(np.float32) +angles_rot_rad = np.deg2rad(data["angles_deg"]) + +acq_geom = tune_acquisition_geometry( + acq_geom, + data=imgs_t, + angles_rot_rad=angles_rot_rad, + params=dict( + theta_deg=np.linspace(-1, 1, 5), + phi_deg=np.linspace(-0.25, 0.25, 5), + u0_pix=np.linspace(-1, 1, 9), + v0_pix=np.linspace(-1, 1, 9), + ), + verbose=True, +) + +print(acq_geom) +``` + +The `tune_acquisition_geometry` expects an initial geometry, projection data (with corresponding rotation angles in radians), and a dictionary listing for each parameter to be optimized the list of values to try. +It will then return the geometry that has the lowest reconstruction residual for the given values. + +::::{grid} +:::{grid-item} +![Phi (detector tilt)](images/alignment_cone-beam_fit-phi.png) +::: +:::{grid-item} +![Theta (detector tilt)](images/alignment_cone-beam_fit-theta.png) +::: +:::: +::::{grid} +:::{grid-item} +![u0 (Horizontal source position)](images/alignment_cone-beam_fit-u0.png) +::: +:::{grid-item} +![v0 (Vertical source position)](images/alignment_cone-beam_fit-v0.png) +::: +:::: + +## Visualization + +After completing the calibration, you can visualize the projection geometry: + +```python +plot_projection_geometry(acq_geom.get_prj_geom(), acq_geom.get_vol_geom()) +``` +![plot-geom](images/alignment_cone-beam_plot-geom.png) diff --git a/doc_sources/images/alignment_cone-beam_ellipse-eta-fit.png b/doc_sources/images/alignment_cone-beam_ellipse-eta-fit.png new file mode 100644 index 0000000..afd3a10 Binary files /dev/null and b/doc_sources/images/alignment_cone-beam_ellipse-eta-fit.png differ diff --git a/doc_sources/images/alignment_cone-beam_fit-phi.png b/doc_sources/images/alignment_cone-beam_fit-phi.png new file mode 100644 index 0000000..3086a9e Binary files /dev/null and b/doc_sources/images/alignment_cone-beam_fit-phi.png differ diff --git a/doc_sources/images/alignment_cone-beam_fit-theta.png b/doc_sources/images/alignment_cone-beam_fit-theta.png new file mode 100644 index 0000000..8741c0b Binary files /dev/null and b/doc_sources/images/alignment_cone-beam_fit-theta.png differ diff --git a/doc_sources/images/alignment_cone-beam_fit-u0.png b/doc_sources/images/alignment_cone-beam_fit-u0.png new file mode 100644 index 0000000..a2587c5 Binary files /dev/null and b/doc_sources/images/alignment_cone-beam_fit-u0.png differ diff --git a/doc_sources/images/alignment_cone-beam_fit-v0.png b/doc_sources/images/alignment_cone-beam_fit-v0.png new file mode 100644 index 0000000..94572a4 Binary files /dev/null and b/doc_sources/images/alignment_cone-beam_fit-v0.png differ diff --git a/doc_sources/images/alignment_cone-beam_plot-geom.png b/doc_sources/images/alignment_cone-beam_plot-geom.png new file mode 100644 index 0000000..d32d2f8 Binary files /dev/null and b/doc_sources/images/alignment_cone-beam_plot-geom.png differ diff --git a/doc_sources/index.rst b/doc_sources/index.rst index 1087c0e..66226a9 100644 --- a/doc_sources/index.rst +++ b/doc_sources/index.rst @@ -17,6 +17,7 @@ Contents: attenuation_tutorial param_tuning alignment_tools + cone_beam_calibration_tutorial apidocs/index changelog diff --git a/examples/alignment/example_02_cone-beam_calibration.py b/examples/alignment/example_02_cone-beam_calibration.py new file mode 100644 index 0000000..86fefc9 --- /dev/null +++ b/examples/alignment/example_02_cone-beam_calibration.py @@ -0,0 +1,68 @@ +""" +Example demonstrating the use of cone-beam geometry calibration routines. + +@author: Nicola VIGANÒ, CEA-IRIG, Grenoble, France +""" + +from pathlib import Path + +import h5py +import numpy as np +from numpy.typing import NDArray + +from corrct.alignment.cone_beam import FitConeBeamGeometry, tune_acquisition_geometry +from corrct.alignment.markers import create_marker_disk, track_marker +from corrct.models import plot_projection_geometry + + +def _get_data(fid: h5py.File, data_path: str) -> NDArray: + dataset = fid[data_path] + if isinstance(dataset, h5py.Dataset): + return np.array(dataset[()]) + else: + raise ValueError(f"Path: {data_path}, is not a h5py.Dataset, but a {dataset} instead") + + +def _load_data(fname: str | Path) -> dict[str, NDArray]: + with h5py.File(fname) as fid: + return {k: _get_data(fid, f"/{k}") for k in fid.keys()} + + +try: + data = _load_data("./data/calibration_scans.h5") +except FileNotFoundError as exc: + raise ValueError("Please download the example dataset from https://doi.org/10.5281/zenodo.20559974") from exc + +prj_size_vu = (data["scan_1"].shape[0], data["scan_1"].shape[2]) + +probe = create_marker_disk(prj_size_vu, 3.5) +pos_l, pos_u = (track_marker(imgs, probe) for imgs in (data["scan_1"], data["scan_2"])) + +pixel_size_um = float(data["pixel_size_um"]) +orbit_radius_pix = float(data["orbit_radius_um"]) / pixel_size_um + +fit_cb_geom = FitConeBeamGeometry( + prj_size_vu, points_ell1=pos_u, points_ell2=pos_l, pix_size_um=pixel_size_um, plot_result=True +) +acq_geom = fit_cb_geom.fit(r=orbit_radius_pix) + +print(acq_geom) + +imgs_t = (data["scan_1"] + data["scan_2"]).astype(np.float32) +angles_rot_rad = np.deg2rad(data["angles_deg"]) + +acq_geom = tune_acquisition_geometry( + acq_geom, + data=imgs_t, + angles_rot_rad=angles_rot_rad, + params=dict( + theta_deg=np.linspace(-1, 1, 5), + phi_deg=np.linspace(-0.25, 0.25, 5), + u0_pix=np.linspace(-1, 1, 9), + v0_pix=np.linspace(-1, 1, 9), + ), + verbose=True, +) + +print(acq_geom) +plot_projection_geometry(acq_geom.get_prj_geom(), acq_geom.get_vol_geom()) diff --git a/src/corrct/alignment/__init__.py b/src/corrct/alignment/__init__.py index 0b1f86f..116b79b 100644 --- a/src/corrct/alignment/__init__.py +++ b/src/corrct/alignment/__init__.py @@ -6,6 +6,8 @@ from . import centering # noqa: F401, F402 from . import fitting # noqa: F401, F402 from . import shifts # noqa: F401, F402 +from . import cone_beam # noqa: F401, F402 +from . import markers # noqa: F401, F402 RecenterVolume = centering.RecenterVolume DetectorShiftsPRE = shifts.DetectorShiftsPRE diff --git a/src/corrct/alignment/cone_beam.py b/src/corrct/alignment/cone_beam.py new file mode 100644 index 0000000..7ebbded --- /dev/null +++ b/src/corrct/alignment/cone_beam.py @@ -0,0 +1,643 @@ +#!/usr/bin/env python3 +""" +Calibrate cone-beam reconstruction geometry. + +@author: Nicola VIGANÒ, ESRF - The European Synchrotron, Grenoble, France, +and CEA-IRIG, Grenoble, France +""" + +import json +from dataclasses import dataclass, replace as dc_replace +from collections.abc import Sequence + +import matplotlib.pyplot as plt +import numpy as np +from numpy.typing import DTypeLike, NDArray +from scipy.spatial.transform import Rotation +from tqdm.auto import tqdm + +from corrct.models import ProjectionGeometry, VolumeGeometry +from corrct.projectors import ProjectorUncorrected +from corrct.solvers import SIRT +from corrct.alignment.fitting import Ellipse, fit_parabola_min, fit_ellipse, fit_ellipse_center + + +def _class_to_json(obj: object) -> str: + return json.dumps(obj, default=lambda o: {o.__class__.__name__: o.__dict__}, sort_keys=True, indent=4) + + +def _get_rot_axis_angle_deg( + center_1_vu: Sequence[float] | NDArray, + center_2_vu: Sequence[float] | NDArray, + decimals: int | None = 4, + dtype: DTypeLike = np.float32, +) -> float: + center_1_vu = np.squeeze(np.array(center_1_vu, dtype=dtype)) + center_2_vu = np.squeeze(np.array(center_2_vu, dtype=dtype)) + + # Check if the arrays are 2D + if center_1_vu.ndim != 1 or center_2_vu.ndim != 1: + raise ValueError("Input arrays must be 1D") + + # Check if the first dimension has length 2 + if center_1_vu.shape[0] != 2 or center_2_vu.shape[0] != 2: + raise ValueError("Input arrays must have length 2") + + diffs_vu = center_1_vu - center_2_vu + angle_rad = np.arctan2(diffs_vu[-1], diffs_vu[-2]) + angle_rad = np.mod(angle_rad, 2 * np.pi) + + angle_deg = np.rad2deg(angle_rad - np.pi) + + if decimals is not None: + angle_deg = np.around(angle_deg, decimals=decimals) + + return float(angle_deg) + + +@dataclass +class ConeBeamGeometry: + """Store the acquisition geometry parameters, used for creating reconstruction geometries. + + A description of the geometry / meaning of the fields can be found here: + - Noo, F., Clackdoyle, R., Mennessier, C., White, T. A. & Roney, T. J. (2000). Phys. Med. Biol. 45, 3489–3508. + doi: 10.1088/0031-9155/45/11/327 + """ + + theta_deg: float = 0.0 + phi_deg: float = 0.0 + eta_deg: float = 0.0 + D_pix: float = 0.0 + R_pix: float = 0.0 + v0_pix: float = 0.0 + u0_pix: float = 0.0 + det_size_v_pix: int = 0 + det_size_u_pix: int = 0 + pix_size_um: float = 0.0 + + def __str__(self) -> str: + """ + Return a human readable representation of the object. + + Returns + ------- + str + The human readable representation of the object. + """ + descr = "AcquisitionGeometry(\n" + for field, value in self.__dict__.items(): + descr += f" {field} = {value}" + if field.lower()[-3:] == "deg": + descr += " [deg]" + elif field.lower()[-3:] == "pix": + descr += " [pix]" + if self.pix_size_um > 0.0: + descr += f" ({value * self.pix_size_um} [um])" + elif field.lower()[-2:] == "um": + descr += " [um]" + descr += ",\n" + return descr + ")" + + def get_prj_geom(self, translate_z_to_center: bool = True) -> ProjectionGeometry: + """ + Create the geometry for reconstruction. + + Returns + ------- + Dict + The geometry to be used for reconstruction. + """ + # Sample base vectors + e_x_xyz = np.array([1, 0, 0]) + + # NOTE: Y coordinate is flipped, with respect to the input!!! + + theta_rad = np.deg2rad(self.theta_deg) + phi_rad = np.deg2rad(self.phi_deg) + + # Rotated detector base vectors + alpha_xyz = np.array([-np.sin(phi_rad), -np.cos(phi_rad), 0]) + beta_xyz = np.array([-np.sin(theta_rad) * np.cos(phi_rad), np.sin(theta_rad) * np.sin(phi_rad), np.cos(theta_rad)]) + + # Detector normal + e_n_xyz = np.array([np.cos(theta_rad) * np.cos(phi_rad), -np.cos(theta_rad) * np.sin(phi_rad), np.sin(theta_rad)]) + + det_pos_xyz = -e_n_xyz * self.D_pix + e_x_xyz * self.R_pix - alpha_xyz * self.u0_pix - beta_xyz * self.v0_pix + src_pos_xyz = e_x_xyz * self.R_pix + + pix2vox_ratio = self.R_pix / self.D_pix * np.abs(np.dot(e_n_xyz, e_x_xyz)) + + if translate_z_to_center: + det_center_xyz = e_n_xyz * self.D_pix + alpha_xyz * self.u0_pix + beta_xyz * self.v0_pix + + translation_z = det_center_xyz[2] / np.abs(det_center_xyz[0]) * self.R_pix + src_pos_xyz[2] += translation_z + det_pos_xyz[2] += translation_z + + rotation = Rotation.from_rotvec(-e_n_xyz * np.deg2rad(self.eta_deg)) + e_u_xyz = rotation.apply(alpha_xyz) + e_v_xyz = rotation.apply(beta_xyz) + + return ProjectionGeometry( + geom_type="cone", + src_pos_xyz=src_pos_xyz, + det_pos_xyz=det_pos_xyz, + det_u_xyz=e_u_xyz, + det_v_xyz=e_v_xyz, + rot_dir_xyz=np.array([0, 0, 1]), + pix2vox_ratio=pix2vox_ratio, + det_shape_vu=np.array((self.det_size_v_pix, self.det_size_u_pix)), + ) + + def get_vol_geom(self, up_sampling: int = 1) -> VolumeGeometry: + """ + Generate volume geometry. + + Returns + ------- + VolumeGeometry + The volume geometry. + """ + return VolumeGeometry( + _vol_shape_xyz=np.array([self.det_size_u_pix, self.det_size_u_pix, self.det_size_v_pix], dtype=int) * up_sampling, + vox_size=1 / up_sampling, + ) + + def update(self, field: str, val: float, is_relative: bool = True, decimals: int | None = 3) -> "ConeBeamGeometry": + """ + Return a copy of the original data, with a replaced field. + + Parameters + ---------- + field : str + The field to replace. + val : float + The new value of the field. + is_relative : bool, optional + Whether the value is relative to the previous. The default is True. + decimals : int | None, optional + The number of decimals (precision) to use for the updated values, by default 3 decimals. + + Returns + ------- + AcquisitionGeometry + The updated geometry. + """ + new_val = getattr(self, field) + val if is_relative else val + if decimals is not None: + new_val = np.around(new_val, decimals=decimals) + return dc_replace(self, **{field: new_val}) + + def get_tuning_params( + self, field: str, val_range: Sequence[float] | NDArray, is_relative: bool = True + ) -> Sequence["ConeBeamGeometry"]: + """ + Generate sequences of acquisition geometries, with a slight variation over a field's value. + + Parameters + ---------- + field : str + The field to tune. + val_range : Sequence[float] | NDArray + The value range. + is_relative : bool, optional + Whether the values are relative. The default is True. + + Returns + ------- + Sequence[AcquisitionGeometry] + The list of new acquisition geometries. + """ + return [self.update(field, val, is_relative) for val in np.array(val_range, ndmin=1)] + + def to_json(self) -> str: + """ + Save instance to JSON. + + Returns + ------- + str + The JSON representation. + """ + return _class_to_json(self) + + def from_json(self, data_json: str) -> None: + """ + Load instance from JSON. + + Parameters + ---------- + data : str + The JSON data to load. + + Raises + ------ + ValueError + In case we were to load more than one instance, or different classes. + """ + data_tree = json.loads(data_json) + if len(data_tree.keys()) > 1: + raise ValueError("Initialization from JSON: More than one class instance passed.") + + class_name = list(data_tree.keys())[0] + if list(data_tree.keys())[0] != self.__class__.__name__: + raise ValueError( + f"Initialization from JSON: expecting {self.__class__.__name__} class instance," + f" but {data_tree.keys()[0]} passed." + ) + + data_dict = data_tree[class_name] + for key in self.__dict__.keys(): + self.__dict__[key] = data_dict[key] + + +class FitConeBeamGeometry: + """Cone-beam geometry calibration object. + + This method is based on the following article: + - Noo, F., Clackdoyle, R., Mennessier, C., White, T. A. & Roney, T. J. (2000). Phys. Med. Biol. 45, 3489–3508. + doi: 10.1088/0031-9155/45/11/327 + """ + + acq_geom: ConeBeamGeometry + + def __init__( + self, + prj_size_vu: Sequence[int] | NDArray, + points_ell1: Sequence[Sequence[float]] | NDArray, + points_ell2: Sequence[Sequence[float]] | NDArray, + points_axis: Sequence[Sequence[float]] | NDArray | None = None, + pix_size_um: float | None = None, + use_l1_norm: bool = False, + verbose: bool = True, + plot_result: bool = False, + ): + """Initialize a cone-beam geometry calibration object. + + Parameters + ---------- + prj_size_vu : Sequence[int] | NDArray + Size of the projections. + points_ell1 : Sequence[Sequence[float]] | NDArray + Points of first ellipse. + points_ell2 : Sequence[Sequence[float]] | NDArray + Points of second ellipse. + points_axis : Sequence[Sequence[float]] | NDArray | None, optional + Points of the rotation axis, by default None + pix_size_um : float | None, optional + The size of the pixel edge in micrometers. Default is None. + use_l1_norm : bool, optional + Whether to use the l1-norm or the least-squares (l2-norm) fit for optimization. + Default is False. + verbose : bool, optional + Whether to produce verbose output, by default True + plot_result : bool, optional + Whether to plot the results of the geometry, by default False + It requires verbose to be True. + """ + prj_size_vu = np.array(prj_size_vu) + if prj_size_vu.ndim != 1 or len(prj_size_vu) != 2: + if len(prj_size_vu) == 1: + prj_size_vu = np.tile(prj_size_vu, 2) + else: + raise ValueError("prj_size_vu must be a 1D array with 2 elements") + self.prj_size_vu = prj_size_vu + + self.center_vu = np.squeeze(self.prj_size_vu)[:, None] / 2 + self.prj_origin_vu = None + + points_ell1 = np.array(points_ell1) + if points_ell1.ndim != 2: + raise ValueError("points_ell1 must be a 2D array") + if points_ell1.shape[0] != 2: + raise ValueError("points_ell1 must have a first dimension equal to 2") + self.points_ell1 = points_ell1 - self.center_vu + + points_ell2 = np.array(points_ell2) + if points_ell2.ndim != 2: + raise ValueError("points_ell2 must be a 2D array") + if points_ell2.shape[0] != 2: + raise ValueError("points_ell2 must have a first dimension equal to 2") + self.points_ell2 = points_ell2 - self.center_vu + + if points_axis is not None: + points_axis = np.array(points_axis) - self.center_vu + self.points_axis = points_axis + + self.acq_geom = ConeBeamGeometry(det_size_v_pix=int(self.prj_size_vu[0]), det_size_u_pix=int(self.prj_size_vu[1])) + if pix_size_um is not None: + self.acq_geom.pix_size_um = pix_size_um + + self.verbose = verbose + self.plot_result = plot_result and verbose + + self._initialize(use_l1_norm=use_l1_norm) + + def _initialize(self, use_l1_norm: bool) -> None: + ell1_fit_prj_c_vu = fit_ellipse_center(self.points_ell1, use_l1_norm=use_l1_norm) + ell2_fit_prj_c_vu = fit_ellipse_center(self.points_ell2, use_l1_norm=use_l1_norm) + fit_eta_deg = _get_rot_axis_angle_deg(ell1_fit_prj_c_vu, ell2_fit_prj_c_vu) + + if self.verbose: + print("Fitted / measured values:") + + if self.points_axis is not None: + # Using measured projected center, whenever available + ell1_acq_prj_c_vu = self.points_axis[:, 0] + ell2_acq_prj_c_vu = self.points_axis[:, 2] + acq_eta_deg = _get_rot_axis_angle_deg(ell1_acq_prj_c_vu, ell2_acq_prj_c_vu) + + if self.verbose: + print(f"- Ellipse 1 center: fitted = {ell1_fit_prj_c_vu} vs acquired = {ell1_acq_prj_c_vu} [pix]") + print(f"- Ellipse 2 center: fitted = {ell2_fit_prj_c_vu} vs acquired = {ell2_acq_prj_c_vu} [pix]") + print("- Detector tilt around its normal (eta):") + print(f" * fitted = {fit_eta_deg:.4} vs acq = {acq_eta_deg:.4} [deg] <= Using acquired!") + + self.ell1_prj_c_vu = ell1_acq_prj_c_vu + self.ell2_prj_c_vu = ell2_acq_prj_c_vu + + self.prj_origin_vu = self.points_axis[:, 1] + self.acq_geom.eta_deg = acq_eta_deg + else: + if self.verbose: + print(f"- Ellipse 1 center: fitted = {ell1_fit_prj_c_vu} [pix]") + print(f"- Ellipse 2 center: fitted = {ell2_fit_prj_c_vu} [pix]") + + self.ell1_prj_c_vu = ell1_fit_prj_c_vu + self.ell2_prj_c_vu = ell2_fit_prj_c_vu + + self.prj_origin_vu = None + self.acq_geom.eta_deg = fit_eta_deg + if self.verbose: + print(f"- Detector tilt around its normal (eta), fitted: {self.acq_geom.eta_deg:.4} [deg]") + + if np.abs(self.acq_geom.eta_deg) > 120: + raise ValueError( + "The order of the ellipses seems to have been inverted." + f" (it suggests an eta of {self.acq_geom.eta_deg}). Please swap them." + ) + + pix_size_um = self.acq_geom.pix_size_um + if self.verbose and self.prj_origin_vu is not None: + print( + f"- Projected origin on the detector: {self.prj_origin_vu} [pix]", + f"({self.prj_origin_vu * pix_size_um} [um])" if pix_size_um > 0.0 else "", + ) + + if np.abs(self.acq_geom.eta_deg) > 0.1: + rot = Rotation.from_rotvec(-np.deg2rad(self.acq_geom.eta_deg) * np.array([0, 0, 1])) + rot_mat = rot.as_matrix()[:2, :2] + + points_ell1_rot = rot_mat.dot(self.points_ell1) + points_ell2_rot = rot_mat.dot(self.points_ell2) + else: + points_ell1_rot = self.points_ell1.copy() + points_ell2_rot = self.points_ell2.copy() + + # Re-instantiate ellipse class, after rotation + self.ell1_rot: Ellipse = fit_ellipse(points_ell1_rot, use_l1_norm=use_l1_norm) + self.ell2_rot: Ellipse = fit_ellipse(points_ell2_rot, use_l1_norm=use_l1_norm) + + if self.plot_result: + fig, axs = plt.subplots() + axs.plot(self.points_ell1[1, :], self.points_ell1[0, :], "C0--", label="Ellipse 1 - Acquired") + axs.plot(self.points_ell2[1, :], self.points_ell2[0, :], "C1--", label="Ellipse 2 - Acquired") + axs.plot(points_ell1_rot[1, :], points_ell1_rot[0, :], "C0", label="Ellipse 1 - Rotated") + axs.plot(points_ell2_rot[1, :], points_ell2_rot[0, :], "C1", label="Ellipse 2 - Rotated") + axs.plot([ell1_fit_prj_c_vu[1], ell2_fit_prj_c_vu[1]], [ell1_fit_prj_c_vu[0], ell2_fit_prj_c_vu[0]], "C2--") + axs.plot([self.ell1_rot.u, self.ell2_rot.u], [self.ell1_rot.v, self.ell2_rot.v], "C2") + if self.points_axis is not None: + axs.scatter(self.points_axis[1], self.points_axis[0], c="C2", marker="*", label="Centers - Acquired") + axs.legend() + axs.grid() + fig.tight_layout() + plt.show(block=False) + + def fit(self, r: float, e: float = 1, meas_D_pix: float | None = None) -> ConeBeamGeometry: + """ + Fit the cone-beam geometry parameters, that will be used for producing the projection geometry. + + Parameters + ---------- + r : float + The radius of the circle performed by the spheres in pixels. + e : float, optional + Either 1 or -1, indicating whether the source is between the circles or not. The default is 1. + meas_D_pix : float, optional + The measured source-detector distance in pixels. This parameter is only necessary when the + computed source-detector distance is invalid or zero. + + Raises + ------ + ValueError + In case of flipped ellipses or invalid computed source-detector distance. + """ + + def get_v0(v: float, b: float, a: float, c: float, D: float, sign_z: float) -> float: + return v - sign_z * np.sqrt(a + (a**2) * (D**2)) / np.sqrt(a * b - (c**2)) + + def get_denom(b: float, a: float, c: float, D: float) -> float: + return np.sqrt(a * b + (a**2) * b * (D**2) - (c**2)) + + def get_rho(b: float, a: float, c: float, D: float) -> float: + return np.sqrt(a * b - (c**2)) / get_denom(b, a, c, D) + + def get_zeta(b: float, a: float, c: float, D: float, sign_zk: float) -> float: + return D * sign_zk * a * np.sqrt(a) / get_denom(b, a, c, D) + + b1, a1, c1, v1, u1 = self.ell1_rot.parameters + b2, a2, c2, v2, u2 = self.ell2_rot.parameters + + if self.verbose: + print("Fitted values from the calibration scan parameters:") + print("- Ellipses' parameters:") + print(f" * upper ({a1 = :.6}, {b1 = :.6}, {c1 = :.6}, {v1 = :.6}, {u1 = :.6})") + print(f" * lower ({a2 = :.6}, {b2 = :.6}, {c2 = :.6}, {v2 = :.6}, {u2 = :.6})") + + pix_size_um = self.acq_geom.pix_size_um + + comp_D_pix = self._fit_distance_det2src(self.ell1_rot, self.ell2_rot, e=e) + if meas_D_pix is None: + if np.isnan(comp_D_pix) or np.isclose(comp_D_pix, 0.0): + raise ValueError( + f"The computed source-detector distance is invalid ({comp_D_pix}), please enter a measured value" + ) + + self.acq_geom.D_pix = comp_D_pix + + if self.verbose: + print( + f"- Computed source-detector distance: {self.acq_geom.D_pix:.6} [pix]", + f"({self.acq_geom.D_pix * pix_size_um:.4e} [um])" if pix_size_um > 0.0 else "", + ) + else: + self.acq_geom.D_pix = meas_D_pix + + if self.verbose: + print("- Source-detector distance (using measured):") + print(f" * Measured = {meas_D_pix:.6} vs computed = {comp_D_pix:.6} [pix]") + print(f" * Measured = {meas_D_pix*pix_size_um:.6} vs computed = {comp_D_pix*pix_size_um:.6} [um]") + + sign_z1 = -1 + sign_z2 = sign_z1 * -e + + v01 = get_v0(v1, b1, a1, c1, self.acq_geom.D_pix, sign_z1) + v02 = get_v0(v2, b2, a2, c2, self.acq_geom.D_pix, sign_z2) + + tilt_ratio1 = c1 / (2 * a1) + tilt_ratio2 = c2 / (2 * a2) + + self.acq_geom.v0_pix = np.array([v01, v02]).mean() + self.acq_geom.u0_pix = ( + (u1 + u2) / 2 + tilt_ratio1 * (v1 - self.acq_geom.v0_pix) + tilt_ratio2 * (v2 - self.acq_geom.v0_pix) + ) + + if self.verbose: + print(f"- Source position over detector: v0 = {self.acq_geom.v0_pix:.6}, u0 = {self.acq_geom.u0_pix:.6}") + print(f" * Separately fitted v (from the two ellipses): {v01 = :.6}, {v02 = :.6}") + + if np.linalg.norm(v01 - v02) > np.linalg.norm(v1 - v2): + raise ValueError( + f"Obtained: {v01 = }, {v02 = }, while {v1 = }, {v2 = }. Probably wrong order of ellipses (please flip them!)" + ) + + rho1 = get_rho(b1, a1, c1, self.acq_geom.D_pix) + rho2 = get_rho(b2, a2, c2, self.acq_geom.D_pix) + + zeta1 = get_zeta(b1, a1, c1, self.acq_geom.D_pix, sign_z1) + zeta2 = get_zeta(b2, a2, c2, self.acq_geom.D_pix, sign_z2) + + sin_phi1 = -tilt_ratio1 * zeta1 + sin_phi2 = -tilt_ratio2 * zeta2 + self.acq_geom.phi_deg = np.rad2deg(np.arcsin(sin_phi1 + sin_phi2)) + + R_e1 = r / rho1 + R_e2 = r / rho2 + + z1 = R_e1 * zeta1 + z2 = R_e2 * zeta2 + + z_full = z1 - z2 + + self.acq_geom.R_pix = (-z2 * R_e1 + z1 * R_e2) / z_full + + if np.isnan(self.acq_geom.R_pix) or np.isclose(self.acq_geom.R_pix, 0.0): + raise ValueError(f"The computed source-origin distance is invalid ({self.acq_geom.R_pix})") + + if self.prj_origin_vu is None: + self.prj_origin_vu = (-z2 * self.ell1_prj_c_vu + z1 * self.ell2_prj_c_vu) / z_full + if self.verbose: + print(f"- Projected origin on the detector: {self.prj_origin_vu} [pix]") + + self.acq_geom.theta_deg = 0.0 + # self.acq_geom.theta_rad = np.arcsin((R1 - R2) / z_full) / np.mean(self.prj_size_vu) + + if self.verbose: + print("- Distances between source and rotation axis:") + print(f" * R1 = {R_e1:.6}, R = {self.acq_geom.R_pix:.6}, R2 = {R_e2:.6} [pix]") + if pix_size_um > 0.0: + print( + f" * R1 = {R_e1 * pix_size_um:.4e},", + f"R = {self.acq_geom.R_pix * pix_size_um:.4e},", + f"R2 = {R_e2 * pix_size_um:.4e} [um]", + ) + print("- Heights of the two ellipses, with respect to the source:") + print(f" * {z1 = :.6}, {z2 = :.6} [pix]") + if pix_size_um > 0.0: + print(f" * z1 = {z1 * pix_size_um:.4e}, z2 = {z2 * pix_size_um:.4e} [um]") + print(f"- Polar angle of the detector (phi deg): {self.acq_geom.phi_deg}") + print(f"- Azimuthal angle of the detector (theta deg): {self.acq_geom.theta_deg}") + + return self.acq_geom + + @staticmethod + def _fit_distance_det2src(ellipse_1: Ellipse, ellipse_2: Ellipse, e: float = 1) -> float: + b1, a1, c1, v1, _ = ellipse_1.parameters + b2, a2, c2, v2, _ = ellipse_2.parameters + + ecc2 = np.sqrt(b2 - c2**2 / a2) + ecc1 = np.sqrt(b1 - c1**2 / a1) + + m0 = (v2 - v1) * ecc2 + m1 = ecc2 / ecc1 + + m0m1 = 2 * m0 * m1 + m1_2 = m1**2 + + n0 = (1 - m0**2 - m1_2) / m0m1 + n1 = (a2 - a1 * m1_2) / m0m1 + + n0n1 = n0 * n1 + n1_2 = n1**2 + + s2d_2 = ((a1 - 2 * n0n1) - e * np.sqrt(a1**2 + 4 * n1_2 - 4 * n0n1 * a1)) / (2 * n1_2) + return np.sqrt(s2d_2) + + +def tune_acquisition_geometry( + acq_geom_init: ConeBeamGeometry, + data: NDArray, + angles_rot_rad: Sequence[float] | NDArray, + params: dict[str, Sequence[float] | NDArray], + data_mask: NDArray | None = None, + verbose: bool = True, +) -> ConeBeamGeometry: + """ + Tune the acquisition geometry, based on calibration data self-consistency. + + Parameters + ---------- + acq_geom : ConeBeamGeometry + The cone-beam geometry to refine. + data : NDArray + The calibration projection data. + angles : Sequence[float] | NDArray + Angles of the projections. + params : dict[str, Sequence[float] | NDArray] + Parameters to tune as a dictionary. + The acquisition parameters to tune are the keys, and their test values are the dictionary values. + data_mask : NDArray | None, optional + Pixel mask of the data, to mask out dead or hot pixels. The default is None. + verbose : bool, optional + Whether to output verbose information or not, by default False. + """ + data = np.array(data) + angles_rot_rad = np.array(angles_rot_rad) + if data_mask is not None: + data_mask = np.array(data_mask) + + solver = SIRT(tolerance=0.0) + + acq_geom_tuned = acq_geom_init + + for par_name, par_vals in params.items(): + par_vals = np.array(par_vals) + desc = f"Tuning '{par_name}': " + residuals = np.atleast_1d(np.empty(len(par_vals))) + for par_ind, acq_geom in enumerate(tqdm(acq_geom_tuned.get_tuning_params(par_name, par_vals), desc=desc)): + vol_geom = acq_geom.get_vol_geom() + prj_geom = acq_geom.get_prj_geom() + with ProjectorUncorrected(vol_geom, angles_rot_rad, prj_geom=prj_geom) as prj: + _, info = solver(prj, data, iterations=100, b_mask=data_mask) + residuals[par_ind] = info.residuals[-1] + + min_par, min_res, fit_info = fit_parabola_min(par_vals, residuals, decimals=6) + + if verbose: + old_par_val = getattr(acq_geom_tuned, par_name) + print(f"Min of {par_name}: {old_par_val}" + f" -> {old_par_val + min_par} (diff: {min_par})\n") + fig, axs = plt.subplots() + axs.plot(par_vals, residuals) + if fit_info is not None: + x = np.linspace(fit_info[1][0], fit_info[1][2]) + y = fit_info[0][0] + x * (fit_info[0][1] + x * fit_info[0][2]) + axs.plot(x, y) + axs.scatter(min_par, min_res, 10, "r") + axs.set_title(par_name) + axs.grid() + fig.tight_layout() + plt.show(block=False) + + acq_geom_tuned = acq_geom_tuned.update(par_name, min_par) + + return acq_geom_tuned diff --git a/src/corrct/alignment/fitting.py b/src/corrct/alignment/fitting.py index 1afd04c..74bd2c3 100644 --- a/src/corrct/alignment/fitting.py +++ b/src/corrct/alignment/fitting.py @@ -8,7 +8,9 @@ and ESRF - The European Synchrotron, Grenoble, France """ +from abc import ABC, abstractmethod from collections.abc import Sequence +from typing import Literal import matplotlib.pyplot as plt import numpy as np @@ -458,7 +460,7 @@ def f(apb: NDArrayFloat) -> float: return float(l1_diff) apb = spopt.minimize(f, np.array([a, p, b])) - (a, p, b) = apb.x + a, p, b = apb.x return a, p, b @@ -719,3 +721,396 @@ def refine_max_position_2d( + f" Input values: {f_vals}" ) return vertex_yx + + +def fit_parabola_min( + fun_x: ArrayLike | NDArray, + fun_vals: ArrayLike | NDArray, + scale: Literal["linear", "log"] = "linear", + decimals: int = 2, +) -> tuple[float, float, tuple[NDArray, NDArray] | None]: + """Parabolic fit local function stationary point. + + Parameters + ---------- + fx : ArrayLike + Parameter values. + f_vals : ArrayLike + Objective function costs of each parameter value. + scale : str, optional + Scale of the fit. Options are: "log" | "linear". The default is "log". + + Returns + ------- + min_fx : float + Expected parameter value of the fitted minimum. + min_f_val : float + Expected objective function cost of the fitted minimum. + """ + fun_x = np.array(fun_x) + fun_vals = np.array(fun_vals) + + if len(fun_x) < 3 or len(fun_vals) < 3 or len(fun_x) != len(fun_vals): + raise ValueError( + "Lengths of the parameter values and function values should be identical and >= 3." + + f"Given: fx={len(fun_x)}, f_vals={len(fun_vals)}" + ) + + if scale.lower() == "log": + + def to_fit(vals: NDArray) -> NDArray: + return np.log10(vals) + + def from_fit(vals: NDArray) -> NDArray: + return 10**vals + + elif scale.lower() == "linear": + + def to_fit(vals: NDArray) -> NDArray: + return vals + + from_fit = to_fit + else: + raise ValueError(f"Parameter 'scale' should be either 'log' or 'linear', given '{scale}' instead") + + min_pos = np.argmin(fun_vals) + if min_pos == 0: + print("WARNING: minimum value at the beginning of the lambda range.") + fx_fit = to_fit(fun_x[:3]) + f_vals_fit = fun_vals[:3] + elif min_pos == (len(fun_vals) - 1): + print("WARNING: minimum value at the end of the lambda range.") + fx_fit = to_fit(fun_x[-3:]) + f_vals_fit = fun_vals[-3:] + else: + fx_fit = to_fit(fun_x[min_pos - 1 : min_pos + 2]) + f_vals_fit = fun_vals[min_pos - 1 : min_pos + 2] + + # using Polynomial.fit, because it is supposed to be more numerically + # stable than previous solutions (according to numpy). + poly = Polynomial.fit(fx_fit, f_vals_fit, deg=2) + coeffs = poly.convert().coef + if coeffs[2] <= 0: + print("WARNING: fitted curve is concave. Returning minimum measured point.") + return fun_x[min_pos], fun_vals[min_pos], None + + # For a 1D parabola `f(x) = c + bx + ax^2`, the vertex position is: + # x_v = -b / 2a. + vertex_pos = -coeffs[1] / (2 * coeffs[2]) + vertex_val = coeffs[0] + vertex_pos * coeffs[1] / 2 + + vertex_pos = np.around(vertex_pos, decimals=decimals) + vertex_val = np.around(vertex_val, decimals=decimals) + + min_fx, min_f_val = from_fit(vertex_pos), vertex_val + if min_fx < fun_x[0] or min_fx > fun_x[-1]: + print( + f"WARNING: fitted stationary point {min_fx} is outside input range [{fun_x[0]}, {fun_x[-1]}]." + + " Returning minimum measured point." + ) + return fun_x[min_pos], fun_vals[min_pos], None + + return min_fx, min_f_val, (coeffs, fx_fit) + + +def fit_ellipse_center( + prj_points_vu: NDArray, rescale: bool = True, use_l1_norm: bool = False, decimals: int | None = 2 +) -> NDArray: + """ + Fit an ellipse center to a set of projected points in VU coordinates. + + The function uses a least-squares approach to fit an ellipse center to the given points. + Optionally, it can use L1 norm for fitting instead of the default L2 norm, and rescale + the points during fitting to have a maximum range of 1 (improve numerical stability). + + Parameters + ---------- + prj_points_vu : NDArray + Projected points in VU coordinates. The expected organization is: + - Last dimension: List of points (each point is a 2D coordinate). + - First dimension: Coordinates (V and U). + rescale : bool, optional + If True, rescale the points to have a maximum range of 1. Default is True. + use_l1_norm : bool, optional + If True, use L1 norm for fitting instead of the default L2 norm. Default is False. + decimals : int | None, optional + The number of decimal places to round the result to. If None, no rounding is performed. + Default is 2. + + Returns + ------- + NDArray + The fitted ellipse center in VU coordinates. + """ + c_vu = np.mean(prj_points_vu, axis=-1) + pos_vu = prj_points_vu - c_vu[:, None] + + if rescale: + scale_vu = np.max(pos_vu, axis=-1) - np.min(pos_vu, axis=-1) + pos_vu /= scale_vu[:, None] + else: + scale_vu = np.ones(2, dtype=pos_vu.dtype) + + num_lines = pos_vu.shape[-1] // 2 + pos1_vu = pos_vu[:, :num_lines] + pos2_vu = pos_vu[:, num_lines : num_lines * 2] + + diffs_vu = pos2_vu - pos1_vu + vandermonde = np.stack([diffs_vu[-1, :], -diffs_vu[-2, :]], axis=-1) + values = np.cross(pos1_vu, pos2_vu, axis=0) + + p_vu = np.linalg.lstsq(vandermonde, values, rcond=None)[0] + if use_l1_norm: + + def _func(params: NDArrayFloat) -> float: + predicted_values = vandermonde.dot(params) + l1_diff = np.linalg.norm(predicted_values - values, ord=1) + return float(l1_diff) + + opt_p_vu = spopt.minimize(_func, p_vu) + p_vu = opt_p_vu.x + + pred_c_vu = p_vu * scale_vu + c_vu + + if decimals is not None: + pred_c_vu = np.around(pred_c_vu, decimals=decimals) + + return pred_c_vu + + +def fit_ellipse_parameters( + prj_points_vu: NDArray, rescale: bool = True, use_l1_norm: bool = False +) -> tuple[float, float, float, float, float]: + """ + Fit ellipse parameters to a set of projected points in VU coordinates. + + Parameters + ---------- + prj_points_vu : NDArray + Projected points in VU coordinates. The expected organization is: + - Last dimension: List of points (each point is a 2D coordinate). + - First dimension: Coordinates (V and U). + rescale : bool, optional + If True, rescale the points to have a maximum range of 1. Default is True. + use_l1_norm : bool, optional + If True, use L1 norm for fitting instead of the default L2 norm. Default is False. + + Returns + ------- + tuple[float, float, float, float, float] + The fitted ellipse parameters: a, b, c, u, v. + """ + # First we fit 5 intermediate variables + p_u: NDArray = prj_points_vu[-1, :] + p_v: NDArray = prj_points_vu[-2, :] + + if rescale: + c_u = float(np.mean(p_u)) + c_v = float(np.mean(p_v)) + p_u = p_u - c_u + p_v = p_v - c_v + + p_u_scaling = float(np.abs(p_u).max()) + p_v_scaling = float(np.abs(p_v).max()) + p_u /= p_u_scaling + p_v /= p_v_scaling + else: + c_u = 0.0 + c_v = 0.0 + p_u_scaling = 1.0 + p_v_scaling = 1.0 + + vandermonde = np.stack([p_u**2, -2 * p_u, -2 * p_v, 2 * p_u * p_v, np.ones_like(p_u)], axis=-1) + values = -(p_v**2) + + coeffs = np.linalg.lstsq(vandermonde, values, rcond=None)[0] + if use_l1_norm: + + def _func(pars: NDArrayFloat) -> float: + predicted_b = vandermonde.dot(pars) + l1_diff = np.linalg.norm(predicted_b - values, ord=1) + return float(l1_diff) + + opt_params = spopt.minimize(_func, coeffs) + coeffs = opt_params.x + + if rescale: + coeffs[0] *= (p_v_scaling**2) / (p_u_scaling**2) + coeffs[1] *= (p_v_scaling**2) / p_u_scaling + coeffs[2] *= p_v_scaling + coeffs[3] *= p_v_scaling / p_u_scaling + coeffs[4] *= p_v_scaling**2 + + u = (coeffs[1] - coeffs[2] * coeffs[3]) / (coeffs[0] - coeffs[3] ** 2) + v = (coeffs[0] * coeffs[2] - coeffs[1] * coeffs[3]) / (coeffs[0] - coeffs[2] * coeffs[3]) + + a = coeffs[0] / (coeffs[0] * u**2 + v**2 + 2 * coeffs[3] * u * v - coeffs[4]) + b = a / coeffs[0] + c = coeffs[3] * b + + u += c_u + v += c_v + + return a, b, c, u, v + + +class Trajectory(ABC): + """Base trajectory class.""" + + @abstractmethod + def __call__(self, uus: Sequence[float] | NDArray) -> Sequence[NDArray]: + """Compute V coordinates, given V coordinates. + + Parameters + ---------- + uus : Sequence[float] | NDArray + The U coordinates + + Returns + ------- + Sequence[NDArray] + Corresponding V coordinates, given the multiplicity of the trajectory + """ + + +class Ellipse(Trajectory): + """Elliptic trajectory class.""" + + a: float + b: float + c: float + u: float + v: float + + c_vu: NDArrayFloat + + def __init__(self, a: float, b: float, c: float, u: float, v: float, c_vu: NDArrayFloat) -> None: + """Initialize ellipse class. + + Ellipse corresponding to the equation: a*(x - u)**2 + b*(y - v)**2 + 2*c*(x - u)*(y - v) = 1 + + Parameters + ---------- + a : float + The semi-major axis of the ellipse. + b : float + The semi-minor axis of the ellipse. + c : float + The rotation angle of the ellipse in radians. + u : float + The center of the ellipse along the x-axis. + v : float + The center of the ellipse along the y-axis. + c_vu : NDArrayFloat + The projected center of the circular orbit generating the ellipse. + """ + self.a = a + self.b = b + self.c = c + self.u = u + self.v = v + self.c_vu = c_vu + + def __repr__(self) -> str: + """Return a string representation of the Ellipse instance.""" + return f"Ellipse(a={self.a}, b={self.b}, c={self.c}, u={self.u}, v={self.v}, c_vu={self.c_vu})" + + @property + def extremes_u(self) -> tuple[float, float]: + """ + Find the most extreme x coordinates (U coordinates) that are still valid for the given ellipse equation. + + Returns + ------- + tuple[float, float] + A tuple containing the most extreme U coordinates (u_min, u_max). + """ + # Calculate the extreme U coordinates + u_min = self.u - np.sqrt(1 / self.a) + u_max = self.u + np.sqrt(1 / self.a) + + return u_min, u_max + + @property + def center_vu(self) -> NDArray: + """Return the fitted ellipse center. + + Returns + ------- + NDArray + The fitted center position. + """ + return self.c_vu + + @property + def parameters(self) -> tuple[float, float, float, float, float]: + """Return the fitted ellipse parameters. + + Returns + ------- + tuple[float, float, float, float, float] + The fitted ellipse parameters: b, a, c, v, u. + """ + return self.b, self.a, self.c, self.v, self.u + + def __call__(self, uus: ArrayLike | NDArray) -> Sequence[NDArray]: + """Predict V coordinates of ellipse from its parameters, and U coordinates. + + Parameters + ---------- + uus : Union[ArrayLike, NDArray] + The U coordinates + + Returns + ------- + tuple[NDArray, NDArray] + The corresponding top and bottom V coordinates + """ + b, a, c, v, u = np.array(self.parameters) + uus = np.array(uus) + + uus_u = uus - u + + a_tilde = b + # b_tilde = 2 * (-b * v + c * uus - c * u) + # c_tilde = -(1 - a * (uus - u) ** 2 - b * v**2 + 2 * c * v * (uus - u)) + b_tilde = 2 * (-b * v + c * uus_u) + c_tilde = -(1 - a * uus_u**2 - b * v**2 + 2 * c * v * uus_u) + delta_tilde = np.sqrt(b_tilde**2 - 4 * a_tilde * c_tilde) + v_1 = (-b_tilde + delta_tilde) / (2 * a_tilde) + v_2 = (-b_tilde - delta_tilde) / (2 * a_tilde) + + return v_1, v_2 + + +def fit_ellipse(prj_points_vu: ArrayLike | NDArray, rescale: bool = True, use_l1_norm: bool = False) -> Ellipse: + """Fit an ellipse to a set of 2D points using either least-squares or l1-norm optimization. + + Parameters + ---------- + prj_points_vu : ArrayLike | NDArray + A list or array of 2D points (shape: Nx2) representing the trajectory to fit. + rescale : bool, optional + Whether to rescale the data within the interval [-1, 1] to improve numerical stability. + Default is True. + use_l1_norm : bool, optional + Whether to use the l1-norm or the least-squares (l2-norm) fit for optimization. + Default is False. + + Returns + ------- + Ellipse + An Ellipse object containing the fitted ellipse parameters (center, axes, and rotation angle). + + Notes + ----- + The function first fits the ellipse parameters (axes and rotation angle) and then fits the center + separately. The optimization method can be switched between least-squares and l1-norm for robustness + against outliers. + """ + prj_points_vu = np.array(prj_points_vu) + + return Ellipse( + *fit_ellipse_parameters(prj_points_vu, rescale, use_l1_norm=use_l1_norm), + fit_ellipse_center(prj_points_vu, rescale, use_l1_norm=use_l1_norm), + ) diff --git a/src/corrct/alignment/markers.py b/src/corrct/alignment/markers.py new file mode 100644 index 0000000..43de374 --- /dev/null +++ b/src/corrct/alignment/markers.py @@ -0,0 +1,215 @@ +#!/usr/bin/env python3 +""" +Fiducial marker tracking routines. + +@author: Nicola VIGANÒ, ESRF - The European Synchrotron, Grenoble, France, +and CEA-IRIG, Grenoble, France +""" + +from collections.abc import Sequence + +import numpy as np +import scipy.ndimage as spimg +from numpy.typing import NDArray + +import matplotlib.pyplot as plt + +from . import fitting + + +def cm2inch(dims: Sequence[float] | NDArray) -> tuple[float]: + """Convert dimensions from centimeters to inches. + + Parameters + ---------- + dims : Sequence[float] | NDArray + The dimensions of the object in centimeters. Can be a sequence of floats or a NumPy array. + + Returns + ------- + tuple[float] + The converted dimensions in inches, as a tuple of floats. + """ + return tuple(np.array(dims) / 2.54) + + +def track_marker(prj_data: NDArray, marker_vu: NDArray, stack_axis: int = -2) -> NDArray: + """Track marker position in a stack of images. + + Parameters + ---------- + prj_data_vwu : NDArray + The projection data. + marker_vu : NDArray + The fiducial marker to track in VU. + stack_axis : int, optional + The axis along which the images are stacked. The default is -2. + + Returns + ------- + NDArray + List of positions for each image. + """ + marker_v1u = np.expand_dims(marker_vu, stack_axis).astype(np.float32) + marker_pos = fitting.fit_shifts_vu_xc(prj_data, marker_v1u, stack_axis=stack_axis, normalize_fourier=False) + return marker_pos + np.array(marker_vu.shape)[:, None] / 2 + + +def create_marker_disk( + data_shape_vu: Sequence[int] | NDArray, radius: float, super_sampling: int = 5, conv: bool = True +) -> NDArray: + """Create a disk-shaped marker for tracking a calibration object's movement. + + Parameters + ---------- + data_shape_vu : Sequence[int] | NDArray + Shape of the images (vertical, horizontal). Can be a sequence of integers or a NumPy array. + radius : float + Radius of the marker in pixels. + super_sampling : int, optional + Super-sampling factor for the coordinates used in marker creation. Default is 5. + conv : bool, optional + Whether to convolve the initial marker with itself. Default is True. + + Returns + ------- + NDArray + An image of the same size as the input projections, with the marker centered. + """ + data_shape_vu = np.array(data_shape_vu, dtype=int) * super_sampling + + # coords = [np.linspace(-(s - 1) / 2, (s - 1) / 2, s, dtype=np.float32) for s in data_shape_vu] + coords = [np.fft.fftfreq(d, 1 / d) for d in data_shape_vu] + coords = np.stack(np.meshgrid(*coords, indexing="ij"), axis=0) + pix_rr = np.sqrt(np.sum(coords**2, axis=0)) + + probe = pix_rr < radius * super_sampling + probe = np.roll(probe, super_sampling // 2, axis=tuple(np.arange(len(data_shape_vu)))) + new_shape = np.stack([data_shape_vu // super_sampling, np.ones_like(data_shape_vu) * super_sampling], axis=1).flatten() + probe = probe.reshape(new_shape) + probe = np.mean(probe, axis=tuple(np.arange(1, len(data_shape_vu) * 2, 2, dtype=int))) + + probe = np.fft.fftshift(probe) + + if conv: + probe = spimg.convolve(probe, probe) + + return probe + + +class MarkerTrackingVisualizer: + """Plotting class to assess the marker tracking quality.""" + + def __init__( + self, + fitted_positions_vu: NDArray, + images: NDArray, + marker: NDArray, + trajectory: fitting.Trajectory | None = None, + ) -> None: + """Initialize the visualization utility for checking the marker position fitting. + + Parameters + ---------- + fitted_positions_vu : NDArray + The fitted positions of the marker + imgs : NDArray + The original images + disk : NDArray + The marker image + trajectory : Union[fitting.Trajectory, None], optional + The trajectory object that the points are supposed to follow, by default None + """ + self.positions_vu = np.array(fitted_positions_vu) + self.images = images + self.marker = marker + self.global_lims = False + + self.trajectory = trajectory + + self.curr_pos = 0 + + if self.trajectory is not None: + uus = np.sort(self.positions_vu[1, :]) + self.vvs = self.trajectory(uus) + + self.fig, self.axs = plt.subplots(1, 3, figsize=cm2inch([36, 12])) # , sharex=True, sharey=True + self.axs[2].imshow(self.marker) + self.axs[0].set_xlim(0, self.images.shape[-1]) + self.axs[0].set_ylim(self.images.shape[-3], 0) + self.fig.tight_layout() + self._update() + + self.fig.canvas.mpl_connect("key_press_event", self._key_event) + self.fig.canvas.mpl_connect("scroll_event", self._scroll_event) + + def _update(self) -> None: + self.curr_pos = self.curr_pos % self.images.shape[-2] + + for img in self.axs[0].get_images(): + img.remove() + x_lims = self.axs[0].get_xlim() + y_lims = self.axs[0].get_ylim() + self.axs[0].cla() + self.axs[0].set_xlim(x_lims[0], x_lims[1]) + self.axs[0].set_ylim(y_lims[0], y_lims[1]) + + for img in self.axs[1].get_images(): + img.remove() + self.axs[1].cla() + + self.axs[0].plot(self.positions_vu[1, :], self.positions_vu[0, :], "bo-", markersize=4) + self.axs[0].scatter(self.positions_vu[1, self.curr_pos], self.positions_vu[0, self.curr_pos], c="r") + + if self.trajectory is not None: + uus = np.sort(self.positions_vu[1, :]) + for vvs in self.vvs: + self.axs[0].plot(uus, vvs, "g") + self.axs[0].grid() + + if self.global_lims: + vmin = self.images.min() + vmax = self.images.max() + else: + vmin = self.images[:, self.curr_pos, :].min() + vmax = self.images[:, self.curr_pos, :].max() + + img = self.axs[1].imshow(self.images[:, self.curr_pos, :], vmin=vmin, vmax=vmax) + self.axs[1].scatter(self.positions_vu[1, self.curr_pos], self.positions_vu[0, self.curr_pos], c="r") + self.axs[1].set_title(f"Range: [{vmin}, {vmax}]") + # plt.colorbar(im, ax=self.axs[1]) + self.fig.canvas.draw() + + def _key_event(self, evnt) -> None: + if evnt.key == "right": + self.curr_pos += 1 + elif evnt.key == "left": + self.curr_pos -= 1 + elif evnt.key == "up": + self.curr_pos += 1 + elif evnt.key == "down": + self.curr_pos -= 1 + elif evnt.key == "pageup": + self.curr_pos += 10 + elif evnt.key == "pagedown": + self.curr_pos -= 10 + elif evnt.key == "escape": + plt.close(self.fig) + elif evnt.key == "ctrl+l": + self.global_lims = not self.global_lims + else: + print(evnt.key) + return + + self._update() + + def _scroll_event(self, evnt) -> None: + if evnt.button == "up": + self.curr_pos += 1 + elif evnt.button == "down": + self.curr_pos -= 1 + else: + print(evnt.key) + return + + self._update() diff --git a/src/corrct/models.py b/src/corrct/models.py index a203cb5..377cca1 100644 --- a/src/corrct/models.py +++ b/src/corrct/models.py @@ -13,8 +13,11 @@ from dataclasses import replace as dc_replace from typing import Any +import matplotlib.pyplot as plt import numpy as np import scipy.spatial.transform as spt +from matplotlib.widgets import CheckButtons +from mpl_toolkits.mplot3d.art3d import Line3DCollection, Poly3DCollection from numpy.typing import ArrayLike, NDArray ROT_DIRS_VALID = ("clockwise", "counter-clockwise") @@ -788,3 +791,136 @@ def get_vol_geom_from_volume(volume: NDArray) -> VolumeGeometry: if len(vol_shape_zxy) < 2: raise ValueError(f"The volume should be at least 2-dimensional, but the following shape was passed: {vol_shape_zxy}") return VolumeGeometry(np.array([vol_shape_zxy[-2], vol_shape_zxy[-1], *np.flip(vol_shape_zxy[:-2])])) + + +def plot_projection_geometry(prj_geom: ProjectionGeometry, vol_geom: VolumeGeometry): + """ + Plot the projection geometry using matplotlib. + + Parameters + ---------- + prj_geom : ProjectionGeometry + ProjectionGeometry class with precomputed geometric properties. + vol_size_xyz : VolumeGeometry + VolumeGeometry class reporting the volume size in number of voxels in X, Y, and Z directions, and their size. + """ + if prj_geom.det_shape_vu is None or len(prj_geom.det_shape_vu) != 2: + raise ValueError(f"prj_geom.det_shape_vu must be a 2-element sequence, got {prj_geom = } instead.") + # Check that the fields have exactly 3 elements + for field in ['det_pos_xyz', 'det_u_xyz', 'det_v_xyz', 'src_pos_xyz']: + field_value = getattr(prj_geom, field) + if field_value.size != 3: + raise ValueError(f"{field} must have exactly 3 elements, got {field_value.size} elements instead.") + if len(vol_geom.shape_xyz) != 3: + raise ValueError(f"vol_geom.shape_xyz must be a 3-element sequence, got {vol_geom = } instead.") + + fig = plt.figure(figsize=(10, 8)) + ax = fig.add_subplot(111, projection='3d') + + src_pos = prj_geom.src_pos_xyz.flatten() + det_pos = prj_geom.det_pos_xyz.flatten() + det_u = prj_geom.det_u_xyz.flatten() + det_v = prj_geom.det_v_xyz.flatten() + + pix_size = vol_geom.vox_size / prj_geom.pix2vox_ratio + + # Plot the detector + det_corners = np.array( + [ + det_pos + det_u * prj_geom.det_shape_vu[1] * pix_size / 2 + det_v * prj_geom.det_shape_vu[0] * pix_size / 2, + det_pos + det_u * prj_geom.det_shape_vu[1] * pix_size / 2 - det_v * prj_geom.det_shape_vu[0] * pix_size / 2, + det_pos - det_u * prj_geom.det_shape_vu[1] * pix_size / 2 - det_v * prj_geom.det_shape_vu[0] * pix_size / 2, + det_pos - det_u * prj_geom.det_shape_vu[1] * pix_size / 2 + det_v * prj_geom.det_shape_vu[0] * pix_size / 2, + ] + ) + + detector = Poly3DCollection([det_corners], alpha=0.5, linewidths=1, edgecolors='k') + detector.set_facecolor('b') + detector.set_label("Detector") + ax.add_collection3d(detector) + + # Plot the volume + x, y, z = vol_geom.shape_xyz * vol_geom.vox_size / 2 + + # Create a cube + cube_vertices = np.array( + [[-x, -y, -z], [x, -y, -z], [x, y, -z], [-x, y, -z], [-x, -y, z], [x, -y, z], [x, y, z], [-x, y, z]] + ) + + # Create the 8 faces of the cube + cube_faces = [ + [cube_vertices[0], cube_vertices[1], cube_vertices[2], cube_vertices[3]], # Bottom face + [cube_vertices[4], cube_vertices[5], cube_vertices[6], cube_vertices[7]], # Top face + [cube_vertices[0], cube_vertices[1], cube_vertices[5], cube_vertices[4]], # Front face + [cube_vertices[2], cube_vertices[3], cube_vertices[7], cube_vertices[6]], # Back face + [cube_vertices[1], cube_vertices[2], cube_vertices[6], cube_vertices[5]], # Right face + [cube_vertices[0], cube_vertices[3], cube_vertices[7], cube_vertices[4]], # Left face + ] + + volume = Poly3DCollection(cube_faces, alpha=0.1, linewidths=1, edgecolors='k') + volume.set_facecolor('g') + volume.set_label("Volume") + ax.add_collection3d(volume) + + # Organize the objects and their labels in a list + objects: list = [detector, volume] + + if prj_geom.geom_type.lower() == "cone": + prj_lines = [np.stack((src_pos, d_c), axis=0) for d_c in det_corners] + prj_lines_c = Line3DCollection(prj_lines, colors="g", linestyles="-.", linewidths=1.0, label="Projection lines") + ax.add_collection(prj_lines_c) + objects.append(prj_lines_c) + + # Plot the source + src_scatter = ax.scatter(*src_pos, color='r', s=100, label='Source Spot') + objects.append(src_scatter) + + # Plot vectors from the origin to the source and detector center + src_quiver = ax.quiver(*([0.0] * 3), *src_pos, color='r', arrow_length_ratio=0.1, label='Source Position') + objects.append(src_quiver) + else: + # plot the projection direction + prj_dir = src_pos / np.linalg.norm(src_pos) * np.linalg.norm(det_pos) + prj_quiver = ax.quiver(*([0.0] * 3), *prj_dir, color='r', arrow_length_ratio=0.1, label='Projection direction') + objects.append(prj_quiver) + + # Plot vectors from the origin to the source and detector center + det_quiver = ax.quiver(*([0.0] * 3), *det_pos, color='b', arrow_length_ratio=0.1, label='Detector Center') + objects.append(det_quiver) + + # Plot vectors from the detector center to the u and v unit vectors + u_quiver = ax.quiver( + *det_pos, *(det_u * prj_geom.det_shape_vu[1] * pix_size / 4), color='m', arrow_length_ratio=0.1, label='U vector' + ) + objects.append(u_quiver) + v_quiver = ax.quiver( + *det_pos, *(det_v * prj_geom.det_shape_vu[0] * pix_size / 4), color='c', arrow_length_ratio=0.1, label='V vector' + ) + objects.append(v_quiver) + + labels = [o.get_label() for o in objects] + colors = [o.get_color() if hasattr(o, "get_color") else o.get_facecolor()[0] for o in objects] + visibility = [o.get_visible() for o in objects] + + # Create a rectangle for the CheckButtons + rax = ax.inset_axes([0.0, 0.0, 0.3, 0.25]) + visibility = [True] * len(objects) + check = CheckButtons(rax, labels, visibility, label_props={'color': colors, "fontsize": [14] * len(objects)}) + + # Define the function to update the visibility of the elements + def func(label): + index = labels.index(label) + objects[index].set_visible(not objects[index].get_visible()) + plt.draw() + + # Register the function with the CheckButtons + check.on_clicked(func) + + # Set labels + ax.set_xlabel('X') + ax.set_ylabel('Y') + ax.set_zlabel('Z') + ax.set_aspect("equal") + fig.tight_layout() + + plt.show() diff --git a/tests/alignment/__init__.py b/tests/alignment/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/alignment/test_cone_beam.py b/tests/alignment/test_cone_beam.py new file mode 100644 index 0000000..f4f570f --- /dev/null +++ b/tests/alignment/test_cone_beam.py @@ -0,0 +1,236 @@ +import numpy as np +import pytest +from numpy.typing import NDArray +from scipy.spatial.transform import Rotation + +from corrct.alignment.cone_beam import FitConeBeamGeometry, ConeBeamGeometry +from corrct.models import ProjectionGeometry, plot_projection_geometry + + +def project_to_uv(points_xyz: NDArray, proj_geom: ProjectionGeometry) -> NDArray: + """ + Project a list of XYZ points to UV coordinates using precomputed geometry. + + Parameters + ---------- + points_xyz : NDArray + Array of (x, y, z) points in the sample volume (XYZ coordinates). + geo : ProjectionGeometry + ProjectionGeometry class with precomputed geometry properties. + + Returns + ------- + NDArray + Array of (u, v) coordinates for each input point. + """ + points_xyz = np.array(points_xyz, dtype=np.float32) + + det_u_xyz = proj_geom.det_u_xyz.flatten() + det_v_xyz = proj_geom.det_v_xyz.flatten() + + det_normal_xyz = np.cross(det_u_xyz, det_v_xyz) + + vecs_src_2_pnts = points_xyz - proj_geom.src_pos_xyz + vecs_src_2_pnts = vecs_src_2_pnts / np.linalg.norm(vecs_src_2_pnts, axis=-1, keepdims=True) + + denominators = np.dot(det_normal_xyz, vecs_src_2_pnts.T) + numerators = (proj_geom.det_pos_xyz - points_xyz).dot(det_normal_xyz) + + invalid_points = np.isclose(denominators, 0.0) + valid_denoms = np.logical_not(invalid_points) + + lams = numerators[list(valid_denoms)] / denominators[list(valid_denoms)] + prj_pnts_xyz = vecs_src_2_pnts[list(valid_denoms), :] * lams[:, None] + points_xyz[list(valid_denoms), :] + + uv_coords = np.empty((len(points_xyz), 2), dtype=np.float32) + uv_coords.fill(np.nan) + uv_coords[list(valid_denoms), :] = np.stack((prj_pnts_xyz.dot(det_u_xyz), prj_pnts_xyz.dot(det_v_xyz)), axis=1) + + return uv_coords + + +def compute_xyz_rotated( + r: float, z: float, ws: NDArray, voxel_size_xyz: tuple[float, float, float] = (1.0, 1.0, 1.0) +) -> NDArray: + """ + Compute XYZ coordinates for an array of (r, z, w) inputs in the rotated volume. + + Parameters + ---------- + r : float + Radial distance from volume center. + z : float + Elevation (Z coordinates). + ws : NDArray + Array of rotation angles (omega) in radians. + voxel_size_xyz : tuple[float, float, float], optional + Voxel sizes in X, Y, Z directions (default: (1.0, 1.0, 1.0)). + + Returns + ------- + NDArray + Array of (x, y, z) coordinates in the rotated volume. + """ + ws = np.squeeze(np.array(ws)) + + # Check if ws is one-dimensional + if ws.ndim != 1: + raise ValueError("Input array ws must be one-dimensional.") + + # At w=0, the point is along Y: (0, r, z) + # Rotation matrix for w (omega) around Z: + # [cos(w), -sin(w), 0] + # [sin(w), cos(w), 0] + # [0, 0, 1] + rot_w = Rotation.from_rotvec(np.array([np.zeros_like(ws), np.zeros_like(ws), ws]).T) + + rot_xyz = rot_w.apply(np.array([0, r, z]).T) + + # Scale by voxel sizes + return rot_xyz * np.array(voxel_size_xyz) + + +@pytest.mark.parametrize("R", [30.0, 50.0]) +@pytest.mark.parametrize("D", [70.0, 80.0]) +def test_cone_beam_ellipse_distances(R: float, D: float): + """Test the cone beam geometry determination for different distance combinations.""" + debug = False + + det_size_vu = 12 + + acq_geom_tst = ConeBeamGeometry(theta_deg=0.0, phi_deg=0.0, eta_deg=0.0, D_pix=D, R_pix=R) + proj_geom_tst = acq_geom_tst.get_prj_geom() + if debug: + print(f"{acq_geom_tst = }") + plot_projection_geometry(proj_geom_tst, acq_geom_tst.get_vol_geom()) + + r = det_size_vu / 2 - 1 + z_u = -2.5 + z_l = 1.5 + ws = np.deg2rad(np.arange(0, 360, 6)) + points_xyz_u = compute_xyz_rotated(r=r, z=z_u, ws=ws) + points_xyz_l = compute_xyz_rotated(r=r, z=z_l, ws=ws) + + # fig, axs = plt.subplots(1, 1) + # axs.plot(points_xyz[:, 0], points_xyz[:, 1]) + # axs.grid() + # fig.tight_layout() + # plt.show() + + prj_uv_coords_u = project_to_uv(points_xyz_u, proj_geom_tst) + prj_uv_coords_l = project_to_uv(points_xyz_l, proj_geom_tst) + + # fig, axs = plt.subplots(1, 1) + # axs.plot(prj_uv_coords_u[:, 0], prj_uv_coords_u[:, 1]) + # axs.plot(prj_uv_coords_l[:, 0], prj_uv_coords_l[:, 1]) + # axs.grid() + # fig.tight_layout() + # plt.show() + + prj_size_vu = np.array((det_size_vu, det_size_vu)) + + points_ell1 = np.flip(prj_uv_coords_u.T, axis=0) + prj_size_vu[:, None] / 2 + points_ell2 = np.flip(prj_uv_coords_l.T, axis=0) + prj_size_vu[:, None] / 2 + + if debug: + print(f"{points_ell1 = }") + print(f"{points_ell2 = }") + + geom_fit = FitConeBeamGeometry(prj_size_vu=prj_size_vu, points_ell1=points_ell1, points_ell2=points_ell2) + geom_fit.fit(r=r, e=1.0) + + if debug: + print(f"{geom_fit.acq_geom = }") + + assert np.isclose( + D, geom_fit.acq_geom.D_pix, rtol=1e-5, atol=1e-3 + ), f"The fitted source-detector distance is wrong ({geom_fit.acq_geom.D_pix:.6} vs expected: {D:.6})" + assert np.isclose( + R, geom_fit.acq_geom.R_pix, rtol=1e-5, atol=1e-3 + ), f"The fitted source-sample distance is wrong ({geom_fit.acq_geom.R_pix:.6} vs expected: {R:.6})" + + +@pytest.mark.parametrize(("p", "n"), [(30.0, 0.0), (5.0, 1.0), (0.0, 15.0)]) +def test_cone_beam_ellipse_angles(p: float, n: float): + """Test the cone beam geometry determination for different angles combinations.""" + debug = False + + det_size_vu = 12 + + acq_geom_tst = ConeBeamGeometry(theta_deg=0.0, phi_deg=p, eta_deg=n, D_pix=80.0, R_pix=50.0) + proj_geom_tst = acq_geom_tst.get_prj_geom() + if debug: + print(f"{acq_geom_tst = }") + print(f"{proj_geom_tst = }") + plot_projection_geometry(proj_geom_tst, acq_geom_tst.get_vol_geom()) + + r = det_size_vu / 2 - 1 + z_u = -2.5 + z_l = 1.5 + ws = np.deg2rad(np.arange(0, 360, 6)) + points_xyz_u = compute_xyz_rotated(r=r, z=z_u, ws=ws) + points_xyz_l = compute_xyz_rotated(r=r, z=z_l, ws=ws) + + # fig, axs = plt.subplots(1, 1) + # axs.plot(points_xyz[:, 0], points_xyz[:, 1]) + # axs.grid() + # fig.tight_layout() + # plt.show() + + prj_uv_coords_u = project_to_uv(points_xyz_u, proj_geom_tst) + prj_uv_coords_l = project_to_uv(points_xyz_l, proj_geom_tst) + + # fig, axs = plt.subplots(1, 1) + # axs.plot(prj_uv_coords_u[:, 0], prj_uv_coords_u[:, 1]) + # axs.plot(prj_uv_coords_l[:, 0], prj_uv_coords_l[:, 1]) + # axs.grid() + # fig.tight_layout() + # plt.show() + + prj_size_vu = np.array((det_size_vu, det_size_vu)) + + points_ell1 = np.flip(prj_uv_coords_u.T, axis=0) + prj_size_vu[:, None] / 2 + points_ell2 = np.flip(prj_uv_coords_l.T, axis=0) + prj_size_vu[:, None] / 2 + + if debug: + print(f"{points_ell1 = }") + print(f"{points_ell2 = }") + + geom_fit = FitConeBeamGeometry(prj_size_vu=prj_size_vu, points_ell1=points_ell1, points_ell2=points_ell2) + geom_fit.fit(r=r, e=1.0) + + if debug: + print(f"{geom_fit.acq_geom = }") + + assert np.isclose( + p, geom_fit.acq_geom.phi_deg, rtol=5e-2, atol=5e-1 + ), f"The fitted phi angle is wrong ({geom_fit.acq_geom.phi_deg:.6} vs expected: {p:.6})" + assert np.isclose( + n, geom_fit.acq_geom.eta_deg, rtol=1e-2, atol=1e-1 + ), f"The fitted yaw angle is wrong ({geom_fit.acq_geom.eta_deg:.6} vs expected: {n:.6})" + + +if __name__ == "__main__": + points_xyz_tst = np.array([[0, 1, 2], [3, 4, 5]]) + + acq_geom = ConeBeamGeometry( + theta_deg=0.0, + phi_deg=0.0, + eta_deg=np.rad2deg(np.pi / 4), + D_pix=260.0, + R_pix=160.0, + det_size_u_pix=100, + det_size_v_pix=100, + ) + print(acq_geom) + + prj_geom = acq_geom.get_prj_geom() + prj_uv_coords = project_to_uv(points_xyz=points_xyz_tst, proj_geom=prj_geom) + print(prj_uv_coords) + + # Compute rotated XYZ for arrays of (r, z, w) + ws_test = np.array([0, np.pi / 4, np.pi / 2]) + xyz_rotated = compute_xyz_rotated(r=2.0, z=1.0, ws=ws_test) + print(xyz_rotated) + + plot_projection_geometry(acq_geom.get_prj_geom(), acq_geom.get_vol_geom()) diff --git a/tests/alignment/test_fitting.py b/tests/alignment/test_fitting.py new file mode 100644 index 0000000..5d88483 --- /dev/null +++ b/tests/alignment/test_fitting.py @@ -0,0 +1,123 @@ +import numpy as np +import pytest +from numpy.typing import NDArray + +from corrct.alignment.fitting import Ellipse, fit_ellipse_center, fit_ellipse_parameters + + +def generate_points(ellipse: Ellipse, num_points: int, noise_std: float = 0.0, outliers: int = 0) -> NDArray: + """Generate points on the ellipse with optional noise and outliers. + + Parameters + ---------- + ellipse : Ellipse + The ellipse to generate points from. + num_points : int + Number of points to generate. + noise : float, optional + Standard deviation of the Gaussian noise to add to the points. Default is 0.0. + outliers : int, optional + Number of outliers to add to the points. Default is 0. + + Returns + ------- + NDArray + Generated points in VU coordinates. + """ + angles = np.linspace(0, np.pi, num_points) + x = np.cos(angles) / np.sqrt(ellipse.a) + ellipse.u + y = ellipse(x) + + points = np.stack([[*y[0], *np.flip(y[1])], [*x, *np.flip(x)]], axis=0) + + if noise_std > 0.0: + points += np.random.normal(scale=noise_std, size=points.shape) + + if outliers > 0: + outlier_indices = np.random.choice(num_points, size=outliers, replace=False) + points[:, outlier_indices] += np.random.normal(scale=5 * noise_std, size=(2, outliers)) + + # import matplotlib.pyplot as plt + + # print(ellipse) + # print(points) + # fig, axs = plt.subplots(1, 1) + # axs.plot(x, y[0]) + # axs.plot(x, y[1]) + # axs.scatter(points[1], points[0], color="C2") + # fig.tight_layout() + # plt.show() + + return points + + +def test_fit_ellipse_center_simple() -> None: + """Test fit_ellipse_center with no noise and 6 points.""" + ellipse = Ellipse(a=2.0, b=1.5, c=0.5, u=1.0, v=1.0, c_vu=np.ones(2)) + points = generate_points(ellipse, num_points=6) + + center = fit_ellipse_center(points, rescale=True, use_l1_norm=False) + + assert np.allclose(center, np.array([[1.0], [1.0]]), atol=1e-2), "Center should be close to [1.0, 1.0]" + + +def test_fit_ellipse_center_noisy() -> None: + """Test fit_ellipse_center with small Gaussian noise and 60 points.""" + ellipse = Ellipse(a=2.0, b=1.5, c=0.5, u=1.0, v=1.0, c_vu=np.ones(2)) + points = generate_points(ellipse, num_points=60, noise_std=0.1) + + center = fit_ellipse_center(points, rescale=True, use_l1_norm=False) + + assert np.allclose(center, np.array([[1.0], [1.0]]), atol=1e-1), "Center should be close to [1.0, 1.0]" + + +def test_fit_ellipse_center_noisy_with_outliers() -> None: + """Test fit_ellipse_center with 60 points and 2 outliers using L1 norm.""" + ellipse = Ellipse(a=2.0, b=1.5, c=0.5, u=1.0, v=1.0, c_vu=np.ones(2)) + points = generate_points(ellipse, num_points=60, noise_std=0.1, outliers=2) + + center = fit_ellipse_center(points, rescale=True, use_l1_norm=True) + + assert np.allclose(center, np.array([[1.0], [1.0]]), atol=1e-1), "Center should be close to [1.0, 1.0]" + + +def test_fit_ellipse_parameters_simple() -> None: + """Test fit_ellipse_parameters with no noise and 6 points.""" + ellipse = Ellipse(a=2.0, b=1.5, c=0.5, u=1.0, v=1.0, c_vu=np.ones(2)) + points = generate_points(ellipse, num_points=6) + + a, b, c, u, v = fit_ellipse_parameters(points, rescale=True, use_l1_norm=False) + + assert np.isclose(a, 2.0, atol=1e-2), "Semi-major axis (a) should be close to 2.0" + assert np.isclose(b, 1.5, atol=1e-2), "Semi-minor axis (b) should be close to 1.5" + assert np.isclose(c, 0.5, atol=1e-2), "Rotation parameter (c) should be close to 0.5" + assert np.isclose(u, 1.0, atol=1e-2), "Center along the x-axis (u) should be close to 1.0" + assert np.isclose(v, 1.0, atol=1e-2), "Center along the y-axis (v) should be close to 1.0" + + +def test_fit_ellipse_parameters_noisy() -> None: + """Test fit_ellipse_parameters with small Gaussian noise and 60 points.""" + ellipse = Ellipse(a=2.0, b=1.5, c=0.5, u=1.0, v=1.0, c_vu=np.ones(2)) + points = generate_points(ellipse, num_points=60, noise_std=0.02) + + a, b, c, u, v = fit_ellipse_parameters(points, rescale=True, use_l1_norm=False) + + assert np.isclose(a, 2.0, atol=1e-1), "Semi-major axis (a) should be close to 2.0" + assert np.isclose(b, 1.5, atol=1e-1), "Semi-minor axis (b) should be close to 1.5" + assert np.isclose(c, 0.5, atol=1e-1), "Rotation parameter (c) should be close to 0.5" + assert np.isclose(u, 1.0, atol=1e-1), "Center along the x-axis (u) should be close to 1.0" + assert np.isclose(v, 1.0, atol=1e-1), "Center along the y-axis (v) should be close to 1.0" + + +def test_fit_ellipse_parameters_noisy_with_outliers() -> None: + """Test fit_ellipse_parameters with 60 points and 2 outliers using L1 norm.""" + ellipse = Ellipse(a=2.0, b=1.5, c=0.5, u=1.0, v=1.0, c_vu=np.ones(2)) + points = generate_points(ellipse, num_points=60, noise_std=0.02, outliers=2) + + a, b, c, u, v = fit_ellipse_parameters(points, rescale=True, use_l1_norm=True) + + assert np.isclose(a, 2.0, atol=1e-1), "Semi-major axis (a) should be close to 2.0" + assert np.isclose(b, 1.5, atol=1e-1), "Semi-minor axis (b) should be close to 1.5" + assert np.isclose(c, 0.5, atol=1e-1), "Rotation parameter (c) should be close to 0.5" + assert np.isclose(u, 1.0, atol=1e-1), "Center along the x-axis (u) should be close to 1.0" + assert np.isclose(v, 1.0, atol=1e-1), "Center along the y-axis (v) should be close to 1.0"