From c1755524de35c34829f6e7bdb2376846d54767ad Mon Sep 17 00:00:00 2001 From: Pierre Raybaut <1311787+PierreRaybaut@users.noreply.github.com> Date: Mon, 26 Oct 2026 07:38:35 +0100 Subject: [PATCH] fix(image): allocate only the occupied bins of radial profiles Since #98, radial profiles bin the pixels in steps of the smallest distance between pixel centers, and get_radial_profile accumulated them with np.bincount on the bin index, which allocates every bin from the center to the largest distance: a tiny spacing (non-uniform coordinates, strongly anisotropic pixels) required arrays of several GB for a few pixels. Coincident pixel centers gave a zero step and an obscure error after divisions by zero. Only the bins holding pixels are now accumulated (np.unique with the inverse index, then np.bincount on it): the memory used never exceeds the size of the image, and the results are bit-identical (same bins, same order of accumulation). A bin width that is not positive and finite, non-finite distances (center or coordinates), and coincident pixel centers are refused with explicit errors; these inputs already failed. The contract of radial_profile is unchanged. Closes #103 Assisted-by: Claude Opus 5.5 --- .../LC_MESSAGES/release_notes/release_1.04.po | 4 +- doc/release_notes/release_1.04.md | 2 +- sigima/proc/image/extraction.py | 20 ++++- .../tests/image/profile_physical_unit_test.py | 79 +++++++++++++++++++ sigima/tools/image/extraction.py | 30 ++++--- 5 files changed, 118 insertions(+), 17 deletions(-) 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 e61b03b..c564dbe 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 @@ -79,8 +79,8 @@ msgstr "**Coordonnées d'ellipse** : `EllipseCoordinates` (`sigima.objects.shape 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 profiles are sampled at regular steps from their start point to their end point included, with one sample per pixel along the axis where the segment spans the most pixels (as in Sigima 1.3), each sample taking the value of the nearest pixel. Radial profiles are binned in steps of the smallest distance between pixel centers, so that they are circular for anisotropic pixels and non-uniform coordinates. Segment and radial profiles 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), line, average and radial profiles are unchanged, and segment profiles take the same pixels for horizontal, vertical and diagonal segments: other oblique segments now take the nearest pixel of each sample instead of truncating its indices. 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 sont échantillonnés à pas réguliers de leur point de départ à leur point d'arrivée inclus, avec un échantillon par pixel le long de l'axe sur lequel le segment couvre le plus de pixels (comme dans Sigima 1.3), chaque échantillon prenant la valeur du pixel le plus proche. Les profils radiaux sont regroupés par pas égaux à la plus petite distance entre centres de pixels, de sorte qu'ils sont circulaires pour des pixels anisotropes et des coordonnées non uniformes. Les profils de segment et radiaux 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 profils ligne, moyen et radial sont inchangés, et les profils de segment prennent les mêmes pixels pour les segments horizontaux, verticaux et diagonaux : les autres segments obliques prennent désormais le pixel le plus proche de chaque échantillon au lieu de tronquer ses indices. Les contrats de `line_profile`, `segment_profile`, `average_profile` et `radial_profile` sont 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 profiles are sampled at regular steps from their start point to their end point included, with one sample per pixel along the axis where the segment spans the most pixels (as in Sigima 1.3), each sample taking the value of the nearest pixel. Radial profiles are binned in steps of the smallest distance between pixel centers, so that they are circular for anisotropic pixels and non-uniform coordinates. Segment and radial profiles are refused when the X and Y units of the image differ (and radial profiles when the centers of neighboring pixels coincide). Coordinates outside the image are refused instead of being clamped. On images without calibration (origin 0, pixel size 1), line, average and radial profiles are unchanged, and segment profiles take the same pixels for horizontal, vertical and diagonal segments: other oblique segments now take the nearest pixel of each sample instead of truncating its indices. 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 sont échantillonnés à pas réguliers de leur point de départ à leur point d'arrivée inclus, avec un échantillon par pixel le long de l'axe sur lequel le segment couvre le plus de pixels (comme dans Sigima 1.3), chaque échantillon prenant la valeur du pixel le plus proche. Les profils radiaux sont regroupés par pas égaux à la plus petite distance entre centres de pixels, de sorte qu'ils sont circulaires pour des pixels anisotropes et des coordonnées non uniformes. Les profils de segment et radiaux sont refusés lorsque les unités X et Y de l'image diffèrent (et les profils radiaux lorsque les centres de pixels voisins sont confondus). 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 profils ligne, moyen et radial sont inchangés, et les profils de segment prennent les mêmes pixels pour les segments horizontaux, verticaux et diagonaux : les autres segments obliques prennent désormais le pixel le plus proche de chaque échantillon au lieu de tronquer ses indices. 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" diff --git a/doc/release_notes/release_1.04.md b/doc/release_notes/release_1.04.md index 0e2abce..d2be646 100644 --- a/doc/release_notes/release_1.04.md +++ b/doc/release_notes/release_1.04.md @@ -27,7 +27,7 @@ The [user guide](../user_guide/overview.rst) shows how readers and connectors at * **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 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 profiles are sampled at regular steps from their start point to their end point included, with one sample per pixel along the axis where the segment spans the most pixels (as in Sigima 1.3), each sample taking the value of the nearest pixel. Radial profiles are binned in steps of the smallest distance between pixel centers, so that they are circular for anisotropic pixels and non-uniform coordinates. Segment and radial profiles 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), line, average and radial profiles are unchanged, and segment profiles take the same pixels for horizontal, vertical and diagonal segments: other oblique segments now take the nearest pixel of each sample instead of truncating its indices. The contracts of `line_profile`, `segment_profile`, `average_profile` and `radial_profile` are 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 profiles are sampled at regular steps from their start point to their end point included, with one sample per pixel along the axis where the segment spans the most pixels (as in Sigima 1.3), each sample taking the value of the nearest pixel. Radial profiles are binned in steps of the smallest distance between pixel centers, so that they are circular for anisotropic pixels and non-uniform coordinates. Segment and radial profiles are refused when the X and Y units of the image differ (and radial profiles when the centers of neighboring pixels coincide). Coordinates outside the image are refused instead of being clamped. On images without calibration (origin 0, pixel size 1), line, average and radial profiles are unchanged, and segment profiles take the same pixels for horizontal, vertical and diagonal segments: other oblique segments now take the nearest pixel of each sample instead of truncating its indices. The contracts of `line_profile`, `segment_profile`, `average_profile` and `radial_profile` are at version 2. ### 🛠️ Bug Fixes since version 1.3.0 ### diff --git a/sigima/proc/image/extraction.py b/sigima/proc/image/extraction.py index d71ca51..c834160 100644 --- a/sigima/proc/image/extraction.py +++ b/sigima/proc/image/extraction.py @@ -378,12 +378,21 @@ def _pixel_spacing(src: ImageObj) -> float: Returns: Smallest pixel spacing along X and Y (uniform or non-uniform coordinates) + + Raises: + ValueError: if the centers of neighboring pixels coincide """ 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))))) + spacing = float(np.min(np.abs(np.diff(centers)))) + if spacing == 0.0: + raise ValueError( + f"Pixel centers coincide along {axis.upper()}: distances " + f"between pixels are undefined" + ) + spacings.append(spacing) 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] @@ -708,8 +717,10 @@ def radial_profile(src: ImageObj, p: RadialProfileParam) -> SignalObj: 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. + pixel sizes and for non-uniform coordinates. Only the distances holding pixels + are kept, so that the memory used never exceeds the size of the image. The X + axis of the profile is the distance to the center (mirrored), in the units of + the image. Args: src: input image object @@ -719,7 +730,8 @@ def radial_profile(src: ImageObj, p: RadialProfileParam) -> SignalObj: Signal object with the radial profile Raises: - ValueError: if the X and Y units of the image differ + ValueError: if the X and Y units of the image differ, if the centers of + neighboring pixels coincide, or if the center is not finite """ _check_same_units(src, "Radial profile") data = src.get_masked_view() diff --git a/sigima/tests/image/profile_physical_unit_test.py b/sigima/tests/image/profile_physical_unit_test.py index c762bc8..4a653f1 100644 --- a/sigima/tests/image/profile_physical_unit_test.py +++ b/sigima/tests/image/profile_physical_unit_test.py @@ -9,12 +9,15 @@ from __future__ import annotations +import tracemalloc + import numpy as np import pytest import sigima.objects import sigima.params import sigima.proc.image as sipi +import sigima.tools.image def _image(data: np.ndarray) -> sigima.objects.ImageObj: @@ -77,6 +80,82 @@ def test_radial_profile_refuses_different_units() -> None: sipi.radial_profile(obj, sigima.params.RadialProfileParam.create()) +def _radial_profile_all_bins(data, center, xcoords, ycoords, step): + """Former implementation of ``get_radial_profile``, allocating every bin from + the center to the largest distance.""" + x, y = np.meshgrid(xcoords, ycoords) + r = np.sqrt((x - center[0]) ** 2 + (y - center[1]) ** 2) + r = np.floor(r / step + 1e-9).astype(int) + tbin = np.bincount(r.ravel(), data.ravel()) + nr = np.bincount(r.ravel()) + bins = np.flatnonzero(nr) + yhalf = tbin[bins] / nr[bins] + xhalf = bins * step + skip = 1 if bins[0] == 0 else 0 + return ( + np.concatenate((-xhalf[::-1], xhalf[skip:])), + np.concatenate((yhalf[::-1], yhalf[skip:])), + ) + + +@pytest.mark.parametrize( + "xcoords, ycoords, center, step", + [ + (np.arange(53.0), np.arange(37.0), (17.3, 20.6), 1.0), # uniform + (np.arange(53.0), np.arange(37.0), (26.0, 18.0), 1.0), # pixel center + (np.arange(53) * 0.1, np.arange(37) * 0.2, (2.6, 3.5), 0.1), # anisotropic + (np.cumsum(np.linspace(0.1, 2.0, 53)), np.arange(37.0), (30.0, 9.5), 0.1), + ], +) +def test_radial_profile_same_as_all_bins(xcoords, ycoords, center, step) -> None: + """Accumulating only the bins holding pixels gives bit-identical profiles.""" + rng = np.random.default_rng(0) + data = rng.normal(100.0, 20.0, (ycoords.size, xcoords.size)) + for values in (data, data.astype(np.uint16), np.ma.masked_less(data, 90.0)): + x, y = sigima.tools.image.get_radial_profile( + values, center, xcoords=xcoords, ycoords=ycoords, step=step + ) + xref, yref = _radial_profile_all_bins(values, center, xcoords, ycoords, step) + assert np.array_equal(x, xref) and np.array_equal(y, yref) + + +def test_radial_profile_tiny_spacing() -> None: + """The memory used does not depend on the pixel spacing (regression: a tiny + spacing allocated one bin per step up to the largest distance, i.e. several + GB for a few pixels).""" + data = np.arange(25.0).reshape(5, 5) + obj = _image(data) + obj.set_coords(np.array([0.0, 1e-8, 1.0, 2.0, 3.0]), np.arange(5.0)) + param = sigima.params.RadialProfileParam.create(center="user", xc=0.0, yc=0.0) + tracemalloc.start() + try: + sig = sipi.radial_profile(obj, param) + _current, peak = tracemalloc.get_traced_memory() + finally: + tracemalloc.stop() + assert peak < 10e6 + assert sig.x.size <= 2 * data.size - 1 + assert np.all(np.isfinite(sig.y)) and np.allclose(sig.x, -sig.x[::-1]) + assert sig.y[sig.x.size // 2] == data[0, 0] + + +def test_radial_profile_refusals() -> None: + """Coincident pixel centers, invalid bin widths and centers are refused + explicitly (regression: obscure error after divisions by zero).""" + obj = _image(np.ones((5, 5))) + obj.set_coords(np.array([0.0, 0.0, 1.0, 2.0, 3.0]), np.arange(5.0)) + param = sigima.params.RadialProfileParam.create(center="user", xc=1.0, yc=1.0) + with pytest.raises(ValueError, match="coincide along X"): + sipi.radial_profile(obj, param) + data = np.ones((5, 5)) + for step in (0.0, -1.0, np.nan, np.inf): + with pytest.raises(ValueError, match="bin width"): + sigima.tools.image.get_radial_profile(data, (2.0, 2.0), step=step) + for center in ((np.nan, 2.0), (2.0, np.inf)): + with pytest.raises(ValueError, match="finite"): + sigima.tools.image.get_radial_profile(data, center) + + @pytest.mark.parametrize( "paramclass, kwargs", [ diff --git a/sigima/tools/image/extraction.py b/sigima/tools/image/extraction.py index 805b0cd..0217fa4 100644 --- a/sigima/tools/image/extraction.py +++ b/sigima/tools/image/extraction.py @@ -28,7 +28,9 @@ def get_radial_profile( """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. + *step*, rounded down), and the values of each bin are averaged. Only the bins + holding pixels are allocated, so that the memory used never exceeds the size + of the image, whatever the step. Args: data: Input data (2D array) @@ -45,7 +47,13 @@ def get_radial_profile( 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) + + Raises: + ValueError: if *step* is not positive and finite, or if the center or the + pixel coordinates are not finite """ + if not (np.isfinite(step) and step > 0): + raise ValueError(f"The bin width must be positive and finite (got {step})") if xcoords is None: xcoords = np.arange(data.shape[1], dtype=float) if ycoords is None: @@ -53,17 +61,19 @@ def get_radial_profile( 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 - # Bin index, robust to rounding noise when distances are multiples of the step - r = np.floor(r / step + 1e-9).astype(int) + if not np.all(np.isfinite(r)): + raise ValueError("The center and the pixel coordinates must be finite") + # Bin of each pixel, robust to rounding noise when distances are multiples of + # the step: only the bins holding pixels are kept (sorted), each pixel being + # given the position of its bin among them + bins, index = np.unique(np.floor(r / step + 1e-9).ravel(), return_inverse=True) - # 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) + # Average over the same distance (distances holding no pixel, possible with + # anisotropic pixels or non-uniform coordinates, are left out of the profile) + tbin = np.bincount(index, data.ravel()) # Sum of pixel values at each distance + nr = np.bincount(index) # Number of pixels at each distance - yhalf = tbin[bins] / nr[bins] # this is the half radial profile + yhalf = tbin / nr # 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