diff --git a/platipy/imaging/dose/dvh.py b/platipy/imaging/dose/dvh.py index 8b397879..8f1d080b 100644 --- a/platipy/imaging/dose/dvh.py +++ b/platipy/imaging/dose/dvh.py @@ -119,7 +119,8 @@ def calculate_d_x(dvh, x, label=None): Args: dvh (pandas.DataFrame): DVH DataFrame as produced by calculate_dvh_for_labels - x (float): The dose which x percent of the volume receives + x (float|list): The dose threshold (or list of dose thresholds) which x percent of the + volume receives label (str, optional): The label to compute the metric for. Computes for all metrics if not set. Defaults to None. @@ -130,15 +131,31 @@ def calculate_d_x(dvh, x, label=None): if label: dvh = dvh[dvh.label == label] + if not isinstance(x, list): + x = [x] + bins = np.array([b for b in dvh.columns if isinstance(b, float)]) values = np.array(dvh[bins]) - i, j = np.where(values >= x / 100) - metrics = [] for idx in range(len(dvh)): d = dvh.iloc[idx] - metrics.append({"label": d.label, "metric": f"D{x}", "value": bins[j][i == idx][-1]}) + + m = {"label": d.label} + + for threshold in x: + value = np.interp(threshold / 100, values[idx][::-1], bins[::-1]) + if values[idx, 0] == np.sum(values[idx]): + value = 0 + + # Interp will return zero when computing D100, do compute this separately + if threshold == 100: + i, j = np.where(values == 1.0) + value = bins[j][i == idx][-1] + + m[f"D{threshold}"] = value + + metrics.append(m) return pd.DataFrame(metrics) @@ -148,7 +165,7 @@ def calculate_v_x(dvh, x, label=None): Args: dvh (pandas.DataFrame): DVH DataFrame as produced by calculate_dvh_for_labels - x (float): The dose to get the volume for. + x (float|list): The dose threshold (or list of dose thresholds) to get the volume for. label (str, optional): The label to compute the metric for. Computes for all metrics if not set. Defaults to None. @@ -159,21 +176,64 @@ def calculate_v_x(dvh, x, label=None): if label: dvh = dvh[dvh.label == label] + if not isinstance(x, list): + x = [x] + bins = np.array([b for b in dvh.columns if isinstance(b, float)]) values = np.array(dvh[bins]) - i = np.where(bins == x) metrics = [] for idx in range(len(dvh)): d = dvh.iloc[idx] - value_idx = values[idx, i] - value = 0.0 - if value_idx.shape[1] > 0: - value = d.cc * values[idx, i][0, 0] - - metric_name = f"V{x}" - if x - int(x) == 0: - metric_name = f"V{int(x)}" - metrics.append({"label": d.label, "metric": metric_name, "value": value}) + + m = {"label": d.label} + + for threshold in x: + value = np.interp(threshold, bins, values[idx]) * d.cc + + metric_name = f"V{threshold}" + if threshold - int(threshold) == 0: + metric_name = f"V{int(threshold)}" + + m[metric_name] = value + + metrics.append(m) + + return pd.DataFrame(metrics) + + +def calculate_d_cc_x(dvh, x, label=None): + """Compute the dose which is received by cc of the volume + + Args: + dvh (pandas.DataFrame): DVH DataFrame as produced by calculate_dvh_for_labels + x (float|list): The cc (or list of cc's) to compute the dose at. + label (str, optional): The label to compute the metric for. Computes for all metrics if not + set. Defaults to None. + + Returns: + pandas.DataFrame: Data frame with a row for each label containing the metric and value. + """ + + if label: + dvh = dvh[dvh.label == label] + + if not isinstance(x, list): + x = [x] + + metrics = [] + for idx in range(len(dvh)): + + d = dvh.iloc[idx] + m = {"label": d.label} + + for threshold in x: + cc_at = (threshold / dvh[dvh.label == d.label].cc.iloc[0]) * 100 + cc_at = min(cc_at, 100) + cc_val = calculate_d_x(dvh[dvh.label == d.label], cc_at)[f"D{cc_at}"].iloc[0] + + m[f"D{threshold}cc"] = cc_val + + metrics.append(m) return pd.DataFrame(metrics) diff --git a/platipy/imaging/dose/metric.py b/platipy/imaging/dose/metric.py new file mode 100644 index 00000000..850600b1 --- /dev/null +++ b/platipy/imaging/dose/metric.py @@ -0,0 +1,183 @@ +# Copyright 2020 University of New South Wales, University of Sydney, Ingham Institute + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import numpy as np +import SimpleITK as sitk +import pandas as pd + + +def calculate_d_mean(dose_grid, label): + """Calculate the mean dose of a structure + + Args: + dose_grid (SimpleITK.Image): The dose grid. + label (SimpleITK.Image): The (binary) label defining a structure. + + Returns: + float: The mean dose in Gy. + """ + + dose_grid = sitk.Resample(dose_grid, label, sitk.Transform(), sitk.sitkLinear) + dose_array = sitk.GetArrayFromImage(dose_grid) + mask_array = sitk.GetArrayFromImage(label) + + return dose_array[mask_array > 0].mean() + + +def calculate_d_max(dose_grid, label): + """Calculate the maximum dose of a structure + + Args: + dose_grid (SimpleITK.Image): The dose grid. + label (SimpleITK.Image): The (binary) label defining a structure. + + Returns: + float: The maximum dose in Gy. + """ + + dose_grid = sitk.Resample(dose_grid, label, sitk.Transform(), sitk.sitkLinear) + dose_array = sitk.GetArrayFromImage(dose_grid) + mask_array = sitk.GetArrayFromImage(label) + + return dose_array[mask_array > 0].max() + + +def calculate_d_to_volume(dose_grid, label, volume, volume_in_cc=False): + """Calculate the dose to a (relative) volume of the label + + Args: + dose_grid (SimpleITK.Image): The dose grid. + label (SimpleITK.Image): The (binary) label defining a structure. + volume (float): The relative volume in %. + volume_in_cc (bool, optional): Whether the volume is in cc (versus percent). + Defaults to False. + + Returns: + float: The dose to volume ratio. + """ + + dose_grid = sitk.Resample(dose_grid, label, sitk.Transform(), sitk.sitkLinear) + dose_array = sitk.GetArrayFromImage(dose_grid) + mask_array = sitk.GetArrayFromImage(label) + + if volume_in_cc: + volume = (volume * 1000 / ((mask_array > 0).sum() * np.product(label.GetSpacing()))) * 100 + + if volume > 100: + volume = 100 + + return np.percentile(dose_array[mask_array > 0], 100 - volume) + + +def calculate_v_receiving_dose(dose_grid, label, dose_threshold, relative=True): + """Calculate the (relative) volume receiving a dose above a threshold + + Args: + dose_grid (SimpleITK.Image): The dose grid. + label (SimpleITK.Image): The (binary) label defining a structure. + dose_threshold (float): The dose threshold in Gy. + relative (bool, optional): If true results will be returned as relative volume, otherwise + as volume in cc. Defaults to True. + + Returns: + float: The (relative) volume receiving a dose above the threshold, as a percent. + """ + + dose_grid = sitk.Resample(dose_grid, label, sitk.Transform(), sitk.sitkLinear) + dose_array = sitk.GetArrayFromImage(dose_grid) + mask_array = sitk.GetArrayFromImage(label) + + dose_array_masked = dose_array[mask_array > 0] + + num_voxels = (mask_array > 0).sum() + + relative_volume = (dose_array_masked >= dose_threshold).sum() / num_voxels * 100 + if relative: + return relative_volume + + total_volume = (mask_array > 0).sum() * np.product(label.GetSpacing()) / 1000 + + return relative_volume * total_volume + + +def calculate_d_to_volume_for_labels(dose_grid, labels, volume, volume_in_cc=False): + """Calculate the dose which x percent of the volume receives for a set of labels + + Args: + dose_grid (SimpleITK.Image): The dose grid. + labels (dict): A Python dictionary containing the label name as key and the SimpleITK.Image + binary mask as value. + volume (float|list): The relative volume (or list of volumes) in %. + volume_in_cc (bool, optional): Whether the volume is in cc (versus percent). + Defaults to False. + + Returns: + pandas.DataFrame: Data frame with a row for each label containing the metric and value. + """ + + if not isinstance(volume, list): + volume = [volume] + + metrics = [] + for label in labels: + + m = {"label": label} + + for v in volume: + col_name = f"D{v}" + if volume_in_cc: + col_name = f"D{v}cc" + + m[col_name] = calculate_d_to_volume( + dose_grid, labels[label], v, volume_in_cc=volume_in_cc + ) + + metrics.append(m) + + return pd.DataFrame(metrics) + + +def calculate_v_receiving_dose_for_labels(dose_grid, labels, dose_threshold, relative=True): + """Get the volume (in cc) which receives x dose for a set of labels + + Args: + dose_grid (SimpleITK.Image): The dose grid. + labels (SimpleITK.Image): The (binary) label defining a structure. + dose_threshold (float|list): The dose threshold (or list of thresholds) in Gy. + relative (bool, optional): If true results will be returned as relative volume, otherwise + as volume in cc. Defaults to True. + + Returns: + pandas.DataFrame: Data frame with a row for each label containing the metric and value. + """ + + if not isinstance(dose_threshold, list): + dose_threshold = [dose_threshold] + + metrics = [] + for label in labels: + + m = {"label": label} + + for dt in dose_threshold: + + metric_name = f"V{dt}" + if dt - int(dt) == 0: + metric_name = f"V{int(dt)}" + + m[metric_name] = calculate_v_receiving_dose(dose_grid, labels[label], dt, relative) + + metrics.append(m) + + return pd.DataFrame(metrics) diff --git a/platipy/imaging/visualisation/visualiser.py b/platipy/imaging/visualisation/visualiser.py index 48ec44ef..2db2ec5f 100644 --- a/platipy/imaging/visualisation/visualiser.py +++ b/platipy/imaging/visualisation/visualiser.py @@ -109,12 +109,13 @@ def clear(self): self.__comparison_overlays = [] self.__vector_overlays = [] - def set_limits_from_label(self, label, expansion=[0, 0, 0]): + def set_limits_from_label(self, label, expansion=2): """Sets the limits of the axes to the bounds of the given label. Args: label (sitk.Image): The label around which to set the limits - expansion (list, optional): Expansion (in mm) around the label. Defaults to [0, 0, 0]. + expansion (list | float, optional): Expansion (in mm) around the label. + Defaults to 2. """ (sag_size, cor_size, ax_size), (sag_0, cor_0, ax_0) = label_to_roi( @@ -496,6 +497,14 @@ def _display_slice(self): sp_plane, _, sp_slice = image.GetSpacing() asp = (1.0 * sp_slice) / sp_plane + # Get the size - this is used for extent + size_sag, size_cor, size_ax = self.__image.GetSize() + extent_dict = { + "x": (0, size_cor, 0, size_ax), + "y": (0, size_sag, 0, size_ax), + "z": (0, size_sag, 0, size_cor), + } + if self.__projection is True: projection = "max" else: @@ -569,6 +578,7 @@ def _display_slice(self): origin={"normal": "upper", "reversed": "lower"}[self.__origin], cmap=self.__colormap, clim=(window[0], window[0] + window[1]), + extent=extent_dict["z"], ) cor_view = ax_cor.imshow( cor_img, @@ -577,6 +587,7 @@ def _display_slice(self): interpolation="none", cmap=self.__colormap, clim=(window[0], window[0] + window[1]), + extent=extent_dict["y"], ) sag_view = ax_sag.imshow( sag_img, @@ -585,6 +596,7 @@ def _display_slice(self): interpolation="none", cmap=self.__colormap, clim=(window[0], window[0] + window[1]), + extent=extent_dict["x"], ) ax_ax.axis("off") @@ -665,6 +677,7 @@ def _display_slice(self): origin=org, cmap=self.__colormap, clim=(window[0], window[0] + window[1]), + extent=extent_dict[self.__axis], ) ax.axis("off") @@ -713,6 +726,14 @@ def _overlay_comparison(self): upper = np.percentile(nda_original, 99) window = (lower, upper - lower) + # Get the size - this is used for extent + size_sag, size_cor, size_ax = self.__image.GetSize() + extent_dict = { + "x": (0, size_cor, 0, size_ax), + "y": (0, size_sag, 0, size_ax), + "z": (0, size_sag, 0, size_cor), + } + if self.__axis == "ortho": figure_size = ( self.__figure_size, @@ -753,6 +774,7 @@ def _overlay_comparison(self): aspect=1.0, origin={"normal": "upper", "reversed": "lower"}[self.__origin], interpolation="none", + extent=extent_dict["z"], ) nda_colormix = generate_comparison_colormix( @@ -767,6 +789,7 @@ def _overlay_comparison(self): origin="lower", aspect=asp, interpolation="none", + extent=extent_dict["y"], ) nda_colormix = generate_comparison_colormix( @@ -781,6 +804,7 @@ def _overlay_comparison(self): origin="lower", aspect=asp, interpolation="none", + extent=extent_dict["x"], ) ax_ax.axis("off") @@ -842,6 +866,7 @@ def _overlay_comparison(self): aspect=asp, interpolation="none", origin=org, + extent=extent_dict[self.__axis], ) ax.axis("off") @@ -889,11 +914,12 @@ def _adjust_view(self): 1 / asp * (cor_orig_1 - cor_orig_0) + (ax_orig_1 - ax_orig_0) ) - if origin == "reversed": - cor_0, cor_1 = cor_1, cor_0 - ax_ax.set_xlim(sag_0, sag_1) - ax_ax.set_ylim(cor_1, cor_0) + + if origin == "reversed": + ax_ax.set_ylim(cor_0, cor_1) + else: + ax_ax.set_ylim(cor_orig_1 - cor_1, cor_orig_1 - cor_0) ax_cor.set_xlim(sag_0, sag_1) ax_cor.set_ylim(ax_0, ax_1) @@ -950,7 +976,7 @@ def _adjust_view(self): x_0, x_1 = sorted([x_0, x_1]) y_0, y_1 = sorted([y_0, y_1]) - if self.__axis == "z": + if self.__axis == "z" and self.__origin == "normal": y_0, y_1 = y_1, y_0 ratio_x = np.abs(x_1 - x_0) / np.abs(x_orig_1 - x_orig_0) @@ -995,6 +1021,14 @@ def _overlay_contours(self): # Test types of axes axes = self.__figure.axes[:4] + # Get the size - this is used for extent + size_sag, size_cor, size_ax = self.__image.GetSize() + extent_dict = { + "x": (0, size_cor, 0, size_ax), + "y": (0, size_sag, 0, size_ax), + "z": (0, size_sag, 0, size_cor), + } + if self.__axis in ["x", "y", "z"]: ax = axes[0] s = return_slice(self.__axis, self.__cut) @@ -1017,6 +1051,11 @@ def _overlay_contours(self): ) contour_disp = sitk.GetArrayFromImage(contour_disp_proj) + if self.__axis == "z": + origin = {"normal": "upper", "reversed": "lower"}[self.__origin] + else: + origin = "lower" + try: ax.contour( contour_disp, @@ -1025,7 +1064,8 @@ def _overlay_contours(self): # alpha=0.8, linewidths=lw_dict[c_name], linestyles=ls_dict[c_name], - origin="lower", + extent=extent_dict[self.__axis], + origin=origin, ) ax.plot( [0], @@ -1091,11 +1131,12 @@ def _overlay_contours(self): ax_ax.contour( contour_ax, - levels=[0.5], + levels=[0], linewidths=lw_dict[c_name], linestyles=ls_dict[c_name], colors=[color_dict[c_name]], - origin="lower", + extent=extent_dict["z"], + origin={"normal": "upper", "reversed": "lower"}[self.__origin], ) ax_ax.plot( [0], @@ -1112,6 +1153,7 @@ def _overlay_contours(self): linewidths=lw_dict[c_name], linestyles=ls_dict[c_name], colors=[color_dict[c_name]], + extent=extent_dict["y"], origin="lower", ) ax_sag.contour( @@ -1120,6 +1162,7 @@ def _overlay_contours(self): linewidths=lw_dict[c_name], linestyles=ls_dict[c_name], colors=[color_dict[c_name]], + extent=extent_dict["x"], origin="lower", ) @@ -1158,10 +1201,20 @@ def _overlay_scalar_field(self): else: norm = None + # Get the spacing sp_plane, _, sp_slice = scalar_image.GetSpacing() + # Get the aspect ratio asp = (1.0 * sp_slice) / sp_plane + # Get the size - this is used for extent + size_sag, size_cor, size_ax = self.__image.GetSize() + extent_dict = { + "x": (0, size_cor, 0, size_ax), + "y": (0, size_sag, 0, size_ax), + "z": (0, size_sag, 0, size_cor), + } + # projection organisation if scalar.projection: projection = scalar.projection @@ -1223,6 +1276,7 @@ def _overlay_scalar_field(self): vmax=s_max, alpha=alpha, norm=norm, + extent=extent_dict["z"], ) cor_view = ax_cor.imshow( @@ -1236,6 +1290,7 @@ def _overlay_scalar_field(self): vmax=s_max, alpha=alpha, norm=norm, + extent=extent_dict["y"], ) sag_view = ax_sag.imshow( @@ -1249,6 +1304,7 @@ def _overlay_scalar_field(self): vmax=s_max, alpha=alpha, norm=norm, + extent=extent_dict["x"], ) # this is for (work-in-progress) dynamic visualisation @@ -1279,18 +1335,24 @@ def _overlay_scalar_field(self): asp = {"x": asp, "y": asp, "z": 1}[self.__axis] + if self.__axis == "z": + origin = {"normal": "upper", "reversed": "lower"}[self.__origin] + else: + origin = "lower" + s = return_slice(self.__axis, self.__cut) ax_view = ax.imshow( disp_img, interpolation="none", cmap=colormap, clim=(s_min, s_max), - origin="lower", + origin=origin, aspect=asp, vmin=s_min, vmax=s_max, alpha=alpha, norm=norm, + extent=extent_dict[self.__axis], ) if scalar.show_colorbar: