{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Features\n\n1. Total Volume of lung - number of pixels except Air\n2. Volume of various HU units range - from -900 to 100, every 100\n3. Height of Lung - first and last segment with lung tissues\n4. Mean, Skew, Kurtosis on -900 to -320 value"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"%matplotlib inline\n\nimport os\nimport cv2\n\nimport pydicom\nimport pandas as pd\nimport numpy as np \nimport torch\nimport scipy\nimport tensorflow as tf \nimport matplotlib.pyplot as plt \nfrom pathlib import Path\nimport scipy.ndimage as ndimage\nfrom skimage import measure, morphology, segmentation\nfrom scipy.ndimage.interpolation import zoom\nfrom PIL import Image \nimport time\n\nfrom tqdm.notebook import tqdm\n%matplotlib inline\n\nfrom skimage import measure, morphology\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_dir = Path('../input/osic-pulmonary-fibrosis-progression')\npatient_paths = list((data_dir/'train').glob('*'))\nsample_patient = Path('../input/osic-pulmonary-fibrosis-progression/train/ID00233637202260580149633')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndef generate_markers(image):\n    \"\"\"\n    Generates markers for a given image.\n    \n    Parameters: image\n    \n    Returns: Internal Marker, External Marker, Watershed Marker\n    \"\"\"\n    \n    #Creation of the internal Marker\n    marker_internal = image < -400\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(image.shape, dtype=np.int)\n    marker_watershed += marker_internal * 255\n    marker_watershed += marker_external * 128\n    \n    return marker_internal, marker_external, marker_watershed\n\n\ndef seperate_lungs(image,iterations = 1):\n    \"\"\"\n    Segments lungs using various techniques.\n    \n    Parameters: image (Scan image), iterations (more iterations, more accurate mask)\n    \n    Returns: \n        - Segmented Lung\n        - Lung Filter\n        - Outline Lung\n        - Watershed Lung\n        - Sobel Gradient\n    \"\"\"\n    \n    # Store the start time\n    # start = time.time()\n    marker_internal, marker_external, marker_watershed = generate_markers(image)\n    \n    \n    '''\n    Creation of Sobel Gradient\n    '''\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    \n    '''\n    Using the watershed algorithm\n    \n    \n    We pass the image convoluted by sobel operator and the watershed marker\n    to morphology.watershed and get a matrix matrix labeled using the \n    watershed segmentation algorithm.\n    '''\n    watershed = morphology.watershed(sobel_gradient, marker_watershed)\n    \n    '''\n    Reducing the image to outlines after Watershed algorithm\n    '''\n    outline = ndimage.morphological_gradient(watershed, size=(3,3))\n    outline = outline.astype(bool)\n    \n    \n    '''\n    Black Top-hat Morphology:\n    \n    The black top hat of an image is defined as its morphological closing\n    minus the original image. This operation returns the dark spots of the\n    image that are smaller than the structuring element. Note that dark \n    spots in the original image are bright spots after the black top hat.\n    '''\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    '''\n    Generate lung filter using internal marker and outline.\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    '''\n    Segment lung using lungfilter and the image.\n    '''\n    segmented = np.where(lungfilter, image, -1000)\n    \n    #return segmented, lungfilter, outline, watershed, sobel_gradient\n    return segmented","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def fix_pxrepr(dcm):\n#     if dcm.PixelRepresentation != 0 or dcm.RescaleIntercept<-100:\n#         return dcm\n    x = dcm.pixel_array + 1000\n    px_mode = 4096\n    x[x>=px_mode] = x[x>=px_mode] - px_mode\n    dcm.PixelData = x.tobytes()\n    #dcm.RescaleIntercept = -1000\n    return dcm\n\ndef load_scan(path):\n    paths = os.listdir(path)\n    paths = sorted(paths, key=lambda x: int(str(x).split('/')[-1].split('.')[0]))\n    slices = [pydicom.read_file(path + '/' + s) for s in paths]\n    #slices = list(map(fix_pxrepr, slices))\n    try:\n        slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n    except:\n        pass\n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        try:\n            slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n        except:\n            slice_thickness = slices[0].SliceThickness\n    for s in slices:\n        s.SliceThickness = slice_thickness\n    return slices","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_pixels_hu(slices):\n    image = np.stack([s.pixel_array for s in slices])\n    # Convert to int16 (from sometimes int16), \n    # should be possible as values should always be low enough (<32k)\n    image = image.astype(np.int16)\n\n    # Set outside-of-scan pixels to 0\n    # The intercept is usually -1024, so air is approximately 0\n    \n    \n    # Convert to Hounsfield units (HU)\n    for slice_number in range(len(slices)):\n        \n        intercept = slices[slice_number].RescaleIntercept\n        slope = slices[slice_number].RescaleSlope\n        \n        if slope != 1:\n            image[slice_number] = slope * image[slice_number].astype(np.float64)\n            image[slice_number] = image[slice_number].astype(np.int16)\n            \n        image[slice_number] += np.int16(intercept)\n    image[image <= -1900] = -1000\n    return np.array(image, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def resample(image, scan, new_spacing=[1,1,1]):\n    st = time.time()\n    slice_thickness = scan[0].SliceThickness\n    spacing = np.array([slice_thickness] + list(scan[0].PixelSpacing), dtype=np.float32)\n    resize_factor = spacing / new_spacing\n    new_real_shape = image.shape * resize_factor\n    new_shape = np.round(new_real_shape)\n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    print('resample factor', time.time() - st, real_resize_factor)\n    image = scipy.ndimage.zoom(image, real_resize_factor, mode='nearest')\n    print('resample time', time.time() - st)\n    return image, new_spacing\n\ndef torch_resample(image, scan, new_spacing=[1,1,1]):\n    st = time.time()\n    slice_thickness = scan[0].SliceThickness\n    spacing = np.array([slice_thickness] + list(scan[0].PixelSpacing), dtype=np.float32)\n    resize_factor = spacing / new_spacing\n    new_real_shape = image.shape * resize_factor\n    new_shape = np.round(new_real_shape)\n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    image = torch.tensor(image).unsqueeze(0).unsqueeze(0)\n    image = torch.nn.functional.interpolate(image, scale_factor=tuple(real_resize_factor), mode='nearest')\n    image = image.squeeze(0).squeeze(0).numpy()\n    return image, new_spacing\n\n\n\n\ndef get_3d_resampled_array(patient_path):\n    start_time = time.time()\n    patient_slices = load_scan(str(patient_path))\n    patient_slices_hu = get_pixels_hu(patient_slices)\n    print('HU loaded', time.time() - start_time)\n    print(patient_slices_hu.shape)\n    #lungmask_3d = np.apply_over_axes(seperate_lungs, patient_slices_hu, 0)\n    idx = np.ndindex(patient_slices_hu.shape[0])\n    patient_slices_hu_masked = np.zeros(patient_slices_hu.shape)\n    for i in idx:\n        patient_slices_hu_masked[i] = seperate_lungs(patient_slices_hu[i])\n        #patient_slices_hu_masked[i, :, :] = np.where(lungmask, patient_slices_hu[i, :, :], -1000)\n    #patient_slices_hu_masked = np.where(lungmask_3d, patient_slices_hu, -1000)\n    \n    print('mask generated', time.time() - start_time)\n    resampled_array, spacing = torch_resample(patient_slices_hu_masked, patient_slices, [1,1,1])\n    print('after resample', time.time() - start_time)\n    return resampled_array, spacing\n\n\ndef get_features_from_3d_array(resampled_array, spacing):\n    features = {}\n    \n    # volume of lungs\n    cube_volume = spacing[0] * spacing[1] * spacing[2]\n    total_lung_volume = (resampled_array[resampled_array > -900].shape[0] * cube_volume)\n    lung_volume_in_liters = total_lung_volume / (1000*1000)\n    features['lung_volume_in_liters'] = lung_volume_in_liters\n    \n    #HU unit binning\n    bins_threshold = (resampled_array <= 300) & (resampled_array >= -900)\n    total_hu_units_bin = resampled_array[bins_threshold].flatten().shape[0]\n    bin_values, bins = np.histogram(resampled_array[bins_threshold].flatten(), bins=range(-900, 400, 100))\n    features['total_hu_units_bin'] = total_hu_units_bin\n    for i, _bin in enumerate(bins[:-1]):\n        features[f'bin_{_bin}'] = bin_values[i] / total_hu_units_bin\n    \n    #mean, skew, kurtosis\n    lung_threshold = (resampled_array <= -320) & (resampled_array >= -900)\n    histogram_values, _ = np.histogram(resampled_array[lung_threshold].flatten(), bins=100)\n    features['lung_mean_hu'] = np.mean(resampled_array[lung_threshold].flatten())\n    features['lung_skew'] = skew(histogram_values)\n    features['lung_kurtosis'] = kurtosis(histogram_values)\n    \n    #height_of_lung\n    n_lung_pixels = lung_threshold.sum(axis=1).sum(axis=1)\n    height_start = np.argwhere(n_lung_pixels > 1000).min()\n    height_end = np.argwhere(n_lung_pixels > 1000).max()\n    features['height_of_lung_cm'] = (height_end - height_start)/10\n    \n    return features","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"start_time = time.time()\nfirst_patient = load_scan(str(sample_patient))\nfirst_patient_pixels = get_pixels_hu(first_patient)\nprint('number of slices', len(first_patient))\n\nfirst_patient_masked_pixels = np.zeros(first_patient_pixels.shape)\nlungmasks = np.zeros(first_patient_pixels.shape)\n\nfor i in range(first_patient_pixels.shape[0]):\n    semented_lung = seperate_lungs(first_patient_pixels[i, :, :])\n    first_patient_masked_pixels[i, :, :] = semented_lung\n\n    \nprint('time taken', time.time() - start_time)\nplt.hist(first_patient_pixels.flatten(), bins=80, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"slice_n = 220\n# Show some slice in the middle\nplt.imshow(first_patient_pixels[slice_n], cmap=plt.cm.gray)\nplt.show()\n\nplt.imshow(first_patient_masked_pixels[slice_n], cmap=plt.cm.gray)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"_ = plt.hist(first_patient_masked_pixels[slice_n].flatten(), bins=80)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"slice_n = 220\nthreshold_filter = (first_patient_masked_pixels[slice_n] >= -900)# & (first_patient_masked_pixels[slice_n] <= -320)\nplt.imshow(np.where(threshold_filter, first_patient_masked_pixels[slice_n], -1000), cmap=plt.cm.gray)\nplt.show()\n\nplt.imshow(first_patient_masked_pixels[slice_n], cmap=plt.cm.gray)\nplt.show()\n\nplt.imshow(first_patient_pixels[slice_n], cmap=plt.cm.gray)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Resampling"},{"metadata":{"trusted":true},"cell_type":"code","source":"start_time = time.time()\npix_resampled, spacing = resample(first_patient_masked_pixels, first_patient, [1,1,1])\nprint('time taken', time.time() - start_time)\nprint(\"Shape before resampling\\t\", first_patient_pixels.shape)\nprint(\"Shape after resampling\\t\", pix_resampled.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"start_time = time.time()\npix_resampled, spacing = torch_resample(first_patient_masked_pixels, first_patient, [1,1,1])\nprint('time taken', time.time() - start_time)\nprint(\"Shape before resampling\\t\", first_patient_pixels.shape)\nprint(\"Shape after resampling\\t\", pix_resampled.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"spacing","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature Engineering"},{"metadata":{"trusted":true},"cell_type":"code","source":"resampled_array, spacing = get_3d_resampled_array(str(patient_paths[6]))\nprint(resampled_array.shape)\nprint(spacing)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cube_volume = spacing[0] * spacing[1] * spacing[2]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"total_lung_volume = (resampled_array[resampled_array > -900].shape[0] * cube_volume)\nlung_volume_in_liters = total_lung_volume / (1000*1000)\nprint(total_lung_volume, lung_volume_in_liters)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"bins_threshold = (resampled_array <= 300) & (resampled_array >= -900)\nbin_values, bins = np.histogram(resampled_array[bins_threshold].flatten(), bins=range(-900, 400, 100))\nprint(bin_values)\nprint(bins)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"list(bins)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"bin_values / sum(bin_values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lung_threshold = (resampled_array <= -320) & (resampled_array >= -900)\nhistogram_values, _ = np.histogram(resampled_array[lung_threshold].flatten(), bins=100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"_ = plt.hist(resampled_array[lung_threshold].flatten(), bins=100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.mean(resampled_array[lung_threshold].flatten())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.stats import kurtosis, skew\n\nprint(skew(histogram_values))\nprint(kurtosis(histogram_values))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(lung_threshold.sum(axis=1).sum(axis=1))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"height_start = np.argwhere(lung_threshold.sum(axis=1).sum(axis=1) > 1000).min()\nheight_end = np.argwhere(lung_threshold.sum(axis=1).sum(axis=1) > 1000).max()\nprint(height_start, height_end)\n\nheight = height_end - height_start\n\nprint(height)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.stats import kurtosis, skew\n\ndef get_features_from_3d_array(resampled_array, spacing):\n    features = {}\n    \n    # volume of lungs\n    cube_volume = spacing[0] * spacing[1] * spacing[2]\n    total_lung_volume = (resampled_array[resampled_array > -900].shape[0] * cube_volume)\n    lung_volume_in_liters = total_lung_volume / (1000*1000)\n    features['lung_volume_in_liters'] = lung_volume_in_liters\n    \n    #HU unit binning\n    bins_threshold = (resampled_array <= 300) & (resampled_array >= -900)\n    total_hu_units_bin = resampled_array[bins_threshold].flatten().shape[0]\n    bin_values, bins = np.histogram(resampled_array[bins_threshold].flatten(), bins=range(-900, 400, 100))\n    features['total_hu_units_bin'] = total_hu_units_bin\n    for i, _bin in enumerate(bins[:-1]):\n        features[f'bin_{_bin}'] = bin_values[i] / total_hu_units_bin\n    \n    #mean, skew, kurtosis\n    lung_threshold = (resampled_array <= -320) & (resampled_array >= -900)\n    histogram_values, _ = np.histogram(resampled_array[lung_threshold].flatten(), bins=100)\n    features['lung_mean_hu'] = np.mean(resampled_array[lung_threshold].flatten())\n    features['lung_skew'] = skew(histogram_values)\n    features['lung_kurtosis'] = kurtosis(histogram_values)\n    \n    #height_of_lung\n    n_lung_pixels = lung_threshold.sum(axis=1).sum(axis=1)\n    height_start = np.argwhere(n_lung_pixels > 1000).min()\n    height_end = np.argwhere(n_lung_pixels > 1000).max()\n    features['height_of_lung_cm'] = (height_end - height_start)/10\n    \n    return features","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.DataFrame(pd.Series(get_features_from_3d_array(resampled_array, spacing)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Generator"},{"metadata":{"trusted":true},"cell_type":"code","source":"import traceback\nimport warnings\n\nwarnings.filterwarnings('ignore')\n\npatients_feature_df = pd.DataFrame()\nst = time.time()\n\n# Remove filter to run for all patients\nfor patient_path in patient_paths[:5]:\n    try:\n        resampled_array, spacing = get_3d_resampled_array(str(patient_path))\n        features = get_features_from_3d_array(resampled_array, spacing)\n        features['patient_id'] = str(patient_path).split('/')[-1]\n        features['missing'] = 0\n    except Exception as e:\n        features = {}\n        features['missing'] = 1\n        print(e)\n    patient_df = pd.DataFrame(pd.Series(features)).T\n    patients_feature_df = pd.concat([patients_feature_df, patient_df], ignore_index=True)\n    patients_feature_df.to_csv('patient_feature2_df.csv', index=False)\n        \nprint('Total Time', time.time() - st)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patients_feature_df","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}