from __future__ import annotations
import numpy as np
from ._utils import _filters, _structures
from ._utils._skimage import skmes
from .specials import PropArray
from .labeledarray import LabeledArray
from ..utils.axesop import complement_axes
from ..utils.deco import check_input_and_output, dims_to_spatial_axes, same_dtype
from ._utils import _misc
from ..collections import DataDict
from .._types import Dims
def _calc_phase_mean(sl, img, periodicity):
a = 2 * np.pi / periodicity
out = np.sum(np.exp(1j*a*img[sl]))
return np.angle(out)/a
def _calc_phase_std(sl, img, periodicity):
a = 2 * np.pi / periodicity
return np.sqrt(-2*np.log(np.abs(np.mean(np.exp(1j*a*img[sl])))))/a
[docs]class PhaseArray(LabeledArray):
additional_props = ["_source", "_metadata", "_name", "unit", "border"]
def __new__(cls, obj, name=None, axes=None, source=None,
metadata=None, dtype=None, unit="rad", border=None) -> PhaseArray:
if dtype is None:
dtype = np.float32
if border is None:
border = {"rad": (0, 2*np.pi), "deg": (0, 360.0)}[unit]
self = super().__new__(cls, obj, name=name, axes=axes, source=source,
metadata=metadata, dtype=dtype)
self.unit = unit
self.border = border
return self
def __add__(self, value) -> PhaseArray:
out = super().__add__(value)
out.fix_border()
return out
def __iadd__(self, value) -> PhaseArray:
out = super().__iadd__(value)
out.fix_border()
return out
def __sub__(self, value) -> PhaseArray:
out = super().__sub__(value)
out.fix_border()
return out
def __isub__(self, value) -> PhaseArray:
out = super().__isub__(value)
out.fix_border()
return out
def __mul__(self, value) -> PhaseArray:
out = super().__mul__(value)
out.fix_border()
return out
def __imul__(self, value) -> PhaseArray:
out = super().__imul__(value)
out.fix_border()
return out
def __truediv__(self, value) -> PhaseArray:
out = super().__truediv__(value)
out.fix_border()
return out
def __itruediv__(self, value) -> PhaseArray:
out = super().__itruediv__(value)
out.fix_border()
return out
@property
def periodicity(self) -> float:
a, b = self.border
return b - a
[docs] def fix_border(self) -> None:
"""
Considering periodic boundary condition, fix the values by `__mod__` method.
"""
self[:] = (self.value - self.border[0]) % self.periodicity + self.border[0]
return None
[docs] def set_border(self, a: float, b: float) -> None:
"""
Set new border safely.
Parameters
----------
a : float
New lower border.
b : float
New higher border.
"""
if not (np.isscalar(a) and np.isscalar(b)):
raise TypeError("Both border values must be scalars.")
if b - a != self.periodicity:
raise ValueError("New border does not match current periodicity.")
self.border = (b, a)
self.fix_border()
return None
[docs] def deg2rad(self) -> PhaseArray:
if self.unit == "rad":
raise ValueError("Array is already in radian.")
np.deg2rad(self, out=self.value[:])
self.unit = "rad"
self.border = tuple(np.deg2rad(self.border))
return self
[docs] def rad2deg(self) -> PhaseArray:
if self.unit == "deg":
raise ValueError("Array is already in degree.")
np.rad2deg(self, out=self.value[:])
self.unit = "deg"
self.border = tuple(np.rad2deg(self.border))
return self
[docs] @check_input_and_output
@dims_to_spatial_axes
@same_dtype(asfloat=True)
def binning(self, binsize: int = 2, *, check_edges: bool = True, dims: Dims = None):
if binsize == 1:
return self
img_to_reshape, shape, scale_ = _misc.adjust_bin(self.value, binsize, check_edges, dims, self.axes)
reshaped_img = img_to_reshape.reshape(shape)
axes_to_reduce = tuple(i*2+1 for i in range(self.ndim))
a = 2 * np.pi / self.periodicity
out = np.sum(np.exp(1j*a*reshaped_img), axis=axes_to_reduce)
out = np.angle(out)/a
out: PhaseArray = out.view(self.__class__)
out._set_info(self)
out.axes = str(self.axes) # _set_info does not pass copy so new axes must be defined here.
out.set_scale({a: self.scale[a]/scale for a, scale in zip(self.axes, scale_)})
return out
[docs] @check_input_and_output
@dims_to_spatial_axes
@same_dtype(asfloat=True)
def mean_filter(self, radius: float = 1, *, dims: Dims = None, update: bool = False) -> PhaseArray:
r"""
Mean filter using phase averaging method:
:math:`\arg{\sum_{k}{e^{i X_k}}}`
Parameters
----------
radius : float, default is 1
Radius of kernel.
dims : str or int, optional
Spatial dimensions.
update : bool, default is False
If update self to filtered image.
Returns
-------
PhaseArray
Filtered image.
"""
disk = _structures.ball_like(radius, len(dims))
a = 2*np.pi/self.periodicity
return self._apply_dask(_filters.phase_mean_filter,
c_axes=complement_axes(dims, self.axes),
args=(disk, a),
dtype=self.dtype)
[docs] def imshow(self, dims="yx", **kwargs):
if "cmap" not in kwargs:
kwargs["cmap"] = "hsv"
return super().imshow(dims=dims, **kwargs)
[docs] def reslice(self, src, dst, *, order: int = 1) -> PropArray:
"""
Measure line profile iteratively for every slice of image. Because input is phase, we can
not apply standard interpolation to calculate intensities on float-coordinates.
Parameters
----------
src : array-like
Source coordinate.
dst : array-like
Destination coordinate.
order : int, default is 1
Spline interpolation order.
dims : int or str, optional
Spatial dimensions.
Returns
-------
PropArray
Line scans.
"""
a = 2*np.pi/self.periodicity
vec_re = np.cos(a*self).view(LabeledArray)
vec_im = np.sin(a*self).view(LabeledArray)
out_re = vec_re.reslice(src, dst, order=order)
out_im = vec_im.reslice(src, dst, order=order)
out = np.arctan2(out_im, out_re)/a
return out
[docs] @check_input_and_output(need_labels=True)
def regionprops(self, properties: tuple[str,...]|str = ("phase_mean",), *,
extra_properties = None) -> DataDict[str, PropArray]:
"""
Run ``skimage``'s ``regionprops()`` function and return the results as PropArray, so
that you can access using flexible slicing. For example, if a tcyx-image is
analyzed with ``properties=("X", "Y")``, then you can get X's time-course profile
of channel 1 at label 3 by ``prop["X"]["p=5;c=1"]`` or ``prop.X["p=5;c=1"]``.
In PhaseArray, instead of mean_intensity you should use "phase_mean". The
phase_mean function is included so that it can be passed in ``properties`` argument.
Parameters
----------
properties : iterable, optional
properties to analyze, see skimage.measure.regionprops.
extra_properties : iterable of callable, optional
extra properties to analyze, see skimage.measure.regionprops.
Returns
-------
DataDict of PropArray
Example
-------
Measure region properties around single molecules.
>>> coords = reference_img.centroid_sm()
>>> img.specify(coords, 3, labeltype="circle")
>>> props = img.regionprops()
"""
def phase_mean(sl, img):
return _calc_phase_mean(sl, img, self.periodicity)
def phase_std(sl, img):
return _calc_phase_std(sl, img, self.periodicity)
additional_props = {"phase_mean": phase_mean, "phase_std": phase_std}
# check arguments
if isinstance(properties, str):
properties = (properties,)
if extra_properties is None:
extra_properties = tuple()
# add extra_properties and move additional properties to extra_properties
properties = properties + tuple(ex.__name__ for ex in extra_properties)
for prop in properties:
if prop in additional_props.keys():
extra_properties = (additional_props[prop],) + extra_properties
if "p" in self.axes:
# this dimension will be label
raise ValueError("axis 'p' is forbidden in regionprops().")
prop_axes = complement_axes(self.labels.axes, self.axes)
shape = self.sizesof(prop_axes)
out = DataDict({p: PropArray(np.empty((self.labels.max(),) + shape, dtype=np.float32),
name=self.name+"-prop",
axes="p"+prop_axes,
source=self.source,
propname=p)
for p in properties})
# calculate property value for each slice
for sl, img in self.iter(prop_axes, exclude=self.labels.axes):
props = skmes.regionprops(self.labels, img, cache=False,
extra_properties=extra_properties)
label_sl = (slice(None),) + sl
for prop_name in properties:
# Both sides have length of p-axis (number of labels) so that values
# can be correctly substituted.
out[prop_name][label_sl] = [getattr(prop, prop_name) for prop in props]
for parr in out.values():
parr.set_scale(self)
return out