def refine_ellipsoid(image3d, params, spacing=1, rad_range=None, maxfit_size=2,
spline_order=3, threshold=0.1):
""" Refines coordinates of a 3D ellipsoid, starting from given parameters.
Interpolates the image along lines perpendicular to the ellipsoid.
The maximum along each line is found using linear regression of the
descrete derivative.
Parameters
----------
image3d : 3d numpy array of numbers
Image indices are interpreted as (z, y, x)
params : tuple
zr, yr, xr, zc, yc, xc
spacing: number
spacing along radial direction
rad_range: tuple of floats
length of the line (distance inwards, distance outwards)
maxfit_size: integer
pixels around maximum pixel that will be used in linear regression
spline_order: integer
interpolation order for edge crossections
threshold: float
a threshold is calculated based on the global maximum
fitregions are rejected if their average value is lower than this
Returns
-------
- zr, yr, xr, zc, yc, xc, skew_y, skew_x
- contour coordinates at z = 0
"""
if not np.all([x > 0 for x in params]):
raise ValueError("All zc, yc, xc, zr, yr, xr params should be positive")
assert image3d.ndim == 3
zr, yr, xr, zc, yc, xc = params
if rad_range is None:
rad_range = (-min(zr, yr, xr) / 2, min(zr, yr, xr) / 2)
steps = np.arange(rad_range[0], rad_range[1] + 1, 1)
pos, normal = ellipsoid_grid((zr, yr, xr), (zc, yc, xc), spacing=spacing)
coords = normal[:, :, np.newaxis] * steps[np.newaxis, np.newaxis, :] + \
pos[:, :, np.newaxis]
# interpolate the image on calculated coordinates
intensity = map_coordinates(image3d, coords, order=spline_order)
# identify the regions around the max value
r_dev = max_linregress(intensity, maxfit_size, threshold)
# calculate new coords
coord_new = pos + (r_dev + rad_range[0])*normal
coord_new = coord_new[:, np.isfinite(coord_new).all(0)]
# fit ellipsoid
radius, center, skew = fit_ellipsoid(coord_new, mode='xy',
return_mode='skew')
return tuple(radius) + tuple(center) + tuple(skew), coord_new.T
评论列表
文章目录
