Airways

Pipfile

@@ -16,8 +16,7 @@ SimpleITK="*"
 tqdm="*"
 fill_voids="*"
 zipfile36="*"
+vtk = "*"
 
 [requires]
 python_version = "3.8"
-
-

airw.py

@@ -0,0 +1,135 @@
+# segment vascular tree and airways
+
+import SimpleITK as sitk
+import os
+import numpy as np
+import argparse
+import vtk
+
+# import itk
+
+
+def connected_threshold(nrrd_image, seeds, tresh):
+    print("reading image")
+    segmentationfilter = sitk.ConnectedThresholdImageFilter()
+    segmentationfilter.SetUpper(tresh[0])
+    segmentationfilter.SetLower(tresh[1])
+    segmentationfilter.SetReplaceValue(1)
+    for i in range(0, len(seeds), 3):
+        seed = seeds[i : i + 3]
+        print(seed)
+        segmentationfilter.AddSeed(seed)
+    print("run filter")
+    result = segmentationfilter.Execute(nrrd_image)
+    return result
+
+
+if __name__ == "__main__":
+    # init arg parser
+
+    parser = argparse.ArgumentParser(description="")
+
+    parser.add_argument(
+        "--image",
+        action="store",
+        default=None,
+        help="the image",
+    )
+
+    parser.add_argument(
+        "--mask",
+        action="store",
+        default=None,
+        help="the mask",
+    )
+
+    parser.add_argument(
+        "--seeds",
+        action="store",
+        nargs="+",
+        type=int,
+        help="array of xyz",
+        default=[],
+    )
+
+    parser.add_argument(
+        "--thresholds",
+        action="store",
+        nargs="+",
+        type=int,
+        help="array of upper and lower tr",
+        default=[],
+    )
+
+    args = parser.parse_args()
+
+print(args.seeds)
+print(args.thresholds)
+
+img = sitk.ReadImage(args.image)
+out_img_path = "./airw.nrrd"
+# tree = connected_threshold(img, args.seeds, args.thresholds)
+
+seg = sitk.ConfidenceConnected(
+    img,
+    seedList=[args.seeds],
+    numberOfIterations=10,
+    multiplier=2.0,
+    initialNeighborhoodRadius=10,
+    replaceValue=1,
+)
+
+# initial sharpening ?
+# erosion, then conn tr ?
+
+# mask = sitk.ReadImage(args.mask)
+# mask = sitk.GetArrayFromImage(mask)
+# mask[mask > 1] = 1
+# tree_arr = sitk.GetArrayFromImage(tree)
+# tree_arr[mask == 1] = 0
+
+# outttt = sitk.GetImageFromArray(tree_arr)
+
+# thresholder = sitk.BinaryThresholdImageFilter()
+# thresholder.SetLowerThreshold(0)
+# # thresholder.SetUpperThreshold()
+# thresholder.SetInsideValue(1)
+# thresholder.SetOutsideValue(0)
+# tree_mask = thresholder.Execute(img)
+
+# # closing
+# closing = sitk.BinaryMorphologicalClosingImageFilter()
+# closing.SetKernelRadius(1)
+# closing.SetForegroundValue(1)
+# closed = closing.Execute(seg)
+
+# # opening
+opening = sitk.BinaryMorphologicalOpeningImageFilter()
+opening.SetKernelRadius(2)
+opening.SetForegroundValue(1)
+opened = opening.Execute(seg)
+
+# # closing
+closing = sitk.BinaryMorphologicalClosingImageFilter()
+closing.SetKernelRadius(1)
+closing.SetForegroundValue(1)
+closed = closing.Execute(opened)
+
+# smooth
+# smoother = sitk.MedianImageFilter()
+# smoother.SetRadius([0, 0, 1])
+# smooth = smoother.Execute(opened)
+
+# gaussian = sitk.SmoothingRecursiveGaussianImageFilter()
+# gaussian.SetSigma(1)
+# smooth = gaussian.Execute(opened)
+
+# surface extraction
+
+
+out = closed
+out.SetSpacing(img.GetSpacing())
+out.SetDirection(img.GetDirection())
+out.SetOrigin(img.GetOrigin())
+print("writing image")
+sitk.WriteImage(out, out_img_path, True)

divide.py

