"""Display functionality, including scan conversion frustrum conversion, etc."""
from functools import partial
from typing import Tuple, Union
import numpy as np
import scipy
from keras import ops
from PIL import Image
from zea import log
from zea.func.tensor import translate
from zea.internal.utils import deprecated, deprecated_alias, deprecated_distance_to_apex
from zea.tools.fit_scan_cone import fit_and_crop_around_scan_cone
# `zea.func.ultrasound.log_compress` maps zeros to 20 * log10(1e-16) = -320 dB. Such pixels
# hold no measurement, so they are excluded when fitting a histogram transform.
EXCLUDE_BELOW: float = -300.0
[docs]
def to_8bit(image, dynamic_range: Union[None, tuple] = None, pillow: bool = True):
"""Convert image to 8 bit image [0, 255]. Clip between dynamic range.
Args:
image (ndarray): Input image(s). Should be in between dynamic range.
dynamic_range (tuple, optional): Dynamic range of input image(s).
pillow (bool, optional): Whether to return PIL image. Defaults to True.
Returns:
image (ndarray): Output 8 bit image(s) [0, 255].
.. note::
If dynamic_range is None, it is assumed that the input image is already in the range
[-60, 0] dB, which is a common range for ultrasound images.
.. note::
NaN values in the input image are replaced with the minimum value of the dynamic range
before scaling, which ensures that they are represented as black (0) in the output image.
+/- inf values are replaced with the min and max values of the dynamic range.
Example:
.. doctest::
>>> import numpy as np
>>> import zea
>>> file_path = (
... "hf://zeahub/camus-sample/val/patient0401/patient0401_4CH_half_sequence.hdf5"
... )
>>> with zea.File(file_path, mode="r") as file:
... data = file.data.image.values[0]
>>> image, _ = zea.display.scan_convert(
... data,
... rlims=(0, 1),
... polar_limits=(-0.78, 0.78),
... fill_value=np.nan,
... )
>>> image = zea.display.to_8bit(image, dynamic_range=(-60, 0))
>>> image.save("image.png") # doctest: +SKIP
"""
if dynamic_range is None:
dynamic_range = (-60, 0)
image = ops.nan_to_num(image, nan=dynamic_range[0])
image = ops.convert_to_numpy(image)
image = np.clip(image, *dynamic_range)
image = translate(image, dynamic_range, (0, 255))
image = image.astype(np.uint8)
if pillow:
image = Image.fromarray(image)
return image
def _masked_pixels(image, roi=None, minimum=2, floor=EXCLUDE_BELOW):
"""Mask out the pixels a transform may not use: unmeasured, or outside ``roi`` if given."""
mask = np.isfinite(image) & (image > floor)
if roi is not None:
# Without this, a mask of a broadcastable shape (a row, a column) selects pixels
# nobody asked for rather than failing.
if roi.shape != image.shape:
raise ValueError(
f"roi shape {roi.shape} does not match the shape {image.shape} of the image "
"it selects in; an roi has to index both the image and the reference."
)
mask &= roi
pixels = image[mask].astype(np.float64)
if pixels.size < minimum:
raise ValueError(
f"Too few valid pixels to fit a transform: {pixels.size}, need at least {minimum}."
)
if pixels.min() == pixels.max():
raise ValueError("Cannot fit a transform: the valid pixels all have the same value.")
return pixels
def _quantile_pairs(x, y, n_bins):
"""One quantile of each distribution per bin, i.e. eq. (13): ``h[i] = j`` such that
``F_X[i] = F_Y[j]``. The pairs are the knots the mapping is interpolated through.
The reference implementation (histmatch.m) bins equispaced in *amplitude*, which on
log-compressed data spends nearly all of its resolution on the tail towards the noise floor.
Binning equispaced in *probability* instead puts equal mass in every bin, so the resolution
follows the data (the paper's own remark in Sec. V that bins should be chosen after
dynamic range compression). Hence the bins are counted out in probability, as the
``n_bins`` levels at which both CDFs are read.
At ``n_bins = x.size`` the quantiles of ``x`` are its order statistics, so the pairs
become the rank transport of Sec. III.C.1: one pixel per bin, with the quantile function
of ``y`` resampled onto the ranks of ``x`` when the two differ in size.
"""
levels = np.linspace(0, 1, n_bins)
def quantiles(v):
# `np.quantile` selects rather than sorts, which is quadratic in the number of
# levels; this is its linear method, interpolating between the order statistics.
v = np.sort(v)
return np.interp(levels, np.linspace(0, 1, v.size), v)
return quantiles(x), quantiles(y)
def _merge_ties(xs, ys):
"""Collapse knots that share an ``x``, which PCHIP cannot accept, onto a single ``y``.
Knots collapse wherever the image distribution has an atom, i.e. a value carrying a
finite fraction of the mass: a pile-up at the edge of the range the image was clipped
to, or a repeated level in a quantized image. A monotone map has to send an atom to one
value, and this takes the lowest of the reference values it spans, which is the CDF read
just below the atom. A clipping floor then stays at the bottom of the matched range
rather than being lifted into the middle of it.
"""
xs, first = np.unique(xs, return_index=True)
return xs, ys[first]
def _end_slopes(xs, ys, frac=0.125, max_ratio=4.0):
"""Slopes for linear extension of the mapping past the first and last knot."""
# Secants between adjacent knots are unusable: an atom collapses knots onto each other
# (see `_merge_ties`), which makes the secant explode. Measure across the outer `frac`
# of the knots instead, which because they are equispaced in probability is the outer
# `frac` of the mass, and cap the result relative to the overall slope.
slope = (ys[-1] - ys[0]) / (xs[-1] - xs[0])
k = max(2, round(frac * xs.size))
def secant(dx, dy):
if dx < 1e-3: # degenerate span, nothing local to learn from
return slope
return float(np.clip(dy / dx, slope / max_ratio, slope * max_ratio))
return (
secant(xs[k - 1] - xs[0], ys[k - 1] - ys[0]),
secant(xs[-1] - xs[-k], ys[-1] - ys[-k]),
)
def _monotone_map(image, xs, ys, spline=True):
"""Monotone interpolation through the knots, one per bin, extended linearly at the ends.
A spline shapes the mapping *between* the knots, which is what a bin pooling many pixels
into one knot leaves undetermined. At one knot per pixel (``spline=False``) there is no
such gap to shape, so the knots are joined by straight lines instead.
"""
# The end slopes are read off the knots before merging, which loses how much mass each
# of them carries and with it the meaning of "the outer `frac` of the distribution".
slope_low, slope_high = _end_slopes(xs, ys)
xs, ys = _merge_ties(xs, ys)
if spline:
out = scipy.interpolate.PchipInterpolator(xs, ys)(image)
else:
out = np.interp(image, xs, ys)
# Cubic extrapolation can fold back on itself and break monotonicity, so replace it by
# a linear extension (Sec. III.C.3): a point target brighter than anything in the ROI
# should stay the brightest instead of piling onto the top knot.
out = np.where(image < xs[0], ys[0] + slope_low * (image - xs[0]), out)
return np.where(image > xs[-1], ys[-1] + slope_high * (image - xs[-1]), out)
[docs]
def histogram_match(
image,
reference,
mode: str = "full",
n_bins: Union[int, str] = 256,
roi=None,
allow_unequal_shapes: bool = False,
exclude_below: Union[float, None] = EXCLUDE_BELOW,
) -> np.ndarray:
"""Match the histogram of an image to a reference image, for fair visual comparison.
Image formation methods apply different dynamic range transformations, which biases
visual as well as quantitative comparison. Matching the histogram of an image to a
reference (typically a conventional B-mode image) puts the two on a common scale.
.. citation:: bottenus2021histogram
Their MATLAB implementation, referred to above as the reference
implementation:
https://github.com/nbottenus/histogram_matching/blob/main/histmatch.m
Args:
image (ndarray): Image to transform, typically log-compressed (dB).
reference (ndarray): Image to match to. Has to be of the shape of ``image``, unless
``allow_unequal_shapes`` says otherwise.
mode (str, optional): Class of transform to match with. Defaults to ``"full"``.
- ``"partial"``: affine match (Sec. III.A), the scale and offset that match
mean and variance. It leaves the *shape* of the image's own distribution
intact, so pick it when that shape is the thing under comparison, when the
fitted region contains structure rather than plain speckle, or when the CNR
must survive the matching.
- ``"full"``: any monotone mapping (Sec. III.B), the one that carries the
histogram of the image onto that of the reference, estimated on ``n_bins``
bins. It matches the shape of the distribution too, hence undoes non-linear
compression, so pick it to compare methods whose transformations differ in
more than scale and offset, or are unknown ("black box" images, Sec. III.C.2).
Together with a speckle ``roi``, the paper's expected common use (Sec. V).
n_bins (int or str, optional): Number of bins the mapping of a ``"full"`` match is
estimated on; ignored by ``"partial"``. Defaults to 256, the number of bins used
in the paper, which it found enough for consistent results.
Pass ``"all"`` for one bin per fitted pixel, the limit of the binning process
(Sec. III.C.1): it reproduces the reference distribution down to its outliers and tails,
which replicates them in the matched image. This is what
`skimage.exposure.match_histograms` does.
roi (ndarray, optional): Boolean mask, of the shape of both images, selecting a
homogeneous speckle region (Sec. III.C.3, the paper's expected common use).
The mapping is fit inside the mask and applied to the whole image, extended
linearly outside the fitted range. Defaults to None (fit on the whole image).
allow_unequal_shapes (bool, optional): Whether ``reference`` may be of a shape other
than that of ``image``. Only the two distributions are compared, so a match is
well defined either way, but images of different shapes are more often a mistake
(a stray transpose, the wrong frame of a sequence) than an intent. Defaults to
False. Pass True to match against a reference deliberately not of the shape of
the image, e.g. a region cut out of a larger acquisition. An ``roi`` has to index
both images regardless, so it cannot be combined with unequal shapes.
exclude_below (float, optional): Value at or below which a pixel of either image
is kept out of the fit and the transform. Defaults to -300.0, which catches
zea's ``log(0)`` sentinel and nothing a log-compressed B-mode measures. Pass `None`
to exclude nothing but the non-finite pixels.
Returns:
ndarray: Transformed image (float64), unclipped: clipping to a display range is
left to the caller (`to_8bit`, or ``vmin``/``vmax``).
.. note::
Which image is matched to which is a choice (Sec. V): matching to B-mode is natural
for a reader used to B-mode, but features outside its dynamic range are lost in the
process, and matching in the other direction then shows them off better.
Example:
.. doctest::
>>> import numpy as np
>>> from zea.display import histogram_match
>>> rng = np.random.default_rng(0)
>>> image = 20 * np.log10(rng.rayleigh(0.1, (64, 64)))
>>> reference = 20 * np.log10(rng.rayleigh(1.0, (64, 64)))
>>> round(image.mean() - reference.mean()) # the image is 20 dB darker
-20
>>> matched = histogram_match(image, reference)
>>> round(matched.mean() - reference.mean())
0
"""
if mode not in ("full", "partial"):
raise ValueError(f"Unknown mode {mode!r}, expected 'full' or 'partial'.")
every_pixel = n_bins == "all"
if not every_pixel and not (isinstance(n_bins, (int, np.integer)) and n_bins >= 2):
raise ValueError(f"Invalid n_bins {n_bins!r}, expected an integer >= 2 or 'all'.")
image = ops.convert_to_numpy(image)
reference = ops.convert_to_numpy(reference)
if not allow_unequal_shapes and image.shape != reference.shape:
raise ValueError(
f"image shape {image.shape} does not match reference shape {reference.shape}; "
"pass allow_unequal_shapes=True to match against a reference of another shape."
)
if roi is not None:
roi = ops.convert_to_numpy(roi).astype(bool)
# A threshold nothing falls at or below excludes nothing, which is what None asks for.
floor = -np.inf if exclude_below is None else exclude_below
# A bin has to be worth estimating: ask for as many pixels as there are bins.
minimum = 2 if mode == "partial" or every_pixel else n_bins
x = _masked_pixels(image, roi, minimum, floor)
y = _masked_pixels(reference, roi, minimum, floor)
image = image.astype(np.float64)
# Pixels that hold no measurement are kept out of the transform as well as out of the fit:
# a mapping they never informed would place them anywhere, and a slope below one is enough
# to lift the sentinel out of the bottom of the range and into the displayed one.
unmeasured = image <= floor
if mode == "partial":
# eq. (8)-(9): the scale and offset that match mean and variance.
slope = y.std() / x.std()
matched = slope * image + (y.mean() - slope * x.mean())
else:
xs, ys = _quantile_pairs(x, y, x.size if every_pixel else n_bins)
matched = _monotone_map(image, xs, ys, spline=not every_pixel)
# Non-finite pixels are not `unmeasured` (a comparison against NaN is False), and pass
# through the transform itself: NaN in, NaN out.
return np.where(unmeasured, image, matched)
[docs]
def overlay_masks(
image,
masks,
alpha: float = 0.5,
colors=None,
):
"""Overlay segmentation masks on top of an image using PIL.
Args:
image (PIL.Image or ndarray): Base image. If grayscale, it is converted
to RGB. If ndarray, it is converted to a PIL Image first.
masks (list of PIL.Image or ndarray): Segmentation masks to overlay.
Each mask should be an 8-bit single-channel image where non-zero
pixels indicate the masked region.
alpha (float, optional): Opacity of the mask overlays in [0, 1].
Defaults to 0.5.
colors (list of tuple, optional): RGB colors for each mask. If None,
a default palette is used. If provided, must contain at least as
many entries as masks (extra entries are ignored).
Returns:
PIL.Image: RGB image with masks overlaid.
"""
# Validate alpha parameter before conversion to uint8
if not (0.0 <= alpha <= 1.0):
raise ValueError(f"alpha must be in the range [0.0, 1.0], got {alpha}")
_DEFAULT_COLORS = [
(255, 0, 0),
(0, 255, 0),
(0, 0, 255),
(255, 255, 0),
(0, 255, 255),
(255, 0, 255),
]
if not isinstance(image, Image.Image):
image = Image.fromarray(np.asarray(image))
if image.mode != "RGB":
image = image.convert("RGB")
# Validate colors list has enough entries if provided
if colors is not None and len(colors) < len(masks):
raise ValueError(
f"colors must have at least as many entries as masks: "
f"got {len(colors)} colors for {len(masks)} masks"
)
result = image.copy()
for i, mask in enumerate(masks):
if not isinstance(mask, Image.Image):
mask = Image.fromarray(np.asarray(mask))
if mask.size != image.size:
raise ValueError(f"Mask {i} size {mask.size} does not match image size {image.size}")
if mask.mode != "L":
mask = mask.convert("L")
color = _DEFAULT_COLORS[i % len(_DEFAULT_COLORS)] if colors is None else colors[i]
# Create a solid color layer the same size as the image
color_layer = Image.new("RGB", image.size, color)
# Build alpha channel from the mask: scale mask values by alpha
mask_np = (np.asarray(mask) > 0).astype(np.uint8)
alpha_channel = Image.fromarray((mask_np * int(alpha * 255)).astype(np.uint8))
result.paste(color_layer, mask=alpha_channel)
return result
def _polar_sampling_coordinates(x_vec, z_vec, rho, theta):
"""Map an explicit Cartesian (x, z) output grid onto fractional (rho, theta) indices
into a polar image, using the frustum convention.
Args:
x_vec, z_vec: 1D Cartesian output axes (x is lateral, z is depth).
rho, theta: 1D polar axes of the input image (``rho`` radial, ``theta`` angular).
Their lengths set the polar image dimensions the indices address.
Returns:
coordinates (Array): ``(2, len(x_vec), len(z_vec))`` stack of ``[rho_idx, theta_idx]``
ready for :func:`map_coordinates`.
"""
# The index mapping is endpoint-based (rho[0]/rho[-1], theta[0]/theta[-1]), so the *order*
# of the axes is preserved: a descending ``theta`` maps the angular axis in reverse. The
# masking limits, by contrast, are sorted so out-of-cone pixels are dropped regardless of
# axis direction.
rho_start, rho_end = rho[0], rho[-1]
theta_start, theta_end = theta[0], theta[-1]
theta_min = ops.minimum(theta_start, theta_end)
theta_max = ops.maximum(theta_start, theta_end)
z_grid, x_grid = ops.meshgrid(z_vec, x_vec)
rho_grid_interp, theta_grid_interp = xz_to_frustum(
x_grid, z_grid, polar_limits=[theta_min, theta_max]
)
n_rho, n_theta = ops.shape(rho)[0], ops.shape(theta)[0]
rho_idx = (rho_grid_interp - rho_start) / (rho_end - rho_start) * (n_rho - 1)
theta_idx = (theta_grid_interp - theta_start) / (theta_end - theta_start) * (n_theta - 1)
return ops.stack([rho_idx, theta_idx], axis=0)
def _deprecated_mirrored_theta_range(polar_limits, theta_range, caller, reverse):
"""Map the deprecated ``theta_range`` of the Cartesian <-> polar matrix helpers onto
``polar_limits``.
These helpers used to measure angles positive toward -x, the mirror image of zea's
``polar_limits``, so the angles change sign. ``reverse`` also swaps them, for the
helpers that laid out their columns from ``theta_range[1]`` to ``theta_range[0]``.
Either way the result gives the same image as ``theta_range`` did.
"""
if theta_range is None:
return polar_limits
if polar_limits is not None:
raise TypeError(
f"{caller}() got both 'polar_limits' and its deprecated alias 'theta_range'."
)
if reverse:
mirrored, rule = (-theta_range[1], -theta_range[0]), "(-b, -a)"
else:
mirrored, rule = (-theta_range[0], -theta_range[1]), "(-a, -b)"
log.deprecated(
f"{caller}(theta_range=...) is deprecated, use polar_limits=... instead. theta_range "
"measured angles positive toward -x, polar_limits measures them toward +x as in the "
f"rest of zea: theta_range=(a, b) is polar_limits={rule}."
)
return mirrored
[docs]
def compute_scan_convert_2d_coordinates(
image_shape,
rlims: Tuple[float, float] | None = None,
polar_limits: Tuple[float, float] | None = None,
resolution: Union[float, None] = None,
dtype: str = "float32",
distance_to_apex=None,
*,
grid_apex=None,
extent=None,
rho_range=None,
theta_range=None,
):
"""Precompute coordinates for 2d scan conversion from polar coordinates.
The polar image has one row per radius and one column per polar angle. The radius is
measured from ``grid_apex``, the point the polar grid fans out from (see
:attr:`zea.Parameters.grid_apex`), and the polar angle is ``atan(x / z)``, from the
z-axis in the x-z plane, positive towards +x. The returned coordinates index the polar
image, so
they live in that apex frame; ``grid_apex`` only moves the output window into the
transducer frame, where the reported ``x_lim`` and ``z_lim`` are given (matching
:attr:`zea.Parameters.extent`).
Args:
image_shape (tuple): Shape of the polar image, ``(..., n_rho, n_theta)``.
rlims (tuple): Radius of the first and last row, measured from ``grid_apex``.
polar_limits (tuple): Angle of the first and last column [rad].
resolution (float, optional): Output pixel size. Derived from the polar image
when not given.
dtype (str, optional): Dtype of the coordinates. Defaults to ``"float32"``.
grid_apex (array-like, optional): The apex ``(x, y, z)`` in the transducer frame.
Defaults to the origin, which makes the two frames coincide.
extent (array-like, optional): The output window in the transducer frame, as
:attr:`zea.Parameters.extent` gives it: ``(xmin, xmax, ..., zmin, zmax)``.
Pixels outside the polar image are filled. Defaults to the bounding box of
the polar image.
distance_to_apex (float, optional): Deprecated, use
``grid_apex=(0, 0, -distance_to_apex)``.
rho_range (tuple, optional): Deprecated, use ``rlims``.
theta_range (tuple, optional): Deprecated, use ``polar_limits``.
Returns:
tuple: The sampling coordinates and a dict with ``resolution``, ``x_lim``,
``z_lim``, ``rlims``, ``polar_limits``, ``dr`` and ``d_polar`` (the radius and angle
steps between rows and columns).
"""
caller = "compute_scan_convert_2d_coordinates"
rlims = deprecated_alias(rlims, rho_range, "rho_range", "rlims", caller)
polar_limits = deprecated_alias(
polar_limits, theta_range, "theta_range", "polar_limits", caller
)
if rlims is None or polar_limits is None:
raise ValueError(f"{caller}() requires rlims and polar_limits.")
assert len(rlims) == 2, "rlims should be a tuple of length 2"
assert len(polar_limits) == 2, "polar_limits should be a tuple of length 2"
assert rlims[0] < rlims[1], "rlims[0] should be less than rlims[1]"
grid_apex = deprecated_distance_to_apex(distance_to_apex, grid_apex, caller)
if grid_apex is None:
apex_x = apex_z = 0.0
else:
grid_apex = ops.cast(ops.convert_to_tensor(grid_apex), dtype)
apex_x, apex_z = grid_apex[0], grid_apex[2]
if extent is not None:
extent = ops.cast(ops.convert_to_tensor(extent), dtype)
rho = ops.linspace(rlims[0], rlims[1], image_shape[-2], dtype=dtype)
theta = ops.linspace(polar_limits[0], polar_limits[1], image_shape[-1], dtype=dtype)
rho_grid, theta_grid = ops.meshgrid(rho, theta, indexing="ij")
x_grid, z_grid = frustum_to_xz(rho_grid, theta_grid)
# Output window in the apex frame
if extent is None:
x_lim = [ops.min(x_grid), ops.max(x_grid)]
z_lim = [ops.min(z_grid), ops.max(z_grid)]
else:
x_lim = [extent[0] - apex_x, extent[1] - apex_x]
z_lim = [extent[-2] - apex_z, extent[-1] - apex_z]
d_rho = rho[1] - rho[0]
d_theta = theta[1] - theta[0]
if resolution is None:
# arc length along the polar direction at 1/4 depth; the angles may run either way
sRT = 0.25 * (rho[0] + rho[-1]) * ops.abs(d_theta)
# average of arc lengths and radial step
resolution = ops.mean([sRT, d_rho]) # mm per pixel
x_vec = ops.arange(x_lim[0], x_lim[1], resolution)
z_vec = ops.arange(z_lim[0], z_lim[1], resolution)
coordinates = _polar_sampling_coordinates(x_vec, z_vec, rho, theta)
parameters = {
"resolution": resolution,
"x_lim": [x_lim[0] + apex_x, x_lim[1] + apex_x],
"z_lim": [z_lim[0] + apex_z, z_lim[1] + apex_z],
"rlims": rlims,
"polar_limits": polar_limits,
"dr": d_rho,
"d_polar": d_theta,
}
return coordinates, parameters
[docs]
def scan_convert_2d(
image,
rlims: Tuple[float, float] | None = None,
polar_limits: Tuple[float, float] | None = None,
resolution: Union[float, None] = None,
coordinates: Union[None, np.ndarray] = None,
fill_value: float = 0.0,
order: int = 1,
distance_to_apex=None,
*,
grid_apex=None,
extent=None,
rho_range=None,
theta_range=None,
**kwargs,
):
"""
Perform scan conversion on a 2D ultrasound image from polar coordinates
(radius, polar angle) to Cartesian coordinates (x, z).
Args:
image (ndarray): The input 2D ultrasound image in polar coordinates.
Has dimensions (n_radius, n_polar) with optional batch.
rlims (tuple): Radius ``(rho_min, rho_max)`` of the first and last row,
measured from ``grid_apex``. In the same unit as the output grid.
polar_limits (tuple): Polar angle ``atan(x / z)`` of the first and last column
[rad], from the z-axis in the x-z plane, positive towards +x.
resolution (float, optional): The resolution for the Cartesian grid.
If None, it is calculated based on the input image. In units of ``rlims``
per pixel.
coordinates (ndarray, optional): Precomputed coordinates for scan conversion.
If provided, it will be used instead of computing new coordinates based on
the input image shape and ranges.
fill_value (float, optional): The value to fill in for coordinates
outside the input image ranges. Defaults to 0.0. When set to NaN,
no interpolation at the edges will happen.
order (int, optional): The order of the spline interpolation. Defaults to 1.
grid_apex (array-like, optional): The point ``(x, y, z)`` that the radius is
measured from, in the transducer frame (see
:attr:`zea.Parameters.grid_apex`). Moves the reported ``x_lim`` and ``z_lim``
into the transducer frame. Defaults to the origin.
extent (array-like, optional): Output window in the transducer frame,
``(xmin, xmax, ..., zmin, zmax)`` as :attr:`zea.Parameters.extent` gives it.
Defaults to the bounding box of the polar image.
distance_to_apex (float, optional): Deprecated, use
``grid_apex=(0, 0, -distance_to_apex)``.
rho_range (tuple, optional): Deprecated, use ``rlims``.
theta_range (tuple, optional): Deprecated, use ``polar_limits``.
Returns:
ndarray: The scan-converted 2D ultrasound image in Cartesian coordinates.
Has dimensions (grid_size_z, grid_size_x). Coordinates outside the input image
ranges are filled with NaNs.
parameters (dict): A dictionary containing information about the scan conversion.
Contains the resolution, x and z limits, ``rlims`` and ``polar_limits``.
Note:
Polar grid is inferred from the input image shape and the supplied
``rlims`` and ``polar_limits``. Cartesian grid is computed based on polar grid
with resolutions specified by resolution parameter.
"""
assert "float" in ops.dtype(image), "Image must be float type"
rlims = deprecated_alias(rlims, rho_range, "rho_range", "rlims", "scan_convert_2d")
polar_limits = deprecated_alias(
polar_limits, theta_range, "theta_range", "polar_limits", "scan_convert_2d"
)
parameters = {}
if coordinates is None:
if rlims is None or polar_limits is None:
raise ValueError(
"rlims and polar_limits are required when coordinates is not provided."
)
coordinates, parameters = compute_scan_convert_2d_coordinates(
image.shape,
rlims,
polar_limits,
resolution,
dtype=image.dtype,
grid_apex=deprecated_distance_to_apex(distance_to_apex, grid_apex, "scan_convert_2d"),
extent=extent,
)
scan_converted = _interpolate_batch(image, coordinates, fill_value, order=order, **kwargs)
# swap axis to match z, x
scan_converted = ops.swapaxes(scan_converted, -1, -2)
return scan_converted, parameters
[docs]
def compute_scan_convert_3d_coordinates(
image_shape,
rlims: Tuple[float, float] | None = None,
polar_limits: Tuple[float, float] | None = None,
elevation_limits: Tuple[float, float] | None = None,
resolution: Union[float, None] = None,
dtype: str = "float32",
*,
rho_range=None,
theta_range=None,
phi_range=None,
):
"""Precompute coordinates for 3d scan conversion from frustum coordinates.
The volume is indexed by radius, the polar angle ``atan(x / z)`` and the elevation
angle ``atan(y / z)`` (see :func:`frustum_to_xyz`). Both angles are
projections onto the x-z and y-z planes, not the spherical ``polar_angles`` and
``azimuth_angles`` of the transmits; at elevation 0 the volume reduces to the 2D sector
of :func:`scan_convert_2d`.
Args:
image_shape (tuple): Shape of the volume, ``(..., n_radius, n_polar, n_elevation)``.
rlims (tuple): Radius of the first and last sample along the radius axis.
polar_limits (tuple): Polar angle ``atan(x / z)`` of the first and last sample
along the polar axis [rad], ascending.
elevation_limits (tuple): Elevation angle ``atan(y / z)`` of the first and last
sample along the elevation axis [rad], ascending.
resolution (float, optional): Output voxel size. Derived from the volume when
not given.
dtype (str, optional): Dtype of the coordinates. Defaults to ``"float32"``.
rho_range (tuple, optional): Deprecated, use ``rlims``.
theta_range (tuple, optional): Deprecated, use ``polar_limits``.
phi_range (tuple, optional): Deprecated, use ``elevation_limits``.
"""
caller = "compute_scan_convert_3d_coordinates"
rlims = deprecated_alias(rlims, rho_range, "rho_range", "rlims", caller)
polar_limits = deprecated_alias(
polar_limits, theta_range, "theta_range", "polar_limits", caller
)
elevation_limits = deprecated_alias(
elevation_limits, phi_range, "phi_range", "elevation_limits", caller
)
if rlims is None or polar_limits is None or elevation_limits is None:
raise ValueError(f"{caller}() requires rlims, polar_limits and elevation_limits.")
assert len(rlims) == 2, "rlims should be a tuple of length 2"
assert len(polar_limits) == 2, "polar_limits should be a tuple of length 2"
assert len(elevation_limits) == 2, "elevation_limits should be a tuple of length 2"
assert rlims[0] < rlims[1], "rlims[0] should be less than rlims[1]"
# Unlike 2D, the 3D masking and indexing assume ascending angles
assert polar_limits[0] < polar_limits[1], "polar_limits[0] should be less than [1]"
assert elevation_limits[0] < elevation_limits[1], (
"elevation_limits[0] should be less than elevation_limits[1]"
)
rho = ops.linspace(rlims[0], rlims[1], image_shape[-3], dtype=dtype)
theta = ops.linspace(polar_limits[0], polar_limits[1], image_shape[-2], dtype=dtype)
phi = ops.linspace(elevation_limits[0], elevation_limits[1], image_shape[-1], dtype=dtype)
rho_grid, theta_grid, phi_grid = ops.meshgrid(rho, theta, phi, indexing="ij")
x_grid, y_grid, z_grid = frustum_to_xyz(rho_grid, theta_grid, phi_grid)
x_lim = [ops.min(x_grid), ops.max(x_grid)]
y_lim = [ops.min(y_grid), ops.max(y_grid)]
z_lim = [ops.min(z_grid), ops.max(z_grid)]
d_rho = rho[1] - rho[0]
d_theta = theta[1] - theta[0]
d_phi = phi[1] - phi[0]
if resolution is None:
# arc length along constant phi at 1/4 depth
sRT = 0.25 * (rho[0] + rho[-1]) * d_theta
# arc length along constant theta at 1/4 depth
sRP = 0.25 * (rho[0] + rho[-1]) * d_phi
# average of arc lengths and radial step
resolution = ops.mean([sRT, sRP, d_rho]) # mm per pixel
z_vec = ops.arange(z_lim[0], z_lim[1], resolution)
x_vec = ops.arange(x_lim[0], x_lim[1], resolution)
y_vec = ops.arange(y_lim[0], y_lim[1], resolution)
z_grid, x_grid, y_grid = ops.meshgrid(z_vec, x_vec, y_vec)
rho_grid_interp, theta_grid_interp, phi_grid_interp = xyz_to_frustum(
x_grid,
y_grid,
z_grid,
polar_limits=[theta[0], theta[-1]],
elevation_limits=[phi[0], phi[-1]],
)
# return volume
rho_min, rho_max = ops.min(rho), ops.max(rho)
theta_min, theta_max = ops.min(theta), ops.max(theta)
phi_min, phi_max = ops.min(phi), ops.max(phi)
rho_idx = (rho_grid_interp - rho_min) / (rho_max - rho_min) * (image_shape[-3] - 1)
theta_idx = (theta_grid_interp - theta_min) / (theta_max - theta_min) * (image_shape[-2] - 1)
phi_idx = (phi_grid_interp - phi_min) / (phi_max - phi_min) * (image_shape[-1] - 1)
# Stack coordinates as required for map_coordinates
coordinates = ops.stack([rho_idx, theta_idx, phi_idx], axis=0)
parameters = {
"resolution": resolution,
"x_lim": x_lim,
"y_lim": y_lim,
"z_lim": z_lim,
"rlims": rlims,
"polar_limits": polar_limits,
"elevation_limits": elevation_limits,
"dr": d_rho,
"d_polar": d_theta,
"d_elevation": d_phi,
}
return coordinates, parameters
[docs]
def scan_convert_3d(
image,
rlims: Tuple[float, float] | None = None,
polar_limits: Tuple[float, float] | None = None,
elevation_limits: Tuple[float, float] | None = None,
resolution: Union[float, None] = None,
coordinates: Union[None, np.ndarray] = None,
fill_value: float = 0.0,
order: int = 1,
*,
rho_range=None,
theta_range=None,
phi_range=None,
):
"""
Perform scan conversion on a 3D ultrasound image from frustum coordinates
(radius, polar angle, elevation angle) to Cartesian coordinates (z, x, y).
The polar angle ``atan(x / z)`` and the elevation angle ``atan(y / z)`` (see
:func:`frustum_to_xyz`) are projections onto the x-z and y-z planes, not the
spherical ``polar_angles`` and ``azimuth_angles`` of the transmits.
Args:
image (ndarray): The input 3D ultrasound image in frustum coordinates.
Has dimensions (n_radius, n_polar, n_elevation) with optional batch.
rlims (tuple): Radius ``(rho_min, rho_max)`` of the first and last sample.
polar_limits (tuple): Polar angle ``atan(x / z)`` of the first and last sample
[rad].
elevation_limits (tuple): Elevation angle ``atan(y / z)`` of the first and last
sample [rad].
resolution (float, optional): The resolution for the Cartesian grid.
If None, it is calculated based on the input image. In units of ``rlims``
per pixel.
coodinates (ndarray, optional): Precomputed coordinates for scan conversion.
If provided, it will be used instead of computing new coordinates based on
the input image shape and ranges.
fill_value (float, optional): The value to fill in for coordinates
outside the input image ranges. Defaults to 0.0. When set to NaN,
no interpolation at the edges will happen.
order (int, optional): The order of the spline interpolation. Defaults to 1.
rho_range (tuple, optional): Deprecated, use ``rlims``.
theta_range (tuple, optional): Deprecated, use ``polar_limits``.
phi_range (tuple, optional): Deprecated, use ``elevation_limits``.
Returns:
ndarray: The scan-converted 3D ultrasound image in Cartesian coordinates.
Has dimensions (grid_size_z, grid_size_x, n_y). Coordinates outside the input image
ranges are filled with NaNs.
parameters (dict): A dictionary containing information about the scan conversion.
Contains the resolution, x, y and z limits, ``rlims``, ``polar_limits`` and
``elevation_limits``.
Note:
Frustum grid is inferred from the input image shape and the supplied limits.
Cartesian grid is computed based on the frustum grid with resolutions specified by
resolution parameter.
"""
assert "float" in ops.dtype(image), "Image must be float type"
caller = "scan_convert_3d"
rlims = deprecated_alias(rlims, rho_range, "rho_range", "rlims", caller)
polar_limits = deprecated_alias(
polar_limits, theta_range, "theta_range", "polar_limits", caller
)
elevation_limits = deprecated_alias(
elevation_limits, phi_range, "phi_range", "elevation_limits", caller
)
parameters = {}
if coordinates is None:
if rlims is None or polar_limits is None or elevation_limits is None:
raise ValueError(
"rlims, polar_limits and elevation_limits are required "
"when coordinates is not provided."
)
coordinates, parameters = compute_scan_convert_3d_coordinates(
image.shape,
rlims,
polar_limits,
elevation_limits,
resolution,
dtype=image.dtype,
)
scan_converted = _interpolate_batch(image, coordinates, fill_value, order=order)
# swap axis to match z, x, y
scan_converted = ops.swapaxes(scan_converted, -2, -3)
return scan_converted, parameters
[docs]
def scan_convert(
image,
rlims: Tuple[float, float] | None = None,
polar_limits: Tuple[float, float] | None = None,
elevation_limits: Tuple[float, float] | None = None,
resolution: Union[float, None] = None,
coordinates: Union[None, np.ndarray] = None,
fill_value: float = 0.0,
order: int = 1,
with_batch_dim: bool = False,
distance_to_apex=None,
*,
grid_apex=None,
extent=None,
rho_range=None,
theta_range=None,
phi_range=None,
):
"""Scan convert image based on number of dimensions.
A 2D image goes to :func:`scan_convert_2d` with ``rlims`` and ``polar_limits``; a 3D
volume goes to :func:`scan_convert_3d`, which adds ``elevation_limits``. The polar
angle ``atan(x / z)`` means the same in both, so a 3D frustum extends the 2D sector in
elevation. ``grid_apex`` and ``extent`` apply to 2D images only.
``rho_range``, ``theta_range`` and ``phi_range`` are deprecated aliases of ``rlims``,
``polar_limits`` and ``elevation_limits``.
"""
rlims = deprecated_alias(rlims, rho_range, "rho_range", "rlims", "scan_convert")
polar_limits = deprecated_alias(
polar_limits, theta_range, "theta_range", "polar_limits", "scan_convert"
)
elevation_limits = deprecated_alias(
elevation_limits, phi_range, "phi_range", "elevation_limits", "scan_convert"
)
if len(image.shape) == 2 + int(with_batch_dim):
return scan_convert_2d(
image,
rlims,
polar_limits,
resolution,
coordinates,
fill_value,
order,
grid_apex=deprecated_distance_to_apex(distance_to_apex, grid_apex, "scan_convert"),
extent=extent,
)
elif len(image.shape) == 3 + int(with_batch_dim):
return scan_convert_3d(
image,
rlims,
polar_limits,
elevation_limits,
resolution,
coordinates,
fill_value,
order,
)
else:
raise ValueError(
f"Image must be 2D or 3D (with optional batch dim). Got shape: {image.shape}"
)
[docs]
def map_coordinates(inputs, coordinates, order, fill_mode="constant", fill_value=0):
"""map_coordinates using keras.ops or scipy.ndimage when order > 1."""
if order > 1:
# Preserve original dtype before conversion
original_dtype = ops.dtype(inputs)
inputs_np = ops.convert_to_numpy(inputs).astype(np.float32)
coordinates_np = ops.convert_to_numpy(coordinates).astype(np.float32)
out = scipy.ndimage.map_coordinates(
inputs_np, coordinates_np, order=order, mode=fill_mode, cval=fill_value
)
return ops.convert_to_tensor(out.astype(original_dtype))
else:
return ops.image.map_coordinates(
inputs,
coordinates,
order=order,
fill_mode=fill_mode,
fill_value=fill_value,
)
def _interpolate_batch(images, coordinates, fill_value=0.0, order=1, vectorize=True):
"""Interpolate a batch of images."""
image_shape = images.shape
num_image_dims = coordinates.shape[0]
batch_dims = images.shape[:-num_image_dims]
images = ops.reshape(images, (-1, *image_shape[-num_image_dims:]))
map_coordinates_fn = partial(
map_coordinates,
coordinates=coordinates,
order=order,
fill_mode="constant",
fill_value=fill_value,
)
if order > 1:
# cpu bound
scan_converted = ops.stack(list(map(map_coordinates_fn, images)))
elif not vectorize:
scan_converted = ops.map(map_coordinates_fn, images)
else:
# gpu bound
scan_converted = ops.vectorized_map(map_coordinates_fn, images)
# ignore batch dim to get image shape
scan_converted_shape = ops.shape(scan_converted)[1:]
scan_converted = ops.reshape(scan_converted, (*batch_dims, *scan_converted_shape))
return scan_converted
[docs]
def cartesian_to_polar_matrix(
cartesian_matrix,
fill_value=0.0,
polar_shape=None,
tip=None,
r_max=None,
theta_range=None,
interpolation_order=1,
*,
polar_limits=None,
):
"""
Convert a Cartesian image matrix to a polar coordinate representation.
Args:
cartesian_matrix (tensor): Input 2D image array in Cartesian coordinates.
fill_value (float): Value to use for points sampled outside the input image.
polar_shape (tuple, optional): Desired shape of the polar output (rows, cols).
Defaults to the shape of the input image.
tip (tuple, optional): (x, y) coordinates of the origin for the polar
transformation (typically the probe tip). Defaults to the center-top of the image.
r_max (float, optional): Maximum radius to consider in the polar transform.
Defaults to the height of the input image.
polar_limits (tuple, optional): Polar angle ``atan(x / z)`` [rad] of the first and
last column of the polar output, positive toward +x (increasing column index of
the Cartesian image), as in :func:`scan_convert_2d`. Unequal magnitudes describe
an asymmetric cone. Defaults to (-45, 45) degrees.
interpolation_order (int): Order of interpolation to use (0 = nearest-neighbor,
1 = linear, 2+ = spline). Matches the convention of `scipy.ndimage.map_coordinates`.
theta_range (tuple, optional): Deprecated, use ``polar_limits``. It measured angles
toward -x, so ``theta_range=(a, b)`` is ``polar_limits=(-b, -a)``.
Returns:
polar_matrix (Array): The image re-sampled in polar coordinates with shape
`polar_shape`.
coordinates (Array): The Cartesian coordinates corresponding to each pixel in
the polar output.
"""
assert "float" in ops.dtype(cartesian_matrix), "Input image must be float type"
polar_limits = _deprecated_mirrored_theta_range(
polar_limits, theta_range, "cartesian_to_polar_matrix", reverse=True
)
if polar_limits is None:
polar_limits = (-np.deg2rad(45), np.deg2rad(45))
# Assume that polar grid is same shape as cartesian grid unless specified
cartesian_rows, cartesian_cols = ops.shape(cartesian_matrix)
if polar_shape is None:
polar_rows, polar_cols = cartesian_rows, cartesian_cols
else:
polar_rows, polar_cols = polar_shape
# assume tip is at center top unless specified
if tip is None:
center_x = cartesian_cols / 2 # center_x can be between two pixels
tip_y = 0
tip = (center_x, tip_y)
# assume r_max is the total height of the input image unless specified
if r_max is None:
r_max = cartesian_rows
center_x, center_y = tip
# Polar sampling grid: rows are radius (0..r_max), columns are angle, from polar_limits[0]
# to polar_limits[1] as in the polar grid of zea.
r = ops.linspace(0, r_max, polar_rows, dtype="float32")
theta = ops.linspace(polar_limits[0], polar_limits[1], polar_cols, dtype="float32")
r_grid, theta_grid = ops.meshgrid(r, theta, indexing="ij")
# Cartesian pixel each (r, theta) samples from. The probe images "downwards" (+depth = +y),
# so depth = center_y + r cos(theta) and lateral = center_x + r sin(theta).
xq = center_x + r_grid * ops.sin(theta_grid)
yq = center_y + r_grid * ops.cos(theta_grid)
coords_for_interp = ops.stack([ops.ravel(yq), ops.ravel(xq)])
polar_values = map_coordinates(
cartesian_matrix,
coords_for_interp,
order=interpolation_order,
fill_mode="constant",
fill_value=fill_value,
)
polar_matrix = ops.reshape(polar_values, (polar_rows, polar_cols))
return polar_matrix
[docs]
def polar_to_cartesian_matrix(
polar_matrix,
cartesian_shape: Tuple[int, int],
tip: Union[Tuple[float, float], None] = None,
r_max: Union[float, None] = None,
theta_range=None,
pitch: float = 1.0,
fill_value: float = 0.0,
order: int = 1,
*,
polar_limits: Union[Tuple[float, float], None] = None,
):
"""Resample a polar image onto a fixed Cartesian canvas with the apex at a chosen pixel.
Faithful inverse of :func:`cartesian_to_polar_matrix`. Where :func:`scan_convert_2d`
fits the cone bounding box into the output (losing absolute position and scale), this
pins the cone apex at pixel ``tip`` on a canvas of shape ``cartesian_shape`` and samples
at ``pitch`` units per pixel, so a forward/inverse round-trip reproduces the original
frame. Both share the :func:`_polar_sampling_coordinates` core; they differ only in how
the output grid is built.
Args:
polar_matrix (tensor): Input polar image of shape ``(n_r, n_polar)``, float type.
cartesian_shape (tuple): Output ``(rows, cols)`` = ``(n_z, n_x)``.
tip (tuple, optional): ``(col, row)`` pixel location of the cone apex in the output.
Defaults to the centre-top ``(cols / 2, 0)``.
r_max (float, optional): Radius spanned by the polar image, in the same units as
``pitch`` (pixels by default). Defaults to ``rows``.
polar_limits (tuple, optional): Polar angle ``atan(x / z)`` [rad] of the first and
last column of ``polar_matrix``, positive toward +x (increasing column index of
the output). Pass the ``polar_limits`` given to :func:`cartesian_to_polar_matrix`
to invert it. Defaults to (-45, 45) degrees.
pitch (float, optional): Output units per pixel (radial units of ``r_max``).
Defaults to 1.0, i.e. ``tip`` and ``r_max`` are in pixels.
fill_value (float, optional): Value for pixels outside the polar domain.
order (int, optional): Spline interpolation order. Defaults to 1.
theta_range (tuple, optional): Deprecated, use ``polar_limits``. It measured angles
toward -x, so ``theta_range=(a, b)`` is ``polar_limits=(-a, -b)``.
Returns:
cartesian_matrix (Array): The polar image resampled onto the Cartesian canvas, with
shape ``cartesian_shape``.
"""
assert "float" in ops.dtype(polar_matrix), "Input image must be float type"
polar_limits = _deprecated_mirrored_theta_range(
polar_limits, theta_range, "polar_to_cartesian_matrix", reverse=False
)
dtype = ops.dtype(polar_matrix)
n_r, n_polar = ops.shape(polar_matrix)[-2], ops.shape(polar_matrix)[-1]
coordinates = polar_to_cartesian_coordinates(
cartesian_shape,
n_r,
n_polar,
tip=tip,
r_max=r_max,
polar_limits=polar_limits,
pitch=pitch,
dtype=dtype,
)
cartesian, _ = scan_convert_2d(
polar_matrix, coordinates=coordinates, fill_value=fill_value, order=order
)
return cartesian
[docs]
def polar_to_cartesian_coordinates(
cartesian_shape: Tuple[int, int],
n_r: int,
n_polar: int,
tip: Union[Tuple[float, float], None] = None,
r_max: Union[float, None] = None,
theta_range=None,
pitch: float = 1.0,
dtype: str = "float32",
*,
polar_limits: Union[Tuple[float, float], None] = None,
):
"""(radius, polar angle) sampling grid that pins a polar cone's apex at pixel ``tip``.
Shared core of :func:`polar_to_cartesian_matrix`: builds the fractional index grid that
resamples an ``(n_r, n_polar)`` polar image onto a ``cartesian_shape`` canvas, with the
apex at pixel ``tip`` and ``pitch`` units per pixel. Pass the result as ``coordinates`` to
:func:`scan_convert_2d` (or the :class:`~zea.ops.ScanConvert` operation) to run the same
resampling inside a pipeline. Unlike :func:`compute_scan_convert_2d_coordinates`, this uses
the ``arctan2`` mapping in :func:`_polar_sampling_coordinates`, so it covers a full circle
rather than only a narrow (polar angle within ``pi/2``) sector.
Args:
cartesian_shape (tuple): Output ``(rows, cols)`` = ``(n_z, n_x)``.
n_r (int): Radial samples (rows) of the polar image the indices address.
n_polar (int): Angular samples (columns) of the polar image the indices address.
tip (tuple, optional): ``(col, row)`` pixel location of the cone apex in the output.
Defaults to the centre-top ``(cols / 2, 0)``.
r_max (float, optional): Radius spanned by the polar image, in the same units as
``pitch`` (pixels by default). Defaults to ``rows``.
polar_limits (tuple, optional): Polar angle ``atan(x / z)`` [rad] of the first and
last column of the polar image, positive toward +x (increasing column index of
the output). Defaults to (-45, 45) degrees.
pitch (float, optional): Output units per pixel (radial units of ``r_max``).
Defaults to 1.0, i.e. ``tip`` and ``r_max`` are in pixels.
dtype (str, optional): Compute dtype for the axes. Defaults to ``"float32"``.
theta_range (tuple, optional): Deprecated, use ``polar_limits``. It measured angles
toward -x, so ``theta_range=(a, b)`` is ``polar_limits=(-a, -b)``.
Returns:
coords (Array): ``(2, cols, rows)`` stack of ``[rho_idx, theta_idx]`` sampling indices.
"""
cart_rows, cart_cols = cartesian_shape
if tip is None:
tip = (cart_cols / 2, 0.0)
if r_max is None:
r_max = cart_rows
polar_limits = _deprecated_mirrored_theta_range(
polar_limits, theta_range, "polar_to_cartesian_coordinates", reverse=False
)
if polar_limits is None:
polar_limits = (-np.deg2rad(45), np.deg2rad(45))
tip_x, tip_y = tip
rho = ops.linspace(0.0, r_max, n_r, dtype=dtype)
theta = ops.linspace(polar_limits[0], polar_limits[1], n_polar, dtype=dtype)
x_vec = (ops.cast(ops.arange(cart_cols), dtype) - tip_x) * pitch
z_vec = (ops.cast(ops.arange(cart_rows), dtype) - tip_y) * pitch
return _polar_sampling_coordinates(x_vec, z_vec, rho, theta)
[docs]
def inverse_scan_convert_2d(
cartesian_image,
fill_value=0.0,
theta_range=None,
output_size=None,
interpolation_order=1,
find_scan_cone=True,
image_range: tuple | None = None,
*,
polar_limits=None,
):
"""
Convert a Cartesian-format ultrasound image to a polar representation.
This function can be used to recover a sector-shaped scan (polar format)
from a Cartesian representation of an image.
Optionally, it can detect and crop around the scan cone before conversion.
Args:
cartesian_image (tensor): 2D image array in Cartesian coordinates of type float.
fill_value (float): Value used to fill regions outside the original image
during interpolation.
polar_limits (tuple, optional): Polar angle ``atan(x / z)`` [rad] of the first and
last column of the polar image, positive toward +x, as in :func:`scan_convert_2d`.
Unequal magnitudes describe an asymmetric cone. Defaults to (-45, 45) degrees.
output_size (tuple, optional): Shape (rows, cols) of the resulting polar image.
If None, the shape of the input image is used.
interpolation_order (int): Order of interpolation used in resampling
(0 = nearest-neighbor, 1 = linear, etc.).
find_scan_cone (bool): If True, automatically detects and crops around the scan cone
in the Cartesian image before polar conversion, ensuring that the scan cone is
centered without padding. Can be set to False if the image is already cropped
and centered.
image_range (tuple, optional): Tuple (vmin, vmax) for display scaling
when detecting the scan cone.
theta_range (tuple, optional): Deprecated, use ``polar_limits``. It measured angles
toward -x, so ``theta_range=(a, b)`` is ``polar_limits=(-b, -a)``.
Returns:
polar_image (Array): 2D image in polar coordinates (sector-shaped scan).
"""
polar_limits = _deprecated_mirrored_theta_range(
polar_limits, theta_range, "inverse_scan_convert_2d", reverse=True
)
if find_scan_cone:
assert image_range is not None, "image_range must be provided when find_scan_cone is True"
cartesian_image, _ = fit_and_crop_around_scan_cone(cartesian_image, image_range)
polar_image = cartesian_to_polar_matrix(
cartesian_image,
fill_value=fill_value,
polar_limits=polar_limits,
polar_shape=output_size,
interpolation_order=interpolation_order,
)
return polar_image
[docs]
def frustum_to_xyz(radius, polar, elevation):
"""Convert frustum coordinates (radius, polar angle, elevation angle) to (x, y, z).
The polar angle ``atan(x / z)`` and the elevation angle ``atan(y / z)`` are projections
onto the x-z and y-z planes, not the spherical ``polar_angles`` and ``azimuth_angles``
of the transmits (see :ref:`angle-conventions`). Angles are defined in radians.
Args:
radius (ndarray): Radial coordinates of the points to convert.
polar (ndarray): Polar angles ``atan(x / z)`` of the points to convert.
elevation (ndarray): Elevation angles ``atan(y / z)`` of the points to convert.
Returns:
x (ndarray): X coordinates of the converted points.
y (ndarray): Y coordinates of the converted points.
z (ndarray): Z coordinates of the converted points.
"""
if ops.size(radius) != ops.size(polar) or ops.size(radius) != ops.size(elevation):
raise ValueError("Number of elements in radius, polar and elevation should be the same")
z = radius / ops.sqrt(1 + ops.tan(polar) ** 2 + ops.tan(elevation) ** 2)
x = z * ops.tan(polar)
y = z * ops.tan(elevation)
return x, y, z
[docs]
def frustum_to_xz(radius, polar):
"""Convert polar coordinates (radius, polar angle) to (x, z): the 2D case of
:func:`frustum_to_xyz`, at elevation 0.
The polar angle ``atan(x / z)`` is signed, positive toward +x. Angles are defined in
radians.
Args:
radius (ndarray): Radial coordinates of the points to convert.
polar (ndarray): Polar angles ``atan(x / z)`` of the points to convert.
Returns:
x (ndarray): X coordinates of the converted points.
z (ndarray): Z coordinates of the converted points.
"""
if ops.size(radius) != ops.size(polar):
raise ValueError("Number of elements in radius and polar should be the same")
z = radius / ops.sqrt(1 + ops.tan(polar) ** 2)
x = z * ops.tan(polar)
return x, z
[docs]
def xz_to_frustum(x, z, polar_limits):
"""Convert (x, z) to polar coordinates (radius, polar angle), the inverse of
:func:`frustum_to_xz`.
The polar angle ``atan(x / z)`` is signed, positive toward +x. Angles are defined in
radians.
Args:
x (ndarray): X coordinates of the points to convert.
z (ndarray): Z coordinates of the points to convert.
polar_limits (list): Polar angle limits of the original image. Any
point that resides outside of these limits is potentially
undefined, and therefore, the radial value for these points is
made to be -1.
Returns:
radius (ndarray): Radial coordinates of the converted points.
polar (ndarray): Polar angles ``atan(x / z)`` of the converted points.
"""
if ops.size(x) != ops.size(z):
raise ValueError("Number of elements in x and z should be the same")
radius = ops.sqrt(x**2 + z**2)
polar = ops.arctan2(x, z)
radius = ops.where(
(radius < 0) | (polar < polar_limits[0]) | (polar > polar_limits[1]),
-1,
radius,
)
return radius, polar
[docs]
def xyz_to_frustum(x, y, z, polar_limits, elevation_limits):
"""Convert (x, y, z) to frustum coordinates (radius, polar angle, elevation angle), the
inverse of :func:`frustum_to_xyz`.
The polar angle ``atan(x / z)`` and the elevation angle ``atan(y / z)`` are projections
onto the x-z and y-z planes (see :ref:`angle-conventions`). Angles are defined in
radians.
Args:
x (ndarray): X coordinates of the points to convert.
y (ndarray): Y coordinates of the points to convert.
z (ndarray): Z coordinates of the points to convert.
polar_limits, elevation_limits:
Polar and elevation angle limits, respectively, of the original volume. Any
point that resides outside of these limits is potentially
undefined, and therefore, the radial value for these points is
made to be -1.
Returns:
radius (ndarray): Radial coordinates of the converted points.
polar (ndarray): Polar angles ``atan(x / z)`` of the converted points.
elevation (ndarray): Elevation angles ``atan(y / z)`` of the converted points.
"""
if ops.size(x) != ops.size(y) or ops.size(x) != ops.size(z):
raise ValueError("Number of elements in x, y, and z should be the same")
radius = ops.sqrt(x**2 + y**2 + z**2)
polar = ops.arctan2(x, z)
elevation = ops.arctan2(y, z)
radius = ops.where(
(radius < 0)
| (polar < polar_limits[0])
| (polar > polar_limits[1])
| (elevation < elevation_limits[0])
| (elevation > elevation_limits[1]),
-1,
radius,
)
return radius, polar, elevation
[docs]
@deprecated(replacement="frustum_to_xyz")
def frustum_convert_rtp2xyz(rho, theta, phi):
"""Deprecated, use :func:`frustum_to_xyz` (``theta`` is its ``polar`` and ``phi`` its
``elevation`` angle)."""
return frustum_to_xyz(rho, theta, phi)
[docs]
@deprecated(replacement="frustum_to_xz")
def frustum_convert_rt2xz(rho, theta):
"""Deprecated, use :func:`frustum_to_xz` (``theta`` is its ``polar`` angle)."""
return frustum_to_xz(rho, theta)
[docs]
@deprecated(replacement="xz_to_frustum")
def frustum_convert_xz2rt(x, z, theta_limits):
"""Deprecated, use :func:`xz_to_frustum` (``theta_limits`` is its ``polar_limits``)."""
return xz_to_frustum(x, z, theta_limits)
[docs]
@deprecated(replacement="xyz_to_frustum")
def frustum_convert_xyz2rtp(x, y, z, theta_limits, phi_limits):
"""Deprecated, use :func:`xyz_to_frustum` (``theta_limits`` is its ``polar_limits`` and
``phi_limits`` its ``elevation_limits``)."""
return xyz_to_frustum(x, y, z, theta_limits, phi_limits)