#!/usr/bin/env python
# -*- coding: utf-8
# Functions processing segmentation data
import math
import platform
import numpy as np
from skimage import measure, transform
import logging
import nibabel
from spinalcordtoolbox.image import Image
from spinalcordtoolbox.aggregate_slicewise import Metric
from spinalcordtoolbox.centerline.core import ParamCenterline, get_centerline
from spinalcordtoolbox.resampling import resample_nib
from spinalcordtoolbox.utils import sct_progress_bar
[docs]def compute_shape(segmentation, angle_correction=True, param_centerline=None, verbose=1):
"""
Compute morphometric measures of the spinal cord in the transverse (axial) plane from the segmentation.
The segmentation could be binary or weighted for partial volume [0,1].
:param segmentation: input segmentation. Could be either an Image or a file name.
:param angle_correction:
:param param_centerline: see centerline.core.ParamCenterline()
:param verbose:
:return metrics: Dict of class Metric(). If a metric cannot be calculated, its value will be nan.
:return fit_results: class centerline.core.FitResults()
"""
# List of properties to output (in the right order)
property_list = ['area',
'angle_AP',
'angle_RL',
'diameter_AP',
'diameter_RL',
'eccentricity',
'orientation',
'solidity',
'length'
]
im_seg = Image(segmentation).change_orientation('RPI')
# Getting image dimensions. x, y and z respectively correspond to RL, PA and IS.
nx, ny, nz, nt, px, py, pz, pt = im_seg.dim
pr = min([px, py])
# Resample to isotropic resolution in the axial plane. Use the minimum pixel dimension as target dimension.
im_segr = resample_nib(im_seg, new_size=[pr, pr, pz], new_size_type='mm', interpolation='linear')
# Update dimensions from resampled image.
nx, ny, nz, nt, px, py, pz, pt = im_segr.dim
# Extract min and max index in Z direction
data_seg = im_segr.data
X, Y, Z = (data_seg > 0).nonzero()
min_z_index, max_z_index = min(Z), max(Z)
# Initialize dictionary of property_list, with 1d array of nan (default value if no property for a given slice).
shape_properties = {key: np.full_like(np.empty(nz), np.nan, dtype=np.double) for key in property_list}
fit_results = None
if angle_correction:
# compute the spinal cord centerline based on the spinal cord segmentation
# here, param_centerline.minmax needs to be False because we need to retrieve the total number of input slices
_, arr_ctl, arr_ctl_der, fit_results = get_centerline(im_segr, param=param_centerline, verbose=verbose)
# Loop across z and compute shape analysis
for iz in sct_progress_bar(range(min_z_index, max_z_index + 1), unit='iter', unit_scale=False, desc="Compute shape analysis",
ascii=True, ncols=80):
# Extract 2D patch
current_patch = im_segr.data[:, :, iz]
if angle_correction:
# Extract tangent vector to the centerline (i.e. its derivative)
tangent_vect = np.array([arr_ctl_der[0][iz - min_z_index] * px,
arr_ctl_der[1][iz - min_z_index] * py,
pz])
# Normalize vector by its L2 norm
tangent_vect = tangent_vect / np.linalg.norm(tangent_vect)
# Compute the angle about AP axis between the centerline and the normal vector to the slice
v0 = [tangent_vect[0], tangent_vect[2]]
v1 = [0, 1]
angle_AP_rad = np.math.atan2(np.linalg.det([v0, v1]), np.dot(v0, v1))
# Compute the angle about RL axis between the centerline and the normal vector to the slice
v0 = [tangent_vect[1], tangent_vect[2]]
v1 = [0, 1]
angle_RL_rad = np.math.atan2(np.linalg.det([v0, v1]), np.dot(v0, v1))
# Apply affine transformation to account for the angle between the centerline and the normal to the patch
tform = transform.AffineTransform(scale=(np.cos(angle_RL_rad), np.cos(angle_AP_rad)))
# Convert to float64, to avoid problems in image indexation causing issues when applying transform.warp
current_patch = current_patch.astype(np.float64)
# TODO: make sure pattern does not go extend outside of image border
current_patch_scaled = transform.warp(current_patch,
tform.inverse,
output_shape=current_patch.shape,
order=1,
)
else:
current_patch_scaled = current_patch
angle_AP_rad, angle_RL_rad = 0.0, 0.0
# compute shape properties on 2D patch
shape_property = _properties2d(current_patch_scaled, [px, py])
if shape_property is not None:
# Add custom fields
shape_property['angle_AP'] = angle_AP_rad * 180.0 / math.pi
shape_property['angle_RL'] = angle_RL_rad * 180.0 / math.pi
shape_property['length'] = pz / (np.cos(angle_AP_rad) * np.cos(angle_RL_rad))
# Loop across properties and assign values for function output
for property_name in property_list:
shape_properties[property_name][iz] = shape_property[property_name]
else:
logging.warning('\nNo properties for slice: {}'.format(iz))
""" DEBUG
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from matplotlib.figure import Figure
fig = Figure()
FigureCanvas(fig)
ax = fig.add_subplot(111)
ax.imshow(current_patch_scaled)
ax.grid()
ax.set_xlabel('y')
ax.set_ylabel('x')
fig.savefig('tmp_fig.png')
"""
metrics = {}
for key, value in shape_properties.items():
# Making sure all entries added to metrics have results
if not value == []:
metrics[key] = Metric(data=np.array(value), label=key)
return metrics, fit_results
def _properties2d(image, dim):
"""
Compute shape property of the input 2D image. Accounts for partial volume information.
:param image: 2D input image in uint8 or float (weighted for partial volume) that has a single object.
:param dim: [px, py]: Physical dimension of the image (in mm). X,Y respectively correspond to AP,RL.
:return:
"""
upscale = 5 # upscale factor for resampling the input image (for better precision)
pad = 3 # padding used for cropping
# Check if slice is empty
if not image.any():
logging.debug('The slice is empty.')
return None
# Normalize between 0 and 1 (also check if slice is empty)
image_norm = (image - image.min()) / (image.max() - image.min())
# Convert to float64
image_norm = image_norm.astype(np.float64)
# Binarize image using threshold at 0. Necessary input for measure.regionprops
image_bin = np.array(image_norm > 0.5, dtype='uint8')
# Get all closed binary regions from the image (normally there is only one)
regions = measure.regionprops(image_bin, intensity_image=image_norm)
# Check number of regions
if len(regions) > 1:
logging.debug('There is more than one object on this slice.')
return None
region = regions[0]
# Get bounding box of the object
minx, miny, maxx, maxy = region.bbox
# Use those bounding box coordinates to crop the image (for faster processing)
image_crop = image_norm[np.clip(minx-pad, 0, image_bin.shape[0]): np.clip(maxx+pad, 0, image_bin.shape[0]),
np.clip(miny-pad, 0, image_bin.shape[1]): np.clip(maxy+pad, 0, image_bin.shape[1])]
# Oversample image to reach sufficient precision when computing shape metrics on the binary mask
image_crop_r = transform.pyramid_expand(image_crop, upscale=upscale, sigma=None, order=1)
# Binarize image using threshold at 0. Necessary input for measure.regionprops
image_crop_r_bin = np.array(image_crop_r > 0.5, dtype='uint8')
# Get all closed binary regions from the image (normally there is only one)
regions = measure.regionprops(image_crop_r_bin, intensity_image=image_crop_r)
region = regions[0]
# Compute area with weighted segmentation and adjust area with physical pixel size
area = np.sum(image_crop_r) * dim[0] * dim[1] / upscale ** 2
# Compute ellipse orientation, modulo pi, in deg, and between [0, 90]
orientation = fix_orientation(region.orientation)
# Find RL and AP diameter based on major/minor axes and cord orientation=
[diameter_AP, diameter_RL] = \
_find_AP_and_RL_diameter(region.major_axis_length, region.minor_axis_length, orientation,
[i / upscale for i in dim])
# TODO: compute major_axis_length/minor_axis_length by summing weighted voxels along axis
# Deal with https://github.com/neuropoly/spinalcordtoolbox/issues/2307
if any(x in platform.platform() for x in ['Darwin-15', 'Darwin-16']):
solidity = np.nan
else:
solidity = region.solidity
# Fill up dictionary
properties = {'area': area,
'diameter_AP': diameter_AP,
'diameter_RL': diameter_RL,
'centroid': region.centroid,
'eccentricity': region.eccentricity,
'orientation': orientation,
'solidity': solidity # convexity measure
}
return properties
[docs]def fix_orientation(orientation):
"""Re-map orientation from skimage.regionprops from [-pi/2,pi/2] to [0,90] and rotate by 90deg because image axis
are inverted"""
orientation_new = orientation * 180.0 / math.pi
if 360 <= abs(orientation_new) <= 540:
orientation_new = 540 - abs(orientation_new)
if 180 <= abs(orientation_new) <= 360:
orientation_new = 360 - abs(orientation_new)
if 90 <= abs(orientation_new) <= 180:
orientation_new = 180 - abs(orientation_new)
return abs(orientation_new)
def _find_AP_and_RL_diameter(major_axis, minor_axis, orientation, dim):
"""
This script checks the orientation of the and assigns the major/minor axis to the appropriate dimension, right-
left (RL) or antero-posterior (AP). It also multiplies by the pixel size in mm.
:param major_axis: major ellipse axis length calculated by regionprops
:param minor_axis: minor ellipse axis length calculated by regionprops
:param orientation: orientation in degree. Ranges between [0, 90]
:param dim: pixel size in mm.
:return: diameter_AP, diameter_RL
"""
if 0 <= orientation < 45.0:
diameter_AP = minor_axis
diameter_RL = major_axis
else:
diameter_AP = major_axis
diameter_RL = minor_axis
# Adjust with pixel size
diameter_AP *= dim[0]
diameter_RL *= dim[1]
return diameter_AP, diameter_RL