{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Generating Features from pixel histograms","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"preprocessing guided by https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm --y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import copy\nimport cv2\nfrom skimage.segmentation import clear_border\nfrom skimage.morphology import ball, disk, dilation, binary_erosion, remove_small_objects, erosion, closing, reconstruction, binary_closing\nfrom skimage.filters import roberts, sobel\nfrom scipy import ndimage as ndi\nimport os\n#import gdcm\nfrom tqdm import tqdm\nfrom skimage import measure, morphology\nfrom datetime import datetime, timedelta\nimport matplotlib.pyplot as plt\nfrom matplotlib import cm\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\nimport pickle\nimport pydicom\nimport pydicom\nfrom scipy.stats import kurtosis\nimport seaborn as sns\nimport scipy\npydicom.config.image_handlers = ['gdcm_handler']\n#pydicom.config.image_handlers = ['pillow_handler']\nfrom sklearn.model_selection import GroupKFold, GroupShuffleSplit\nfrom torch.utils.data import Dataset\nfrom tqdm import tqdm\nimport torch\nfrom torch.utils.data import DataLoader, Subset\nfrom torch import optim\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom tqdm import trange\nfrom time import time\nimport warnings\nfrom scipy.ndimage.interpolation import zoom\nfrom enum import Enum\nfrom torchvision import transforms\nfrom skimage.measure import label, regionprops\nfrom skimage.segmentation import clear_border\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train=pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntest=pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\nsubmission=pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['base_Weeks']=train.groupby(['Patient'])['Weeks'].transform('min')\nbase=train[train.Weeks==train.base_Weeks]\nbase = base.rename(columns={'FVC': 'base_FVC','Percent': 'base_Percent'})\nbase.drop_duplicates(subset=['Patient', 'Weeks'], keep='first',inplace=True)\ntrain=train.merge(base[['Patient','base_FVC','base_Percent']],on='Patient',how='left')\ntrain['Week_passed'] = train['Weeks'] - train['base_Weeks']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = test.rename(columns={'Weeks': 'base_Weeks', 'FVC': 'base_FVC','Percent': 'base_Percent'})\n\n# Adding Sample Submission\nsubmission = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/sample_submission.csv\")\n\n# In submisison file, format: ID_'week', using lambda to split the ID\nsubmission['Patient'] = submission['Patient_Week'].apply(lambda x:x.split('_')[0])\n\n# In submisison file, format: ID_'week', using lambda to split the Week\nsubmission['Weeks'] = submission['Patient_Week'].apply(lambda x:x.split('_')[1]).astype(int)\n\ntest = submission.drop(columns = [\"FVC\", \"Confidence\"]).merge(test, on = 'Patient')\n\ntest['Week_passed'] = test['Weeks'] - test['base_Weeks']\n\ntest=test[train.columns.drop(['FVC','Percent'])]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load the scans in given folder path\ndef load_scan(path):\n\n    #slices = [pydicom.read_file(path / s) for s in os.listdir(path)]\n    slices = [pydicom.read_file(path / s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n    if slice_thickness==0:\n        slice_thickness=slices[0].SliceThickness\n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_pixels_hu(slices):\n    image = np.stack([np.array(s.pixel_array,dtype=np.int16) 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    image[image == -2000] = 0\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    \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    # Determine current pixel spacing\n    #spacing = np.array([scan[0].SliceThickness] + scan[0].PixelSpacing, dtype=np.float32)\n    spacing = np.array([scan[0].SliceThickness] + 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    \n    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')\n    return image, new_spacing","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_segmented_lungs(im, plot=False):\n    \n    '''\n    This funtion segments the lungs from the given 2D slice.\n    '''\n    if plot == True:\n        f, plots = plt.subplots(8, 1, figsize=(5, 40))\n    '''\n    Step 1: Convert into a binary image. \n    '''\n    binary = im < -200\n    if plot == True:\n        plots[0].axis('off')\n        plots[0].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 2: Remove the blobs connected to the border of the image.\n    '''\n    cleared = clear_border(binary)\n    if plot == True:\n        plots[1].axis('off')\n        plots[1].imshow(cleared, cmap=plt.cm.bone) \n    '''\n    Step 3: Label the image.\n    '''\n    label_image = label(cleared)\n    if plot == True:\n        plots[2].axis('off')\n        plots[2].imshow(label_image, cmap=plt.cm.bone) \n    '''\n    Step 4: Keep the labels with 2 largest areas.\n    '''\n    areas = [r.area for r in regionprops(label_image)]\n    areas.sort()\n    if len(areas) > 2:\n        for region in regionprops(label_image):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:                \n                       label_image[coordinates[0], coordinates[1]] = 0\n    binary = label_image > 0\n    if plot == True:\n        plots[3].axis('off')\n        plots[3].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 5: Erosion operation with a disk of radius 2. This operation is \n    seperate the lung nodules attached to the blood vessels.\n    '''\n    selem = disk(2)\n    binary = binary_erosion(binary, selem)\n    if plot == True:\n        plots[4].axis('off')\n        plots[4].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 6: Closure operation with a disk of radius 10. This operation is \n    to keep nodules attached to the lung wall.\n    '''\n    selem = disk(10)\n    binary = binary_closing(binary, selem)\n    if plot == True:\n        plots[5].axis('off')\n        plots[5].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 7: Fill in the small holes inside the binary mask of lungs.\n    '''\n    edges = roberts(binary)\n    binary = ndi.binary_fill_holes(edges)\n    if plot == True:\n        plots[6].axis('off')\n        plots[6].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 8: Superimpose the binary mask on the input image.\n    '''\n    get_high_vals = binary == 0\n    im[get_high_vals] = 0\n    if plot == True:\n        plots[7].axis('off')\n        plots[7].imshow(im, cmap=plt.cm.bone) \n        \n    return im","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import random\nroot_dir = Path('/kaggle/input/osic-pulmonary-fibrosis-progression')\nctscans_dir=root_dir/'train'\ncache_dir = Path('/kaggle/input/osic-cache/cache')\nlatent_dir = Path('/kaggle/working/latent')\nids=train.Patient.unique()\n#index = np.argwhere(ids=='ID00011637202177653955184')\n#ids = list(np.delete(ids, index))\n#random.shuffle(ids)\nids=np.array(ids)\ntest_ids=test.Patient.unique()\ntrain_ids,val_ids=np.split(ids, [int(round(0.9 * len(ids), 0))])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_kurtosis_stats(ids):\n    kurt=[]\n    std=[]\n    fvc=[]\n    mean=[]\n    median=[]\n    for i in ids:\n        print(i)\n        try:\n            patient_path= ctscans_dir / i\n            scan = load_scan(patient_path)\n            image=get_pixels_hu(scan)\n            image, new_spacing = resample(image, scan, new_spacing=[2,2,2])\n            image=np.asarray([get_segmented_lungs(slice) for slice in image])\n            kurt_i=kurtosis(image.ravel()[image.ravel() < -200])\n            std_i=image.ravel()[image.ravel() < -200].std()\n            fvc_i=train.base_FVC[train.Patient==i].values[0]\n            mean_i=image.ravel()[image.ravel() < -200].mean()\n            median_i=np.median(image.ravel()[image.ravel() < -200])\n            print('Kurtosis: ', kurt_i)\n            print('Standard Deviation: ', std_i)\n            print('FVC: ', fvc_i)\n            kurt.append(kurt_i)\n            std.append(std_i)\n            fvc.append(fvc_i)\n            mean.append(mean_i)\n            median.append(median_i)\n            ax=sns.kdeplot(image.ravel()[(image.ravel() < 0)&(image.ravel() > -1200)], bw=0.5)\n            ax.set(xlabel='HU', ylabel='% voxels',title='Histogram of voxel characteristics')\n            plt.show()\n            plt.imshow(image[round(image.shape[0]/2),:,:])\n            plt.show()\n        except:\n            print('error')\n            kurt.append(np.nan)\n            std.append(np.nan)\n            fvc.append(np.nan)\n            mean.append(np.nan)\n            median.append(np.nan)\n    return kurt,std,fvc,mean,median\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_ct_scan(scan):\n    f, plots = plt.subplots(int(scan.shape[0] / 20) + 1, 4, figsize=(25, 25))\n    for i in range(0, scan.shape[0], 5):\n        plots[int(i / 20), int((i % 20) / 5)].axis('off')\n        plots[int(i / 20), int((i % 20) / 5)].imshow(scan[i], cmap=plt.cm.bone)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def check():\n    patient_path= ctscans_dir / ids[8]\n    scan = load_scan(patient_path)\n    image=get_pixels_hu(scan)\n    image, new_spacing = resample(image, scan, new_spacing=[2,2,2])\n    image=np.asarray([get_segmented_lungs(slice) for slice in image])\n    #plt.imshow(image[50,:,:])\n    return image\n\nplot_ct_scan(check())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"k,s,f,m,me=get_kurtosis_stats(ids)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(k,f)\nplt.title('Scatter Plot of base FVC against pixel histogram kurtosis')\nplt.xlabel('Kurtosis')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(s,f)\nplt.title('Scatter Plot of base FVC against pixel histogram standard deviation')\nplt.xlabel('Standard Deviation')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pixel_stats=train.copy()\npixel_stats=pixel_stats.drop_duplicates(subset=['Patient'])\npixel_stats['kurtosis']=np.array(k)\npixel_stats['std']=np.array(s)\npixel_stats['mean']=np.array(m)\npixel_stats['median']=np.array(me)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train=train.merge(pixel_stats[['Patient','kurtosis','std','mean','median']],how='left',on='Patient')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.to_csv('train_pixel_stats.csv')","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}