{"cells":[{"metadata":{},"cell_type":"markdown","source":"Starting point is from the excellent notebook \"OSIC - Comprehensive EDA\", but only to get things set up. \n\nHere I asked myself the question: is the image data of any use at all? I have avoided the \"throw it at a huge deep learning network and see what happens\" reflex. Instead I took patches from a single slice of the CT scan and tried to identify small features that matter. Rather than attempt to improve prediction overall, I instead look at an intermediate objective. Can we predict the final volume accurately? I look at the mean of the last three measurements. \n\nThe patch features are clustered using k-means. It seems that around 30 clusters is optimal. Then I simply count the number of occurrences of cluster features in each image. Yes, this is very, very rough. But if one of the common features is important to the final outcome then it should show up here. \n\nThe results? Not awful, but not good. I will go on from here to feed these predictions into a table based notebook to see if I can improve my overall score. \n","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt # Plotting\nimport seaborn as sns # statistical data visualization\nimport plotly.express as px\nimport plotly.graph_objects as go # interactive plots\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\n#for dirname, _, filenames in os.walk('/kaggle/input'):\n#    for filename in filenames:\n#        print(os.path.join(dirname, filename))\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train_df = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This dataset will contain a single entry for each patient, leaving only their personal information, to have a clearer and more accurate representation of the distribution of these categorical variables.\nIt will also include how many times a patient appears in the dataset, to observe how many visits each one has done.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"unique_patient_df = train_df.drop(['Weeks', 'FVC', 'Percent'], axis=1).drop_duplicates().reset_index(drop=True)\nunique_patient_df['# visits'] = [train_df['Patient'].value_counts().loc[pid] for pid in unique_patient_df['Patient']]\n\nprint('Number of data points: ' + str(len(unique_patient_df)))\nprint('----------------------')\n\nfor col in unique_patient_df.columns:\n    print('{} : {} unique values, {} missing.'.format(col, \n                                                          str(len(unique_patient_df[col].unique())), \n                                                          str(unique_patient_df[col].isna().sum())))\nunique_patient_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Expected FVC**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df['Expected FVC'] = train_df['FVC'] + (train_df['Percent']/100)*train_df['FVC']\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Add final lung volume as percentage of initial lung volume to unique_df**\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import statistics\n\nvolume_path = []\nfor patient in unique_patient_df['Patient']:\n    visits = train_df.loc[train_df['Patient'] == patient]   \n    volumes = visits['FVC'].to_numpy()\n    descent = 100*(statistics.mean(volumes[-3:])/volumes[0])\n    volume_path.append([patient,descent])\n    \nsc = pd.DataFrame(volume_path,columns=['Patient','Final_FVC'])\n    \nunique_patient_df['Final_FVC'] = sc['Final_FVC']\n\nunique_patient_df.tail()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3 - CT Scan Features ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import pydicom\nfrom glob import glob\nimport scipy.ndimage\nfrom skimage import morphology\nfrom skimage import measure\nfrom skimage.filters import threshold_otsu, median\nfrom scipy.ndimage import binary_fill_holes\nfrom skimage.segmentation import clear_border\nfrom scipy.stats import describe","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_single_slice(path):\n    slice = pydicom.dcmread(path) \n    return slice","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Look for image features.\n\nMakes use of: \n\nhttps://scikit-learn.org/stable/modules/generated/sklearn.feature_extraction.image.extract_patches_2d.html#sklearn.feature_extraction.image.extract_patches_2d","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"\n\nslice_no = 20\nimages = []\n\n\nfor patient_id in unique_patient_df['Patient']:\n    patient_dir = '../input/osic-pulmonary-fibrosis-progression/train/' + patient_id + '/'\n    \n    fname = str(slice_no) + \".dcm\"\n    path = patient_dir + fname\n    if (os.path.exists(path)):\n        patient_img = load_single_slice(path)\n        images.append(patient_img)\n        \n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Clustering","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"\nimport time\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\n\nfrom sklearn import datasets\nfrom sklearn.cluster import MiniBatchKMeans\nfrom sklearn.feature_extraction.image import extract_patches_2d\n\nfrom collections import Counter\n\n# #############################################################################\n# Learn the dictionary of images\n\nno_clusters = 10\n\nprint('Learning the dictionary... ')\nrng = np.random.RandomState(0)\nkmeans = MiniBatchKMeans(n_clusters=no_clusters, random_state=rng, verbose=False)\npatch_size = (10,10)\n\nbuffer = []\nt0 = time.time()\nno_patches = 200\n\n# cycle over the whole dataset 3 times\n\nfor _ in range(3):\n    index = 0\n    for img in images:\n        #\n        # abandoning images that fail extraction\n        #\n        try:\n            imgf = (img.pixel_array).astype(np.float64)\n            data = extract_patches_2d(imgf, patch_size, max_patches=no_patches,random_state=rng)\n            data = np.reshape(data, (len(data), -1))\n            buffer.append(data)\n            index += 1\n            if (index % 10 == 0):\n                data = np.concatenate(buffer, axis=0)\n                data -= np.mean(data, axis=0)\n                data /= np.std(data, axis=0)\n                kmeans.partial_fit(data)\n                buffer = []\n                #print(\"Partial fit of %4i out of %i\" % (index, 6 * len(images))) \n                \n        except Exception as e: print(e)\n           \n\ndt = time.time() - t0\nprint('done in %.2fs.' % dt)\n\nlabels = kmeans.labels_\ncluster_count = Counter(kmeans.labels_)\n\nprint(cluster_count)\n\n# #############################################################################\n# Plot the results\n\n\nplt.figure(figsize=(4.2, 4))\nfor i, patch in enumerate(kmeans.cluster_centers_):\n    plt.subplot(9, 9, i + 1)\n    plt.imshow(patch.reshape(patch_size), cmap=plt.cm.gray, interpolation='nearest')\n    plt.xticks(())\n    plt.yticks(())\n\n\nplt.suptitle('Patches of CT slice\\nTrain time %.1fs on %d patches' % (dt, 8 * len(images)), fontsize=16)\nplt.subplots_adjust(0.08, 0.02, 0.92, 0.85, 0.08, 0.23)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Radically compress the image to a list of valid clusters that appear in the image**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_images = []\npatient_clusters = []\n\nfor patient_id in unique_patient_df['Patient']:\n    patient_dir = '../input/osic-pulmonary-fibrosis-progression/train/' + patient_id + '/'\n    \n    fname = str(slice_no) + \".dcm\"\n    path = patient_dir + fname\n    if (os.path.exists(path)):\n        patient_img = load_single_slice(path)\n        patient_images.append([patient_id,patient_img])\n           \n            \n            \nfor (patient_id,img) in patient_images:\n    cluster_list = []\n    row=[patient_id]\n    #\n    try:\n        imgf = (img.pixel_array).astype(np.float64)\n        data = extract_patches_2d(imgf, patch_size, max_patches=300,random_state=rng)\n        data = np.reshape(data, (len(data), -1))\n        closest = list(kmeans.predict(data)) \n                \n    except Exception as e: print(e) \n     \n    cls = [0]*no_clusters\n    for i in range(no_clusters):\n        cls[i] = closest.count(i)       \n    row = row + cls    \n    patient_clusters.append(row)\n    \nhd = ['Patient'] + [str(x) for x in range(no_clusters)]\n\ncluster_df = pd.DataFrame(patient_clusters,columns=hd)\n\ncluster_df.head()\n  \n           ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cluster_df = pd.merge(cluster_df,unique_patient_df,on='Patient')\n\ncluster_df = cluster_df.drop(['Sex','SmokingStatus','# visits'],axis=1)\n\ncluster_df.head()\n\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = cluster_df.drop(['Patient','Final_FVC'],axis=1)\ny = cluster_df['Final_FVC']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Random forest predictor of final lung volume ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)\n\nfrom sklearn.ensemble import RandomForestRegressor\nregressor = RandomForestRegressor(n_estimators = 1000, random_state = 42)\nregressor.fit(X_train, y_train)\n\ny_pred = regressor.predict(X_test)\n\ndf=pd.DataFrame({'Actual':y_test, 'Predicted':y_pred})\n\n\nfrom sklearn import metrics\nprint('Mean Absolute Error:', metrics.mean_absolute_error(y_test, y_pred))\nprint('Mean Squared Error:', metrics.mean_squared_error(y_test, y_pred))\nprint('Root Mean Squared Error:', np.sqrt(metrics.mean_squared_error(y_test, y_pred)))\n\nprint(df)","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}