{"cells":[{"metadata":{},"cell_type":"markdown","source":"This kernel provides an environment for experiments with DICOM images.\n\ncurrent pipeline:\n* crop if size > 512\n* resize to (512, 512)\n* convert to Hounsfield Units\n* clip values (-1000, 500)\n* lung masks\n* norm\n* grayscale png 512\n\nThe result dataset link is available in the last cell.\n\nreferences:\n* https://www.kaggle.com/carlossouza/osic-autoencoder-training\n* https://www.kaggle.com/aadhavvignesh/lung-segmentation-by-marker-controlled-watershed\n* https://www.kaggle.com/gunesevitan/osic-pulmonary-fibrosis-progression-eda","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-output":true,"collapsed":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm -y\nimport gdcm # required for ID00011637202177653955184","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nimport matplotlib.animation as animation\nfrom IPython.display import HTML\n\nimport pydicom\nimport cv2\nfrom skimage import measure, morphology, segmentation\nimport scipy.ndimage as ndimage\n\nfrom multiprocessing import Pool\nfrom tqdm.notebook import tqdm\n\nDCOM_DIR = '../input/osic-pulmonary-fibrosis-progression/train'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntrain = train.drop_duplicates(keep=False, subset=['Patient','Weeks'])\ntrain","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# copypaste https://www.kaggle.com/carlossouza/osic-autoencoder-training\n\ndef seperate_lungs(image, min_hu, iterations):\n    h, w = image.shape[0], image.shape[1]\n\n    marker_internal, marker_external, marker_watershed = generate_markers(image)\n\n    # Sobel-Gradient\n    sobel_filtered_dx = ndimage.sobel(image, 1)\n    sobel_filtered_dy = ndimage.sobel(image, 0)\n    sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)\n    sobel_gradient *= 255.0 / np.max(sobel_gradient)\n\n    watershed = morphology.watershed(sobel_gradient, marker_watershed)\n\n    outline = ndimage.morphological_gradient(watershed, size=(3,3))\n    outline = outline.astype(bool)\n\n    # Structuring element used for the filter\n    blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [0, 0, 1, 1, 1, 0, 0]]\n\n    blackhat_struct = ndimage.iterate_structure(blackhat_struct, iterations)\n\n    # Perform Black Top-hat filter\n    outline += ndimage.black_tophat(outline, structure=blackhat_struct)\n\n    lungfilter = np.bitwise_or(marker_internal, outline)\n    lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)\n\n    segmented = np.where(lungfilter == 1, image, min_hu * np.ones((h, w)))\n\n    return segmented  #, lungfilter, outline, watershed, sobel_gradient\n\ndef generate_markers(image, threshold=-400):\n    h, w = image.shape[0], image.shape[1]\n\n    marker_internal = image < threshold\n    marker_internal = segmentation.clear_border(marker_internal)\n    marker_internal_labels = measure.label(marker_internal)\n\n    areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n    areas.sort()\n\n    if len(areas) > 2:\n        for region in measure.regionprops(marker_internal_labels):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:\n                    marker_internal_labels[coordinates[0], coordinates[1]] = 0\n\n    marker_internal = marker_internal_labels > 0\n\n    # Creation of the External Marker\n    external_a = ndimage.binary_dilation(marker_internal, iterations=10)\n    external_b = ndimage.binary_dilation(marker_internal, iterations=55)\n    marker_external = external_b ^ external_a\n\n    # Creation of the Watershed Marker\n    marker_watershed = np.zeros((h, w), dtype=np.int)\n    marker_watershed += marker_internal * 255\n    marker_watershed += marker_external * 128\n\n    return marker_internal, marker_external, marker_watershed","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def _parts(pid:int, dcom_dir=DCOM_DIR)->list:\n    return sorted([int(i.split('.')[0]) for i in os.listdir(os.path.join(dcom_dir, pid))])\n\ndef _crop(s, size=512):\n    if (s.shape[0]<=size):\n        return s\n    s_cropped = s[~np.all(s == 0, axis=1)]\n    s_cropped = s_cropped[:, ~np.all(s_cropped == 0, axis=0)]\n    return s_cropped\n\ndef _resize(s, size=512):\n    if (s.shape[0]<=size):\n        return s\n    s_resized = cv2.resize(s, (0,0), fx=size/s.shape[0], fy=size/s.shape[1], interpolation=cv2.INTER_AREA)        \n    return s_resized\n\ndef _clip(s, lo=-1000, hi=500):\n    s[s<lo]=lo\n    s[s>hi]=hi\n    return s\n\ndef _hu(s, slope:float, intercept:int):     \n    s_hu = (s * slope + intercept).astype(np.int16)\n    return s_hu\n\ndef _norm(s, lo=-1000, hi=500):\n    s_normed= s.astype(np.float)\n    s_normed = (s_normed-lo)*255.0/(hi-lo)\n    return s_normed.astype(np.int16)\n\ndef _meta(dcom):\n    return {\n        #'Modality': dcom.Modality, # const\n        'Manufacturer' : dcom.Manufacturer,\n        'ManufacturerModelName' : dcom.ManufacturerModelName,\n        #'BodyPartExamined': dcom.BodyPartExamined, # const\n        'SliceThickness' : dcom.SliceThickness,\n        'KVP' : dcom.KVP,\n        'TableHeight' : dcom.TableHeight,\n        #'RotationDirection' : dcom.RotationDirection, # const       \n        'ConvolutionKernel' : dcom.ConvolutionKernel,\n        'PatientPosition' : dcom.PatientPosition,             \n        #'PhotometricInterpretation': dcom.PhotometricInterpretation, # const\n        #'SamplesPerPixel': dcom.SamplesPerPixel, # const\n        #'BitsAllocated': dcom.BitsAllocated, # const\n        'BitsStored': dcom.BitsStored,\n        'HighBit': dcom.HighBit,\n        'PixelRepresentation': dcom.PixelRepresentation,\n        'PixelSpacing0': dcom.PixelSpacing[0],\n        'PixelSpacing1': dcom.PixelSpacing[1],\n        'WindowCenter': dcom.WindowCenter,\n        'WindowWidth': dcom.WindowWidth,\n        'RescaleIntercept': dcom.RescaleIntercept,\n        #'RescaleSlope': dcom.RescaleSlope # const\n    }","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\ndef dicom_cube(pid):\n    arr = []\n    for part in _parts(pid)[200:250]:\n        dicom = pydicom.dcmread(os.path.join(DCOM_DIR, pid, f'{part}.dcm'))\n        img = dicom.pixel_array        \n        meta = _meta(dicom)\n        img = _crop(img, size=512)\n        img = _hu(img, slope=1., intercept=meta['RescaleIntercept'])\n        img = _resize(img, size=512)\n        img = _clip(img, lo=-1000, hi=500)\n        img = seperate_lungs(img, min_hu=-1000, iterations=1) # time consuming operation\n        img = _norm(img, lo=-1000, hi=500)\n        arr.append(img)\n    return np.array(arr)\n\ncube = dicom_cube('ID00067637202189903532242')\n\nfig = plt.figure(figsize=(8,8))\n\nims = []\nfor i in cube:\n    im = plt.imshow(i, animated=True, cmap='gray')\n    plt.axis('off')\n    ims.append([im])\n\nani = animation.ArtistAnimation(fig, ims, interval=50)\nplt.close()\n\nHTML(ani.to_jshtml())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"OUT_DIR = '/tmp/osic'\n\ntrash = [\n    os.path.join(DCOM_DIR, 'ID00105637202208831864134', '1.dcm'), # black\n    os.path.join(DCOM_DIR, 'ID00052637202186188008618', '4.dcm'), # corrupted\n]\n\ndef process_one(pid):\n    os.makedirs(os.path.join(OUT_DIR, pid), exist_ok=True)    \n    metas = []\n    parts = _parts(pid)\n    for part in parts:\n        path = os.path.join(DCOM_DIR, pid, f'{part}.dcm')\n        if path in trash:\n            continue\n        \n        dicom = pydicom.dcmread(path)\n        img = dicom.pixel_array        \n        meta = _meta(dicom)\n        img = _crop(img, size=512)\n        img = _hu(img, slope=1., intercept=meta['RescaleIntercept'])\n        img = _resize(img, size=512)\n        img = _clip(img, lo=-1000, hi=500)\n        img = seperate_lungs(img, min_hu=-1000, iterations=1)\n        img = _norm(img, lo=-1000, hi=500)     \n        \n        meta['Patient'] = pid\n        meta['Part'] = part\n        meta['Min'] = np.min(img)\n        meta['Max'] = np.max(img)\n        meta['Mean'] = np.mean(img)\n        meta['Std'] = np.std(img)\n        metas.append(meta)\n        \n        res = cv2.imwrite(os.path.join(OUT_DIR, pid, f'{part}.png'), img)\n        \n    return metas\n\n#metas = process_one('ID00067637202189903532242')\n#metas[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"batch = train.Patient.unique()#[:5]\nwith Pool(processes=4) as pool:\n    res = list(\n        tqdm(pool.imap(process_one, list(batch)), total = len(batch))\n    )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"stack = []\nfor i in range(len(res)):\n    for j in range(len(res[i])):\n        stack.append(res[i][j])\n\ndf = pd.DataFrame(stack)\ndf.loc[df.Manufacturer == 'GE MEDICAL SYSTEMS', 'Manufacturer'] = 'GE'\ndf.loc[df.Manufacturer == 'Hitachi Medical Corporation', 'Manufacturer'] = 'Hitachi'\ndf.ManufacturerModelName.fillna('TOSHIBA', inplace=True)\ndf","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"mykaggleapi is a private dataset. https://github.com/Kaggle/kaggle-api for more.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"!mkdir -p /root/.kaggle/\n!cp ../input/mykaggleapi/kaggle.json /root/.kaggle/\n!chmod 600 /root/.kaggle/kaggle.json","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DS_DIR = '/tmp/dataset'\n!mkdir -p {DS_DIR}\n!tar czf {DS_DIR}/osic.tar.gz -C {OUT_DIR} .\ndf.to_csv(os.path.join(DS_DIR, f'CT.csv'), index=False)\n!kaggle datasets init -p {DS_DIR}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import json\n\nwith open(f'{DS_DIR}/dataset-metadata.json', 'r+') as f:\n    data = json.load(f)\n    data['title'] = f'osic sandbox'\n    data['id'] = f'dlarionov/osic-sandbox'\n    f.seek(0)\n    json.dump(data, f, indent=4)\n    f.truncate()\n\n!cat {DS_DIR}/dataset-metadata.json","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!kaggle datasets create -p {DS_DIR} -q -r tar\n!rm -rf {OUT_DIR}\n!rm -rf {DS_DIR}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}