{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Positioning of Cervical Spines pt 2","metadata":{}},{"cell_type":"markdown","source":"### Imports","metadata":{}},{"cell_type":"code","source":"! pip install python-gdcm\n! pip install pylibjpeg pylibjpeg-libjpeg pydicom","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-08-30T22:24:27.207465Z","iopub.execute_input":"2022-08-30T22:24:27.208371Z","iopub.status.idle":"2022-08-30T22:24:55.662459Z","shell.execute_reply.started":"2022-08-30T22:24:27.208278Z","shell.execute_reply":"2022-08-30T22:24:55.66136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%matplotlib inline\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\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\n\nimport os\nimport scipy.ndimage\nimport matplotlib.pyplot as plt\nfrom glob import glob\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nimport seaborn as sns\n\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 *\nimport gc\nimport nibabel as nib\nfrom tqdm import tqdm\nimport gdcm\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-30T22:24:55.664705Z","iopub.execute_input":"2022-08-30T22:24:55.665412Z","iopub.status.idle":"2022-08-30T22:24:57.754724Z","shell.execute_reply.started":"2022-08-30T22:24:55.665354Z","shell.execute_reply":"2022-08-30T22:24:57.753351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Helper Functions\n\nload_scan to load the dicom scans and get_pixels_hu to transform the image data into [Hounsfield Units](https://en.wikipedia.org/wiki/Hounsfield_scale)","metadata":{}},{"cell_type":"code","source":"def load_scan(path):\n    \n    try: \n        slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n        \n    except: \n        reader = gdcm.ImageReader()        \n        ret = reader.Read()\n        slices = [reader.SetFileName(path + '/' + s).Read() for s in os.listdir(path)]\n    \n    \n        \n    slices.sort(key = lambda x: int(x.InstanceNumber))\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        \n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices\n\ndef get_pixels_hu(scans):\n    image = np.stack([s.pixel_array for s in scans])\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 1\n    # The intercept is usually -1024, so air is approximately 0\n    image[image == -2000] = 0\n    \n    # Convert to Hounsfield units (HU)\n    intercept = scans[0].RescaleIntercept\n    slope = scans[0].RescaleSlope\n    \n    if slope != 1:\n        image = slope * image.astype(np.float64)\n        image = image.astype(np.int16)\n        \n    image += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)\n\n","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:24:57.75636Z","iopub.execute_input":"2022-08-30T22:24:57.757343Z","iopub.status.idle":"2022-08-30T22:24:57.769531Z","shell.execute_reply.started":"2022-08-30T22:24:57.757304Z","shell.execute_reply":"2022-08-30T22:24:57.768123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Analysing a sample","metadata":{}},{"cell_type":"code","source":"data_path = '../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.10001'\npatient = load_scan(data_path)\nimgs = get_pixels_hu(patient)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:24:57.77207Z","iopub.execute_input":"2022-08-30T22:24:57.772827Z","iopub.status.idle":"2022-08-30T22:25:01.915629Z","shell.execute_reply.started":"2022-08-30T22:24:57.77279Z","shell.execute_reply":"2022-08-30T22:25:01.914424Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Hounsfield Units","metadata":{}},{"cell_type":"markdown","source":"<table>\n<tbody><tr>\n<th>Substance</th>\n<th>HU</th>\n</tr>\n<tr>\n<td>Air</td>\n<td>−1000</td>\n</tr>\n<tr>\n<td>Lung</td>\n<td>−500</td>\n</tr>\n<tr>\n<td>Fat</td>\n<td>−100 to −50</td>\n</tr>\n<tr>\n<td>Water</td>\n<td>0</td>\n</tr>\n<tr>\n<td>Blood</td>\n<td>+30 to +70</td>\n</tr>\n<tr>\n<td>Muscle</td>\n<td>+10 to +40</td>\n</tr>\n<tr>\n<td>Liver</td>\n<td>+40 to +60</td>\n</tr>\n<tr>\n<td>Bone</td>\n<td>+700 (cancellous bone) to +3000 (cortical bone)</td>\n</tr>\n</tbody>\n</table>","metadata":{}},{"cell_type":"code","source":"plt.hist(imgs.flatten(), bins=50, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:25:01.917161Z","iopub.execute_input":"2022-08-30T22:25:01.917675Z","iopub.status.idle":"2022-08-30T22:25:03.777116Z","shell.execute_reply.started":"2022-08-30T22:25:01.917626Z","shell.execute_reply":"2022-08-30T22:25:03.775878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Displaying sample","metadata":{}},{"cell_type":"code","source":"def sample_stack(stack, rows=6, cols=6, start_with=10, show_every=1):\n    fig,ax = plt.subplots(rows,cols,figsize=[12,12])\n    for i in range(rows*cols):\n        ind = start_with + i*show_every\n        ax[int(i/rows),int(i % rows)].set_title('slice %d' % ind)\n        ax[int(i/rows),int(i % rows)].imshow(stack[ind],cmap='gray')\n        ax[int(i/rows),int(i % rows)].axis('off')\n    plt.show()\n\nsample_stack(imgs)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:25:03.778657Z","iopub.execute_input":"2022-08-30T22:25:03.779642Z","iopub.status.idle":"2022-08-30T22:25:06.688354Z","shell.execute_reply.started":"2022-08-30T22:25:03.779591Z","shell.execute_reply":"2022-08-30T22:25:06.687307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3D Plot","metadata":{}},{"cell_type":"markdown","source":"Now we want to see how will our data look if we transform all our scans into a single 3D image. We will set a threshold to separate Bone from other structures in the scans based on the hounsfield units","metadata":{}},{"cell_type":"code","source":"def make_mesh(image, threshold=-300, step_size=1):\n\n    print(\"Transposing surface\")\n    p = image.transpose(2,1,0)\n    \n    print (\"Calculating surface\")\n    verts, faces, norm, val = measure.marching_cubes(p, threshold, step_size=step_size, allow_degenerate=True) \n    return verts, faces\n\n\n","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:25:06.690213Z","iopub.execute_input":"2022-08-30T22:25:06.690592Z","iopub.status.idle":"2022-08-30T22:25:06.698035Z","shell.execute_reply.started":"2022-08-30T22:25:06.690557Z","shell.execute_reply":"2022-08-30T22:25:06.696808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"400 is our threshold, pixels below that threshold will become 0 valued","metadata":{}},{"cell_type":"code","source":"v, f = make_mesh(imgs, 400,2)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:25:06.699903Z","iopub.execute_input":"2022-08-30T22:25:06.700713Z","iopub.status.idle":"2022-08-30T22:25:08.50376Z","shell.execute_reply.started":"2022-08-30T22:25:06.700671Z","shell.execute_reply":"2022-08-30T22:25:08.500556Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plt_3d(verts, faces):\n    print (\"Drawing\")\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,edgecolor=(0,0,0))\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_facecolor((0.7, 0.7, 0.7))\n    plt.show()\n    \nplt_3d(v,f)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:25:08.508929Z","iopub.execute_input":"2022-08-30T22:25:08.510434Z","iopub.status.idle":"2022-08-30T22:27:08.461283Z","shell.execute_reply.started":"2022-08-30T22:25:08.510209Z","shell.execute_reply":"2022-08-30T22:27:08.459753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analysing DICOM positioning and spacing in the Dataset","metadata":{}},{"cell_type":"code","source":"def read_dicom(path, voi_lut = True, fix_monochrome = True):\n    '''ref: https://www.kaggle.com/code/raddar/convert-dicom-to-np-array-the-correct-way\n    '''\n    dicom = pydicom.read_file(path)\n    \n    # VOI LUT (if available by DICOM device) is used to transform raw DICOM data to \"human-friendly\" view\n    if voi_lut:\n        data = apply_voi_lut(dicom.pixel_array, dicom)\n    else:\n        data = dicom.pixel_array\n               \n    # depending on this value, X-ray may look inverted - fix that:\n    if fix_monochrome and dicom.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n        \n    data = data - np.min(data)\n    data = data / np.max(data)\n    data = (data * 255).astype(np.uint8)\n    return data, dicom","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:08.466616Z","iopub.execute_input":"2022-08-30T22:27:08.467056Z","iopub.status.idle":"2022-08-30T22:27:08.476756Z","shell.execute_reply.started":"2022-08-30T22:27:08.467019Z","shell.execute_reply":"2022-08-30T22:27:08.475127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"coordinates_df = pd.read_csv('../input/coordinates-rsna2022/coordinates.csv')\n\nx=coordinates_df['X']\ny=coordinates_df['Y']\nz=coordinates_df['Z']\n\n\n    ","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:08.479004Z","iopub.execute_input":"2022-08-30T22:27:08.479571Z","iopub.status.idle":"2022-08-30T22:27:08.973055Z","shell.execute_reply.started":"2022-08-30T22:27:08.479511Z","shell.execute_reply":"2022-08-30T22:27:08.971616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nplt.figure(figsize=(15,15))\nsns.histplot(x=x)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:08.974938Z","iopub.execute_input":"2022-08-30T22:27:08.975367Z","iopub.status.idle":"2022-08-30T22:27:09.988584Z","shell.execute_reply.started":"2022-08-30T22:27:08.975327Z","shell.execute_reply":"2022-08-30T22:27:09.986841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15,15))\nsns.histplot(x=y)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:09.990595Z","iopub.execute_input":"2022-08-30T22:27:09.990974Z","iopub.status.idle":"2022-08-30T22:27:11.064336Z","shell.execute_reply.started":"2022-08-30T22:27:09.990942Z","shell.execute_reply":"2022-08-30T22:27:11.063002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15,15))\nsns.histplot(x=z)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:11.066068Z","iopub.execute_input":"2022-08-30T22:27:11.06723Z","iopub.status.idle":"2022-08-30T22:27:12.845478Z","shell.execute_reply.started":"2022-08-30T22:27:11.067174Z","shell.execute_reply":"2022-08-30T22:27:12.844005Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It looks like our Z axis is not well distributed as we thought, according to [this website](https://dicom.innolitics.com/ciods/rt-dose/image-plane/00200032) on DICOM data:\n\n        Sxyz The three values of Image Position (Patient) (0020,0032). It is the location in mm from the origin of the RCS.\n\nFrom the **origin** of the RCS, so this may vary from machine to machine, or something like that. But all axis are in milimeters, so maybe if we get the difference of the top Z value to the bottom Z value we will have the distribution we are looking for","metadata":{}},{"cell_type":"code","source":"z_diffs = pd.read_csv('../input/coordinates-rsna2022/Z_diffs.csv')['Z_diffs']\n\n    ","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:12.847302Z","iopub.execute_input":"2022-08-30T22:27:12.848693Z","iopub.status.idle":"2022-08-30T22:27:12.862822Z","shell.execute_reply.started":"2022-08-30T22:27:12.848644Z","shell.execute_reply":"2022-08-30T22:27:12.861187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.histplot(x=z_diffs)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:12.864795Z","iopub.execute_input":"2022-08-30T22:27:12.865971Z","iopub.status.idle":"2022-08-30T22:27:13.191684Z","shell.execute_reply.started":"2022-08-30T22:27:12.865911Z","shell.execute_reply":"2022-08-30T22:27:13.190373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Much better, for now on this will be our metric for the Z axis, the distance of the slice to the top slice.","metadata":{}},{"cell_type":"markdown","source":"Here are some patients which Z axis ranges above 500","metadata":{}},{"cell_type":"markdown","source":"Let us look into the distribution of the axis values in one of them","metadata":{}},{"cell_type":"code","source":"images = glob('../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.3850/*')\nx0 = []\ny0 = []\nz0 = []\nfor image in images:\n    data = read_dicom(image)[1]\n    \n    x0.append(data[0x0020, 0x0032][0])\n    y0.append(data[0x0020, 0x0032][1])\n    z0.append(data[0x0020, 0x0032][2])\n    \n","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:13.193956Z","iopub.execute_input":"2022-08-30T22:27:13.194899Z","iopub.status.idle":"2022-08-30T22:27:25.219694Z","shell.execute_reply.started":"2022-08-30T22:27:13.194843Z","shell.execute_reply":"2022-08-30T22:27:25.217887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.histplot(x=x0)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:25.221703Z","iopub.execute_input":"2022-08-30T22:27:25.222121Z","iopub.status.idle":"2022-08-30T22:27:25.423828Z","shell.execute_reply.started":"2022-08-30T22:27:25.222081Z","shell.execute_reply":"2022-08-30T22:27:25.422929Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.histplot(x=y0)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:25.424962Z","iopub.execute_input":"2022-08-30T22:27:25.426195Z","iopub.status.idle":"2022-08-30T22:27:25.670423Z","shell.execute_reply.started":"2022-08-30T22:27:25.42614Z","shell.execute_reply":"2022-08-30T22:27:25.669215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.histplot(x=z0)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:25.672212Z","iopub.execute_input":"2022-08-30T22:27:25.672619Z","iopub.status.idle":"2022-08-30T22:27:25.942156Z","shell.execute_reply.started":"2022-08-30T22:27:25.672582Z","shell.execute_reply":"2022-08-30T22:27:25.940254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Lets look at how many slices we have for this patient","metadata":{"execution":{"iopub.status.busy":"2022-08-11T00:55:06.204648Z","iopub.execute_input":"2022-08-11T00:55:06.205232Z","iopub.status.idle":"2022-08-11T00:58:57.768225Z","shell.execute_reply.started":"2022-08-11T00:55:06.205198Z","shell.execute_reply":"2022-08-11T00:58:57.766962Z"}}},{"cell_type":"code","source":"print(len(z0))","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:25.944669Z","iopub.execute_input":"2022-08-30T22:27:25.945578Z","iopub.status.idle":"2022-08-30T22:27:25.953004Z","shell.execute_reply.started":"2022-08-30T22:27:25.94551Z","shell.execute_reply":"2022-08-30T22:27:25.951875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Checking the distance in milimeters from top Z to bottom Z","metadata":{}},{"cell_type":"code","source":"max(z0) - min(z0)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:25.955167Z","iopub.execute_input":"2022-08-30T22:27:25.955565Z","iopub.status.idle":"2022-08-30T22:27:25.970189Z","shell.execute_reply.started":"2022-08-30T22:27:25.955529Z","shell.execute_reply":"2022-08-30T22:27:25.968931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let us check if the slice thickness of our dicoms vary for this patient","metadata":{}},{"cell_type":"code","source":"images = glob('../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.3850/*')\n\nslice_thickness=[]\n\nfor image in images:\n    data = read_dicom(image)[1]\n    slice_thickness.append(data.SliceThickness)\n    \n    \n    \n    ","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:25.971783Z","iopub.execute_input":"2022-08-30T22:27:25.972144Z","iopub.status.idle":"2022-08-30T22:27:31.088738Z","shell.execute_reply.started":"2022-08-30T22:27:25.972102Z","shell.execute_reply":"2022-08-30T22:27:31.087407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.unique(np.array(slice_thickness))","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:31.091256Z","iopub.execute_input":"2022-08-30T22:27:31.091652Z","iopub.status.idle":"2022-08-30T22:27:31.101083Z","shell.execute_reply.started":"2022-08-30T22:27:31.091619Z","shell.execute_reply":"2022-08-30T22:27:31.099982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So one thing to note, is that 688 times 0.5 equals 344,that should be the distance between top Z and bottom Z.","metadata":{}},{"cell_type":"markdown","source":"Our `load_scan` function does the following operation to our scans:\n\n`slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])`","metadata":{}},{"cell_type":"markdown","source":"This makes it so that our slice thickness actually reflects on the milimeter positioning of the scan, while the data in the dicom file doesn't seem to reflect that as it should.","metadata":{}},{"cell_type":"code","source":"images = '../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.3850'\n\nslice_thickness=[]\ndata = load_scan(images)\nprint(data[1].SliceThickness)\n","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:31.102632Z","iopub.execute_input":"2022-08-30T22:27:31.102965Z","iopub.status.idle":"2022-08-30T22:27:32.488723Z","shell.execute_reply.started":"2022-08-30T22:27:31.102935Z","shell.execute_reply":"2022-08-30T22:27:32.487445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now lets check if we get the right Z axis range","metadata":{}},{"cell_type":"code","source":"0.40000000000009095 * 688","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:32.490804Z","iopub.execute_input":"2022-08-30T22:27:32.491176Z","iopub.status.idle":"2022-08-30T22:27:32.499658Z","shell.execute_reply.started":"2022-08-30T22:27:32.491142Z","shell.execute_reply":"2022-08-30T22:27:32.497862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Not exactly the same, but much better.","metadata":{}},{"cell_type":"markdown","source":"Lets see how this 27cm length scan looks like","metadata":{}},{"cell_type":"code","source":"import gc\ngc.collect()\ndata_path = '../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.3850'\npatient = load_scan(data_path)\nimgs = get_pixels_hu(patient)\nv, f = make_mesh(imgs, 400,2)\nplt_3d(v,f)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:27:32.501333Z","iopub.execute_input":"2022-08-30T22:27:32.502173Z","iopub.status.idle":"2022-08-30T22:31:24.574542Z","shell.execute_reply.started":"2022-08-30T22:27:32.502127Z","shell.execute_reply":"2022-08-30T22:31:24.573193Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It looks like this patient scan reaches all the way down to the chest","metadata":{}},{"cell_type":"code","source":"\ndata_path = '../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.3850'\npatient = load_scan(data_path)\nimgs = get_pixels_hu(patient)\nsample_stack(imgs, rows=10, cols=10,show_every=3,start_with=350)\n    ","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:31:24.585372Z","iopub.execute_input":"2022-08-30T22:31:24.586454Z","iopub.status.idle":"2022-08-30T22:31:41.08717Z","shell.execute_reply.started":"2022-08-30T22:31:24.58638Z","shell.execute_reply":"2022-08-30T22:31:41.085911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(imgs[458],cmap='gray')","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:31:41.088973Z","iopub.execute_input":"2022-08-30T22:31:41.089737Z","iopub.status.idle":"2022-08-30T22:31:41.288883Z","shell.execute_reply.started":"2022-08-30T22:31:41.089687Z","shell.execute_reply":"2022-08-30T22:31:41.287279Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, we have lung tissue in this scan. We don't want this part of the scan as our targets are the Cervical spine not the Thorax","metadata":{}},{"cell_type":"markdown","source":"## The Problem","metadata":{}},{"cell_type":"markdown","source":"The submission file asks us to predict the probability of a fracture for each cervical vertebrae (C1-C7) and also predict in overall probability of a cervical fracture in the patient, i.e. a fracute in any of the cervical vertebrae.\n\nWhat we need to do is normalize the data in a way that we know where each vertebrae is and deal with each vertebrae accordingly, we can possibly do this by using our Z axis values. After we have done this, then we can proceed to predictions and modeling.\n\nWe also need to remove the thorax from our data and leave only our targets, the cervical spine. Also in my case, I will be aiming at using a 3D convolutional neural net, so we need to normalize the data into a fixed amount of slices for the model. \n\nThe way we can achieve both of these is by analysing the distance of each cervical vertebrae to the head, and maybe come up with a threshold with that data. C1 for example should have a distance closer to 0, C7 should have a further distance, but not so far, as at some point we will reach the thorax if we go further down the spine.","metadata":{}},{"cell_type":"markdown","source":"The way we will do this is by using the labeled Nii images. So below we will store all the Patient IDs that have a labeled nii file","metadata":{}},{"cell_type":"code","source":"from glob import glob\nnii_files = glob('../input/rsna-2022-cervical-spine-fracture-detection/segmentations/*')\npatients_with_labels =[]\nfor patient in nii_files:\n    patient_id =  patient.split('/')[-1][:-4]\n    patients_with_labels.append(patient_id)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:31:41.290457Z","iopub.execute_input":"2022-08-30T22:31:41.290825Z","iopub.status.idle":"2022-08-30T22:31:41.327203Z","shell.execute_reply.started":"2022-08-30T22:31:41.290791Z","shell.execute_reply":"2022-08-30T22:31:41.325839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we will load a single sample, and transpose it to the Axial view. After that we will check if the shape matches with the DICOM files for the same patient","metadata":{}},{"cell_type":"code","source":"nii_example = nib.load(f'../input/rsna-2022-cervical-spine-fracture-detection/segmentations/{patients_with_labels[0]}.nii').get_fdata()[:, ::-1, ::-1].transpose(2, 1, 0)\nnii_example.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:31:41.329534Z","iopub.execute_input":"2022-08-30T22:31:41.329935Z","iopub.status.idle":"2022-08-30T22:31:42.155991Z","shell.execute_reply.started":"2022-08-30T22:31:41.32989Z","shell.execute_reply":"2022-08-30T22:31:42.154571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_path = f'../input/rsna-2022-cervical-spine-fracture-detection/train_images/{patients_with_labels[0]}'\nimgs = np.stack([s.pixel_array for s in load_scan(data_path)])\nimgs.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:31:42.157885Z","iopub.execute_input":"2022-08-30T22:31:42.15833Z","iopub.status.idle":"2022-08-30T22:31:45.077544Z","shell.execute_reply.started":"2022-08-30T22:31:42.158294Z","shell.execute_reply":"2022-08-30T22:31:45.076288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So there is a nii image for each slice in the DICOM files. Lets take a look at what the nii images contain inside, so you can understand where the label comes from\n\n","metadata":{}},{"cell_type":"code","source":"np.unique(nii_example)","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:31:45.079216Z","iopub.execute_input":"2022-08-30T22:31:45.080587Z","iopub.status.idle":"2022-08-30T22:31:47.019114Z","shell.execute_reply.started":"2022-08-30T22:31:45.08054Z","shell.execute_reply":"2022-08-30T22:31:47.017752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As you can see we have numbers from 0 to 9. 0 in this case is just the background, as the files serves as masks as well. But 1 to 9 are pixels that label the vertebrae, being 1-7 the cervical vertebrae and beyond that are the thoraxic vertebrae. \n\nNow since there is a nii slice for each dicom slice, we can get the label for each dicom slice. This is how we can start to figure out which bone is in a dicom slice","metadata":{}},{"cell_type":"markdown","source":"Here we iterate through all the patients that have labels, load the Nii and Dicom files, sorted and shaped as we need them. Then for each dicom slice, we will get its **position in the Z axis**, subtract that from the **top Z axis** for the **whole patient scan**. Then we get the unique values of the respective nii slice, which are the labels, and put both the **Z axis difference** and the **nii labels** in a **list** ","metadata":{}},{"cell_type":"code","source":"bone_labels=[]\n\nfor patient in tqdm(patients_with_labels):\n    nii = nib.load(f'../input/rsna-2022-cervical-spine-fracture-detection/segmentations/{patient}.nii').get_fdata()[:, ::-1, ::-1].transpose(2, 1, 0)\n    \n    path = f'../input/rsna-2022-cervical-spine-fracture-detection/train_images/{patient}'\n    patient = load_scan(path)\n    zmax = patient[0][0x0020, 0x0032][-1]\n\n    for i, _ in enumerate(nii):\n        \n        \n        respective_nii = np.unique(nii[i])\n        data = patient[i]\n        z = np.absolute(zmax - data[0x0020, 0x0032][-1])\n        bone_labels.append([respective_nii , z])\n        \n\n    \n    del patient\n    del nii\n    \n    ","metadata":{"execution":{"iopub.status.busy":"2022-08-30T22:31:47.020932Z","iopub.execute_input":"2022-08-30T22:31:47.021445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we will create a list for each cervical vertebrae where we **append the Z axis difference to the head** for **each nii image** that has its label on it.","metadata":{}},{"cell_type":"code","source":"c1 = []\nc2 = []\nc3 = []\nc4 = []\nc5 = []\nc6 = []\nc7 = []\n\n\n  \n\nfor data in bone_labels:\n\n    \n       \n    if len(data[0])==1:continue\n    \n    elif 1.0 in data[0]:\n        \n        c1.append(data[1])    \n    elif 2.0 in data[0]:\n\n        c2.append(data[1])       \n    elif 3.0 in data[0]:\n\n        c3.append(data[1]) \n    elif 4.0 in data[0]:\n\n        c4.append(data[1])\n    elif 5.0 in data[0]:\n\n        c5.append(data[1])\n    elif 6.0 in data[0]:\n\n        c6.append(data[1])\n    elif 7.0 in data[0]:\n\n        c7.append(data[1])\n    \n\n   \n    ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# BINGO!!","metadata":{}},{"cell_type":"markdown","source":"As we can see below, there is a nice distribution for each of the cervical vertebrae, and they follow a clear trend in distance as we thought. ","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(20,20))\n\n\nsns.histplot(x=c1,label='C1', color='red')\n\nsns.histplot(x=c2,label='C2', color='pink')\n\nsns.histplot(x=c3,label='C3', color='blue')\n\nsns.histplot(x=c4,label='C4', color='yellow')\n\nsns.histplot(x=c5,label='C5', color='green')\n\nsns.histplot(x=c6,label='C6', color='purple')\n\nsns.histplot(x=c7,label='C7', color='black')\n\nplt.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are a few outliers but it's probably no biggie, we can deal with them easily","metadata":{}},{"cell_type":"markdown","source":"Now a one more variables we need to worry about:\n\n* Slice Thickness\n\n\nThe reason we need to check these variables, specially the slice thickness, is that for each slice we should expect a change on our **position values** because of these variables. Lets say a patients' scan has a slice thickness of 0.5, we should expect a 0.5 variation in the Z axis for each dicom slice we open form that scan. \n\nSo lets say, from the plot above, we want all slices in a 200mm range from the top of the scan, so we can get only our cervical spine from the scan. If the slice thickness of that scan is 0.5, we can expect to get 400 slices, if the slice thickness is 1, we can expect 200 slices. So if the slice thickness varies, our number of slices vary as well, and that's not good if we want to train a 3D model on the data.\n\nPixel spacing on the other hand influences on the Sagittal and Coronal plane, i.e. the X and Y axis. So let us check the distribution of these variables.","metadata":{}},{"cell_type":"code","source":"\ndf = pd.read_csv('../input/slicethickness-rsna2022/slice_thickness.csv')\n\n\nslice_thickness=df['SliceThickness']\nspacing0=df['Spacing0']\nspacing1=df['Spacing1']\n\n\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.histplot(x=slice_thickness)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.histplot(x=spacing0)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.histplot(x=spacing1)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our slice thickness varies alot, this is not good. Luckily there is a way to actually alter this, and theoretically change the slice thickness of a scan\n\nLets take a scan that has a slice thickness of 1mm","metadata":{}},{"cell_type":"markdown","source":"This is our function to change our slice thickness","metadata":{}},{"cell_type":"code","source":"def resample(image, scan, new_spacing=1):\n    # Determine current pixel spacing\n    new_spacing = [new_spacing,scan[0].PixelSpacing[0], scan[0].PixelSpacing[1]]\n    spacing = np.array([scan[0].SliceThickness, scan[0].PixelSpacing[0], scan[0].PixelSpacing[1]], 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)\n    \n    return image, new_spacing\n\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Code to resample the images with a slice thickness of 0.5","metadata":{}},{"cell_type":"code","source":"data_path = '../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.28580'\npatient = load_scan(data_path)\nimgs_before = get_pixels_hu(patient)\nimgs, spacing= resample(imgs_before, patient, new_spacing=0.5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Old slice thickness\n- 3D Shape vefore\n- New resampled 3D shape","metadata":{}},{"cell_type":"code","source":"print(patient[0].SliceThickness)\nprint(imgs_before.shape)\nprint(imgs.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Expected distance in mm from top to bottom scan","metadata":{}},{"cell_type":"code","source":"478 * 0.5","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Real distance in mm from top to bottom","metadata":{}},{"cell_type":"code","source":"images = glob('../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.28580/*')\nz0 = []\nfor image in images:\n    data = read_dicom(image)[1]\n    z0.append(data[0x0020, 0x0032][2])\n    \nprint(max(z0)-min(z0))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plot resampled images","metadata":{}},{"cell_type":"code","source":"gc.collect()\n\nv, f = make_mesh(imgs, 400,2)\nplt_3d(v,f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Changing samples","metadata":{}},{"cell_type":"code","source":"data_path = '../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.1480'\npatient = load_scan(data_path)\nimgs = get_pixels_hu(patient)\nimgs, spacing= resample(imgs, patient, new_spacing=0.5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Lets check for some outliers in our axes positioning with the 3-sigma rule","metadata":{}},{"cell_type":"code","source":"import math\nmu  = sum(c1)/len(c1)\nstd = np.std(c1)\noutliers = [c for c in c1 if (abs(c-mu) > 3*std )]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.histplot(x=outliers)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see above, none of these fit the expected distance from the top of the scan for a c1 vertebrae.\n\nLet's create a function that removes these outliers from our data","metadata":{}},{"cell_type":"code","source":"\ndef remove_outliers(lst:list):\n\n    mu  = sum(lst)/len(lst)\n    std = np.std(lst)\n    outlier_clean = [c for c in lst if not(abs(c-mu) > 3*std )]\n    \n    return outlier_clean\n\nc1 = remove_outliers(c1)\nc2 = remove_outliers(c2)\nc3 = remove_outliers(c3)\nc4 = remove_outliers(c4)\nc5 = remove_outliers(c5)\nc6 = remove_outliers(c6)\nc7 = remove_outliers(c7)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now to plot the data without the outliers","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\n\n\nsns.histplot(x=c1,label='C1', color='red')\n\nsns.histplot(x=c2,label='C2', color='pink')\n\nsns.histplot(x=c3,label='C3', color='blue')\n\nsns.histplot(x=c4,label='C4', color='yellow')\n\nsns.histplot(x=c5,label='C5', color='green')\n\nsns.histplot(x=c6,label='C6', color='purple')\n\nsns.histplot(x=c7,label='C7', color='black')\n\nplt.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data Augmentation","metadata":{}},{"cell_type":"markdown","source":"Lets start transforming our data in the X and y values for training a classifier model","metadata":{}},{"cell_type":"code","source":"from sklearn import svm\n\nX = []\ny = []\n\nfor c in c1:\n    X.append(c)\n    y.append(1)     \nfor c in c2:\n    X.append(c)\n    y.append(2) \nfor c in c3:\n    X.append(c)\n    y.append(3) \nfor c in c4:\n    X.append(c)\n    y.append(4)\nfor c in c5:\n    X.append(c)\n    y.append(5) \nfor c in c6:\n    X.append(c)\n    y.append(6) \nfor c in c7:\n    X.append(c)\n    y.append(7) \n    \n    \n    \n\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is the distribution of labels that we have in the dataset","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.histplot(x=y)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see the distribution of labels is not balanced. We will be using **np.random.normal** to create random numbers inside the normal distribution of each label, so we have our classes balanced and our data distributions still normalized","metadata":{}},{"cell_type":"code","source":"c2_diff= len(c1)-len(c2)\nc3_diff= len(c1)-len(c3)\nc4_diff= len(c1)-len(c4)\nc5_diff= len(c1)-len(c5)\nc6_diff= len(c1)-len(c6)\nc7_diff= len(c1)-len(c7)\n\nc2_train = np.concatenate((np.array(c2), np.random.normal(np.mean(c2), np.std(c2), c2_diff)))\nc3_train = np.concatenate((np.array(c3), np.random.normal(np.mean(c3), np.std(c3), c3_diff)))\nc4_train = np.concatenate((np.array(c4), np.random.normal(np.mean(c4), np.std(c4), c4_diff)))\nc5_train = np.concatenate((np.array(c5), np.random.normal(np.mean(c5), np.std(c5), c5_diff)))\nc6_train = np.concatenate((np.array(c6), np.random.normal(np.mean(c6), np.std(c6), c6_diff)))\nc7_train = np.concatenate((np.array(c7), np.random.normal(np.mean(c7), np.std(c7), c7_diff)))\n\nc1_train=c1","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we want to transform our data in ways we can help our model distinguish between each vertebrae's distance to the top of the scan","metadata":{}},{"cell_type":"markdown","source":"Below we will shift the data on the X axis behind 0 *(0 - x)*, then multiply all the data points to the power of 3","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\n\n\n\n\ndef transform_list(list):\n    list = [abs(0-c) for c in list]\n    \n    \n    return [c**3 for c in list]\n\n\n\nc1_t1 = transform_list(c1_train)\nc2_t1 = transform_list(c2_train)\nc3_t1 = transform_list(c3_train)\nc4_t1 = transform_list(c4_train)\nc5_t1 = transform_list(c5_train)\nc6_t1 = transform_list(c6_train)\nc7_t1 = transform_list(c7_train)\n\n\nsns.histplot(x=c1_t1,label='C1', color='red')\n\nsns.histplot(x=c2_t1,label='C2', color='pink')\n\nsns.histplot(x=c3_t1,label='C3', color='blue')\n\nsns.histplot(x=c4_t1,label='C4', color='yellow')\n\nsns.histplot(x=c5_t1,label='C5', color='green')\n\nsns.histplot(x=c6_t1,label='C6', color='purple')\n\nsns.histplot(x=c7_t1,label='C7', color='black')\n\nplt.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we will subtract all data points from the max distance value *(200 - x)* and multiply the data points to the power of 3","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\n\n\n\n\ndef transform_list(list):\n    list = [abs(200-c) for c in list]\n    \n    return [c**3 for c in list]\n\n\n\nc1_t = transform_list(c1_train)\nc2_t = transform_list(c2_train)\nc3_t = transform_list(c3_train)\nc4_t = transform_list(c4_train)\nc5_t = transform_list(c5_train)\nc6_t = transform_list(c6_train)\nc7_t = transform_list(c7_train)\n\n\nsns.histplot(x=c1_t,label='C1', color='red')\n\nsns.histplot(x=c2_t,label='C2', color='pink')\n\nsns.histplot(x=c3_t,label='C3', color='blue')\n\nsns.histplot(x=c4_t,label='C4', color='yellow')\n\nsns.histplot(x=c5_t,label='C5', color='green')\n\nsns.histplot(x=c6_t,label='C6', color='purple')\n\nsns.histplot(x=c7_t,label='C7', color='black')\n\nplt.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we append all these features and labels","metadata":{}},{"cell_type":"code","source":"from sklearn import svm\n\nX = []\ny = []\n\n\nfor i, c in enumerate(c1_train):\n    X.append([c , c1_t[i] , c1_t1[i]])\n    y.append(1)     \nfor i, c in enumerate(c2_train):\n    X.append([c , c2_t[i] , c2_t1[i]])\n    y.append(2)  \nfor i, c in enumerate(c3_train):\n    X.append([c , c3_t[i] , c3_t1[i]])\n    y.append(3)  \nfor i, c in enumerate(c4_train):\n    X.append([c , c4_t[i] , c4_t1[i]])\n    y.append(4)  \nfor i, c in enumerate(c5_train):\n    X.append([c , c5_t[i] , c5_t1[i]])\n    y.append(5)  \nfor i, c in enumerate(c6_train):\n    X.append([c , c6_t[i] , c6_t1[i]])\n    y.append(6)  \nfor i, c in enumerate(c7_train):\n    X.append([c , c7_t[i] , c7_t1[i]])\n    y.append(7)   \n    \n    \n    \n\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And check the distribution of classes again","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.histplot(x=y)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Building the Model","metadata":{}},{"cell_type":"markdown","source":"Creating train and test dataset","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\n\nX = np.array(X)\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=0)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Building an SVC One-vs-All classifier with the data","metadata":{}},{"cell_type":"code","source":"\nfrom numpy import mean\nfrom numpy import std\nfrom sklearn.experimental import enable_hist_gradient_boosting\nfrom sklearn.ensemble import HistGradientBoostingClassifier\nfrom sklearn.model_selection import RepeatedStratifiedKFold\nfrom sklearn.model_selection import cross_val_score\nimport lightgbm as lgb\n\n\nc1_classifier = svm.SVC(decision_function_shape='ovo', probability=True, C=0.1,tol=1e-4)#One vs Rest\n\n\nc1_classifier.fit(X_train, y_train)\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Checking the ROC score for each label classification","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import roc_auc_score\n\n\n\n\n\nprobs = c1_classifier.predict_proba(X_test)  \n\nfor i in range(probs.shape[1]):\n    new_y=[]\n    for y in y_test:\n        if y==i+1: new_y.append(1)\n        else: new_y.append(0)\n    \n    roc_auc = roc_auc_score(new_y,probs[:,i])\n    print(f'c{i+1} roc auc score:  {roc_auc}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And here we have the overall ROC score for the model","metadata":{}},{"cell_type":"code","source":"probs = c1_classifier.predict_proba(X_test)  \nroc_auc = roc_auc_score(y_test,probs,multi_class='ovo')\nroc_auc","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Predicted X Real labels","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import plot_confusion_matrix\nplot_confusion_matrix(c1_classifier, X_test, y_test)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In this code, we grab the probabilistic output (p) of the model for each class. Then we create a threshold that starts at 0.02. Every output **p** below the threshold is classified as 0 or False, every output **p** above the threshold is classified as 1 or True.\n\nThen, we loop this code, adding to the threshold, and plot the ROC score for each threshold, for each class","metadata":{}},{"cell_type":"code","source":"thresholds=[]\nscores=[]\n\nfor i in range(probs.shape[1]):\n    \n    threshold=0.02\n    class_scores = []\n    class_thresholds=[]\n    for _ in range(50):\n        \n        new_y=[]\n        for y in y_test:\n            if y==i+1: new_y.append(1)\n            else: new_y.append(0)\n        \n        new_x = []\n        for x in probs[:,i]:\n            if x>threshold: new_x.append(1)\n            else: new_x.append(0)\n        \n        \n        roc_auc = roc_auc_score(new_y,new_x)\n        class_scores.append(roc_auc)\n        class_thresholds.append(threshold)\n        threshold+=0.02\n        \n    scores.append(class_scores)\n    thresholds.append(class_thresholds)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This way, we can see that the model is best at classifying the C1 and C7 vertebrae, all the rest lose a lot of performance as the threshold goes up, indicating the model doesn't have much \"confidence\" predicting them","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.lineplot(y=scores[0],x=thresholds[0], label='ROC AUC x C1 thresholds')\nsns.lineplot(y=scores[1],x=thresholds[1], label='ROC AUC x C2 thresholds')\nsns.lineplot(y=scores[2],x=thresholds[2], label='ROC AUC x C3 thresholds')\nsns.lineplot(y=scores[3],x=thresholds[3], label='ROC AUC x C4 thresholds')\nsns.lineplot(y=scores[4],x=thresholds[4], label='ROC AUC x C5 thresholds')\nsns.lineplot(y=scores[5],x=thresholds[5], label='ROC AUC x C6 thresholds')\nsns.lineplot(y=scores[6],x=thresholds[6], label='ROC AUC x C7 thresholds')\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Saving the model for later use","metadata":{}},{"cell_type":"code","source":"import pickle\nfilename = 'Z_axis_classifier.sav'\npickle.dump(c1_classifier, open(filename, 'wb'))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## MLP","metadata":{}},{"cell_type":"markdown","source":"We will try using a Multi-layer perceptron and see how it does on the data","metadata":{}},{"cell_type":"code","source":"import tensorflow\nimport keras","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import Dense, Dropout, BatchNormalization\n\nmodel = Sequential()\nmodel.add(keras.Input(shape=(3)))\nmodel.add(Dense(9, activation='relu'))\nmodel.add(BatchNormalization())\nmodel.add(Dense(81, activation='relu'))\nmodel.add(Dropout(0.3))\nmodel.add(BatchNormalization())\nmodel.add(Dense(243, activation='relu'))\nmodel.add(Dropout(0.3))\nmodel.add(BatchNormalization())\nmodel.add(Dense(7, activation='softmax'))\nmodel.compile(optimizer='adam', loss='categorical_crossentropy', metrics=['AUC','accuracy'])","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.summary()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = []\ny = []\n\n\nfor i, c in enumerate(c1_train):\n    X.append([c , c1_t[i] , c1_t1[i]])\n    y.append(0)     \nfor i, c in enumerate(c2_train):\n    X.append([c , c2_t[i] , c2_t1[i]])\n    y.append(1)  \nfor i, c in enumerate(c3_train):\n    X.append([c , c3_t[i] , c3_t1[i]])\n    y.append(2)  \nfor i, c in enumerate(c4_train):\n    X.append([c , c4_t[i] , c4_t1[i]])\n    y.append(3)  \nfor i, c in enumerate(c5_train):\n    X.append([c , c5_t[i] , c5_t1[i]])\n    y.append(4)  \nfor i, c in enumerate(c6_train):\n    X.append([c , c6_t[i] , c6_t1[i]])\n    y.append(5)  \nfor i, c in enumerate(c7_train):\n    X.append([c , c7_t[i] , c7_t1[i]])\n    y.append(6)   \n    \n    ","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tensorflow.keras.utils import to_categorical\n\ny=to_categorical(y)\n\n\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=0)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train = np.array(X_train)\nX_test = np.array(X_test)\ny_train = np.array(y_train)\ny_test = np.array(y_test)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.fit(X_train, y_train, validation_data=(X_test,y_test), epochs=15)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Lets check the performance metrics","metadata":{}},{"cell_type":"code","source":"probs = model.predict(X_test)  \n\n\n    \n    \n    \nroc_auc = roc_auc_score(y_test,probs)\nprint(f'roc auc score:  {roc_auc}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.metrics import confusion_matrix\n\npred = model.predict(X_test)\n\ncm = confusion_matrix(y_test.argmax(axis=1), pred.argmax(axis=1))\n\n\nax= plt.subplot()\nsns.heatmap(cm, annot=True, fmt='g', ax=ax);  #annot=True to annotate cells, ftm='g' to disable scientific notation\n\n\n# labels, title and ticks\nax.set_xlabel('Predicted labels')\nax.set_ylabel('True labels')\nax.set_title('Confusion Matrix')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Saving model for later","metadata":{}},{"cell_type":"code","source":"model.save('Z_axis_perceptron.h5')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now what we want is to use both models to see if the classification performance increases","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import confusion_matrix\n\npred = model.predict(X_test)\npred += c1_classifier.predict_proba(X_test)\n\ncm = confusion_matrix(y_test.argmax(axis=1), pred.argmax(axis=1))\n\n\nax= plt.subplot()\nsns.heatmap(cm, annot=True, fmt='g', ax=ax);  #annot=True to annotate cells, ftm='g' to disable scientific notation\n\n\n# labels, title and ticks\nax.set_xlabel('Predicted labels')\nax.set_ylabel('True labels')\nax.set_title('Confusion Matrix')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we repeat the threshold process","metadata":{}},{"cell_type":"code","source":"thresholds=[]\nscores=[]\n\nfor i in range(pred.shape[1]):\n    \n    threshold=0.02\n    class_scores = []\n    class_thresholds=[]\n    for _ in range(50):\n        \n        new_y=[]\n        for y in y_test[:,i]:\n            \n            if y==1: new_y.append(1)\n            else: new_y.append(0)\n        \n        new_x = []\n        for x in probs[:,i]:\n            if x>threshold: new_x.append(1)\n            else: new_x.append(0)\n        \n        \n        roc_auc = roc_auc_score(new_y,new_x)\n        class_scores.append(roc_auc)\n        class_thresholds.append(threshold)\n        threshold+=0.02\n        \n    scores.append(class_scores)\n    thresholds.append(class_thresholds)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nsns.lineplot(y=scores[0],x=thresholds[0], label='ROC AUC x C1 thresholds')\nsns.lineplot(y=scores[1],x=thresholds[1], label='ROC AUC x C2 thresholds')\nsns.lineplot(y=scores[2],x=thresholds[2], label='ROC AUC x C3 thresholds')\nsns.lineplot(y=scores[3],x=thresholds[3], label='ROC AUC x C4 thresholds')\nsns.lineplot(y=scores[4],x=thresholds[4], label='ROC AUC x C5 thresholds')\nsns.lineplot(y=scores[5],x=thresholds[5], label='ROC AUC x C6 thresholds')\nsns.lineplot(y=scores[6],x=thresholds[6], label='ROC AUC x C7 thresholds')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now lets just show a simple prediction of the model","metadata":{}},{"cell_type":"code","source":"def transform_list(list):\n    list = [abs(200-c) for c in list]\n    \n    return [c**3 for c in list]\n\ndef transform_list1(list):\n    list = [abs(0-c) for c in list]\n    \n    \n    return [c**3 for c in list]\n\n\nz_distance = [100]\nz_t = transform_list(z_distance)\nz_t1 = transform_list1(z_distance)\n\nx = np.array([z_distance, z_t, z_t1]).reshape(1,-1)\npred = c1_classifier.predict_proba(x)\npred += model.predict(x)\n\nprint('Predicted vertebrae 100mm from top Z position:  C'+str(pred.argmax()+1))\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Stacking slices\n\nHere we will classify each of the vertebrae images and piece each vertebrae in a 3D image.\n\nWe will be adding adding neighbor vertebraes to our stacked images for 2 reasons:\n\n- Some images contain parts of a neighbor vertebrae in them, even though the image may focus on a specific vertebrae\n\n- Not stacking neighbor vertebrae leads to small 3D images and inbalance in Z axis size.\n\nThere is a correlation between vertebraes being fractured and the neighbor vertebrae being fractured as well, so I believe this approach is justified.","metadata":{}},{"cell_type":"code","source":"c1_imgs = []\nc2_imgs = []\nc3_imgs = []\nc4_imgs = []\nc5_imgs = []\nc6_imgs = []\nc7_imgs = []\n\n\npath = f'../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.10001'\npatient = load_scan(path)\nzmax = patient[0][0x0020, 0x0032][-1]\nimgs = get_pixels_hu(patient)\n\nfor i, _ in enumerate(patient):\n\n\n\n    data = patient[i]\n    \n    z = np.absolute(zmax - data[0x0020, 0x0032][-1])\n    z_t = transform_list([z])\n    z_t1 = transform_list1([z])\n    \n    x = np.array([z_distance, z_t, z_t1]).reshape(1,-1)\n    \n    pred = c1_classifier.predict(x)\n    \n    if 1 in pred:\n        c1_imgs.append(imgs[i])\n        c2_imgs.append(imgs[i])\n        \n        \n    if 2 in pred:\n        c2_imgs.append(imgs[i])\n        c3_imgs.append(imgs[i])\n        \n    if 3 in pred:\n        c3_imgs.append(imgs[i])\n        c4_imgs.append(imgs[i])\n        \n    if 4 in pred:\n        c4_imgs.append(imgs[i])\n        c5_imgs.append(imgs[i])\n        \n    if 5 in pred:\n        c5_imgs.append(imgs[i])\n        c6_imgs.append(imgs[i])\n        \n    if 6 in pred:\n        c6_imgs.append(imgs[i])\n        c7_imgs.append(imgs[i])\n\n    if 7 in pred:\n        c7_imgs.append(imgs[i])\n    \n        ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Function to Stack Images","metadata":{}},{"cell_type":"code","source":"import cv2\n\ndef stack(imgs):\n    thresh_imgs = []\n    for img in imgs:\n        _,img = cv2.threshold(img,400,img.max(),cv2.THRESH_TOZERO)\n        img=np.expand_dims(img, axis=0)\n        thresh_imgs.append(img)\n        del img\n\n\n    stack = np.vstack(thresh_imgs)\n    del thresh_imgs\n    return stack\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Stacking images of each vertebrae","metadata":{}},{"cell_type":"code","source":"c1_3D = stack(c1_imgs)\nc2_3D = stack(c2_imgs)\nc3_3D = stack(c3_imgs)\nc4_3D = stack(c4_imgs)\nc5_3D = stack(c5_imgs)\nc6_3D = stack(c6_imgs)\nc7_3D = stack(c7_imgs)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Resampling the slices","metadata":{}},{"cell_type":"code","source":"c1_3D,_ = resample(c1_3D, patient, new_spacing=0.5)\nc2_3D,_ = resample(c2_3D, patient, new_spacing=0.5)\nc3_3D,_  = resample(c3_3D, patient, new_spacing=0.5)\nc4_3D,_  = resample(c4_3D, patient, new_spacing=0.5)\nc5_3D,_  = resample(c5_3D, patient, new_spacing=0.5)\nc6_3D,_  = resample(c6_3D, patient, new_spacing=0.5)\nc7_3D,_  = resample(c7_3D, patient, new_spacing=0.5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plotting","metadata":{}},{"cell_type":"markdown","source":"## C1","metadata":{}},{"cell_type":"code","source":"\nv, f = make_mesh(c1_3D, 400,1)\nplt_3d(v,f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## C2","metadata":{}},{"cell_type":"code","source":"v, f = make_mesh(c2_3D, 300,1)\nplt_3d(v,f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## C3","metadata":{}},{"cell_type":"code","source":"v, f = make_mesh(c3_3D, 300,1)\nplt_3d(v,f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## C4","metadata":{}},{"cell_type":"code","source":"v, f = make_mesh(c4_3D, 300,1)\nplt_3d(v,f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## C5","metadata":{}},{"cell_type":"code","source":"v, f = make_mesh(c5_3D, 300,1)\nplt_3d(v,f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## C6","metadata":{}},{"cell_type":"code","source":"v, f = make_mesh(c6_3D, 300,1)\nplt_3d(v,f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## C7","metadata":{}},{"cell_type":"code","source":"v, f = make_mesh(c7_3D, 300,1)\nplt_3d(v,f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Conclusion\n\nThis classification model is not all great, but should be very fitting working with a image classification model to determine which cervical vertebrae is presented in a certain image. The normalization process for the images made in this notebook should also help further in this competition","metadata":{}}]}