save distance maps to NiFti format

Author: rmsandu

calc_margin.py

@@ -3,7 +3,7 @@
 import numpy as np
 import pandas as pd
 
-
+import scripts.plot_ablation_margin_hist as pm
 from customradiomics.margin import compute_distances
 from utils.niftireader import load_image
 
@@ -36,33 +36,46 @@ def get_args():
         ap.add_argument("-l", "--liver", required=True, help="path to the liver segmentation")
         ap.add_argument("-i", "--lesion-id", required=True, help="lesion id")
         ap.add_argument("-p", "--patient-id", required=True, help="patient id from study")
-        ap.add_argument("-o", "--OUTPUT",  required=True, help="output file (csv)")
+        ap.add_argument("-o", "--OUTPUT", required=True, help="output file (csv)")
         args = vars(ap.parse_args())
     return args
 
 
 if __name__ == '__main__':
     args = get_args()
-
+    # -t
+    # "data\B04\01\B04_L01_Tumor.nii.gz" - a
+    # "data\B04\01\B04_L01_Ablation.nii.gz" - l
+    # "data\B04\01\B04_L01_Liver.nii.gz" - i
+    # 1 - p
+    # "B04" - o
+    # "output111.csv"
     tumor_file = args['tumor']
     ablation_file = args['ablation']
     liver_file = args['liver']
     patient_id = args['patient_id']
     lesion_id = args['lesion_id']
     output_file = args['OUTPUT']
-
+    rootdir = r"C:\develop\segmentation-eval\figures"
     tumor, tumor_np = load_image(tumor_file)
     ablation, ablation_np = load_image(ablation_file)
     liver, liver_np = load_image(liver_file)
+    affine_transform = tumor.affine
 
     has_liver_segmented = np.sum(liver_np.astype(np.uint8)) > 0
 
     pixdim = liver.header['pixdim']
     spacing = (pixdim[1], pixdim[2], pixdim[3])
