import math
import numpy as np
from scipy import ndimage as ndi
from scipy.spatial import distance
from skimage import draw, segmentation, measure, morphology
import fibsem.conversions
def _d_points_to_line(points, point_on_line, direction):
points = np.asanyarray(points)
point_on_line = np.asanyarray(point_on_line)
direction = np.asanyarray(direction)
points_centered = points - point_on_line
unit_direction = direction / np.linalg.norm(direction)
line_vector = points_centered @ unit_direction
d_vector = points_centered - line_vector[:, np.newaxis] * unit_direction
return line_vector, np.linalg.norm(d_vector, axis=1)
def _edge_coordinates(mask_image, direction_vector):
"""Finds the outermost edges of a mask, given a direction from the center.
Parameters
----------
mask_image : array of bool
The image containing the object. Pixels with value True are considered
part of the object.
direction_vector : 2 element numpy array
Orientation vector describing a direction.
Returns
-------
src, dst : tuple of int
Source and destination points along the edge of the mask.
"""
boundary = segmentation.find_boundaries(mask_image,
connectivity=mask_image.ndim,
mode='inner')
boundary_points = np.argwhere(boundary)
centroid_coords = centroid(mask_image)
along, distances = _d_points_to_line(boundary_points,
centroid_coords,
direction_vector)
pts_by_dist = np.argsort(distances)
pt1 = boundary_points[pts_by_dist[0]]
distances_to_pt1 = distance.cdist(boundary_points, [pt1])[:, 0]
std_distance = np.std(along)
for pt in pts_by_dist[1:]:
if distance.euclidean(boundary_points[pt], pt1) > 2 * std_distance:
return pt1, boundary_points[pt]
def _orientation_vector(image):
"""Orientation vector of the major axis."""
T = measure.inertia_tensor(image)
v, V = np.linalg.eig(T)
smallest = np.argmin(np.abs(v))
return np.flip(V[:, smallest], axis=0)
def _orthogonal_orientation_vector(image):
"""Orientation vector of minor axis.
Orthogonal to the major axis orientation vector.
"""
u, v = _orientation_vector(image)
return np.array([v, -u])
def _calculate_roi_height(coord_realspace, sample):
"""[summary]
Parameters
----------
coord_realspace : [type]
Coordinates in row, column format.
sample : [type]
[description]
Returns
-------
[type]
[description]
"""
mask = fibsem.conversions._convert_to_numpy(sample.mask)
orthog_theta = -rotation_angle(sample.minor_axis)
max_line_length = np.ceil(fibsem.conversions.meters_to_pixels(
line_length(sample.minor_axis),
sample.pixel_size))
coord_pixels = fibsem.conversions._realspace_coord_to_pixels(
coord_realspace, sample.pixel_size, mask.shape)
temp_coord = _point_given_angle_and_distance(
coord_pixels, orthog_theta, max_line_length)
rr, cc = draw.line(int(coord_pixels[0]),
int(coord_pixels[1]),
int(temp_coord[0]),
int(temp_coord[1]))
rr = np.clip(rr, 0, sample.mask.shape[0] - 1)
cc = np.clip(cc, 0, sample.mask.shape[1] - 1)
temp_img = np.zeros(sample.mask.shape)
temp_img[rr, cc] = 1
roi_height = fibsem.conversions.pixels_to_meters(
np.sum(mask * temp_img) * 2, sample.pixel_size)
return roi_height
[docs]def centroid(mask_image):
"""Centroid coordinate of a masked object region.
Parameters
----------
mask_image : array of bool
The image containing the object. Pixels with value True are considered
part of the object.
Returns
-------
centroid
Centroid pixel coordinates in row, column format. Origin at top left.
"""
return np.mean(np.argwhere(mask_image), axis=0)
[docs]def major_minor_axes_pixels(mask_image):
"""Major and minor axes for object in mask image, in pixel coordinates.
Each axis contains two coordinates, the start and end of the axis line
at the outer edges of the object border.
Coordinates are in row, column format.
Parameters
----------
mask_image : array of bool
The image containing the object. Pixels with value True are considered
part of the object.
Returns
-------
major_axis, minor_axis : tuple of numpy arrays
Start and end coordinates where axis intersects with the object edge.
Coordinates are in row, column format, units are pixels.
"""
mask_image = fibsem.conversions._convert_to_numpy(mask_image)
major_direction_vector = _orientation_vector(mask_image)
minor_direction_vector = _orthogonal_orientation_vector(mask_image)
major_axis = _edge_coordinates(mask_image, major_direction_vector)
minor_axis = _edge_coordinates(mask_image, minor_direction_vector)
return major_axis, minor_axis
[docs]def major_minor_axes_realspace(mask_image, pixelsize, reference_coord=None):
"""Major and minor axes for object in mask image, real space coordinates.
Parameters
----------
mask_image : array of bool
The image containing the object. Pixels with value True are considered
part of the object.
pixelsize : Pixel size object, with attributes x and y. Units of meters.
Returns
-------
major_axis, minor_axis :
Start and end coordinates where axis intersects with the object edge.
Units are in meters, zero coord is at the center of the field of view.
Coordinates are in row, column format.
"""
height, width = fibsem.conversions._convert_to_numpy(mask_image).shape
major_axis_px, minor_axis_px = major_minor_axes_pixels(mask_image)
# major axis
major_axis_x0 = (+major_axis_px[0][1] - (width / 2)) * pixelsize.x
major_axis_y0 = (-major_axis_px[0][0] + (height / 2)) * pixelsize.y
major_axis_x1 = (+major_axis_px[1][1] - (width / 2)) * pixelsize.x
major_axis_y1 = (-major_axis_px[1][0] + (height / 2)) * pixelsize.y
major_axis = (np.array([major_axis_y0, major_axis_x0]),
np.array([major_axis_y1, major_axis_x1]))
# minor axis
minor_axis_x0 = (+minor_axis_px[0][1] - (width / 2)) * pixelsize.x
minor_axis_y0 = (-minor_axis_px[0][0] + (height / 2)) * pixelsize.y
minor_axis_x1 = (+minor_axis_px[1][1] - (width / 2)) * pixelsize.x
minor_axis_y1 = (-minor_axis_px[1][0] + (height / 2)) * pixelsize.y
minor_axis = (np.array([minor_axis_y0, minor_axis_x0]),
np.array([minor_axis_y1, minor_axis_x1]))
# Align with sample head, if position is specified
if reference_coord:
major_axis, minor_axis = _align_with_head(
major_axis, minor_axis, reference_coord)
return major_axis, minor_axis
[docs]def multi_axis(label_image):
"""Like :func:`major_minor_axes`, but does so for each label in an image.
Parameters
----------
label_image : array of int
The image containing one integer label per object.
Returns
-------
lines : dict of tuple of tuple of int
A dictionary mapping labels (int) to start and end points for the
longest axis line. Coordinates are in row, column format.
"""
results = {}
for label in np.unique(label_image)[1:]: # ignore 0 as background
mask_image = label_image == label
results[label] = major_minor_axes_pixels(mask_image)
return results
[docs]def line_length(axis_coords):
"""Length in pixels of line connecting two points.
Parameters
----------
pt1, pt2 : Source and destination points along the edge of the object.
Coordinates are in row, column format.
Returns
-------
np.hypot(p1, p2) : line length in units of pixels
"""
pt2, pt1 = np.array(axis_coords)
rise = pt2[0] - pt1[0]
run = pt2[1] - pt1[1]
hypotenuse = np.sqrt(rise**2 + run**2)
return hypotenuse
[docs]def rotation_angle(axis_coords):
"""Find the rotation angle for Autoscript.
Autoscript requires: Positive direction is clockwise, value in radians.
Parameters
----------
axis_coords : Start and end coordinates of axis.
Coordinates are in row, column format.
Returns
-------
theta : in radians, positive direction is clockwise.
"""
dy, dx = axis_coords[1] - axis_coords[0]
theta = -math.atan2(dy, dx) # negative because of autoscript convention
return theta
def _align_with_head(major_axis, minor_axis, reference_coord):
length_1 = line_length((major_axis[0], reference_coord))
length_2 = line_length((major_axis[1], reference_coord))
if length_1 > length_2:
major_axis = np.flip(major_axis, axis=0)
minor_axis = np.flip(minor_axis, axis=0)
return major_axis, minor_axis
def _point_given_angle_and_distance(reference_coord, theta, distance):
theta = -theta # Autoscript is positive clockwise, everything else is not
dx = distance * np.cos(theta)
dy = distance * np.sin(theta)
new_x = reference_coord[1] + dx
new_y = reference_coord[0] + dy
new_coordinate = np.array([new_y, new_x])
return new_coordinate