"""Analysis for patellar cartilage.
Attributes:
BOUNDS (dict): Upper bounds for quantitative values.
"""
import itertools
import os
import warnings
from copy import deepcopy
import numpy as np
import pandas as pd
import scipy.ndimage as sni
from dosma.defaults import preferences
from dosma.quant_vals import QuantitativeValueType
from dosma.tissues.tissue import Tissue, largest_cc
from dosma.utils import io_utils
import matplotlib.pyplot as plt
# milliseconds
BOUNDS = {
QuantitativeValueType.T2: 60.0,
QuantitativeValueType.T1_RHO: 100.0,
QuantitativeValueType.T2_STAR: 50.0,
}
__all__ = ["PatellarCartilage"]
[docs]class PatellarCartilage(Tissue):
"""Handles analysis and visualization for patellar cartilage."""
ID = 3
STR_ID = "pc"
FULL_NAME = "patellar cartilage"
# Expected quantitative values
T1_EXPECTED = 1000 # milliseconds
# Region Keys
_ANTERIOR_KEY = 0
_POSTERIOR_KEY = 1
_CORONAL_KEYS = [_ANTERIOR_KEY, _POSTERIOR_KEY]
_MEDIAL_KEY = 0
_LATERAL_KEY = 1
_SAGITTAL_KEYS = [_MEDIAL_KEY, _LATERAL_KEY]
_REGION_DEEP_KEY = 0
_REGION_SUPERFICIAL_KEY = 1
_TOTAL_AXIAL_KEY = -1
[docs] def __init__(self, weights_dir: str = None, medial_to_lateral: bool = None):
super().__init__(weights_dir=weights_dir, medial_to_lateral=medial_to_lateral)
self.regions_mask = None
[docs] def unroll_coronal(self, quant_map: np.ndarray):
"""Unroll patellar cartilage in the coronal direction.
Because patellar cartilage is flat, "unrolling" is projecting the patellar
cartilage onto the coronal axis.
Args:
quant_map (np.ndarray):
"""
mask = self.__mask__.volume
assert (
self.regions_mask is not None
), "region_mask not initialized. Should be initialized when mask is set"
region_mask_deep_superficial = self.regions_mask[..., 0]
superficial = (
(region_mask_deep_superficial == self._REGION_SUPERFICIAL_KEY) * mask * quant_map
)
superficial[superficial == 0] = np.nan
superficial = np.nanmean(superficial, axis=2)
deep = (region_mask_deep_superficial == self._REGION_DEEP_KEY) * mask * quant_map
deep[deep == 0] = np.nan
deep = np.nanmean(deep, axis=2)
total = mask * quant_map
total[total == 0] = np.nan
total = np.nanmean(total, axis=2)
return total, superficial, deep
[docs] def split_regions(self, base_map):
"""Split patellar cartilage into deep/superficial regions.
For patellar cartilage, the superficial/deep transition occurs in
the anterior/posterior (A/P) direction. The boundary is determined
for each non-zero 1D column spanning independently by the local
center-of-mass (COM). The medial/lateral (M/L) plane is computed
using the global COM.
Args:
base_map (ndarray): Binary 3D mask with orientation (SI, AP, ML/LM).
If `self.medial_to_lateral`, last dimension should be ML.
"""
if np.sum(base_map) == 0:
warnings.warn("No mask for `%s` was found." % self.FULL_NAME)
# Superficial/Deep (A/P)
locs = base_map.sum(axis=1).nonzero()
voxels = base_map[locs[0], :, locs[1]]
com_sup_inf = np.asarray(
[
int(np.ceil(sni.measurements.center_of_mass(voxels[i, :])[0]))
for i in range(voxels.shape[0])
]
)
region_mask_sup_deep = np.full(base_map.shape, self._REGION_DEEP_KEY)
for i in range(len(com_sup_inf)):
region_mask_sup_deep[
locs[0][i], : com_sup_inf[i], locs[1][i]
] = self._REGION_SUPERFICIAL_KEY
# M/L
cum_ml = np.nonzero(base_map.sum(axis=(0, 1)))[0] # noqa: F841
# midpoint_ml = int(np.ceil((np.min(cum_ml) + np.max(cum_ml)) / 2))
midpoint_ml = int(np.ceil(sni.measurements.center_of_mass(base_map)[2]))
region_mask_med_lat = np.full(base_map.shape, self._LATERAL_KEY)
medial_span = slice(0, midpoint_ml) if self.medial_to_lateral else slice(midpoint_ml, None)
region_mask_med_lat[:, :, medial_span] = self._MEDIAL_KEY
self.regions_mask = np.stack([region_mask_sup_deep, region_mask_med_lat], axis=-1)
def __calc_quant_vals__(self, quant_map, map_type):
subject_pid = self.pid
super().__calc_quant_vals__(quant_map, map_type)
assert (
self.regions_mask is not None
), "region_mask not initialized. Should be initialized when mask is set"
quant_map_volume = quant_map.volume
mask = self.__mask__.volume
quant_map_volume = mask * quant_map_volume
deep_superficial_map = self.regions_mask[..., 0]
med_lat_map = self.regions_mask[..., 1]
axial_names = ["superficial", "deep", "total"]
sagittal_names = ["medial", "lateral"]
pd_header = ["Subject", "Location", "Condyle", "Mean", "Std", "Median"]
pd_list = []
regions = itertools.product(
[self._REGION_SUPERFICIAL_KEY, self._REGION_DEEP_KEY, self._TOTAL_AXIAL_KEY],
[self._MEDIAL_KEY, self._LATERAL_KEY],
)
for axial, sagittal in regions:
if axial == self._TOTAL_AXIAL_KEY:
axial_map = np.asarray(
deep_superficial_map == self._REGION_SUPERFICIAL_KEY, dtype=np.float32
) + np.asarray(deep_superficial_map == self._REGION_DEEP_KEY, dtype=np.float32)
axial_map = np.asarray(axial_map, dtype=np.bool)
else:
axial_map = deep_superficial_map == axial
sagittal_map = med_lat_map == sagittal
curr_region_mask = quant_map_volume * axial_map * sagittal_map
curr_region_mask[curr_region_mask == 0] = np.nan
# discard all values that are 0
c_mean = np.nanmean(curr_region_mask)
c_std = np.nanstd(curr_region_mask)
c_median = np.nanmedian(curr_region_mask)
row_info = [
subject_pid,
axial_names[axial],
sagittal_names[sagittal],
c_mean,
c_std,
c_median,
]
pd_list.append(row_info)
# Generate 2D unrolled matrix
total, superficial, deep = self.unroll_coronal(quant_map.volume)
df = pd.DataFrame(pd_list, columns=pd_header)
qv_name = map_type.name
maps = [
{
"title": "%s superficial" % qv_name,
"data": superficial,
"xlabel": "Slice",
"ylabel": "Angle (binned)",
"filename": "%s_superficial" % qv_name,
"raw_data_filename": "%s_superficial.data" % qv_name,
},
{
"title": "%s deep" % qv_name,
"data": deep,
"xlabel": "Slice",
"ylabel": "Angle (binned)",
"filename": "%s_deep" % qv_name,
"raw_data_filename": "%s_deep.data" % qv_name,
},
{
"title": "%s total" % qv_name,
"data": total,
"xlabel": "Slice",
"ylabel": "Angle (binned)",
"filename": "%s_total" % qv_name,
"raw_data_filename": "%s_total.data" % qv_name,
},
]
self.__store_quant_vals__(maps, df, map_type)
[docs] def set_mask(self, mask):
msk = np.asarray(largest_cc(mask.volume), dtype=np.uint8)
mask_copy = deepcopy(mask)
mask_copy.volume = msk
super().set_mask(mask_copy)
self.split_regions(self.__mask__.volume)
def __save_quant_data__(self, dirpath):
"""Save quantitative data and 2D visualizations of patellar cartilage
Check which quantitative values (T2, T1rho, etc) are defined for
patellar cartilage and analyze these:
1. Save 2D total, superficial, and deep visualization maps
2. Save {'medial', 'lateral'}, {'anterior', 'posterior'},
{'superior', 'inferior', 'total'} data to excel file
:param dirpath: base filepath to save data
"""
q_names = []
dfs = []
for quant_val in QuantitativeValueType:
if quant_val.name not in self.quant_vals.keys():
continue
q_names.append(quant_val.name)
q_val = self.quant_vals[quant_val.name]
dfs.append(q_val[1])
q_name_dirpath = io_utils.mkdirs(os.path.join(dirpath, quant_val.name.lower()))
for q_map_data in q_val[0]:
filepath = os.path.join(q_name_dirpath, q_map_data["filename"])
xlabel = ""
ylabel = ""
title = q_map_data["title"]
data_map = q_map_data["data"]
axs_bounds = self.__get_axis_bounds__(data_map, leave_buffer=True)
plt.clf()
upper_bound = BOUNDS[quant_val]
if preferences.visualization_use_vmax:
# Hard bounds - clipping
plt.imshow(data_map, cmap="jet", vmin=0.0, vmax=BOUNDS[quant_val])
else:
# Try to use a soft bounds
if np.sum(data_map <= upper_bound) == 0:
plt.imshow(data_map, cmap="jet", vmin=0.0, vmax=BOUNDS[quant_val])
else:
warnings.warn(
"%s: Pixel value exceeded upper bound (%0.1f). Using normalized scale."
% (quant_val.name, upper_bound)
)
plt.imshow(data_map, cmap="jet")
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.title(title)
plt.ylim(axs_bounds[0])
plt.gca().invert_yaxis()
plt.xlim(axs_bounds[1])
# plt.axis('tight')
clb = plt.colorbar()
clb.ax.set_ylabel("(ms)")
plt.savefig(filepath)
# Save data
raw_data_filepath = os.path.join(
q_name_dirpath, "raw_data", q_map_data["raw_data_filename"]
)
io_utils.save_pik(raw_data_filepath, data_map)
if len(dfs) > 0:
io_utils.save_tables(os.path.join(dirpath, "data.xlsx"), dfs, q_names)