-    surface_distance = compute_distances(tumor_np, ablation_np,
+    # compute_distances(mask_gt, mask_pred, exclusion_zone, spacing_mm, connectivity=1, crop=True, exclusion_distance=5)
+    surface_distance = compute_distances(tumor_np, ablation_np, affine_transform,
                                          exclusion_zone=liver_np if has_liver_segmented else None,
                                          spacing_mm=spacing, connectivity=1)
 
+    pm.plot_histogram_surface_distances(patient_id, lesion_id, rootdir, distance_map=surface_distance['distances_gt_to_pred'],
+                                        num_voxels=len(surface_distance['distances_gt_to_pred']),
+                                        title='SurfaceDistance_Histogram_B04_L1', ablation_date="20160101",
+                                        flag_to_plot=True)
+
     df = pd.DataFrame(data={
         'Patient': [patient_id] * len(surface_distance['distances_gt_to_pred']),
         'Lesion': [lesion_id] * len(surface_distance['distances_gt_to_pred']),

customradiomics/margin.py

@@ -1,6 +1,7 @@
 import numpy as np
 from scipy import ndimage
-
+import matplotlib.pyplot as plt
+import nibabel as nib
 
 def compute_bounding_box(mask_gt, mask_pred, exclusion_zone):
     mask_all = mask_gt | mask_pred
@@ -27,25 +28,26 @@ def compute_bounding_box(mask_gt, mask_pred, exclusion_zone):
     bbox_min[2] = np.min(idx_nonzero_2)
     bbox_max[2] = np.max(idx_nonzero_2)
 
-    # crop the processing subvolume.
-    # we need to zeropad the cropped region with 1 voxel at the lower,
-    # the right and the back side. This is required to obtain the "full"
-    # convolution result with the 2x2x2 kernel
-    # cropmask = np.zeros((bbox_max - bbox_min)+2, np.uint8)
 
     return bbox_min, bbox_max
 
 
 def crop_mask(mask, bbox_min, bbox_max):
+
     cropmask = np.zeros((bbox_max - bbox_min) + 2)
 
     cropmask[0:-1, 0:-1, 0:-1] = mask[bbox_min[0]:bbox_max[0] + 1,
                                  bbox_min[1]:bbox_max[1] + 1,
                                  bbox_min[2]:bbox_max[2] + 1]
+    # TODO:  This is correct only if the object is interior to the bounding box. Must perform testing.
+    # crop the processing subvolume.
+    # we need to zeropad the cropped region with 1 voxel at the lower,
+    # the right and the back side. This is required to obtain the "full"
+    # convolution result with the 2x2x2 kernel
     return cropmask
 
 
-def compute_distances(mask_gt, mask_pred, exclusion_zone, spacing_mm, connectivity=1, crop=True, exclusion_distance=5):
+def compute_distances(mask_gt, mask_pred, affine_transform, exclusion_zone, spacing_mm, connectivity=1, crop=True, exclusion_distance=5):
     if crop:
         if exclusion_zone is not None:
             bbox_min, bbox_max = compute_bounding_box(mask_gt, mask_pred, exclusion_zone)
@@ -58,16 +60,24 @@ def compute_distances(mask_gt, mask_pred, exclusion_zone, spacing_mm, connectivi
             mask_pred = crop_mask(mask_pred, bbox_min, bbox_max).astype(np.bool)
 
     border_inside = ndimage.binary_erosion(mask_gt, structure=ndimage.generate_binary_structure(3, connectivity))
-    borders_gt = mask_gt ^ border_inside
+    ni_img = nib.Nifti1Image(border_inside.astype(np.uint8), affine_transform)
+    nib.save(ni_img, 'border_inside_tumor.nii.gz')
+    borders_gt = mask_gt ^ border_inside  # contours
+    ni_img = nib.Nifti1Image(borders_gt.astype(np.uint8), affine_transform)
+    nib.save(ni_img, 'contour_tumor.nii.gz')
     border_inside = ndimage.binary_erosion(mask_pred, structure=ndimage.generate_binary_structure(3, connectivity))
-    borders_pred = mask_pred ^ border_inside
-
-    # compute the distance transform (closest distance of each voxel to the
-    # surface voxels)
+    borders_pred = mask_pred ^ border_inside  # contours
+    ni_img = nib.Nifti1Image(borders_pred.astype(np.uint8), affine_transform)
+    nib.save(ni_img, 'contour_ablation.nii.gz')
+    # compute the distance transform (closest distance of each voxel to the surface voxels)
     if borders_gt.any():
         distmap_gt = ndimage.morphology.distance_transform_edt(
             ~borders_gt, sampling=spacing_mm)
+        ni_img = nib.Nifti1Image(distmap_gt.astype(np.uint8), affine_transform)
+        nib.save(ni_img, 'distmap_tumor.nii.gz')
         distmask_gt = mask_gt.astype(np.int8)
+        # todo: visualize with matplotlib
+        # mark outside zones with -1
         distmask_gt[distmask_gt == 0] = -1
         distmap_gt *= distmask_gt
     else:
@@ -76,7 +86,10 @@ def compute_distances(mask_gt, mask_pred, exclusion_zone, spacing_mm, connectivi
     if borders_pred.any():
         distmap_pred = ndimage.morphology.distance_transform_edt(
             ~borders_pred, sampling=spacing_mm)
+        ni_img = nib.Nifti1Image(distmap_pred, affine_transform)
+        nib.save(ni_img, 'distmap_ablation.nii.gz')
         distmask_pred = mask_pred.astype(np.int8)
+        # mark outside zones with -1
         distmask_pred[distmask_pred == 0] = -1
         distmap_pred *= distmask_pred
     else:
@@ -90,6 +103,7 @@ def compute_distances(mask_gt, mask_pred, exclusion_zone, spacing_mm, connectivi
             ~borders_exclusion, sampling=spacing_mm)
 
         distmask_exclusion = exclusion_zone.astype(np.int8)
+        # mark outside zones with -1
         distmask_exclusion[distmask_exclusion == 0] = -1
         distmap_exclusion *= distmask_exclusion