From efb80bbb2e6e6c438b78ae79e2d0d4d1fd3341c8 Mon Sep 17 00:00:00 2001 From: rnfinnegan Date: Tue, 25 May 2021 10:40:21 +1000 Subject: [PATCH] updates --- platipy/dicom/io/crawl.py | 10 +- platipy/imaging/generation/image.py | 13 ++- platipy/imaging/label/iar.py | 2 +- platipy/imaging/projects/cardiac/run.py | 59 ++++++++++-- platipy/imaging/utils/valve.py | 122 ++++++++++++++++++++++++ 5 files changed, 186 insertions(+), 20 deletions(-) create mode 100644 platipy/imaging/utils/valve.py diff --git a/platipy/dicom/io/crawl.py b/platipy/dicom/io/crawl.py index 55d9c61b..9f904186 100644 --- a/platipy/dicom/io/crawl.py +++ b/platipy/dicom/io/crawl.py @@ -667,7 +667,7 @@ def write_output_data_to_disk( # This will depend on the name format chosen # If there is a list, we append an index as we write to disk - if hasattr(field_list, "__iter__"): + if isinstance(field_list, (tuple, list)): # Iterate for suffix, file_to_write in enumerate(field_list): field_filename = field_filename_base + f"_{suffix}" @@ -680,10 +680,6 @@ def write_output_data_to_disk( field_filename = field_filename[:-1] # Save image! - """ - ! TO DO - Use pathlib, and perform some checks so we don"t overwrite anything! - """ output_name = ( pathlib.Path(output_directory) / parent_sorting_data @@ -882,7 +878,9 @@ def process_dicom_directory( try: study_uid_index = max(study_uid_dict.values()) + 1 except AttributeError: - study_uid_index = 0 + study_uid_index = 0 # Study UID dict might not exist + except ValueError: + study_uid_index = 0 # Study UID dict might be empty logger.info(f" Setting study instance UID index: {study_uid_index}") diff --git a/platipy/imaging/generation/image.py b/platipy/imaging/generation/image.py index 060e6b23..975729d1 100644 --- a/platipy/imaging/generation/image.py +++ b/platipy/imaging/generation/image.py @@ -39,9 +39,9 @@ def insert_sphere(arr, sp_radius=4, sp_centre=(0, 0, 0)): sp_radius_x, sp_radius_y, sp_radius_z = sp_radius arr_copy[ - ((x - sp_centre[0]) / sp_radius_x) ** 2 - + ((y - sp_centre[1]) / sp_radius_y) ** 2 - + ((z - sp_centre[2]) / sp_radius_z) ** 2 + ((x - sp_centre[0]) / sp_radius_x) ** 2.0 + + ((y - sp_centre[1]) / sp_radius_y) ** 2.0 + + ((z - sp_centre[2]) / sp_radius_z) ** 2.0 <= 1 ] = 1 @@ -53,7 +53,7 @@ def insert_sphere_image(image, sp_radius, sp_centre): Args: image (sitk.Image): Image in which to insert sphere - sp_radius (int, optional): The radius of the sphere. Defaults to 4. + sp_radius (int | list, optional): The radius of the sphere. Can also be defined as a vector. Defaults to 4. sp_centre (tuple, optional): The position at which the sphere should be inserted. Defaults to (0, 0, 0). @@ -61,7 +61,10 @@ def insert_sphere_image(image, sp_radius, sp_centre): np.array: An array with the sphere inserted """ - sp_radius_image = [sp_radius * i for i in image.GetSpacing()] + if not hasattr(sp_radius, "__iter__"): + sp_radius = [sp_radius] * 3 + + sp_radius_image = [i / j for i, j in zip(sp_radius, image.GetSpacing()[::-1])] arr = sitk.GetArrayFromImage(image) diff --git a/platipy/imaging/label/iar.py b/platipy/imaging/label/iar.py index 4b17d9e4..df450499 100644 --- a/platipy/imaging/label/iar.py +++ b/platipy/imaging/label/iar.py @@ -220,7 +220,7 @@ def run_iar( z_ideal = gaussian_curve(bin_centers, *popt) z_diff = np.abs(z_density - z_ideal) - except RuntimeError: + except (RuntimeError, ValueError): logger.debug("IAR couldnt fit curve, estimating with sampled statistics.") z_ideal = gaussian_curve(bin_centers, a=1, m=z_density.mean(), s=z_density.std()) z_diff = np.abs(z_density - z_ideal) diff --git a/platipy/imaging/projects/cardiac/run.py b/platipy/imaging/projects/cardiac/run.py index 6b324e06..e6c93653 100644 --- a/platipy/imaging/projects/cardiac/run.py +++ b/platipy/imaging/projects/cardiac/run.py @@ -118,6 +118,7 @@ "stop_condition_value_dict": {"LANTDESCARTERY_SPLINE": 1}, }, "return_as_cropped": False, + "return_additional_probability_dict": True, } @@ -135,6 +136,10 @@ def run_cardiac_segmentation(img, guide_structure=None, settings=CARDIAC_SETTING results = {} return_as_cropped = settings["return_as_cropped"] + return_additional_probability_dict = settings["return_additional_probability_dict"] + + if return_additional_probability_dict: + results_prob = {} """ Initialisation - Read in atlases @@ -269,12 +274,12 @@ def run_cardiac_segmentation(img, guide_structure=None, settings=CARDIAC_SETTING img_crop = crop_to_roi(img, crop_box_size, crop_box_index) - logger.info( - f"Calculated crop box\n\ - {crop_box_index}\n\ - {crop_box_size}\n\n\ - Volume reduced by factor {np.product(img.GetSize())/np.product(crop_box_size)}" - ) + logger.info( + f"Calculated crop box\n\ + {crop_box_index}\n\ + {crop_box_size}\n\n\ + Volume reduced by factor {np.product(img.GetSize())/np.product(crop_box_size)}" + ) """ Step 2 - Rigid registration of target images @@ -501,6 +506,8 @@ def run_cardiac_segmentation(img, guide_structure=None, settings=CARDIAC_SETTING """ template_img_binary = sitk.Cast((img * 0), sitk.sitkUInt8) + if return_additional_probability_dict: + template_img_prob = sitk.Cast((img * 0), sitk.sitkFloat64) vote_structures = settings["label_fusion_settings"]["optimal_threshold"].keys() @@ -515,8 +522,11 @@ def run_cardiac_segmentation(img, guide_structure=None, settings=CARDIAC_SETTING if return_as_cropped: results[structure_name] = binary_struct + if return_additional_probability_dict: + results_prob[structure_name] = probability_map + else: - paste_binary_img = sitk.Paste( + paste_img_binary = sitk.Paste( template_img_binary, binary_struct, binary_struct.GetSize(), @@ -524,7 +534,17 @@ def run_cardiac_segmentation(img, guide_structure=None, settings=CARDIAC_SETTING crop_box_index, ) - results[structure_name] = paste_binary_img + results[structure_name] = paste_img_binary + + if return_additional_probability_dict: + paste_prob_img = sitk.Paste( + template_img_prob, + probability_map, + probability_map.GetSize(), + (0, 0, 0), + crop_box_index, + ) + results_prob[structure_name] = paste_prob_img for structure_name in vessel_spline_settings["vessel_name_list"]: binary_struct = segmented_vessel_dict[structure_name] @@ -532,6 +552,12 @@ def run_cardiac_segmentation(img, guide_structure=None, settings=CARDIAC_SETTING if return_as_cropped: results[structure_name] = binary_struct + if return_additional_probability_dict: + results_prob[structure_name] = [ + atlas_set[atlas_id]["DIR"][structure_name] + for atlas_id in list(atlas_set.keys()) + ] + else: paste_img_binary = sitk.Paste( template_img_binary, @@ -543,7 +569,24 @@ def run_cardiac_segmentation(img, guide_structure=None, settings=CARDIAC_SETTING results[structure_name] = paste_img_binary + if return_additional_probability_dict: + vessel_list = [] + for atlas_id in list(atlas_set.keys()): + paste_img_binary = sitk.Paste( + template_img_binary, + atlas_set[atlas_id]["DIR"][structure_name], + atlas_set[atlas_id]["DIR"][structure_name].GetSize(), + (0, 0, 0), + crop_box_index, + ) + vessel_list.append(paste_img_binary) + + results_prob[structure_name] = vessel_list + if return_as_cropped: results["CROP_IMAGE"] = img_crop + if return_additional_probability_dict: + return results, results_prob + return results diff --git a/platipy/imaging/utils/valve.py b/platipy/imaging/utils/valve.py new file mode 100644 index 00000000..b6b713c7 --- /dev/null +++ b/platipy/imaging/utils/valve.py @@ -0,0 +1,122 @@ +# 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 + +from platipy.imaging.label.utils import get_com + +from platipy.imaging.generation.image import insert_sphere_image + +""" +Generate valves +First example - meeting point of LA/LV (mitral valve) +""" + + +def generate_valves_from_chambers(img_chamber_1, img_chamber_2, radius_mm=10): + + # Find the mid-point of the chambers + com_1 = np.array(get_com(img_chamber_1)) + com_2 = np.array(get_com(img_chamber_2)) + + midpoint = com_1 + 0.5 * (com_2 - com_1) + + # Define the overlap region (using binary dilation) + # Increment overlap to make sure we have enough voxels + dilation = 1 + overlap_vol = 0 + while overlap_vol <= 2000: + overlap = sitk.BinaryDilate(img_chamber_1, (dilation,) * 3) & sitk.BinaryDilate( + img_chamber_2, (dilation,) * 3 + ) + overlap_vol = np.sum(sitk.GetArrayFromImage(overlap) * np.product(overlap.GetSpacing())) + dilation += 1 + + print(f"Sufficient overlap found at dilation = {dilation} [V = {overlap_vol/1000:.2f} cm^3]") + + # Find the point in the overlap region closest to the mid-point + separation_vector_pixels = ( + np.stack(np.where(sitk.GetArrayFromImage(overlap))) - midpoint[:, None] + ) ** 2 + spacing = np.array(img_chamber_1.GetSpacing()) + separation_vector_mm = separation_vector_pixels / spacing[:, None] + + separation_mm = np.sum(separation_vector_mm, axis=0) + closest_overlap_point = np.argmin(separation_mm) + + com_valve = np.stack(np.where(sitk.GetArrayFromImage(overlap)))[:, closest_overlap_point] + + # Define the valve as a sphere + auto_valve = insert_sphere_image(0 * overlap, sp_radius=radius_mm, sp_centre=com_valve) + + return auto_valve + + +def generate_valves_from_vessels(vessel, thickness_mm=4, erosion_mm=2): + + # Thickness can be defined by (inferior_thickness, superior_thickness), + # or just total_thickness + if hasattr(thickness_mm, "__iter__"): + thickness_inferior, thickness_superior = thickness_mm + else: + thickness_inferior, thickness_superior = thickness_mm * 0.5, thickness_mm * 0.5 + + # Get the most inferior slice + arr_vessel = sitk.GetArrayFromImage(vessel) + inferior_slice = np.where(arr_vessel)[0].min() + + # Get interior and superior limits + filled_superior_slice = np.ceil( + np.where(arr_vessel)[0].min() + (thickness_superior / vessel.GetSpacing()[2]) + ).astype(int) + filled_inferior_slice = np.floor( + np.where(arr_vessel)[0].min() - (thickness_inferior / vessel.GetSpacing()[2]) + ).astype(int) + + # Create vessel interior + vessel_interior = sitk.BinaryErode(vessel, (erosion_mm, erosion_mm, 0)) + arr = sitk.GetArrayFromImage(vessel_interior) + + # Erase upper slices + arr[filled_superior_slice:, :, :] = 0 + + # Copy down (using mirroring about inferior slice) + for s_in, s_out in zip( + range(filled_inferior_slice, inferior_slice), + range(inferior_slice, filled_superior_slice)[::-1], + ): + arr[s_in, :, :] = arr[s_out, :, :] + + # Define the valve + auto_valve = sitk.GetImageFromArray(arr) + auto_valve.CopyInformation(vessel) + + # Post-processing + # 1. Extend the actual vessel downwards (continued) + arr_vessel[:inferior_slice, :, :] = arr_vessel[inferior_slice, :, :] + continued_vessel = sitk.GetImageFromArray(arr_vessel) + continued_vessel.CopyInformation(vessel) + + # 2. Erode this continued vessel + continued_vessel = sitk.BinaryErode(continued_vessel, (erosion_mm, erosion_mm, 0)) + + # 3. Mask + auto_valve = sitk.Mask(auto_valve, continued_vessel) + + # 4. Fill small holes + auto_valve = sitk.BinaryMorphologicalClosing(auto_valve, (1, 1, 1)) + + return auto_valve \ No newline at end of file