@@ -0,0 +1,126 @@
+import SimpleITK as sitk
+import numpy as np
+import argparse
+import os
+
+# divide image with same serie id but different acquisition number in the same folder
+# write img_1 and img_2 in output
+
+if __name__ == "__main__":
+    # init arg parser
+
+    parser = argparse.ArgumentParser(
+        description=""
+    )
+
+    parser.add_argument(
+        "--image",
+        action="store",
+        default=None,
+        help="the target image",
+    )
+
+    parser.add_argument(
+        "--dicomdir",
+        action="store",
+        default=None,
+        help="the dicom folder of target image",
+    )
+
+    parser.add_argument(
+        "--mask",
+        action="store",
+        default=None,
+        help="the mask image",
+    )
+
+    args = parser.parse_args()
+
+    # image = sitk.ReadImage(args.image)
+    # mask = sitk.ReadImage(args.mask)
+
+    # crop = sitk.CropImageFilter()
+    # crop.SetUpperBoundaryCropSize([0,0,23])
+    # crop.SetLowerBoundaryCropSize([0,0,0])
+    # res = crop.Execute(image)
+
+
+    # res = image[:, :, 200:423]
+
+    # sitk.WriteImage(res, "./res.nrrd", True)
+
+    # img = sitk.GetImageFromArray(image)
+    # msk = sitk.GetImageFromArray(mask)
+
+    # print('done conversion')
+
+    # out_img = img[msk]
+
+    # print(out_img)
+
+    def getImageSeriesId(file_name, series_list, desc_list):
+        print("Reading image...")
+        # A file name that belongs to the series we want to read
+
+        # Read the file's meta-information without reading bulk pixel data
+        # print('Reading image...')
+        file_reader = sitk.ImageFileReader()
+        file_reader.SetFileName(file_name)
+
+        try:
+            file_reader.ReadImageInformation()
+        except:
+            print("ERROR while reading: ", file_name)
+            print("SKIP file")
+            return
+
+        # Get the sorted file names, opens all files in the directory and reads the meta-information
+        # without reading the bulk pixel data
+        series_ID = file_reader.GetMetaData("0020|000e")
+        description = file_reader.GetMetaData("0008|103e")
+        # print('seriesId', series_ID, '\t\t descr', description)
+
+        if series_ID not in series_list:
+            series_list.append(series_ID)
+            desc_list.append(description)
+
+        return series_ID
+
+    def divideByAcqNumber(file_names):
+        seriesObject = [[],[]]
+
+        for file_name in file_names:
+            file_reader = sitk.ImageFileReader()
+            file_reader.SetFileName(file_name)
+
+            try:
+                file_reader.ReadImageInformation()
+            except:
+                print("ERROR while reading: ", file_name)
+                print("SKIP file")
+                return
+
+            acqNumber = file_reader.GetMetaData("0020|0012")
+            print(acqNumber, file_name)
+            seriesObject[int(acqNumber)-1].append(file_name)
+
+        return seriesObject
+
+    path_image = args.dicomdir
+
+    for (root, dirs, files) in os.walk(path_image):
+        series_id = getImageSeriesId(os.path.join(root, files[0]), [], [])
+
+        sorted_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
+            path_image, series_id
+        )
+
+    obj = divideByAcqNumber(sorted_file_names)
+
+    image1 = sitk.ReadImage(obj[0])
+    image2 = sitk.ReadImage(obj[1])
+
+    sitk.WriteImage(image1, "./image1.nrrd", True)
+    sitk.WriteImage(image2, "./image2.nrrd", True)
+
+

fillMask.py

@@ -0,0 +1,87 @@
+import SimpleITK as sitk
+import numpy as np
+import argparse
+
+if __name__ == "__main__":
+    # init arg parser
+
+    parser = argparse.ArgumentParser(description="")
+
+    parser.add_argument(
+        "--mask",
+        action="store",
+        default=None,
+        help="the binary lung mask",
+    )
+
+    parser.add_argument(
+        "--image",
+        action="store",
+        default=None,
+        help="the original image",
+    )
+
+    args = parser.parse_args()
+
+    mask = sitk.ReadImage(args.mask)
+    mask_arr = sitk.GetArrayFromImage(mask)
+
+    out_mask = mask
+
+    dims = mask.GetSize()
+
+    def fill_line(k, j, arr_in, arr_out):
+        fill = False
+        for i in range(arr_in.shape[2]):
+            px0 = arr_in[k, j, i]
+            px1 = arr_in[k, j, i + 1]
+            if px0 == 1 and px1 == 0:
+                fill = True
+            if fill:
+                arr_out[k, j, i] = 3
+            if fill and px0 == 0 and px1 == 2:  # right lung value is 2
+                break
+
+    llv = 1
+    rlv = 2
+
+    out_arr = np.copy(mask_arr)
+
+    for k in range(mask_arr.shape[0]):
+        for j in range(mask_arr.shape[1]):
+            row = mask_arr[k, j, :]
+            if (llv in row) and (rlv in row):
+                print(k, j)
+                fill_line(k, j, mask_arr, out_arr)
+
+    # binarize!
+    out_arr[out_arr == 2] = 1
+    out_arr[out_arr == 3] = 1
+
+    spacing = mask.GetSpacing()
+    direction = mask.GetDirection()
+    origin = mask.GetOrigin()
+
+    out_mask = sitk.GetImageFromArray(out_arr)
+    out_mask.SetSpacing(spacing)
+    out_mask.SetDirection(direction)
+    out_mask.SetOrigin(origin)
+
+    # close little holes
+    filler = sitk.BinaryFillholeImageFilter()
+    filler.SetForegroundValue(1)
+    filled_mask = filler.Execute(out_mask)
+
+    sitk.WriteImage(filled_mask, "./filled_mask.nrrd")
+
+    if args.image:
+        img = sitk.ReadImage(args.image)
+        img_arr = sitk.GetArrayFromImage(img)
+        o_arr = np.zeros(img_arr.shape)
+        o_arr = o_arr - 1500  # background value
+        o_arr[out_arr == 1] = img_arr[out_arr == 1]
+        out_img = sitk.GetImageFromArray(o_arr)
+        out_img.SetDirection(img.GetDirection())
+        out_img.SetSpacing(img.GetSpacing())
+        out_img.SetOrigin(img.GetOrigin())
+        sitk.WriteImage(out_img, "./filled_lungs.nrrd")