{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport pandas as pd\nfrom pydicom import dcmread\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport plotly.express as px\nimport scipy\nfrom scipy import ndimage\n\nfrom skimage import measure\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.linear_model import LinearRegression\n\nimport numpy as np\n\nfrom glob import glob\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport scipy.ndimage\nfrom skimage import morphology\nfrom skimage import measure\nfrom skimage.transform import resize\nfrom sklearn.cluster import KMeans\nfrom plotly import __version__\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nfrom plotly.tools import FigureFactory as FF\nfrom plotly.graph_objs import *\n\nfrom scipy.ndimage import zoom\nfrom scipy.stats import kurtosis\nfrom scipy.stats import skew\n\nfrom tqdm import tqdm\n\nimport tensorflow as tf\nimport tensorflow_hub as hub\nfrom tensorflow.keras.layers.experimental import preprocessing\n\nimport keras\n\nos.getcwd()\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_path='../input/osic-pulmonary-fibrosis-progression'\ntraindf=pd.read_csv(f'{data_path}/train.csv')\n\ntraindf","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"coredf=pd.DataFrame({'pid': traindf['Patient'].unique()})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_weeks=[]\nweeks=[]\nfvc=[]\npercent=[]\nage=[]\nsex=[]\nsmoking=[]\n\n\nfor pid in coredf['pid']:\n    tempdf=traindf[traindf['Patient']==pid]\n    \n    n_weeks.append(len(tempdf))\n    weeks.append(list(tempdf['Weeks'].values))\n    fvc.append(list(tempdf['FVC'].values))\n    percent.append(list(tempdf['Percent'].values))\n    age.append(tempdf['Age'].values[0])\n    sex.append(tempdf['Sex'].values[0])\n    smoking.append(tempdf['SmokingStatus'].values[0])\n    \n\ncoredf['n_weeks']=n_weeks\ncoredf['weeks']=weeks\ncoredf['fvc']=fvc\ncoredf['percent']=percent\ncoredf['age']=age\ncoredf['sex']=sex\ncoredf['smoking']=smoking\n\ndisplay(coredf)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dicoms=[]\nn_dicoms=[]\nbad_ids=[]\n\nfor pid in coredf['pid']:\n    try:\n        files=os.listdir(f\"{data_path}/train/{pid}\")\n\n        for i in range(len(files)):\n            files[i]=files[i][:-4]\n\n        dicoms.append(files)\n        n_dicoms.append(len(files))\n    \n    except:\n        bad_ids.append(pid)\n\ncoredf['dicoms']=dicoms\ncoredf['n_dicoms']=n_dicoms\n\ndisplay(coredf)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"coredf.query('pid == \"ID00052637202186188008618\"')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_scans(path):\n    dicoms=[]\n\n    for filename in os.listdir(path):\n        dicoms.append(dcmread(os.path.join(path,filename)))\n\n    return dicoms\n\ndef get_images(scans):\n  temp=[]\n\n  for i in range(len(scans)):\n    try:\n        temp.append(scans[i].pixel_array)\n    except:\n        return None\n  \n  return np.asarray(temp,dtype=np.int16)\n\ndef crop_image(image):\n    if image.shape != (512,512):\n        left=((image.shape[0]-512)//2)\n        right=image.shape[0]-left\n        top=((image.shape[1]-512)//2)\n        bottom=image.shape[1]-top\n        \n        return image[left-1:right-1,top-1:bottom-1]\n    else:\n        return image\n    \ndef crop_images(images):\n    if images[0].shape != (512,512):\n        temp=[]\n        for image in images:\n            temp.append(crop_image(image))\n        \n        return np.asarray(temp)\n    else:\n        return images\n\ndef resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    spacing=[scan[0].SliceThickness,scan[0].PixelSpacing[0],scan[0].PixelSpacing[1]]\n    spacing = np.array(list(spacing))\n\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)\n    \n    return image\n\ndef normalize(images):\n    '''\n      Returns normalized images\n\n      Inputs:\n      imageslices:  NumpyArray  Imageslices extracted from the DICOM files using fnc get_images_from_slices\n\n    '''\n    images=(images-images.min())/(images.max()-images.min())\n    \n    return images\n\ndef dimension_scaling(images):\n  '''\n  Returns images with channels. e.g. converts 512x512 images to 512x512x3\n\n  Inputs:\n  imageslices:  NumpyArray  Imageslices extracted from the DICOM files using fnc get_images_from_slices\n\n  '''\n\n  images=np.expand_dims(images,axis=-1)\n  images=np.concatenate((images,)*3, axis=-1)\n\n  return np.asarray(images)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_mesh(image, threshold=-300, step_size=1):\n    p = image.transpose(2,1,0)\n\n    verts, faces, norm, val = measure.marching_cubes(p, threshold, step_size=step_size, allow_degenerate=True) \n    return verts, faces\n\ndef plotly_3d(verts, faces):\n    x,y,z = zip(*verts) \n    \n    # Make the colormap single color since the axes are positional not intensity. \n#    colormap=['rgb(255,105,180)','rgb(255,255,51)','rgb(0,191,255)']\n    colormap=['rgb(236, 236, 212)','rgb(236, 236, 212)']\n    \n    fig = FF.create_trisurf(x=x,\n                        y=y, \n                        z=z, \n                        plot_edges=False,\n                        colormap=colormap,\n                        simplices=faces,\n                        backgroundcolor='rgb(64, 64, 64)',\n                        title=\"Interactive Visualization\")\n    iplot(fig)\n\ndef plt_3d(verts, faces):\n    x,y,z = zip(*verts) \n    fig = plt.figure(figsize=(10, 10))\n    ax = fig.add_subplot(111, projection='3d')\n\n    # Fancy indexing: `verts[faces]` to generate a collection of triangles\n    mesh = Poly3DCollection(verts[faces], linewidths=0.05, alpha=1)\n    face_color = [1, 1, 0.9]\n    mesh.set_facecolor(face_color)\n    ax.add_collection3d(mesh)\n\n    ax.set_xlim(0, max(x))\n    ax.set_ylim(0, max(y))\n    ax.set_zlim(0, max(z))\n    ax.set_axis_bgcolor((0.7, 0.7, 0.7))\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Standardize the pixel values\ndef make_lungmask(img, display=False):\n    row_size= img.shape[0]\n    col_size = img.shape[1]\n    \n    mean = np.mean(img)\n    std = np.std(img)\n    img = img-mean\n    img = img/std\n    # Find the average pixel value near the lungs\n    # to renormalize washed out images\n    middle = img[int(col_size/5):int(col_size/5*4),int(row_size/5):int(row_size/5*4)] \n    mean = np.mean(middle)  \n    max = np.max(img)\n    min = np.min(img)\n    # To improve threshold finding, I'm moving the \n    # underflow and overflow on the pixel spectrum\n    img[img==max]=mean\n    img[img==min]=mean\n    #\n    # Using Kmeans to separate foreground (soft tissue / bone) and background (lung/air)\n    #\n    kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))\n    centers = sorted(kmeans.cluster_centers_.flatten())\n    threshold = np.mean(centers)\n    thresh_img = np.where(img<threshold,1.0,0.0)  # threshold the image\n\n    # First erode away the finer elements, then dilate to include some of the pixels surrounding the lung.  \n    # We don't want to accidentally clip the lung.\n\n    eroded = morphology.erosion(thresh_img,np.ones([3,3]))\n    dilation = morphology.dilation(eroded,np.ones([8,8]))\n\n    labels = measure.label(dilation) # Different labels are displayed in different colors\n    label_vals = np.unique(labels)\n    regions = measure.regionprops(labels)\n    good_labels = []\n    for prop in regions:\n        B = prop.bbox\n        if B[2]-B[0]<row_size/10*9 and B[3]-B[1]<col_size/10*9 and B[0]>row_size/5 and B[2]<col_size/5*4:\n            good_labels.append(prop.label)\n    mask = np.ndarray([row_size,col_size],dtype=np.int8)\n    mask[:] = 0\n\n    #\n    #  After just the lungs are left, we do another large dilation\n    #  in order to fill in and out the lung mask \n    #\n    for N in good_labels:\n        mask = mask + np.where(labels==N,1,0)\n    mask = morphology.dilation(mask,np.ones([10,10])) # one last dilation\n\n    if (display):\n        fig, ax = plt.subplots(3, 2, figsize=[12, 12])\n        ax[0, 0].set_title(\"Original\")\n        ax[0, 0].imshow(img, cmap='gray')\n        ax[0, 0].axis('off')\n        ax[0, 1].set_title(\"Threshold\")\n        ax[0, 1].imshow(thresh_img, cmap='gray')\n        ax[0, 1].axis('off')\n        ax[1, 0].set_title(\"After Erosion and Dilation\")\n        ax[1, 0].imshow(dilation, cmap='gray')\n        ax[1, 0].axis('off')\n        ax[1, 1].set_title(\"Color Labels\")\n        ax[1, 1].imshow(labels)\n        ax[1, 1].axis('off')\n        ax[2, 0].set_title(\"Final Mask\")\n        ax[2, 0].imshow(mask, cmap='gray')\n        ax[2, 0].axis('off')\n        ax[2, 1].set_title(\"Apply Mask on Original\")\n        ax[2, 1].imshow(mask*img, cmap='gray')\n        ax[2, 1].axis('off')\n        \n        plt.show()\n        \n    return mask","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"make_lungmask(get_images(get_scans(f'{data_path}/train/ID00007637202177411956430'))[0],display=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"coredf","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"m = tf.keras.Sequential([\n    hub.KerasLayer(\"https://tfhub.dev/tensorflow/efficientnet/b6/feature-vector/1\",\n                   trainable=False),  # Can be True, see below.\n    tf.keras.layers.Dense(7, activation='relu')\n])\nm.build([None, 512, 512, 3])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"eff_features=[]\n\n\nfor i in tqdm(range(1,len(coredf))):\n    pid=coredf['pid'][i]\n    \n    images=get_images(get_scans(f\"{data_path}/train/{pid}\"))\n    \n    if type(images)!=type(None):\n        images=crop_images(images)\n        images=normalize(images)\n\n        for i in tqdm(range(len(images)), leave=True, position=0,desc=f\"{pid}: {i}\"):\n            try:\n                images[i]=images[i]*make_lungmask(images[i],False)\n            except:\n                images[i]=images[i]*np.zeros((512,512))\n        \n        \n\n        \n        images=dimension_scaling(images)\n        \n        features=np.mean(m.predict(images),axis=0)\n        \n        eff_features.append(features)\n \n    else:\n        eff_features.append(np.zeros(10,np.float32))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"age_list=[]\nsex_list=[]\nsmoking_list=[]\nini_week=[]\nini_fvc=[]\nnew_week=[]\nweeks_passed=[]\nfinal_week=[]\nfvc_passed=[]\n\nfor i in tqdm(range(len(coredf))):\n    weeks=coredf['weeks'][i]\n    fvc=coredf['fvc'][i]\n    \n    age=coredf['age'][i]\n    \n    if coredf['sex'][i]=='Male':\n        sex=0\n    else:\n        sex=1\n\n    if coredf['smoking'][i] == 'Ex-smoker':\n        smoking=0\n    elif coredf['smoking'][i] == 'Never smoked':\n        smoking=1\n    else:\n        smoking=2\n    \n    for i in range(len(weeks)):            \n        for j in range(len(weeks)):\n            age_list.append(age)\n            sex_list.append(sex)\n            smoking_list.append(smoking)\n        \n            ini_week.append(weeks[i])\n            ini_fvc.append(fvc[i])\n            final_week.append(weeks[j])\n            \n            weeks_passed.append(weeks[j]-weeks[i])\n            fvc_passed.append(fvc[i]-fvc[j])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"modeldf=pd.DataFrame({'age':age_list,\n                      'sex':sex_list,\n                      'smoking':smoking_list,\n                       'ini_week':ini_week,\n                      'ini_fvc':ini_fvc,\n                      'weeks_passed':weeks_passed,\n                      'final_week':final_week,\n                      'fvc_change':fvc_passed\n                     })","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"modeldf","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def score(y_true, y_pred):\n    tf.dtypes.cast(y_true, tf.float32)\n    tf.dtypes.cast(y_pred, tf.float32)\n    sigma = y_pred[:, 2] - y_pred[:, 0]\n    fvc_pred = y_pred[:, 1]\n    \n    #sigma_clip = sigma + C1\n    sigma_clip = tf.maximum(sigma, C1)\n    delta = tf.abs(y_true[:, 0] - fvc_pred)\n    delta = tf.minimum(delta, C2)\n    sq2 = tf.sqrt( tf.dtypes.cast(2, dtype=tf.float32) )\n    metric = (delta / sigma_clip)*sq2 + tf.math.log(sigma_clip* sq2)\n    \n    return K.mean(metric)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"px.imshow(modeldf.corr())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"trainx=modeldf[['age','sex','smoking','ini_week','ini_fvc','weeks_passed','final_week']]\ntrainy=modeldf['fvc_change']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"initializer = tf.keras.initializers.GlorotNormal()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ini_input=keras.Input(shape=(7,))\n\ny = keras.layers.LSTM(25)(tf.expand_dims(ini_input, axis=-1))\ny = keras.layers.LSTM(25)(y)\ny = keras.layers.LSTM(25)(y)\ny = keras.layers.Reshape([25*7])(y)\nx = keras.layers.Dense(100,kernel_initializer=initializer, activation=\"relu\", name=\"d1\")(y)\nx = keras.layers.Dense(100,kernel_initializer=initializer, activation=\"relu\", name=\"d2\")(x)\np1 = keras.layers.Dense(3, activation=\"linear\", name=\"p1\")(x)\np2 = keras.layers.Dense(3, activation=\"relu\", name=\"p2\")(x)\npreds = keras.layers.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), \n                     name=\"preds\")([p1, p2])\nmodel=keras.Model(inputs=ini_input,outputs=p2)\n\nmodel.compile(loss=\"mean_absolute_error\", optimizer=tf.keras.optimizers.Adam(lr=0.1, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.01, amsgrad=False),metrics=['accuracy'])\nmodel.fit(normalize(trainx),trainy,epochs=10,validation_split = 0.20,shuffle=True)","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}