{"cells":[{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm -y","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np \nimport pandas as pd \nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport pydicom\nimport scipy.ndimage\nimport gdcm\n\nfrom skimage import measure \nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nfrom skimage.morphology import disk, opening, closing\nfrom tqdm import tqdm\n\nfrom IPython.display import HTML\nfrom PIL import Image\n\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)\n\nfrom os import listdir, mkdir","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Table of contents\n\n1. [Qu'est-ce que la fibrose pulmonaire ?](#fibrosis)\n2. [References](#references)\n * [Data Science Bowl 2017 - Preprocessing Tutorial by Guido Zuidhof](#bowl_2017)\n * [Papers](#papers)\n2. [Data paths](#prepare)\n3. [Travailler avec des fichiers dicom](#dicom)\n * [Chargement des CT-scans par patient](#ct_scans)\n * [Transformation en unités Hounsfield](#hunits)\n * [La taille des voxels](#voxel)\n * [Surface et volume de la tranche du CT-scan - EDA](#scan_eda)\n * [Reconstruction 3D des scanners](#reconstruction)\n * [Segmentation des tissus](#segmentation)\n4. [Génération d'un ensemble de données pour les fichiers prétraités](#datagenerator)\n * [Lien vers le data-set](https://www.kaggle.com/allunia/osic-pulmonary-fibrosis-progression-huscans)"},{"metadata":{},"cell_type":"markdown","source":"# Qu'est-ce que la fibrose pulmonaire ?"},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"\n#### **Quand parle-t-on de fibrose pulmonaire ?**\n\n On parle de fibrose pulmonaire lorsque se développe dans le poumon du tissu fibreux qui remplace peu à peu le tissu normal. Certaines fibroses pulmonaires résultent d’une toxicité de médicaments ; d’autres sont associées à des maladies auto-immunes, c’est-à-dire des maladies dans lesquelles notre système immunitaire attaque nos organes.\n Parfois, on ne retrouve pas de cause et on parle de « fibrose pulmonaire idiopathique » (FPI). La FPI est une maladie rare (1 personne sur 2 500 à 1 sur 7 000 personnes), d’origine inexpliquée, qu’on rencontre plus volontiers après la soixantaine.\n\n#### **Comment se manifeste la fibrose pulmonaire ?**\n\n La fibrose pulmonaire se traduit par un essoufflement progressif et une toux sèche. Les patients peuvent aussi présenter un amaigrissement, une perte d’appétit, une fatigue importante. Dans un cas sur deux, les doigts revêtent un aspect caractéristique en baguette de tambour, avec des ongles bombés.\n\n\n#### **Comment diagnostique-t-on la fibrose pulmonaire ?**\n\n Le diagnostic de fibrose pulmonaire est souvent difficile. À l’auscultation, le médecin peut entendre des bruits pulmonaires évocateurs.\n\n La radiographie pulmonaire peut être normale au début. Mais, le scanner visualise, dans un cas sur deux, les zones de fibrose sous forme d’un aspect « en rayon de miel ».\n\n Pour écarter d’autres maladies (maladie liée à l’amiante, silicose des mineurs), on pourra effectuer un lavage des bronches et des alvéoles, ce qui nécessite une fibroscopie ; ce lavage permet de recueillir des cellules pulmonaires pour les analyser.\n\n Parfois, une biopsie de poumon sera demandée. Cet examen exige un geste chirurgical.\n Afin d’évaluer le degré de handicap respiratoire, on a recours à un test de marche (capacité à l’effort), des épreuves fonctionnelles respiratoires (mesure du souffle), une mesure des gaz sanguins (oxygène, gaz carbonique). \n"},{"metadata":{},"cell_type":"markdown","source":"# References "},{"metadata":{},"cell_type":"markdown","source":"\n\n* [Intrinsic dependencies of CT radiomic features on voxel size and number of gray levels](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5462462/)"},{"metadata":{},"cell_type":"markdown","source":"# Data paths "},{"metadata":{"trusted":true},"cell_type":"code","source":"listdir(\"../input/\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#basepath = \"../input/osic-pulmonary-fibrosis-progression/\"\n# or if you are taking part in RSNA pulmonary embolism detection:\nbasepath = \"../input/rsna-str-pulmonary-embolism-detection/\"\nlistdir(basepath)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's load the csv-files:"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(basepath + \"train.csv\")\ntest = pd.read_csv(basepath + \"test.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if basepath == \"../input/osic-pulmonary-fibrosis-progression/\":\n train[\"dcm_path\"] = basepath + \"train/\" + train.Patient + \"/\"\nelse:\n train[\"dcm_path\"] = basepath + \"train/\" + train.StudyInstanceUID + \"/\" + train.SeriesInstanceUID ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Travailler avec des fichiers dicom "},{"metadata":{},"cell_type":"markdown","source":"## Chargement des CT-scans par patient \n\n* Pour charger le scan 3D complet, nous devons commander les fichiers/tranches dicom uniques par \"ImagePositionPatient\": "},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"def load_scans(dcm_path):\n if basepath == \"../input/osic-pulmonary-fibrosis-progression/\":\n # in this competition we have missing values in ImagePosition, this is why we are sorting by filename number\n files = listdir(dcm_path)\n file_nums = [np.int(file.split(\".\")[0]) for file in files]\n sorted_file_nums = np.sort(file_nums)[::-1]\n slices = [pydicom.dcmread(dcm_path + \"/\" + str(file_num) + \".dcm\" ) for file_num in sorted_file_nums]\n else:\n # otherwise we sort by ImagePositionPatient (z-coordinate) or by SliceLocation\n slices = [pydicom.dcmread(dcm_path + \"/\" + file) for file in listdir(dcm_path)]\n slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n return slices","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"example = train.dcm_path.values[0]\nscans = load_scans(example)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Examinons le premier fichier dicom de notre exemple de patient :"},{"metadata":{"trusted":true},"cell_type":"code","source":"scans[0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### information\n\n1. Le CT-scan capture des informations sur la radiodensité d'un objet ou d'un tissu exposé aux rayons X. Une tranche transversale d'un scan est reconstituée après avoir pris des mesures dans plusieurs directions différentes.\n2. Nous devons passer à des unités Hounsfield car la composition spectrale des rayons X dépend des paramètres de mesure tels que les paramètres d'acquisition et la tension du tube. En se normalisant aux valeurs de l'eau et de l'air (l'eau a un HU 0 et l'air -1000), les images des différentes mesures deviennent comparables.\n3. Un ct-scanner donne environ 4000 valeurs de gris qui ne peuvent pas être capturées par nos yeux. C'est pourquoi on procède à un fenêtrage. De cette façon, l'image est affichée dans une plage de HU qui correspond le mieux à la région d'intérêt. "},{"metadata":{},"cell_type":"markdown","source":"## Transformation en unités Hounsfield "},{"metadata":{},"cell_type":"markdown","source":"Avant de commencer, traçons la distribution des pixels de certains fichiers dicom pour avoir une impression des données brutes:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nfor n in range(10):\n image = scans[n].pixel_array.flatten()\n rescaled_image = image * scans[n].RescaleSlope + scans[n].RescaleIntercept\n sns.distplot(image.flatten(), ax=ax[0]);\n sns.distplot(rescaled_image.flatten(), ax=ax[1])\nax[0].set_title(\"Raw pixel array distributions for 10 examples\")\nax[1].set_title(\"HU unit distributions for 10 examples\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Pour quelques exemples, nous pouvons voir qu'il existe des valeurs brutes à -2000. Elles correspondent à des images ayant une limite circulaire à l'intérieur de l'image. L'\"extérieur\" de cette valeur circulaire est souvent fixé par défaut à -2000 (ou dans d'autres concours, j'ai également trouvé -3000)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def transform_to_hu(slices):\n images = np.stack([file.pixel_array for file in slices])\n images = images.astype(np.int16)\n\n # convert ouside pixel-values to air:\n # I'm using <= -1000 to be sure that other defaults are captured as well\n images[images <= -1000] = 0\n \n # convert to HU\n for n in range(len(slices)):\n \n intercept = slices[n].RescaleIntercept\n slope = slices[n].RescaleSlope\n \n if slope != 1:\n images[n] = slope * images[n].astype(np.float64)\n images[n] = images[n].astype(np.int16)\n \n images[n] += np.int16(intercept)\n \n return np.array(images, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"hu_scans = transform_to_hu(scans)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,4,figsize=(20,3))\nax[0].set_title(\"Original CT-scan\")\nax[0].imshow(scans[0].pixel_array, cmap=\"bone\")\nax[1].set_title(\"Pixelarray distribution\");\nsns.distplot(scans[0].pixel_array.flatten(), ax=ax[1]);\n\nax[2].set_title(\"CT-scan in HU\")\nax[2].imshow(hu_scans[0], cmap=\"bone\")\nax[3].set_title(\"HU values distribution\");\nsns.distplot(hu_scans[0].flatten(), ax=ax[3]);\n\nfor m in [0,2]:\n ax[m].grid(False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Tri de nos tranches et création de GIF"},{"metadata":{"trusted":true},"cell_type":"code","source":"first_patient = load_slice('../input/rsna-str-pulmonary-embolism-detection/train/0003b3d648eb/d2b2960c2bbf')\nfirst_patient_pixels = transform_to_hu(first_patient)\n\ndef sample_stack(stack, rows=6, cols=6, start_with=10, show_every=5):\n fig,ax = plt.subplots(rows,cols,figsize=[18,20])\n for i in range(rows*cols):\n ind = start_with + i*show_every\n ax[int(i/rows),int(i % rows)].set_title(f'slice {ind}')\n ax[int(i/rows),int(i % rows)].imshow(stack[ind],cmap='bone')\n ax[int(i/rows),int(i % rows)].axis('off')\n plt.show()\n\nsample_stack(first_patient_pixels)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Maintenant, toutes les valeurs brutes par tranche sont mises à l'échelle des H-units."},{"metadata":{},"cell_type":"markdown","source":"## La taille des voxels \n\nLe voxel représente le pixel 3D qui est donné dans un scanner. Pour autant que je sache, il est couvert par le plan 2D de l'attribut d'espacement des pixels dans les directions x et y et par l'épaisseur de la tranche dans la direction z."},{"metadata":{},"cell_type":"markdown","source":"### Pixelspacing\n\n* L'attribut \"pixelspacing\" que vous pouvez trouver dans les fichiers dicom est important. Il nous indique la distance physique parcourue par un pixel. Vous pouvez voir qu'il n'y a que 2 valeurs qui décrivent les directions x et y dans le plan d'une tranche transversale. \n* Pour un patient, cet espacement des pixels est généralement le même pour toutes les tranches.\n* Mais entre les patients, l'espacement des pixels peut varier en raison des préférences personnelles ou institutionnelles des médecins et de la clinique, et il dépend également du type de scanner. Par conséquent, si vous comparez deux images dans la taille des poumons, cela ne signifie pas automatiquement que la plus grande est vraiment plus grande dans la taille physique de l'organe !\n\nExaminons les distributions des largeurs et des hauteurs d'espacement des pixels des patients"},{"metadata":{"trusted":true},"cell_type":"code","source":"N = 100","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Pour accélérer le calcul, j'ai sélectionné N patients à prendre en considération. Utilisez N = train.shape [0] pour le faire pour tous les patients de l'ensemble de données"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def get_window_value(feature):\n if type(feature) == pydicom.multival.MultiValue:\n return np.int(feature[0])\n else:\n return np.int(feature)\n\npixelspacing_r = []\npixelspacing_c = []\nslice_thicknesses = []\npatient_id = []\npatient_pth = []\nrow_values = []\ncolumn_values = []\nwindow_widths = []\nwindow_levels = []\n\nif basepath == \"../input/osic-pulmonary-fibrosis-progression/\":\n patients = train.Patient.unique()[0:N]\nelse:\n patients = train.SeriesInstanceUID.unique()[0:N]\n\nfor patient in patients:\n patient_id.append(patient)\n if basepath == \"../input/osic-pulmonary-fibrosis-progression/\":\n path = train[train.Patient == patient].dcm_path.values[0]\n else:\n path = train[train.SeriesInstanceUID == patient].dcm_path.values[0]\n example_dcm = listdir(path)[0]\n patient_pth.append(path)\n dataset = pydicom.dcmread(path + \"/\" + example_dcm)\n \n window_widths.append(get_window_value(dataset.WindowWidth))\n window_levels.append(get_window_value(dataset.WindowCenter))\n \n spacing = dataset.PixelSpacing\n slice_thicknesses.append(dataset.SliceThickness)\n \n row_values.append(dataset.Rows)\n column_values.append(dataset.Columns)\n pixelspacing_r.append(spacing[0])\n pixelspacing_c.append(spacing[1])\n \nscan_properties = pd.DataFrame(data=patient_id, columns=[\"patient\"])\nscan_properties.loc[:, \"rows\"] = row_values\nscan_properties.loc[:, \"columns\"] = column_values\nscan_properties.loc[:, \"area\"] = scan_properties[\"rows\"] * scan_properties[\"columns\"]\nscan_properties.loc[:, \"pixelspacing_r\"] = pixelspacing_r\nscan_properties.loc[:, \"pixelspacing_c\"] = pixelspacing_c\nscan_properties.loc[:, \"pixelspacing_area\"] = scan_properties.pixelspacing_r * scan_properties.pixelspacing_c\nscan_properties.loc[:, \"slice_thickness\"] = slice_thicknesses\nscan_properties.loc[:, \"patient_pth\"] = patient_pth\nscan_properties.loc[:, \"window_width\"] = window_widths\nscan_properties.loc[:, \"window_level\"] = window_levels\nscan_properties.head()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(pixelspacing_r, ax=ax[0], color=\"Limegreen\", kde=False)\nax[0].set_title(\"Pixel spacing distribution \\n in row direction \")\nax[0].set_ylabel(\"Counts in train\")\nax[0].set_xlabel(\"mm\")\nsns.distplot(pixelspacing_c, ax=ax[1], color=\"Mediumseagreen\", kde=False)\nax[1].set_title(\"Pixel spacing distribution \\n in column direction\");\nax[1].set_ylabel(\"Counts in train\");\nax[1].set_xlabel(\"mm\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"On voit que les valeurs varient vraiment beaucoup d'un patient à l'autre ! Comme elles sont données en mm et que les scans ct couvrent généralement 512 valeurs de lignes et de colonnes... **Nous pouvons calculer la distance minimale et maximale couverte par les images"},{"metadata":{},"cell_type":"markdown","source":"### Épaisseur de la tranche et surface des pixels\n\nL'épaisseur de la tranche nous indique la distance parcourue par une tranche dans la direction Z. Traçons également la distribution de celle-ci. En outre, le tableau de pixels des valeurs brutes couvre une zone spécifique donnée par les valeurs des lignes et des colonnes. Examinons-le également"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"counts = scan_properties.groupby([\"rows\", \"columns\"]).size()\ncounts = counts.unstack()\ncounts.fillna(0, inplace=True)\n\n\nfig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(slice_thicknesses, color=\"orangered\", kde=False, ax=ax[0])\nax[0].set_title(\"Slice thicknesses of all patients\");\nax[0].set_xlabel(\"Slice thickness in mm\")\nax[0].set_ylabel(\"Counts in train\");\n\nfor n in counts.index.values:\n for m in counts.columns.values:\n ax[1].scatter(n, m, s=counts.loc[n,m], c=\"midnightblue\")\nax[1].set_xlabel(\"rows\")\nax[1].set_ylabel(\"columns\")\nax[1].set_title(\"Pixel area of ct-scan per patient\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Des tranches très fines permettent de montrer plus de détails. En revanche, les tranches épaisses contiennent moins de bruit mais sont plus sujettes aux artefacts. Hmm... Je suis très excité de voir quelques exemples ici aussi. \n* Même s'il est courant d'avoir des zones de taille 512x512 pixels, on peut voir que ce n'est pas toujours vrai ! On peut trouver beaucoup d'exceptions et même une ou quelques très grandes zones de pixels (1300x1300) ! !! (OSIC)\n* Un prétraitement correct de ces scans pourrait être très important... nous devons le vérifier."},{"metadata":{},"cell_type":"markdown","source":"## Zone physique et volume de la tranche couverts par un seul ct-scan"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Maintenant, nous connaissons des quantités importantes pour calculer la distance physique couverte par un ct-scan !"},{"metadata":{"trusted":true},"cell_type":"code","source":"scan_properties[\"r_distance\"] = scan_properties.pixelspacing_r * scan_properties.rows\nscan_properties[\"c_distance\"] = scan_properties.pixelspacing_c * scan_properties[\"columns\"]\nscan_properties[\"area_cm2\"] = 0.1* scan_properties[\"r_distance\"] * 0.1*scan_properties[\"c_distance\"]\nscan_properties[\"slice_volume_cm3\"] = 0.1*scan_properties.slice_thickness * scan_properties.area_cm2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(scan_properties.area_cm2, ax=ax[0], color=\"purple\")\nsns.distplot(scan_properties.slice_volume_cm3, ax=ax[1], color=\"magenta\")\nax[0].set_title(\"CT-slice area in $cm^{2}$\")\nax[1].set_title(\"CT-slice volume in $cm^{3}$\")\nax[0].set_xlabel(\"$cm^{2}$\")\nax[1].set_xlabel(\"$cm^{3}$\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Nous avons des images avec des zones et des volumes de sliches extrêmement larges"},{"metadata":{},"cell_type":"markdown","source":"## Surface et volume de la tranche du CT-scan - EDA "},{"metadata":{},"cell_type":"markdown","source":"### La plus petite et la plus grande zone de coupe CT"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"max_path = scan_properties[\n scan_properties.area_cm2 == scan_properties.area_cm2.max()].patient_pth.values[0]\nmin_path = scan_properties[\n scan_properties.area_cm2 == scan_properties.area_cm2.min()].patient_pth.values[0]\n\nmin_scans = load_scans(min_path)\nmin_hu_scans = transform_to_hu(min_scans)\n\nmax_scans = load_scans(max_path)\nmax_hu_scans = transform_to_hu(max_scans)\n\nbackground_water_hu_scans = max_hu_scans.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def set_manual_window(hu_image, custom_center, custom_width):\n w_image = hu_image.copy()\n min_value = custom_center - (custom_width/2)\n max_value = custom_center + (custom_width/2)\n w_image[w_image < min_value] = min_value\n w_image[w_image > max_value] = max_value\n return w_image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,10))\nax[0].imshow(set_manual_window(min_hu_scans[np.int(len(min_hu_scans)/2)], -500, 1000), cmap=\"YlGnBu\")\nax[1].imshow(set_manual_window(max_hu_scans[np.int(len(max_hu_scans)/2)], -500, 1000), cmap=\"YlGnBu\");\nax[0].set_title(\"CT-scan with small slice area\")\nax[1].set_title(\"CT-scan with large slice area\");\nfor n in range(2):\n ax[n].axis(\"off\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Perspectives\n\n* En regardant une tranche de scan avec la plus petite et la plus grande surface, on peut voir que la grande tranche a beaucoup de région inutile couverte. Nous pourrions la recadrer.\n* Bizarre... dans la deuxième image avec la grande surface, la région extérieure du tube du scanner n'est pas réglée sur la valeur de l'air mais plutôt sur une valeur située au milieu de la plage de -1000 à 1000."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(max_hu_scans[np.int(len(max_hu_scans)/2)].flatten(), kde=False, ax=ax[1])\nax[1].set_title(\"Large area image\")\nsns.distplot(min_hu_scans[np.int(len(min_hu_scans)/2)].flatten(), kde=False, ax=ax[0])\nax[0].set_title(\"Small area image\")\nax[0].set_xlabel(\"HU values\")\nax[1].set_xlabel(\"HU values\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Pour l'exemple de l'OSIC, nous pouvons trouver : Ihh... il était réglé sur l'eau par défaut dans la grande image... pourquoi ? ! C'est mauvais ! Nous devons trouver une stratégie pour régler ce problème. Ce n'est pas bon que nous ayons parfois des régions extérieures qui ressemblent à de l'\"eau\" et parfois des régions qui ressemblent à de l'\"air\"."},{"metadata":{},"cell_type":"markdown","source":"### Le plus petit et le plus grand volume de CT-slice"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"max_path = scan_properties[\n scan_properties.slice_volume_cm3 == scan_properties.slice_volume_cm3.max()].patient_pth.values[0]\nmin_path = scan_properties[\n scan_properties.slice_volume_cm3 == scan_properties.slice_volume_cm3.min()].patient_pth.values[0]\n\nmin_scans = load_scans(min_path)\nmin_hu_scans = transform_to_hu(min_scans)\n\nmax_scans = load_scans(max_path)\nmax_hu_scans = transform_to_hu(max_scans)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,10))\nax[0].imshow(set_manual_window(min_hu_scans[np.int(len(min_hu_scans)/2)], -500, 1000), cmap=\"YlGnBu\")\nax[1].imshow(set_manual_window(max_hu_scans[np.int(len(max_hu_scans)/2)], -500, 1000), cmap=\"YlGnBu\");\nax[0].set_title(\"CT-scan with small slice volume\")\nax[1].set_title(\"CT-scan with large slice volume\");\nfor n in range(2):\n ax[n].axis(\"off\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"je ne vois pas une grande différence. Peut-être que celle avec le grand volume de la tranche semble un peu plus floue. Mais comme ci-dessus, il y a une région du scanner extérieur qui a été réglée sur la valeur de l'eau (valeur HU de 0) au lieu de celle de l'air."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(max_hu_scans[np.int(len(max_hu_scans)/2)].flatten(), kde=False, ax=ax[1])\nax[1].set_title(\"Large slice volume\")\nsns.distplot(min_hu_scans[np.int(len(min_hu_scans)/2)].flatten(), kde=False, ax=ax[0])\nax[0].set_title(\"Small slice volume\")\nax[0].set_xlabel(\"HU values\")\nax[1].set_xlabel(\"HU values\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Reconstruction 3D du CT-scans "},{"metadata":{},"cell_type":"markdown","source":"Le plot_3d fonctionne bien dans le Data Science Bowl 2017, mais dans notre cas, les résultats ne sont pas aussi bons. Cela dépend du seuil, mais jusqu'à présent, je ne sais pas pourquoi nos reconstructions ont souvent l'air floues ou montrent aussi des régions du tube"},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_3d(image, threshold=700, color=\"navy\"):\n \n # Position the scan upright, \n # so the head of the patient would be at the top facing the camera\n p = image.transpose(2,1,0)\n \n verts, faces,_,_ = measure.marching_cubes_lewiner(p, threshold)\n\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], alpha=0.2)\n mesh.set_facecolor(color)\n ax.add_collection3d(mesh)\n\n ax.set_xlim(0, p.shape[0])\n ax.set_ylim(0, p.shape[1])\n ax.set_zlim(0, p.shape[2])\n\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(max_hu_scans)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Vous pouvez voir que celui-ci est très différent de celui-là"},{"metadata":{"trusted":true},"cell_type":"code","source":"old_distribution = max_hu_scans.flatten()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"example = train.dcm_path.values[0]\nscans = load_scans(example)\nhu_scans = transform_to_hu(scans)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(hu_scans)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Par rapport à la version précédente, celle-ci est bien plus belle. Traçons les distributions. Peut-être pouvons-nous comprendre ce qui ne va pas en les regardant :"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nsns.distplot(old_distribution, label=\"weak 3d plot\", kde=False)\nsns.distplot(hu_scans.flatten(), label=\"strong 3d plot\", kde=False)\nplt.title(\"HU value distribution\")\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(len(max_hu_scans), len(hu_scans))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Je pense que nous devons comprendre l'algorithme de marching_cubes_lewiner pour comprendre pourquoi l'intrigue fonctionne parfois bien et parfois pas. Mais je pense que ce n'est pas vraiment important pour la compétition elle-même. Pour l'instant, je n'aime pas passer plus de temps sur ce sujet. Il est peut-être plus important de garder à l'esprit que les distributions globales peuvent être différentes."},{"metadata":{},"cell_type":"markdown","source":"## Rééchantillonnage de la taille des voxels"},{"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] + list(scan[0].PixelSpacing), dtype=np.float32)\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, mode='nearest')\n \n return image, new_spacing","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img_resampled, spacing = resample(max_hu_scans, scans, [1,1,1])\nprint(\"Shape before resampling\\t\", hu_scans.shape)\nprint(\"Shape after resampling\\t\", img_resampled.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Tissue segmentation "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def largest_label_volume(im, bg=-1):\n vals, counts = np.unique(im, return_counts=True)\n\n counts = counts[vals != bg]\n vals = vals[vals != bg]\n\n if len(counts) > 0:\n return vals[np.argmax(counts)]\n else:\n return None\n \ndef fill_lungs(binary_image):\n image = binary_image.copy()\n # For every slice we determine the largest solid structure\n for i, axial_slice in enumerate(image):\n axial_slice = axial_slice - 1\n labeling = measure.label(axial_slice)\n l_max = largest_label_volume(labeling, bg=0)\n\n if l_max is not None: #This slice contains some lung\n image[i][labeling != l_max] = 1\n return image\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def segment_lung_mask(image):\n segmented = np.zeros(image.shape) \n \n for n in range(image.shape[0]):\n binary_image = np.array(image[n] > -320, dtype=np.int8)+1\n labels = measure.label(binary_image)\n \n background_label_1 = labels[0,0]\n background_label_2 = labels[0,-1]\n background_label_3 = labels[-1,0]\n background_label_4 = labels[-1,-1]\n \n #Fill the air around the person\n binary_image[background_label_1 == labels] = 2\n binary_image[background_label_2 == labels] = 2\n binary_image[background_label_3 == labels] = 2\n binary_image[background_label_4 == labels] = 2\n \n #We have a lot of remaining small signals outside of the lungs that need to be removed. \n #In our competition closing is superior to fill_lungs \n selem = disk(4)\n binary_image = closing(binary_image, selem)\n \n binary_image -= 1 #Make the image actual binary\n binary_image = 1-binary_image # Invert it, lungs are now 1\n \n segmented[n] = binary_image.copy() * image[n]\n \n return segmented","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Comprendre la segmentation étape par étape:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nsns.distplot(hu_scans[20], kde=False)\nplt.title(\"Example HU value distribution\");\nplt.xlabel(\"HU-value\")\nplt.ylabel(\"count\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Avec -320, nous séparons les poumons (-700) / l'air (-1000) et les tissus avec des valeurs proches de l'eau (0)."},{"metadata":{"trusted":true},"cell_type":"code","source":"binary_image = np.array((hu_scans[20]>-320), dtype=np.int8) + 1\nnp.unique(binary_image)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Les poumons ont des valeurs de 1 ainsi que les milieux atmosphériques. En revanche, les milieux par défaut semblables à ceux de l'eau et de nombreux autres tissus ou fluides organiques ont des valeurs de 2. Comme nous aimons seulement segmenter les poumons, nous devons éliminer le milieu. Dans le cas de l'air et des valeurs par défaut similaires à l'air, nous devons définir manuellement leurs valeurs à 2. Nous pouvons le faire en étiquetant les régions connectées dans l'image binaire et en extrayant les étiquettes de chaque région qui correspond aux coins de la tranche d'image 2D "},{"metadata":{"trusted":true},"cell_type":"code","source":"labels = measure.label(binary_image)\n\nbackground_label_1 = labels[0,0]\nbackground_label_2 = labels[0,-1]\nbackground_label_3 = labels[-1,0]\nbackground_label_4 = labels[-1,-1]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Nous pouvons maintenant définir toutes les régions étiquetées de l'image binaire qui correspondent à ces étiquettes d'angle à la valeur \"non poumon\" 2 "},{"metadata":{"trusted":true},"cell_type":"code","source":"binary_image_2 = binary_image.copy()\nbinary_image_2[background_label_1 == labels] = 2\nbinary_image_2[background_label_2 == labels] = 2\nbinary_image_2[background_label_3 == labels] = 2\nbinary_image_2[background_label_4 == labels] = 2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Le résultat de ces étapes se présente comme suit "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,3,figsize=(20,7))\nax[0].imshow(binary_image, cmap=\"binary\", interpolation='nearest')\nax[1].imshow(labels, cmap=\"jet\", interpolation='nearest')\nax[2].imshow(binary_image_2, cmap=\"binary\", interpolation='nearest')\n\nax[0].set_title(\"Binary image\")\nax[1].set_title(\"Labelled image\");\nax[2].set_title(\"Binary image - background removed\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Perspectives\n\n1. La première image montre le binaire brut. Dans ce cas, nous trouvons l'air comme arrière-plan et nous devons le régler sur la valeur \"non poumon\" de 2.\n2. Pour cela, nous étiquetons toutes les régions connectées dans l'image binaire. Il existe de nombreuses régions étiquetées, mais les seules qui nous intéressent sont les 4 régions d'angle à (0,0), (0,500), (500,0) et (500,500).\n3. La connaissance des étiquettes correspondantes nous aide à régler manuellement le fond à la valeur 2 (noir).\n4. Au final, nous pouvons voir que les poumons sont blancs (1), mais nous trouvons encore beaucoup de signaux restants qui correspondent à des tissus corporels qu'il nous faut encore enlever."},{"metadata":{},"cell_type":"markdown","source":"Pour supprimer les signaux, nous pouvons utiliser la fermeture morphologique. Si vous aimez jouer avec la valeur du disque. Elle doit être suffisamment grande pour annuler les signaux du corps mais suffisamment petite pour garder suffisamment de détails à l'intérieur des poumons"},{"metadata":{"trusted":true},"cell_type":"code","source":"selem = disk(4)\nclosed_binary_2 = closing(binary_image_2, selem)\n\nclosed_binary_2 -= 1 #Make the image actual binary\nclosed_binary_2 = 1-closed_binary_2 # Invert it, lungs are now 1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Comparons avec la méthode de remplissage des poumons de Guidos:"},{"metadata":{"trusted":true},"cell_type":"code","source":"filled_lungs_binary = fill_lungs(binary_image_2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Et avec la suppression de sa poche d'air "},{"metadata":{"trusted":true},"cell_type":"code","source":"air_pocket_binary = closed_binary_2.copy()\n# Remove other air pockets insided body\nlabels_2 = measure.label(air_pocket_binary, background=0)\nl_max = largest_label_volume(labels_2, bg=0)\nif l_max is not None: # There are air pockets\n air_pocket_binary[labels_2 != l_max] = 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,3,figsize=(20,7))\n\nax[0].imshow(closed_binary_2, cmap=\"binary\", interpolation='nearest')\nax[1].imshow(filled_lungs_binary, cmap=\"binary\", interpolation='nearest')\nax[2].imshow(air_pocket_binary, cmap=\"binary\", interpolation='nearest')\n\n\nax[0].set_title(\"Morphological closing\");\nax[1].set_title(\"Guidos filling lung structures\");\nax[2].set_title(\"Guidos air pocket removal\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Perspectives\n\n* La fermeture morphologique a mieux fonctionné que la méthode Guidos. Mais il nous manque souvent beaucoup d'informations à l'intérieur des poumons, juste pour éliminer ces signaux corporels. Il y a certainement place pour des améliorations ! ;-)\n* En revanche, la méthode du fill-lung a des problèmes avec les signaux restants du corps et ne donne pas ce que nous aimons obtenir.\n* En outre, l'élimination des poches d'air ne fonctionne pas bien non plus.\n\n"},{"metadata":{},"cell_type":"markdown","source":"### solution Finale"},{"metadata":{},"cell_type":"markdown","source":"Et voici à quoi cela ressemble si nous masquons l'image originale avec les poumons binaires segmentés dans notre exemple en 2D"},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented = segment_lung_mask(np.array([hu_scans[20]]))\n\nfig, ax = plt.subplots(1,2,figsize=(20,10))\nax[0].imshow(hu_scans[20], cmap=\"Blues_r\")\nax[1].imshow(segmented[0], cmap=\"Blues_r\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Et nous pouvons également vérifier l'aspect de l'affaire 3D"},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented_lungs = segment_lung_mask(hu_scans)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(6,5, figsize=(20,20))\nfor n in range(6):\n for m in range(5):\n ax[n,m].imshow(segmented_lungs[n*5+m], cmap=\"Blues_r\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Il est intéressant de noter que nous avons encore quelques signaux en dehors des poumons"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(segmented_lungs, threshold=-600)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"améliorer la segmentation"},{"metadata":{},"cell_type":"markdown","source":"# Génération d'un ensemble de données pour les fichiers prétraités \n\n\n## Traitement des différentes tailles d'images \n\nPour générer les données, nous devons réexaminer les différentes tailles d'images : Pour l'OSIC, nous avons deux grands groupes de taille et quelques petites valeurs aberrantes. Par exemple, nous pourrions redimensionner ou recadrer manuellement les valeurs aberrantes et trouver une stratégie pour les deux grands groupes. Examinons à nouveau les tailles :"},{"metadata":{"trusted":true},"cell_type":"code","source":"image_sizes = scan_properties.groupby([\"rows\", \"columns\"]).size().sort_values(ascending=False)\nimage_sizes","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Pour l'OSIC, il est intéressant de constater que nous avons deux types de modèles différents. Dans la compétition RSNA, toutes les images d'entraînement sont de forme 512, 512. Vous n'avez donc pas besoin d'explorer davantage les tailles des images. Mais peut-être est-il encore utile de parcourir les images pour trouver des idées de bonnes augmentations."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plt.figure(figsize=(8,8))\nfor n in counts.index.values:\n for m in counts.columns.values:\n plt.scatter(n, m, s=counts.loc[n,m], c=\"dodgerblue\", alpha=0.7)\nplt.xlabel(\"rows\")\nplt.ylabel(\"columns\")\nplt.title(\"Pixel area of ct-scan per patient\");\nplt.plot(np.arange(0,1400), '-.', c=\"purple\", label=\"squared\")\nplt.plot(888 * np.ones(1400), '-.', c=\"crimson\", label=\"888 rows\");\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ce genre de modèles évidents vaut toujours la peine d'être examiné. Peut-être pouvons-nous trouver une sorte de règle qui nous permette de créer une bonne stratégie de redimensionnement. Voici un petit observateur qui vous permet de parcourir les fichiers en exécutant à nouveau la cellule de code :"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"class ImageObserver:\n \n def __init__(self, scan_properties, batch_size):\n self.scan_properties = scan_properties\n self.batch_size = batch_size\n \n def select_group(self, group=(512,512)):\n self.group = group\n self.name = \"rows {}, columns {}\".format(group[0], group[1])\n self.batch_shape = (self.batch_size, group[0], group[1])\n self.selection = self.scan_properties[\n (self.scan_properties[\"rows\"]==group[0]) & (self.scan_properties[\"columns\"]==group[1])\n ].copy()\n self.patient_pths = self.selection.patient_pth.unique()\n \n \n def get_loader(self):\n \n idx=0\n images = np.zeros(self.batch_shape)\n \n for path in self.patient_pths:\n \n scans = load_scans(path)\n hu_scans = transform_to_hu(scans)\n images[idx,:,:] = hu_scans[0]\n \n idx += 1\n if idx == self.batch_shape[0]:\n yield images\n images = np.zeros(self.batch_shape)\n idx = 0\n if idx > 0:\n yield images","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"my_choice = image_sizes.index.values[0]\nprint(my_choice)\nto_display = 4","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"observer = ImageObserver(scan_properties, to_display)\nobserver.select_group(my_choice)\nobserver_iterator = observer.get_loader()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Il suffit de relancer la cellule suivante, pour observer le prochain lot d'images :"},{"metadata":{"trusted":true},"cell_type":"code","source":"images = next(observer_iterator)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,to_display,figsize=(20,5))\n\n\nfor m in range(to_display):\n image = images[m]\n ax[m].imshow(set_manual_window(image, -500, 1000), cmap=\"YlGnBu\")\n ax[m].set_title(observer.name)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* Les grandes tailles d'images carrées présentent souvent des résolutions plus élevées qu'avec 512 lignes et 512 colonnes. Le redimensionnement en fonction des deux grands groupes (512, 512) ou (768, 768) est ici logique.\n* Dans les cas non carrés, les cultures centrales devraient fonctionner au mieux car elles n'ont que des valeurs de fond plus importantes mais appartiennent toujours aux grands groupes de la région du scanner interne."},{"metadata":{},"cell_type":"markdown","source":"### pour en savoir plus sur les unités Hounsfield"},{"metadata":{"trusted":true},"cell_type":"code","source":"from IPython.display import HTML\nHTML('