diff --git a/doc/examples/use_cases/fabry_perot.py b/doc/examples/use_cases/fabry_perot.py
index 904246f6..58a1f9ed 100644
--- a/doc/examples/use_cases/fabry_perot.py
+++ b/doc/examples/use_cases/fabry_perot.py
@@ -116,11 +116,11 @@
# Configure line profile extraction
profile_param = sigima.proc.image.LineProfileParam()
profile_param.direction = "horizontal"
-profile_param.row = center_y # Extract profile through image center
+profile_param.y = center_y # Extract profile through image center (Y coordinate)
print("\n✓ Horizontal profile configured:")
print(f"Direction: {profile_param.direction}")
-print(f"Row: {profile_param.row} (image center)")
+print(f"Y: {profile_param.y} (image center)")
# Extract intensity profile
profile_signal1 = sigima.proc.image.line_profile(img1, profile_param)
diff --git a/doc/examples/use_cases/laser_beam.py b/doc/examples/use_cases/laser_beam.py
index c7e43360..a39287de 100644
--- a/doc/examples/use_cases/laser_beam.py
+++ b/doc/examples/use_cases/laser_beam.py
@@ -167,16 +167,17 @@ def load_laser_beam_images():
line_profiles = []
for i, (img, centroid_coords) in enumerate(zip(clipped_images, centroids)):
- # Create line profile parameters for horizontal line through centroid
+ # Create line profile parameters for horizontal line through centroid (the
+ # centroid and the profile coordinates are in the physical units of the image)
line_param = sigima.proc.image.LineProfileParam()
line_param.direction = "horizontal"
- line_param.row = int(centroid_coords[1]) # Use centroid y-coordinate as row
+ line_param.y = centroid_coords[1] # Use centroid y-coordinate
# Extract line profile
profile = sigima.proc.image.line_profile(img, line_param)
profile.title = f"Line_profile_{img.title}"
line_profiles.append(profile)
- print(f" ✓ Extracted line profile for {img.title} at row {line_param.row}")
+ print(f" ✓ Extracted line profile for {img.title} at y = {line_param.y:.1f}")
print(f"\n✓ Generated {len(line_profiles)} line profiles")
diff --git a/doc/locale/fr/LC_MESSAGES/release_notes/release_1.04.po b/doc/locale/fr/LC_MESSAGES/release_notes/release_1.04.po
index 65b58a38..872510a0 100644
--- a/doc/locale/fr/LC_MESSAGES/release_notes/release_1.04.po
+++ b/doc/locale/fr/LC_MESSAGES/release_notes/release_1.04.po
@@ -76,12 +76,18 @@ msgstr "**Centre de l'image** : `ImageObj.xc` et `ImageObj.yc` sont désormais l
msgid "**Ellipse coordinates**: `EllipseCoordinates` (`sigima.objects.shape`) now holds the 5 values `(x, y, a, b, theta)` of ellipse geometry results, which can thus be transformed."
msgstr "**Coordonnées d'ellipse** : `EllipseCoordinates` (`sigima.objects.shape`) contient désormais les 5 valeurs `(x, y, a, b, theta)` des résultats géométriques d'ellipse, qui peuvent ainsi être transformés."
-msgid "**Operation contract versions**: The operation contracts keep version 1 in this release, although the result of some operations changed since 1.3.0 (for example `resampling`, `extract_roi` and the product of signals): the contracts had just been published with Sigima 1.3.0 and were not used yet. This one-off exception does not apply to later releases: from now on, the version of a contract is increased whenever a change of the computation changes its result, whatever the releases of the applications using Sigima, and calls recorded with an earlier version are refused instead of being replayed (see [Operation contracts and versions](../user_guide/contracts.rst))."
-msgstr "**Versions des contrats d'opération** : les contrats d'opération conservent la version 1 dans cette version, bien que le résultat de certaines opérations ait changé depuis la 1.3.0 (par exemple `resampling`, `extract_roi` et le produit de signaux) : les contrats venaient d'être publiés avec Sigima 1.3.0 et n'étaient pas encore exploités. Cette exception ponctuelle ne s'applique pas aux versions ultérieures : désormais, la version d'un contrat est augmentée dès qu'une modification du calcul change son résultat, indépendamment des publications des applications utilisant Sigima, et les appels enregistrés avec une version antérieure sont refusés au lieu d'être rejoués (voir [Contrats d'opération et versions](../user_guide/contracts.rst))."
+msgid "**Operation contract versions**: The operation contracts keep version 1 in this release, although the result of some operations changed since 1.3.0 (for example `resampling`, `extract_roi` and the product of signals): the contracts had just been published with Sigima 1.3.0 and were not used yet. This one-off exception does not apply to later changes: from now on, the version of a contract is increased whenever a change of the computation changes its result, whatever the releases of the applications using Sigima, and calls recorded with an earlier version are refused instead of being replayed (see [Operation contracts and versions](../user_guide/contracts.rst)). The image profiles, changed afterwards (see below), are thus at version 2."
+msgstr "**Versions des contrats d'opération** : les contrats d'opération conservent la version 1 dans cette version, bien que le résultat de certaines opérations ait changé depuis la 1.3.0 (par exemple `resampling`, `extract_roi` et le produit de signaux) : les contrats venaient d'être publiés avec Sigima 1.3.0 et n'étaient pas encore exploités. Cette exception ponctuelle ne s'applique pas aux modifications ultérieures : désormais, la version d'un contrat est augmentée dès qu'une modification du calcul change son résultat, indépendamment des publications des applications utilisant Sigima, et les appels enregistrés avec une version antérieure sont refusés au lieu d'être rejoués (voir [Contrats d'opération et versions](../user_guide/contracts.rst)). Les profils d'image, modifiés ensuite (voir ci-dessous), sont donc en version 2."
+
+msgid "**Image profiles in physical coordinates**: The line, segment, average and radial profiles of images are now defined and returned in the physical coordinates of the image, like regions of interest and analysis results, instead of pixel indices. Their parameters were renamed, so that scripts using the former names fail explicitly instead of silently using another meaning: `LineProfileParam.row`/`col` became `y`/`x` (the row or column whose pixel centers are the nearest), `SegmentProfileParam.row1`/`col1`/`row2`/`col2` and `AverageProfileParam.row1`/`col1`/`row2`/`col2` became `x1`/`y1`/`x2`/`y2` (for the average profile, `None` stands for the edge of the image instead of `-1`, and the pixels are those whose centers lie in the rectangle, like rectangular ROIs), and `RadialProfileParam.x0`/`y0` became `xc`/`yc`. The X axis of the profiles holds physical coordinates or distances, with the unit and label of the corresponding axis of the image. Segment and radial profiles are sampled in steps of the smallest distance between pixel centers, so that the radial profile is circular for anisotropic pixels and non-uniform coordinates; they are refused when the X and Y units of the image differ. Coordinates outside the image are refused instead of being clamped. On images without calibration (origin 0, pixel size 1), results are unchanged. The contracts of `line_profile`, `segment_profile`, `average_profile` and `radial_profile` are at version 2."
+msgstr "**Profils d'image en coordonnées physiques** : les profils ligne, segment, moyen et radial des images sont désormais définis et renvoyés dans les coordonnées physiques de l'image, comme les régions d'intérêt et les résultats d'analyse, au lieu d'indices de pixels. Leurs paramètres ont été renommés, afin que les scripts utilisant les anciens noms échouent explicitement au lieu d'utiliser silencieusement un autre sens : `LineProfileParam.row`/`col` deviennent `y`/`x` (la ligne ou la colonne dont les centres de pixels sont les plus proches), `SegmentProfileParam.row1`/`col1`/`row2`/`col2` et `AverageProfileParam.row1`/`col1`/`row2`/`col2` deviennent `x1`/`y1`/`x2`/`y2` (pour le profil moyen, `None` désigne le bord de l'image au lieu de `-1`, et les pixels retenus sont ceux dont le centre est dans le rectangle, comme pour les ROI rectangulaires), et `RadialProfileParam.x0`/`y0` deviennent `xc`/`yc`. L'axe X des profils contient des coordonnées ou des distances physiques, avec l'unité et le libellé de l'axe correspondant de l'image. Les profils de segment et radiaux sont échantillonnés par pas égaux à la plus petite distance entre centres de pixels, de sorte que le profil radial est circulaire pour des pixels anisotropes et des coordonnées non uniformes ; ils sont refusés lorsque les unités X et Y de l'image diffèrent. Les coordonnées hors de l'image sont refusées au lieu d'être ramenées au bord. Sur les images sans calibration (origine 0, taille de pixel 1), les résultats sont inchangés. Les contrats de `line_profile`, `segment_profile`, `average_profile` et `radial_profile` sont en version 2."
msgid "🛠️ Bug Fixes since version 1.3.0"
msgstr "🛠️ Corrections de bugs depuis la version 1.3.0"
+msgid "**Image profiles**: The radial profile centered on the image center was computed around a wrong point on calibrated images (the physical center was used as pixel indices), and could even contain undefined values. Vertical profiles were labeled with the unit of the X axis instead of the Y axis, and all profiles with the unit of the image although their X axis was in pixels. Anti-diagonal segment profiles (e.g. from the top-right to the bottom-left corner) were extracted along the other diagonal."
+msgstr "**Profils d'image** : le profil radial centré sur le centre de l'image était calculé autour d'un mauvais point sur les images calibrées (le centre physique était utilisé comme indices de pixels), et pouvait même contenir des valeurs indéfinies. Les profils verticaux portaient l'unité de l'axe X au lieu de l'axe Y, et tous les profils l'unité de l'image alors que leur axe X était en pixels. Les profils de segment anti-diagonaux (par exemple du coin supérieur droit au coin inférieur gauche) étaient extraits le long de l'autre diagonale."
+
msgid "**Exact integer import**: Integer data that does not fit any integer type supported by the destination object (for example `uint32` counts above 2 147 483 647 in an image, or large integers in a signal) is now converted to a floating-point type that represents every value exactly. Such data was previously converted to `float32` when the source used 4-byte integers, which rounded counts above 16 777 216. A warning is now emitted when no supported floating-point type is exact (integers beyond 2⁵³)."
msgstr "**Import exact des entiers** : les données entières qui ne tiennent dans aucun type entier pris en charge par l'objet de destination (par exemple des comptages `uint32` supérieurs à 2 147 483 647 dans une image, ou de grands entiers dans un signal) sont désormais converties vers un type flottant représentant chaque valeur exactement. Ces données étaient auparavant converties en `float32` lorsque la source utilisait des entiers sur 4 octets, ce qui arrondissait les comptages au-delà de 16 777 216. Un avertissement est désormais émis lorsqu'aucun type flottant pris en charge n'est exact (entiers au-delà de 2⁵³)."
diff --git a/doc/release_notes/release_1.04.md b/doc/release_notes/release_1.04.md
index 305839be..0c5443d4 100644
--- a/doc/release_notes/release_1.04.md
+++ b/doc/release_notes/release_1.04.md
@@ -26,10 +26,12 @@ The [user guide](../user_guide/overview.rst) shows how readers and connectors at
* **Replay of processing in image ROIs**: The selection of image pixels by regions of interest is now recorded with each call on images carrying regions of interest (rule `sigima.image.roi.pixel_centers`, `sigima.proc.contracts.IMAGE_ROI_CONVENTION`, in the execution context described by `describe_call` and `prepare_execution`). Applications that replay recorded processing (e.g. DataLab capsules) refuse to replay such calls recorded with earlier versions, which selected other pixels, instead of reporting a different result. Calls on images without regions of interest are not concerned.
* **Image center**: `ImageObj.xc` and `ImageObj.yc` are now the center of the pixel grid (halfway between the centers of the first and last pixels), for uniform and non-uniform coordinates: they were half a pixel off for uniform coordinates. `width` and `height` are unchanged.
* **Ellipse coordinates**: `EllipseCoordinates` (`sigima.objects.shape`) now holds the 5 values `(x, y, a, b, theta)` of ellipse geometry results, which can thus be transformed.
-* **Operation contract versions**: The operation contracts keep version 1 in this release, although the result of some operations changed since 1.3.0 (for example `resampling`, `extract_roi` and the product of signals): the contracts had just been published with Sigima 1.3.0 and were not used yet. This one-off exception does not apply to later releases: from now on, the version of a contract is increased whenever a change of the computation changes its result, whatever the releases of the applications using Sigima, and calls recorded with an earlier version are refused instead of being replayed (see [Operation contracts and versions](../user_guide/contracts.rst)).
+* **Operation contract versions**: The operation contracts keep version 1 in this release, although the result of some operations changed since 1.3.0 (for example `resampling`, `extract_roi` and the product of signals): the contracts had just been published with Sigima 1.3.0 and were not used yet. This one-off exception does not apply to later changes: from now on, the version of a contract is increased whenever a change of the computation changes its result, whatever the releases of the applications using Sigima, and calls recorded with an earlier version are refused instead of being replayed (see [Operation contracts and versions](../user_guide/contracts.rst)). The image profiles, changed afterwards (see below), are thus at version 2.
+* **Image profiles in physical coordinates**: The line, segment, average and radial profiles of images are now defined and returned in the physical coordinates of the image, like regions of interest and analysis results, instead of pixel indices. Their parameters were renamed, so that scripts using the former names fail explicitly instead of silently using another meaning: `LineProfileParam.row`/`col` became `y`/`x` (the row or column whose pixel centers are the nearest), `SegmentProfileParam.row1`/`col1`/`row2`/`col2` and `AverageProfileParam.row1`/`col1`/`row2`/`col2` became `x1`/`y1`/`x2`/`y2` (for the average profile, `None` stands for the edge of the image instead of `-1`, and the pixels are those whose centers lie in the rectangle, like rectangular ROIs), and `RadialProfileParam.x0`/`y0` became `xc`/`yc`. The X axis of the profiles holds physical coordinates or distances, with the unit and label of the corresponding axis of the image. Segment and radial profiles are sampled in steps of the smallest distance between pixel centers, so that the radial profile is circular for anisotropic pixels and non-uniform coordinates; they are refused when the X and Y units of the image differ. Coordinates outside the image are refused instead of being clamped. On images without calibration (origin 0, pixel size 1), results are unchanged. The contracts of `line_profile`, `segment_profile`, `average_profile` and `radial_profile` are at version 2.
### 🛠️ Bug Fixes since version 1.3.0 ###
+* **Image profiles**: The radial profile centered on the image center was computed around a wrong point on calibrated images (the physical center was used as pixel indices), and could even contain undefined values. Vertical profiles were labeled with the unit of the X axis instead of the Y axis, and all profiles with the unit of the image although their X axis was in pixels. Anti-diagonal segment profiles (e.g. from the top-right to the bottom-left corner) were extracted along the other diagonal.
* **Exact integer import**: Integer data that does not fit any integer type supported by the destination object (for example `uint32` counts above 2 147 483 647 in an image, or large integers in a signal) is now converted to a floating-point type that represents every value exactly. Such data was previously converted to `float32` when the source used 4-byte integers, which rounded counts above 16 777 216. A warning is now emitted when no supported floating-point type is exact (integers beyond 2⁵³).
* **Image resampling**: Resampling an image onto its own grid now gives back the same image. The output pixels were sampled half a pixel away from the requested positions, so that every value was shifted and the last row and column were filled with NaN, while the output origin claimed the requested position. The `xmin`/`xmax` and `ymin`/`ymax` parameters of `resampling` are now the coordinates of the centers of the first and last output pixels, like the `xmin`/`xmax` of an image, and default to the source grid. In pixel size mode, the requested pixel size is now kept exactly (the last pixel is the last one not beyond `xmax`/`ymax`).
* **Real and complex signals**: The sum, average, product and standard deviation of several signals, and their conversion into an image, no longer depend on the order of the signals: when a real signal came first, the imaginary parts of the complex signals were discarded.
diff --git a/sigima/locale/fr/LC_MESSAGES/sigima.po b/sigima/locale/fr/LC_MESSAGES/sigima.po
index 64ca2857..f9ba69ed 100644
--- a/sigima/locale/fr/LC_MESSAGES/sigima.po
+++ b/sigima/locale/fr/LC_MESSAGES/sigima.po
@@ -1812,3 +1812,17 @@ msgstr "Barycentre"
msgid "Plot dialog"
msgstr "Fenêtre de tracé"
+msgid "Y coordinate of the row (nearest pixel center)"
+msgstr "Coordonnée Y de la ligne (centre de pixel le plus proche)"
+
+msgid "X coordinate of the column (nearest pixel center)"
+msgstr "Coordonnée X de la colonne (centre de pixel le plus proche)"
+
+msgid "Start point"
+msgstr "Point de départ"
+
+msgid "End point"
+msgstr "Point d'arrivée"
+
+msgid "None: edge of the image"
+msgstr "Aucune valeur : bord de l'image"
diff --git a/sigima/params.py b/sigima/params.py
index 0488ed66..2922a149 100644
--- a/sigima/params.py
+++ b/sigima/params.py
@@ -80,7 +80,9 @@
- :class:`Resampling2DParam`: Updates bounds based on signal range
- :class:`ResizeParam`: Updates bounds based on image dimensions
- :class:`TranslateParam`: Updates bounds based on image dimensions
-- :class:`LineProfileParam`: Updates line coordinates based on image dimensions
+- :class:`LineProfileParam`, :class:`SegmentProfileParam`,
+ :class:`RadialProfileParam`: Set profile coordinates inside the image (physical
+ units)
- :class:`BandPassFilterParam`, :class:`BandStopFilterParam`,
:class:`HighPassFilterParam`, :class:`LowPassFilterParam`: Update frequency bounds
diff --git a/sigima/proc/image/extraction.py b/sigima/proc/image/extraction.py
index 48b18657..3484ca08 100644
--- a/sigima/proc/image/extraction.py
+++ b/sigima/proc/image/extraction.py
@@ -38,9 +38,10 @@
from sigima.objects.image import ImageObj, ImageROI, RectangularROI, ROI2DParam
from sigima.objects.image.object import PIXEL_INDEX_TOLERANCE
from sigima.objects.signal import SignalObj
+from sigima.proc.base import new_signal_result
from sigima.proc.contracts import EXACT, ReplayQualification
from sigima.proc.decorator import computation_function
-from sigima.proc.image.base import dst_1_to_1, dst_1_to_1_signal
+from sigima.proc.image.base import dst_1_to_1
# NOTE: Only parameter classes DEFINED in this module should be included in __all__.
# Parameter classes imported from other modules (like sigima.proc.base) should NOT
@@ -325,170 +326,328 @@ def generate_image_grid_roi(src: ImageObj, p: ROIGridParam) -> ImageROI:
return roi
+def _pixel_centers(src: ImageObj, axis: str) -> np.ndarray:
+ """Return the physical coordinates of the pixel centers along an axis.
+
+ Args:
+ src: image object
+ axis: ``"x"`` (columns) or ``"y"`` (rows)
+
+ Returns:
+ Coordinates of the pixel centers
+ """
+ size = src.data.shape[1 if axis == "x" else 0]
+ return src.axis_to_physical(np.arange(size, dtype=float), axis)
+
+
+def _nearest_pixels(
+ src: ImageObj, values: np.ndarray | list[float], axis: str, name: str
+) -> np.ndarray:
+ """Return the indices of the pixels whose centers are the nearest to coordinates.
+
+ Args:
+ src: image object
+ values: physical coordinates along the axis
+ axis: ``"x"`` (columns) or ``"y"`` (rows)
+ name: name of the coordinates, for error messages
+
+ Returns:
+ Pixel indices (ties go to the next pixel)
+
+ Raises:
+ ValueError: if a coordinate lies outside the image (beyond the outer edges
+ of its pixels)
+ """
+ size = src.data.shape[1 if axis == "x" else 0]
+ findices = np.atleast_1d(src.axis_to_float_indices(values, axis))
+ outside = (findices < -0.5 - PIXEL_INDEX_TOLERANCE) | (
+ findices > size - 0.5 + PIXEL_INDEX_TOLERANCE
+ )
+ if np.any(outside):
+ raise ValueError(
+ f"{name} = {np.atleast_1d(values)[outside][0]:g} lies outside the image"
+ )
+ return np.clip(np.floor(findices + 0.5), 0, size - 1).astype(int)
+
+
+def _pixel_spacing(src: ImageObj) -> float:
+ """Return the smallest distance between the centers of neighboring pixels.
+
+ Args:
+ src: image object
+
+ Returns:
+ Smallest pixel spacing along X and Y (uniform or non-uniform coordinates)
+ """
+ spacings = []
+ for axis in ("x", "y"):
+ centers = _pixel_centers(src, axis)
+ if centers.size > 1:
+ spacings.append(float(np.min(np.abs(np.diff(centers)))))
+ if not spacings:
+ # Single pixel: any positive step gives the same profile
+ spacings = [abs(float(src.dx)) if src.is_uniform_coords else 1.0]
+ return min(spacings)
+
+
+def _check_same_units(src: ImageObj, name: str) -> None:
+ """Check that distances can be computed in the image plane.
+
+ Args:
+ src: image object
+ name: name of the computation, for error messages
+
+ Raises:
+ ValueError: if the X and Y units of the image differ
+ """
+ if (src.xunit or "") != (src.yunit or ""):
+ raise ValueError(
+ f"{name}: the X and Y units of the image differ "
+ f"({src.xunit!r} and {src.yunit!r}), so that distances are undefined"
+ )
+
+
+def _axis_units(src: ImageObj, axis: str) -> tuple[tuple[str, str], tuple[str, str]]:
+ """Return the units and labels of a profile along an axis of the image."""
+ if axis == "x":
+ return (src.xunit, src.zunit), (src.xlabel, src.zlabel)
+ return (src.yunit, src.zunit), (src.ylabel, src.zlabel)
+
+
class LineProfileParam(gds.DataSet):
- """Horizontal or vertical profile parameters"""
+ """Horizontal or vertical profile parameters
+
+ The profile is extracted along the row (or column) of pixels whose center is
+ the nearest to the given coordinate, in the physical units of the image.
+ """
_prop = gds.GetAttrProp("direction")
_directions = (("horizontal", _("horizontal")), ("vertical", _("vertical")))
direction = gds.ChoiceItem(_("Direction"), _directions, radio=True).set_prop(
"display", store=_prop
)
- row = gds.IntItem(_("Row"), default=0, min=0).set_prop(
- "display", active=gds.FuncProp(_prop, lambda x: x == "horizontal")
- )
- col = gds.IntItem(_("Column"), default=0, min=0).set_prop(
- "display", active=gds.FuncProp(_prop, lambda x: x == "vertical")
- )
+ y = gds.FloatItem(
+ "Y", default=0.0, help=_("Y coordinate of the row (nearest pixel center)")
+ ).set_prop("display", active=gds.FuncProp(_prop, lambda x: x == "horizontal"))
+ x = gds.FloatItem(
+ "X", default=0.0, help=_("X coordinate of the column (nearest pixel center)")
+ ).set_prop("display", active=gds.FuncProp(_prop, lambda x: x == "vertical"))
+
+ def update_from_obj(self, obj: ImageObj) -> None:
+ """Update parameters from image: profiles through the center of the image"""
+ self.x, self.y = obj.xc, obj.yc
@computation_function(
operation_id="sigima.image.line_profile",
- contract_version=1,
- qualification=ReplayQualification(contract_version=1, comparison=EXACT),
+ contract_version=2,
+ qualification=ReplayQualification(contract_version=2, comparison=EXACT),
)
def line_profile(src: ImageObj, p: LineProfileParam) -> SignalObj:
"""Compute horizontal or vertical profile
+ The row (or column) is the one whose pixel centers are the nearest to the
+ coordinate given in the physical units of the image. The X axis of the profile
+ holds the coordinates of the pixel centers, with the unit and label of the
+ corresponding axis of the image.
+
Args:
src: input image object
p: parameters
Returns:
Signal object with the profile
+
+ Raises:
+ ValueError: if the coordinate lies outside the image
"""
data = src.get_masked_view()
- row = min(p.row, data.shape[0] - 1)
- col = min(p.col, data.shape[1] - 1)
if p.direction == "horizontal":
- suffix, shape_index, pdata = f"row={row}", 1, data[row, :]
+ (row,) = _nearest_pixels(src, [p.y], "y", "Y")
+ suffix, axis, pdata = f"y={p.y:g}", "x", data[row, :]
else:
- suffix, shape_index, pdata = f"col={col}", 0, data[:, col]
- pdata: ma.MaskedArray
- x = np.arange(data.shape[shape_index])[~pdata.mask]
- y = np.array(pdata, dtype=float)[~pdata.mask]
- dst = dst_1_to_1_signal(src, "profile", suffix)
+ (col,) = _nearest_pixels(src, [p.x], "x", "X")
+ suffix, axis, pdata = f"x={p.x:g}", "y", data[:, col]
+ keep = ~ma.getmaskarray(pdata)
+ x = _pixel_centers(src, axis)[keep]
+ y = np.array(pdata, dtype=float)[keep]
+ units, labels = _axis_units(src, axis)
+ dst = new_signal_result(src, "profile", suffix, units, labels)
dst.set_xydata(x, y)
return dst
class SegmentProfileParam(gds.DataSet):
- """Segment profile parameters"""
+ """Segment profile parameters
+
+ The segment goes from (X1, Y1) to (X2, Y2), in the physical units of the image.
+ """
+
+ x1 = gds.FloatItem("X1", default=0.0, help=_("Start point"))
+ y1 = gds.FloatItem("Y1", default=0.0, help=_("Start point")).set_pos(
+ col=1
+ )
+ x2 = gds.FloatItem("X2", default=0.0, help=_("End point"))
+ y2 = gds.FloatItem("Y2", default=0.0, help=_("End point")).set_pos(col=1)
- row1 = gds.IntItem(_("Start row"), default=0, min=0)
- col1 = gds.IntItem(_("Start column"), default=0, min=0)
- row2 = gds.IntItem(_("End row"), default=0, min=0)
- col2 = gds.IntItem(_("End column"), default=0, min=0)
+ def update_from_obj(self, obj: ImageObj) -> None:
+ """Update parameters from image: segment along the middle row"""
+ x_centers = _pixel_centers(obj, "x")
+ self.x1, self.x2 = float(x_centers[0]), float(x_centers[-1])
+ self.y1 = self.y2 = obj.yc
def csline(data: np.ndarray, row0, col0, row1, col1) -> tuple[np.ndarray, np.ndarray]:
- """Return intensity profile of data along a line
+ """Return intensity profile of data along a line, in pixel indices
Args:
data: 2D array
row0, col0: start point
row1, col1: end point
+
+ Returns:
+ Sample indices and values along the line, from the start point to the end
+ point (nearest pixels, one sample per pixel along the longest direction)
"""
# Keep coordinates inside the image
row0 = max(0, min(row0, data.shape[0] - 1))
col0 = max(0, min(col0, data.shape[1] - 1))
row1 = max(0, min(row1, data.shape[0] - 1))
col1 = max(0, min(col1, data.shape[1] - 1))
- # Keep coordinates in the right order
- row0, row1 = min(row0, row1), max(row0, row1)
- col0, col1 = min(col0, col1), max(col0, col1)
- # Extract the line
- line = np.zeros((2, max(abs(row1 - row0), abs(col1 - col0)) + 1), dtype=int)
- line[0, :] = np.linspace(row0, row1, line.shape[1]).astype(int)
- line[1, :] = np.linspace(col0, col1, line.shape[1]).astype(int)
- # Interpolate the line
- y = np.ma.array(data[line[0], line[1]], float).filled(np.nan)
+ # Sample from the start point to the end point (rows and columns are not
+ # reordered independently, which would mirror anti-diagonal segments)
+ count = max(abs(row1 - row0), abs(col1 - col0)) + 1
+ rows = np.rint(np.linspace(row0, row1, count)).astype(int)
+ cols = np.rint(np.linspace(col0, col1, count)).astype(int)
+ y = np.ma.array(data[rows, cols], float).filled(np.nan)
x = np.arange(y.size)
return x, y
@computation_function(
operation_id="sigima.image.segment_profile",
- contract_version=1,
- qualification=ReplayQualification(contract_version=1, comparison=EXACT),
+ contract_version=2,
+ qualification=ReplayQualification(contract_version=2, comparison=EXACT),
)
def segment_profile(src: ImageObj, p: SegmentProfileParam) -> SignalObj:
"""Compute segment profile
+ The segment is sampled at regular steps, equal to the smallest distance between
+ the centers of neighboring pixels, from (X1, Y1) to (X2, Y2): each sample takes
+ the value of the pixel whose center is the nearest (no interpolation). The X
+ axis of the profile is the distance from the start point, in the units of the
+ image.
+
Args:
src: input image object
p: parameters
Returns:
Signal object with the segment profile
+
+ Raises:
+ ValueError: if an end point lies outside the image, or if the X and Y units
+ of the image differ
"""
- data = src.get_masked_view()
- row1 = min(p.row1, data.shape[0] - 1)
- col1 = min(p.col1, data.shape[1] - 1)
- row2 = min(p.row2, data.shape[0] - 1)
- col2 = min(p.col2, data.shape[1] - 1)
- suffix = f"({row1}, {col1})-({row2}, {col2})"
- x, y = csline(data, row1, col1, row2, col2)
- x, y = x[~np.isnan(y)], y[~np.isnan(y)] # Remove NaN values
- dst = dst_1_to_1_signal(src, "segment_profile", suffix)
- dst.set_xydata(np.array(x, dtype=float), np.array(y, dtype=float))
+ _check_same_units(src, "Segment profile")
+ _nearest_pixels(src, [p.x1, p.x2], "x", "X")
+ _nearest_pixels(src, [p.y1, p.y2], "y", "Y")
+ length = float(np.hypot(p.x2 - p.x1, p.y2 - p.y1))
+ step = _pixel_spacing(src)
+ count = int(np.floor(length / step + PIXEL_INDEX_TOLERANCE)) + 1
+ dist = np.arange(count) * step
+ fraction = dist / length if length > 0 else np.zeros(count)
+ cols = _nearest_pixels(src, p.x1 + (p.x2 - p.x1) * fraction, "x", "X")
+ rows = _nearest_pixels(src, p.y1 + (p.y2 - p.y1) * fraction, "y", "Y")
+ values = src.get_masked_view()[rows, cols]
+ keep = ~ma.getmaskarray(values)
+ suffix = f"({p.x1:g}, {p.y1:g})-({p.x2:g}, {p.y2:g})"
+ units, labels = _axis_units(src, "x")
+ dst = new_signal_result(src, "segment_profile", suffix, units, labels)
+ dst.set_xydata(dist[keep], np.array(values, dtype=float)[keep])
return dst
class AverageProfileParam(gds.DataSet):
- """Average horizontal or vertical profile parameters"""
+ """Average horizontal or vertical profile parameters
+
+ The rectangular area goes from (X1, Y1) to (X2, Y2), in the physical units of
+ the image (None: edge of the image). It holds the pixels whose centers lie
+ inside it, bounds included.
+ """
_directions = (("horizontal", _("horizontal")), ("vertical", _("vertical")))
direction = gds.ChoiceItem(_("Direction"), _directions, radio=True)
_hgroup_begin = gds.BeginGroup(_("Profile rectangular area"))
- row1 = gds.IntItem(_("Row 1"), default=0, min=0)
- row2 = gds.IntItem(_("Row 2"), default=-1, min=-1)
- col1 = gds.IntItem(_("Column 1"), default=0, min=0)
- col2 = gds.IntItem(_("Column 2"), default=-1, min=-1)
+ _none_help = _("None: edge of the image")
+ x1 = gds.FloatItem("X1", default=None, allow_none=True, help=_none_help)
+ y1 = gds.FloatItem(
+ "Y1", default=None, allow_none=True, help=_none_help
+ ).set_pos(col=1)
+ x2 = gds.FloatItem("X2", default=None, allow_none=True, help=_none_help)
+ y2 = gds.FloatItem(
+ "Y2", default=None, allow_none=True, help=_none_help
+ ).set_pos(col=1)
_hgroup_end = gds.EndGroup(_("Profile rectangular area"))
@computation_function(
operation_id="sigima.image.average_profile",
- contract_version=1,
- qualification=ReplayQualification(contract_version=1, comparison=EXACT),
+ contract_version=2,
+ qualification=ReplayQualification(contract_version=2, comparison=EXACT),
)
def average_profile(src: ImageObj, p: AverageProfileParam) -> SignalObj:
"""Compute horizontal or vertical average profile
+ The pixels whose centers lie in the rectangle (bounds included, like
+ rectangular ROIs) are averaged along the columns (horizontal profile) or the
+ rows (vertical profile). The X axis of the profile holds the coordinates of the
+ pixel centers, with the unit and label of the corresponding axis of the image.
+
Args:
src: input image object
p: parameters
Returns:
Signal object with the average profile
+
+ Raises:
+ ValueError: if no pixel center lies in the rectangle
"""
data = src.get_masked_view()
- row1, row2, col1, col2 = p.row1, p.row2, p.col1, p.col2
- if row2 == -1:
- row2 = data.shape[0] - 1
- if col2 == -1:
- col2 = data.shape[1] - 1
- if row1 > row2:
- row1, row2 = row2, row1
- if col1 > col2:
- col1, col2 = col2, col1
- row1 = min(row1, data.shape[0] - 1)
- row2 = min(row2, data.shape[0] - 1)
- col1 = min(col1, data.shape[1] - 1)
- col2 = min(col2, data.shape[1] - 1)
- suffix = f"{p.direction}, rows=[{row1}, {row2}], cols=[{col1}, {col2}]"
+ nrows, ncols = data.shape
+ left, top, right, bottom = src.pixel_range_to_physical([0, 0, ncols, nrows])
+ x1 = left if p.x1 is None else p.x1
+ y1 = top if p.y1 is None else p.y1
+ x2 = right if p.x2 is None else p.x2
+ y2 = bottom if p.y2 is None else p.y2
+ ix0, iy0, ix1, iy1 = src.physical_to_pixel_range([x1, y1, x2, y2])
+ ix0, iy0 = max(0, ix0), max(0, iy0)
+ ix1, iy1 = min(ncols, ix1), min(nrows, iy1)
+ if ix1 <= ix0 or iy1 <= iy0:
+ raise ValueError(
+ f"Average profile: no pixel center lies in the rectangle "
+ f"({x1:g}, {y1:g})-({x2:g}, {y2:g})"
+ )
+ suffix = f"{p.direction}, x=[{x1:g}, {x2:g}], y=[{y1:g}, {y2:g}]"
+ area = data[iy0:iy1, ix0:ix1]
if p.direction == "horizontal":
- x, axis = np.arange(col1, col2 + 1), 0
+ axis = "x"
+ x, y = _pixel_centers(src, "x")[ix0:ix1], ma.mean(area, axis=0)
else:
- x, axis = np.arange(row1, row2 + 1), 1
- y = ma.mean(data[row1 : row2 + 1, col1 : col2 + 1], axis=axis)
- dst = dst_1_to_1_signal(src, "average_profile", suffix)
+ axis = "y"
+ x, y = _pixel_centers(src, "y")[iy0:iy1], ma.mean(area, axis=1)
+ units, labels = _axis_units(src, axis)
+ dst = new_signal_result(src, "average_profile", suffix, units, labels)
dst.set_xydata(x, y)
return dst
class RadialProfileParam(gds.DataSet):
- """Radial profile parameters"""
+ """Radial profile parameters
+
+ The center is given in the physical units of the image.
+ """
def __init__(self, *args, **kwargs) -> None:
super().__init__(*args, **kwargs)
@@ -497,19 +656,17 @@ def __init__(self, *args, **kwargs) -> None:
def update_from_obj(self, obj: ImageObj) -> None:
"""Update parameters from image"""
self.__obj = obj
- self.x0 = obj.xc
- self.y0 = obj.yc
+ self.xc = obj.xc
+ self.yc = obj.yc
def choice_callback(self, item, value): # pylint: disable=unused-argument
"""Callback for choice item"""
if self.__obj is None:
return
if value == "centroid":
- self.y0, self.x0 = sigima.tools.image.get_centroid_fourier(
- self.__obj.get_masked_view()
- )
+ self.xc, self.yc = _centroid(self.__obj, self.__obj.get_masked_view())
elif value == "center":
- self.x0, self.y0 = self.__obj.xc, self.__obj.yc
+ self.xc, self.yc = self.__obj.xc, self.__obj.yc
_prop = gds.GetAttrProp("center")
center = gds.ChoiceItem(
@@ -524,39 +681,61 @@ def choice_callback(self, item, value): # pylint: disable=unused-argument
_func_prop = gds.FuncProp(_prop, lambda x: x == "user")
_xyl = "" + _("Center") + ""
- x0 = gds.FloatItem(f"X{_xyl}", default=0.0, unit="pixel").set_prop(
- "display", active=_func_prop
- )
- y0 = gds.FloatItem(f"Y{_xyl}", default=0.0, unit="pixel").set_prop(
- "display", active=_func_prop
+ xc = gds.FloatItem(f"X{_xyl}", default=0.0).set_prop("display", active=_func_prop)
+ yc = gds.FloatItem(f"Y{_xyl}", default=0.0).set_prop("display", active=_func_prop)
+
+
+def _centroid(src: ImageObj, data: np.ndarray) -> tuple[float, float]:
+ """Return the centroid of the image, in physical coordinates."""
+ yi, xi = sigima.tools.image.get_centroid_fourier(data)
+ return (
+ float(src.axis_to_physical([xi], "x")[0]),
+ float(src.axis_to_physical([yi], "y")[0]),
)
@computation_function(
operation_id="sigima.image.radial_profile",
- contract_version=1,
- qualification=ReplayQualification(contract_version=1, comparison=EXACT),
+ contract_version=2,
+ qualification=ReplayQualification(contract_version=2, comparison=EXACT),
)
def radial_profile(src: ImageObj, p: RadialProfileParam) -> SignalObj:
- """Compute radial profile around the centroid
+ """Compute radial profile around the centroid, the image center or a given point
with :py:func:`sigima.tools.image.get_radial_profile`
+ The center and the distances are physical: the pixels are binned by their
+ distance to the center, in steps equal to the smallest distance between the
+ centers of neighboring pixels, so that the profile is circular whatever the
+ pixel sizes and for non-uniform coordinates. The X axis of the profile is the
+ distance to the center (mirrored), in the units of the image.
+
Args:
src: input image object
p: parameters
Returns:
Signal object with the radial profile
+
+ Raises:
+ ValueError: if the X and Y units of the image differ
"""
+ _check_same_units(src, "Radial profile")
data = src.get_masked_view()
if p.center == "centroid":
- y0, x0 = sigima.tools.image.get_centroid_fourier(data)
+ xc, yc = _centroid(src, data)
elif p.center == "center":
- x0, y0 = src.xc, src.yc
+ xc, yc = src.xc, src.yc
else:
- x0, y0 = p.x0, p.y0
- suffix = f"center=({x0:.3f}, {y0:.3f})"
- dst = dst_1_to_1_signal(src, "radial_profile", suffix)
- x, y = sigima.tools.image.get_radial_profile(data, (x0, y0))
+ xc, yc = p.xc, p.yc
+ suffix = f"center=({xc:.3f}, {yc:.3f})"
+ units, labels = _axis_units(src, "x")
+ dst = new_signal_result(src, "radial_profile", suffix, units, labels)
+ x, y = sigima.tools.image.get_radial_profile(
+ data,
+ (xc, yc),
+ xcoords=_pixel_centers(src, "x"),
+ ycoords=_pixel_centers(src, "y"),
+ step=_pixel_spacing(src),
+ )
dst.set_xydata(x, y)
return dst
diff --git a/sigima/tests/image/proc_extraction_extra_unit_test.py b/sigima/tests/image/proc_extraction_extra_unit_test.py
index cacb7adb..6bd45b03 100644
--- a/sigima/tests/image/proc_extraction_extra_unit_test.py
+++ b/sigima/tests/image/proc_extraction_extra_unit_test.py
@@ -5,9 +5,8 @@
Covers:
-* ``average_profile`` boundary handling: swapped row/col bounds and
- ``-1`` default sentinel values for both vertical and horizontal
- directions.
+* ``average_profile`` boundary handling: swapped corners and default
+ (``None``) bounds for both vertical and horizontal directions.
* ``line_profile`` vertical / horizontal direction branches.
* ``segment_profile`` basic execution.
"""
@@ -42,45 +41,40 @@ def _make_image(shape: tuple[int, int] = (8, 10)) -> ImageObj:
@pytest.mark.parametrize("direction", ["horizontal", "vertical"])
def test_average_profile_swapped_bounds(direction: str) -> None:
- """``average_profile`` must transparently swap inverted bounds (e.g.
- ``row1=5, row2=2``) so the user can pass them in either order without
- getting an empty profile, regardless of direction."""
+ """``average_profile`` accepts the corners of its rectangle in any order, so
+ that the user can pass them in either order without getting an empty profile,
+ regardless of direction."""
img = _make_image(shape=(8, 10))
- p = AverageProfileParam.create(direction=direction, row1=5, row2=2, col1=8, col2=1)
+ p = AverageProfileParam.create(direction=direction, x1=8, y1=5, x2=1, y2=2)
sig = average_profile(img, p)
- ordered = AverageProfileParam.create(
- direction=direction, row1=2, row2=5, col1=1, col2=8
- )
+ ordered = AverageProfileParam.create(direction=direction, x1=1, y1=2, x2=8, y2=5)
expected = average_profile(img, ordered)
assert np.array_equal(sig.x, expected.x) and np.array_equal(sig.y, expected.y)
# The parameters are not modified: they apply unchanged to other images
- assert (p.row1, p.row2, p.col1, p.col2) == (5, 2, 8, 1)
+ assert (p.x1, p.y1, p.x2, p.y2) == (8, 5, 1, 2)
@pytest.mark.parametrize(
"direction, expected_axis",
[("horizontal", 1), ("vertical", 0)],
)
-def test_average_profile_default_bounds_minus_one(
+def test_average_profile_default_bounds_none(
direction: str, expected_axis: int
) -> None:
- """The sentinel value ``-1`` for ``row2`` / ``col2`` is resolved to the
- last valid index of the image, covering the whole array by default,
- for both directions."""
+ """Bounds left to None are the edges of the image, covering the whole array by
+ default, for both directions."""
img = _make_image(shape=(4, 5))
- p = AverageProfileParam.create(
- direction=direction, row1=0, row2=-1, col1=0, col2=-1
- )
+ p = AverageProfileParam.create(direction=direction)
sig = average_profile(img, p)
- assert (p.row2, p.col2) == (-1, -1)
+ assert (p.x1, p.y1, p.x2, p.y2) == (None, None, None, None)
assert sig.y.size == img.data.shape[expected_axis]
@pytest.mark.parametrize(
"direction, kwargs, expected_axis",
[
- ("vertical", {"col": 2}, 0),
- ("horizontal", {"row": 3}, 1),
+ ("vertical", {"x": 2}, 0),
+ ("horizontal", {"y": 3}, 1),
],
)
def test_line_profile_directions(
@@ -98,7 +92,7 @@ def test_segment_profile_basic() -> None:
"""Smoke test: ``segment_profile`` returns a non-empty signal for an
arbitrary diagonal segment inside the image."""
img = _make_image()
- p = SegmentProfileParam.create(row1=0, col1=0, row2=5, col2=5)
+ p = SegmentProfileParam.create(x1=0, y1=0, x2=5, y2=5)
sig = segment_profile(img, p)
assert sig.y.size > 0
@@ -111,19 +105,18 @@ def _values(p) -> dict:
@pytest.mark.parametrize(
"func, p",
[
- (line_profile, LineProfileParam.create(direction="horizontal", row=6)),
- (line_profile, LineProfileParam.create(direction="vertical", col=7)),
- (segment_profile, SegmentProfileParam.create(row1=1, col1=1, row2=7, col2=9)),
+ (line_profile, LineProfileParam.create(direction="horizontal", y=3)),
+ (line_profile, LineProfileParam.create(direction="vertical", x=4)),
+ (segment_profile, SegmentProfileParam.create(x1=1, y1=1, x2=4, y2=3)),
(
average_profile,
- AverageProfileParam.create(
- direction="horizontal", row1=6, row2=2, col1=0, col2=-1
- ),
+ AverageProfileParam.create(direction="horizontal", y1=3, y2=1),
),
],
)
def test_profile_parameters_apply_unchanged_to_other_images(func, p) -> None:
- """Coordinates clamped to a small image are not clamped for the next one."""
+ """Parameters are not modified by a computation: they apply unchanged to
+ another image."""
small, large = _make_image(shape=(4, 5)), _make_image(shape=(8, 10))
before = _values(p)
expected = func(large, copy.deepcopy(p))
diff --git a/sigima/tests/image/profile_physical_unit_test.py b/sigima/tests/image/profile_physical_unit_test.py
new file mode 100644
index 00000000..d361c654
--- /dev/null
+++ b/sigima/tests/image/profile_physical_unit_test.py
@@ -0,0 +1,109 @@
+# Copyright (c) DataLab Platform Developers, BSD 3-Clause license, see LICENSE file.
+
+"""
+Image profiles in physical coordinates: anisotropic pixels, non-uniform
+coordinates, refusals and renamed parameters.
+"""
+
+# pylint: disable=invalid-name # Allows short reference names like x, y, ...
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+import sigima.objects
+import sigima.params
+import sigima.proc.image as sipi
+
+
+def _image(data: np.ndarray) -> sigima.objects.ImageObj:
+ obj = sigima.objects.create_image("i", data)
+ obj.xunit = obj.yunit = "mm"
+ return obj
+
+
+def test_radial_profile_anisotropic_pixels() -> None:
+ """With dx != dy, distances are physical: a disk gives a step profile."""
+ xc, yc, radius = 0.0, 0.0, 3.0
+ xcoords = np.arange(-40, 41) * 0.1 # dx = 0.1
+ ycoords = np.arange(-20, 21) * 0.2 # dy = 0.2
+ xx, yy = np.meshgrid(xcoords, ycoords)
+ disk = (np.hypot(xx - xc, yy - yc) <= radius).astype(float)
+ obj = _image(disk)
+ obj.set_uniform_coords(0.1, 0.2, xcoords[0], ycoords[0])
+ param = sigima.params.RadialProfileParam.create(center="user", xc=xc, yc=yc)
+ sig = sipi.radial_profile(obj, param)
+ assert np.allclose(np.diff(sig.x), 0.1) # Step: smallest pixel spacing
+ half = sig.x >= 0
+ assert np.all(sig.y[half][sig.x[half] < radius - 0.1] == 1.0)
+ assert np.all(sig.y[half][sig.x[half] > radius + 0.1] == 0.0)
+
+
+def test_profiles_non_uniform_coordinates() -> None:
+ """Profiles of images with non-uniform coordinates use the pixel centers."""
+ data = np.arange(20.0).reshape(4, 5)
+ obj = _image(data)
+ xcoords = np.array([0.0, 1.0, 3.0, 6.0, 10.0])
+ ycoords = np.array([0.0, 0.5, 1.5, 3.0])
+ obj.set_coords(xcoords, ycoords)
+ # Line profile: nearest row center to y = 1.2 is row 2 (1.5)
+ sig = sipi.line_profile(obj, sigima.params.LineProfileParam.create(y=1.2))
+ assert np.array_equal(sig.x, xcoords) and np.array_equal(sig.y, data[2])
+ # Vertical: nearest column center to x = 5 is column 3 (6.0)
+ sig = sipi.line_profile(
+ obj, sigima.params.LineProfileParam.create(direction="vertical", x=5.0)
+ )
+ assert np.array_equal(sig.x, ycoords) and np.array_equal(sig.y, data[:, 3])
+ # Average profile: centers in [1, 6] x [0.5, 3]
+ param = sigima.params.AverageProfileParam.create(x1=1, y1=0.5, x2=6, y2=3)
+ sig = sipi.average_profile(obj, param)
+ assert np.array_equal(sig.x, xcoords[1:4])
+ assert np.allclose(sig.y, data[1:4, 1:4].mean(axis=0))
+ # Radial profile: steps of the smallest spacing (0.5); distances holding no
+ # pixel are left out (no NaN)
+ param = sigima.params.RadialProfileParam.create(center="user", xc=0, yc=0)
+ sig = sipi.radial_profile(obj, param)
+ assert np.all(np.isfinite(sig.y))
+ assert np.allclose(sig.x / 0.5, np.rint(sig.x / 0.5))
+ assert sig.x[len(sig.x) // 2] == 0.0 and np.allclose(sig.x, -sig.x[::-1])
+
+
+def test_radial_profile_refuses_different_units() -> None:
+ """Distances are undefined when the X and Y units differ."""
+ obj = _image(np.ones((5, 5)))
+ obj.yunit = "s"
+ with pytest.raises(ValueError, match="units"):
+ sipi.radial_profile(obj, sigima.params.RadialProfileParam.create())
+
+
+@pytest.mark.parametrize(
+ "paramclass, kwargs",
+ [
+ (sigima.params.LineProfileParam, {"row": 1}),
+ (sigima.params.SegmentProfileParam, {"row1": 1}),
+ (sigima.params.AverageProfileParam, {"col2": -1}),
+ (sigima.params.RadialProfileParam, {"x0": 1.0}),
+ ],
+)
+def test_former_parameter_names_are_refused(paramclass, kwargs) -> None:
+ """Parameters in pixel indices were renamed: former names fail explicitly."""
+ with pytest.raises((AttributeError, TypeError)):
+ paramclass.create(**kwargs)
+
+
+def test_default_parameters_from_image() -> None:
+ """Defaults set from an image lie inside it (center, middle row)."""
+ obj = _image(np.arange(30.0).reshape(5, 6))
+ obj.set_uniform_coords(2.0, 3.0, 100.0, -50.0)
+ line = sigima.params.LineProfileParam()
+ line.update_from_obj(obj)
+ assert (line.x, line.y) == (obj.xc, obj.yc)
+ sipi.line_profile(obj, line)
+ segment = sigima.params.SegmentProfileParam()
+ segment.update_from_obj(obj)
+ sig = sipi.segment_profile(obj, segment)
+ assert np.array_equal(sig.y, obj.data[2])
+ radial = sigima.params.RadialProfileParam()
+ radial.update_from_obj(obj)
+ assert (radial.xc, radial.yc) == (obj.xc, obj.yc)
diff --git a/sigima/tests/image/profile_unit_test.py b/sigima/tests/image/profile_unit_test.py
index 556d9e93..22e73005 100644
--- a/sigima/tests/image/profile_unit_test.py
+++ b/sigima/tests/image/profile_unit_test.py
@@ -2,6 +2,10 @@
"""
Profile extraction unit test
+
+Profiles are defined and returned in the physical coordinates of the image: the
+tests use calibrated images (origin, anisotropic pixel sizes, non-uniform
+coordinates), whose expected profiles are given by the pixel indices.
"""
# pylint: disable=invalid-name # Allows short reference names like x, y, ...
@@ -15,87 +19,127 @@
from sigima.tests.data import create_sincos_image
from sigima.tests.helpers import check_array_result
+WIDTH, HEIGHT = 256, 128
+# Calibration: X = X0 + DX * column, Y = Y0 + DY * row
+X0, Y0, DX, DY = -10.0, 5.0, 0.5, 2.0
-@pytest.mark.validation
-def test_line_profile() -> None:
- """Test line profile computation"""
- width, height = 256, 128
- dtype = sigima.objects.ImageDatatypes.UINT16
+
+def _image(calibrated: bool = True) -> sigima.objects.ImageObj:
+ """Return the test image, calibrated with anisotropic pixels by default."""
newparam = sigima.objects.NewImageParam.create(
- dtype=dtype, height=height, width=width
+ dtype=sigima.objects.ImageDatatypes.UINT16, height=HEIGHT, width=WIDTH
)
ima = create_sincos_image(newparam)
+ if calibrated:
+ ima.set_uniform_coords(DX, DY, X0, Y0)
+ ima.xunit = ima.yunit = "mm"
+ ima.xlabel, ima.ylabel = "Position X", "Position Y"
+ return ima
+
+
+@pytest.mark.validation
+def test_line_profile() -> None:
+ """Test line profile computation"""
+ ima = _image()
- # Test horizontal line profile
- row = 100
- param = sigima.params.LineProfileParam.create(row=row, direction="horizontal")
+ # Horizontal profile: row whose center is the nearest to Y (row 100 at 205)
+ param = sigima.params.LineProfileParam.create(
+ y=Y0 + 100.4 * DY, direction="horizontal"
+ )
sig = sigima.proc.image.line_profile(ima, param)
- assert sig is not None
- assert len(sig.y) == width
- exp = np.array(ima.data[row, :], dtype=float)
- check_array_result("Horizontal line profile", sig.y, exp)
+ check_array_result("Horizontal line profile", sig.y, ima.data[100, :].astype(float))
+ check_array_result("Horizontal X axis", sig.x, X0 + DX * np.arange(WIDTH))
+ assert (sig.xunit, sig.xlabel) == ("mm", "Position X")
- # Test vertical line profile
- col = 50
- param = sigima.params.LineProfileParam.create(col=col, direction="vertical")
+ # Vertical profile: column 50; X axis along Y, with the Y unit and label
+ param = sigima.params.LineProfileParam.create(
+ x=X0 + 49.6 * DX, direction="vertical"
+ )
sig = sigima.proc.image.line_profile(ima, param)
- assert sig is not None
- assert len(sig.y) == height
- exp = np.array(ima.data[:, col], dtype=float)
- check_array_result("Vertical line profile", sig.y, exp)
+ check_array_result("Vertical line profile", sig.y, ima.data[:, 50].astype(float))
+ check_array_result("Vertical X axis", sig.x, Y0 + DY * np.arange(HEIGHT))
+ assert sig.xlabel == "Position Y"
+
+ # Coordinates outside the image are refused (beyond the outer pixel edges)
+ for direction, kwargs in (
+ ("horizontal", {"y": Y0 - 0.51 * DY}),
+ ("vertical", {"x": X0 + (WIDTH - 0.49) * DX}),
+ ):
+ param = sigima.params.LineProfileParam.create(direction=direction, **kwargs)
+ with pytest.raises(ValueError, match="outside the image"):
+ sigima.proc.image.line_profile(ima, param)
@pytest.mark.validation
def test_segment_profile() -> None:
"""Test segment profile computation"""
- width, height = 256, 128
- dtype = sigima.objects.ImageDatatypes.UINT16
- newparam = sigima.objects.NewImageParam.create(
- dtype=dtype, height=height, width=width
+ ima = _image(calibrated=False)
+ ima.set_uniform_coords(0.5, 0.5, X0, Y0) # Isotropic: one sample per pixel
+
+ # Vertical segment along column 20, from row 10 to row 120
+ param = sigima.params.SegmentProfileParam.create(
+ x1=X0 + 20 * 0.5, y1=Y0 + 10 * 0.5, x2=X0 + 20 * 0.5, y2=Y0 + 120 * 0.5
)
- ima = create_sincos_image(newparam)
+ sig = sigima.proc.image.segment_profile(ima, param)
+ check_array_result("Segment profile", sig.y, ima.data[10:121, 20].astype(float))
+ check_array_result("Segment distances", sig.x, 0.5 * np.arange(111))
+ assert sig.xunit == "mm"
- # Test segment profile
- row1, col1, row2, col2 = 10, 20, 200, 20
+ # Anti-diagonal segment: sampled from its start point (row 0, column 30) to its
+ # end point (row 30, column 0), each sample taking the nearest pixel
param = sigima.params.SegmentProfileParam.create(
- row1=row1, col1=col1, row2=row2, col2=col2
+ x1=X0 + 30 * 0.5, y1=Y0, x2=X0, y2=Y0 + 30 * 0.5
)
sig = sigima.proc.image.segment_profile(ima, param)
- assert sig is not None
- assert len(sig.y) == min(row2, height - 1) - max(row1, 0) + 1
- exp = np.array(ima.data[10:200, 20], dtype=float)
- check_array_result("Segment profile", sig.y, exp)
+ fraction = np.arange(43) * 0.5 / np.hypot(15.0, 15.0)
+ rows = np.floor(30 * fraction + 0.5).astype(int)
+ cols = np.floor(30 - 30 * fraction + 0.5).astype(int)
+ check_array_result(
+ "Anti-diagonal segment", sig.y, ima.data[rows, cols].astype(float)
+ )
+ assert sig.y[0] == ima.data[0, 30]
+
+ # End points outside the image and different X/Y units are refused
+ param = sigima.params.SegmentProfileParam.create(x1=X0, y1=Y0, x2=X0, y2=Y0 + 1e3)
+ with pytest.raises(ValueError, match="outside the image"):
+ sigima.proc.image.segment_profile(ima, param)
+ ima.yunit = "s"
+ param = sigima.params.SegmentProfileParam.create(x1=X0, y1=Y0, x2=X0, y2=Y0 + 1)
+ with pytest.raises(ValueError, match="units"):
+ sigima.proc.image.segment_profile(ima, param)
@pytest.mark.validation
def test_average_profile() -> None:
"""Test average profile computation"""
- width, height = 256, 128
- dtype = sigima.objects.ImageDatatypes.UINT16
- newparam = sigima.objects.NewImageParam.create(
- dtype=dtype, height=height, width=width
- )
- ima = create_sincos_image(newparam)
- row1, col1, row2, col2 = 10, 20, 200, 230
+ ima = _image()
+ # Pixels whose centers lie in the rectangle: rows 10 to 100, columns 20 to 230
param = sigima.params.AverageProfileParam.create(
- row1=row1, col1=col1, row2=row2, col2=col2
+ x1=X0 + 19.6 * DX, y1=Y0 + 10 * DY, x2=X0 + 230.4 * DX, y2=Y0 + 100.4 * DY
)
+ area = ima.data[10:101, 20:231].astype(float)
- # Test horizontal average profile
param.direction = "horizontal"
sig = sigima.proc.image.average_profile(ima, param)
- assert sig is not None
- assert len(sig.y) == col2 - col1 + 1
- exp = np.array(ima.data[row1 : row2 + 1, col1 : col2 + 1].mean(axis=0), dtype=float)
- check_array_result("Horizontal average profile", sig.y, exp)
+ check_array_result("Horizontal average profile", sig.y, area.mean(axis=0))
+ check_array_result("Horizontal X axis", sig.x, X0 + DX * np.arange(20, 231))
- # Test vertical average profile
param.direction = "vertical"
sig = sigima.proc.image.average_profile(ima, param)
- assert sig is not None
- assert len(sig.y) == min(row2, height - 1) - max(row1, 0) + 1
- exp = np.array(ima.data[row1 : row2 + 1, col1 : col2 + 1].mean(axis=1), dtype=float)
- check_array_result("Vertical average profile", sig.y, exp)
+ check_array_result("Vertical average profile", sig.y, area.mean(axis=1))
+ check_array_result("Vertical X axis", sig.x, Y0 + DY * np.arange(10, 101))
+ assert sig.xlabel == "Position Y"
+
+ # Default bounds: the whole image
+ sig = sigima.proc.image.average_profile(
+ ima, sigima.params.AverageProfileParam.create(direction="horizontal")
+ )
+ check_array_result("Whole image", sig.y, ima.data.astype(float).mean(axis=0))
+
+ # A rectangle holding no pixel center is refused
+ param = sigima.params.AverageProfileParam.create(x1=X0 + 0.1 * DX, x2=X0 + 0.4 * DX)
+ with pytest.raises(ValueError, match="no pixel center"):
+ sigima.proc.image.average_profile(ima, param)
def __test_radial_profile_center(
@@ -123,12 +167,7 @@ def __test_radial_profile_center(
@pytest.mark.validation
def test_radial_profile() -> None:
"""Test radial profile computation"""
- width, height = 256, 128
- dtype = sigima.objects.ImageDatatypes.UINT16
- newparam = sigima.objects.NewImageParam.create(
- dtype=dtype, height=height, width=width
- )
- ima = create_sincos_image(newparam)
+ ima = _image(calibrated=False)
# Test radial profile with centroid center
param = sigima.params.RadialProfileParam.create(center="centroid")
@@ -141,8 +180,8 @@ def test_radial_profile() -> None:
__test_radial_profile_center(ima, param)
# Test radial profile with user-defined center
- x0, y0 = width // 2, height // 2
- param = sigima.params.RadialProfileParam.create(center="user", x0=x0, y0=y0)
+ xc, yc = WIDTH // 2, HEIGHT // 2
+ param = sigima.params.RadialProfileParam.create(center="user", xc=xc, yc=yc)
sig = __test_radial_profile_center(ima, param)
# Test that the x-axis represents distance from center (symmetric around 0)
@@ -150,8 +189,22 @@ def test_radial_profile() -> None:
assert sig.x[-1] > 0 # Last element should be positive
assert sig.x[len(sig.x) // 2] == 0 # Center should be at distance 0
- # Additional validation using the helper function for the last profile
- check_array_result("Radial profile", sig.y, sig.y)
+ # Calibrated image (isotropic pixels): same profile, physical distances
+ calibrated = ima.copy()
+ calibrated.set_uniform_coords(0.25, 0.25, X0, Y0)
+ param = sigima.params.RadialProfileParam.create(
+ center="user", xc=X0 + 0.25 * xc, yc=Y0 + 0.25 * yc
+ )
+ csig = __test_radial_profile_center(calibrated, param)
+ check_array_result("Calibrated radial profile", csig.y, sig.y)
+ check_array_result("Calibrated distances", csig.x, 0.25 * sig.x)
+ for center in ("center", "centroid"):
+ param = sigima.params.RadialProfileParam.create(center=center)
+ expected = sigima.proc.image.radial_profile(ima, param)
+ result = sigima.proc.image.radial_profile(calibrated, param)
+ check_array_result(
+ f"Calibrated radial profile ({center})", result.y, expected.y
+ )
if __name__ == "__main__":
diff --git a/sigima/tools/image/extraction.py b/sigima/tools/image/extraction.py
index 59a47130..805b0cd4 100644
--- a/sigima/tools/image/extraction.py
+++ b/sigima/tools/image/extraction.py
@@ -18,31 +18,57 @@
@check_2d_array
def get_radial_profile(
- data: np.ndarray, center: tuple[int, int]
+ data: np.ndarray,
+ center: tuple[float, float],
+ *,
+ xcoords: np.ndarray | None = None,
+ ycoords: np.ndarray | None = None,
+ step: float = 1.0,
) -> tuple[np.ndarray, np.ndarray]:
"""Return radial profile of image data
+ The pixels are binned by their distance to the center (distance divided by
+ *step*, rounded down), and the values of each bin are averaged.
+
Args:
data: Input data (2D array)
- center: Coordinates of the center of the profile (x, y)
+ center: Coordinates of the center of the profile (x, y), in the same units
+ as *xcoords* and *ycoords* (pixel indices by default)
+ xcoords: Coordinates of the pixel centers along X (columns); by default,
+ the column indices
+ ycoords: Coordinates of the pixel centers along Y (rows); by default, the
+ row indices
+ step: Width of the distance bins, in the same units
Returns:
- Radial profile (X, Y) where X is the distance from the center (1D array)
- and Y is the average value of pixels at this distance (1D array)
+ Radial profile (X, Y) where X is the distance from the center (1D array,
+ mirrored: from minus to plus the largest distance, in steps of *step*;
+ distances holding no pixel are left out) and Y is the average value of
+ pixels at this distance (1D array)
"""
- y, x = np.indices((data.shape)) # Get the indices of pixels
+ if xcoords is None:
+ xcoords = np.arange(data.shape[1], dtype=float)
+ if ycoords is None:
+ ycoords = np.arange(data.shape[0], dtype=float)
+ x, y = np.meshgrid(np.asarray(xcoords, float), np.asarray(ycoords, float))
x0, y0 = center
r = np.sqrt((x - x0) ** 2 + (y - y0) ** 2) # Calculate distance to the center
- r = r.astype(int)
+ # Bin index, robust to rounding noise when distances are multiples of the step
+ r = np.floor(r / step + 1e-9).astype(int)
# Average over the same distance
tbin = np.bincount(r.ravel(), data.ravel()) # Sum of pixel values at each distance
nr = np.bincount(r.ravel()) # Number of pixels at each distance
+ # Distances holding no pixel (possible with anisotropic pixels or non-uniform
+ # coordinates) are left out of the profile
+ bins = np.flatnonzero(nr)
- yprofile = tbin / nr # this is the half radial profile
- # Let's mirror it to get the full radial profile (the first element is the center)
- yprofile = np.concatenate((yprofile[::-1], yprofile[1:]))
+ yhalf = tbin[bins] / nr[bins] # this is the half radial profile
+ xhalf = bins * step
+ # Mirror it to get the full radial profile (the first element is the center)
+ skip = 1 if bins[0] == 0 else 0
+ yprofile = np.concatenate((yhalf[::-1], yhalf[skip:]))
# The x axis is the distance from the center (0 is the center)
- xprofile = np.arange(len(yprofile)) - len(yprofile) // 2
+ xprofile = np.concatenate((-xhalf[::-1], xhalf[skip:]))
return xprofile, yprofile