From de462a80995e72609a705741a96e0b61151b47eb Mon Sep 17 00:00:00 2001 From: Malte Date: Wed, 31 Jan 2024 02:02:00 +0100 Subject: [PATCH] added agatson option to aortic calcium --- bin/C2C | 12 ++- comp2comp/aortic_calcium/aortic_calcium.py | 79 +++++++++++++++---- .../aortic_calcium_visualization.py | 54 ++++++++++++- comp2comp/inference_pipeline.py | 3 +- 4 files changed, 127 insertions(+), 21 deletions(-) diff --git a/bin/C2C b/bin/C2C index 06233d1..82b2a00 100755 --- a/bin/C2C +++ b/bin/C2C @@ -128,7 +128,8 @@ def AorticCalciumPipelineBuilder(path, args): aortic_calcium.AorticCalciumMetrics(), aortic_calcium_visualization.AorticCalciumVisualizer(), aortic_calcium_visualization.AorticCalciumPrinter(), - ] + ], + args=args ) return pipeline @@ -201,8 +202,13 @@ def argument_parser(): "liver_spleen_pancreas", parents=[base_parser] ) - aortic_calcium = subparsers.add_parser("aortic_calcium", parents=[base_parser]) - + aortic_calcium = subparsers.add_parser( + "aortic_calcium", parents=[base_parser]) + + aortic_calcium.add_argument( + "--threshold", default="adaptive", type=str + ) + contrast_phase_parser = subparsers.add_parser( "contrast_phase", parents=[base_parser] ) diff --git a/comp2comp/aortic_calcium/aortic_calcium.py b/comp2comp/aortic_calcium/aortic_calcium.py index abfda49..bfe9722 100644 --- a/comp2comp/aortic_calcium/aortic_calcium.py +++ b/comp2comp/aortic_calcium/aortic_calcium.py @@ -126,7 +126,8 @@ def __call__(self, inference_pipeline): spine_mask = inference_pipeline.spine_segmentation.get_fdata() > 0 calcification_results = self.detectCalcifications( - ct, aorta_mask, exclude_mask=spine_mask, remove_size=3, return_dilated_mask=True + ct, aorta_mask, exclude_mask=spine_mask, remove_size=3, return_dilated_mask=True, + threshold=inference_pipeline.args.threshold ) inference_pipeline.calc_mask = calcification_results['calc_mask'] @@ -182,6 +183,7 @@ def detectCalcifications( exclude_center_aorta=True, return_eroded_aorta=False, aorta_erode_iteration=6, + threshold = 'adaptive', ): """ Function that takes in a CT image and aorta segmentation (and optionally volumes to use @@ -222,6 +224,9 @@ def detectCalcifications( Return the eroded center aorta. Defaults to False. aorta_erode_iteration (int, optional): Number of iterations for the strcturing element. Defaults to 6. + threshold: (str, int): + Can either be 'adaptive', 'agatson', or int. Choosing 'agatson' + Will mean a threshold of 130 HU. Returns: results: array of only the mask is returned, or dict if other volumes are also returned. @@ -268,16 +273,25 @@ def slicedDilationOrErosion(input_mask, struct, num_iteration, operation): biggest_idx = np.argmax(aorta_vols) + 1 aorta_mask[labelled_aorta != biggest_idx] = 0 - - # Get aortic CT point to set adaptive threshold - aorta_ct_points = ct[aorta_mask == 1] - - # equal to one standard deviation to the left of the curve - quant = 0.158 - quantile_median_dist = np.median(aorta_ct_points) - np.quantile( - aorta_ct_points, q=quant - ) - calc_thres = np.median(aorta_ct_points) + quantile_median_dist * num_std + + ### Choose the threshold ### + if threshold == 'adaptive': + # Get aortic CT point to set adaptive threshold + aorta_ct_points = ct[aorta_mask == 1] + + # equal to one standard deviation to the left of the curve + quant = 0.158 + quantile_median_dist = np.median(aorta_ct_points) - np.quantile( + aorta_ct_points, q=quant + ) + calc_thres = np.median(aorta_ct_points) + quantile_median_dist * num_std + elif threshold == 'agatson': + calc_thres = 130 + elif isinstance(threshold, int): + calc_thres = threshold + else: + raise ValueError('Error in threshold value for aortic calcium segmentaiton. \ + Should be \'adaptive\', \'agatson\' or int, but got: ' + str(threshold)) t0 = time.time() @@ -448,15 +462,17 @@ def __call__(self, inference_pipeline): ct_full = inference_pipeline.medical_volume.get_fdata() - for i in range(2): + for i in range(len(region_names)): # count statistics for individual calcifications if i == 0: - labelled_calc, num_lesions = ndimage.label(calc_mask[:,:,:sep_plane]) + calc_mask_region = calc_mask[:,:,:sep_plane] ct = ct_full[:,:,:sep_plane] elif i == 1: - labelled_calc, num_lesions = ndimage.label(calc_mask[:,:,sep_plane:]) + calc_mask_region = calc_mask[:,:,sep_plane:] ct = ct_full[:,:,sep_plane:] + labelled_calc, num_lesions = ndimage.label(calc_mask_region) + metrics = { "volume": [], "mean_hu": [], @@ -481,9 +497,42 @@ def __call__(self, inference_pipeline): metrics["volume_total"] = calc_vol metrics["num_calc"] = len(metrics["volume"]) - + + if inference_pipeline.args.threshold == 'agatson': + metrics["agatson_score"] = self.CalculateAgatsonScore(calc_mask_region, ct, inference_pipeline.pix_dims) + all_regions[region_names[i]] = metrics inference_pipeline.metrics = all_regions return {} + + def CalculateAgatsonScore(self, calc_mask_region, ct, pix_dims): + def get_hu_factor(max_hu): + # if max_hu >< + if max_hu < 200: + factor = 1 + elif 200 <= max_hu < 300: + factor = 2 + elif 300 <= max_hu < 400: + factor = 3 + elif max_hu > 400: + factor = 4 + + return factor + + # dims are in mm here + area_per_pixel = pix_dims[0] * pix_dims[1] + agatson = 0 + + for i in range(calc_mask_region.shape[2]): + tmp_slice = calc_mask_region[:,:,i] + tmp_ct_slice = ct[:,:,i] + + labelled_calc, num_lesions = ndimage.label(tmp_slice) + + for j in range(1, num_lesions + 1): + tmp_mask = labelled_calc == j + agatson += tmp_mask.sum() * area_per_pixel * get_hu_factor(tmp_ct_slice[tmp_mask].max()) + + return agatson \ No newline at end of file diff --git a/comp2comp/aortic_calcium/aortic_calcium_visualization.py b/comp2comp/aortic_calcium/aortic_calcium_visualization.py index ba11ec4..7d89be6 100644 --- a/comp2comp/aortic_calcium/aortic_calcium_visualization.py +++ b/comp2comp/aortic_calcium/aortic_calcium_visualization.py @@ -57,6 +57,7 @@ def __call__(self, inference_pipeline): ): f.write("{},{:.1f},{:.1f},{:.1f}\n".format(vol, mean, median, max)) + # Write total results with open( os.path.join( @@ -68,8 +69,49 @@ def __call__(self, inference_pipeline): f.write(region + ',\n') f.write("Total number,{}\n".format(metrics["num_calc"])) - f.write("Total volume (cm^3),{}\n".format(metrics["volume_total"])) - f.write("Threshold (HU),{}\n".format(inference_pipeline.calcium_threshold)) + f.write("Total volume (cm^3),{:.1f}\n".format(metrics["volume_total"])) + f.write("Threshold (HU),{:.1f}\n".format(inference_pipeline.calcium_threshold)) + + f.write("{},{:.1f}+/-{:.1f}\n".format( + "Mean HU", + np.mean(metrics["mean_hu"]), + np.std(metrics["mean_hu"]), + ) + ) + f.write("{},{:.1f}+/-{:.1f}\n".format( + "Median HU", + np.mean(metrics["median_hu"]), + np.std(metrics["median_hu"]), + ) + ) + f.write("{},{:.1f}+/-{:.1f}\n".format( + "Max HU", + np.mean(metrics["max_hu"]), + np.std(metrics["max_hu"]), + ) + ) + f.write("{},{:.1f}+/-{:.1f}\n".format( + "Mean volume (cm³):", + np.mean(metrics["volume"]), + np.std(metrics["volume"]), + ) + ) + f.write("{},{:.1f}\n".format( + "Median volume (cm³)", np.median(metrics["volume"]) + ) + ) + f.write("{},{:.1f}\n".format( + "Max volume (cm³)", np.max(metrics["volume"]) + ) + ) + f.write("{},{:.1f}\n".format( + "Min volume (cm³):", np.min(metrics["volume"]) + ) + ) + + if inference_pipeline.args.threshold == 'agatson': + f.write("Agatson score,{:.1f}\n".format(metrics["agatson_score"])) + distance = 25 print("\n") @@ -139,7 +181,15 @@ def __call__(self, inference_pipeline): "Threshold (HU):", distance, inference_pipeline.calcium_threshold ) ) + if inference_pipeline.args.threshold == 'agatson': + print( + "{:<{}}{:.1f}".format( + "Agatson score:", distance, metrics["agatson_score"] + ) + ) + print("\n") return {} + diff --git a/comp2comp/inference_pipeline.py b/comp2comp/inference_pipeline.py index a70175b..c424b62 100644 --- a/comp2comp/inference_pipeline.py +++ b/comp2comp/inference_pipeline.py @@ -14,8 +14,9 @@ class InferencePipeline(InferenceClass): """Inference pipeline.""" - def __init__(self, inference_classes: List = None, config: Dict = None): + def __init__(self, inference_classes: List = None, config: Dict = None, args: Dict = None): self.config = config + self.args = args # assign values from config to attributes if self.config is not None: for key, value in self.config.items():