{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":71549,"databundleVersionId":8561470,"sourceType":"competition"},{"sourceId":188023704,"sourceType":"kernelVersion"},{"sourceId":206417594,"sourceType":"kernelVersion"},{"sourceId":206408250,"sourceType":"kernelVersion"}],"dockerImageVersionId":30786,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"### Setting up Model and Dataset. \nThe dataset consists of about 2000*.8 images of patient's spine annotated for the L3-L4. 20% of the dataset was left of the model to valid it's performance on unseen images. The objective for the object detection model is to identify the l3-l4 region on the spine so that we can easily diagnosis spinal stenosis.\n\nThe exact training process can be found here:\nhttps://www.kaggle.com/code/junglebeastds/spinal-stenosis-detection-demo\nNote that model was not optimized at all for performance. So in practice, we can probably get training and detection much faster.","metadata":{}},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport sys\nimport os\nfrom matplotlib import pyplot as plt\nimport cv2 as cv\nfrom PIL import Image\nimport ctypes\nimport gc\n\ndef clean_memory():\n    gc.collect()\n    ctypes.CDLL(\"libc.so.6\").malloc_trim(0)\n    torch.cuda.empty_cache()\nimport transformers\nprint(transformers.__version__)\n\nimport torch\n\nfrom PIL import Image\nfrom transformers import RTDetrForObjectDetection, RTDetrImageProcessor,RTDetrConfig\nfrom transformers.image_transforms import center_to_corners_format\nfrom torch.utils.data import DataLoader\n\n\n#image_processor = RTDetrImageProcessor.from_pretrained(\"PekingU/rtdetr_r50vd_coco_o365\")\nimage_processor = RTDetrImageProcessor.from_pretrained(\"/kaggle/input/spinal-stenosis-detection-demo/image_processor\")\nconfig = RTDetrConfig.from_pretrained(\"/kaggle/input/spinal-stenosis-detection-demo/epoch9\")\nconfig.num_labels = 1\n# label_to_id = {k:v for k,v in zip(['left_neural_foraminal_narrowing_l1_l2',\n#        'left_neural_foraminal_narrowing_l2_l3',\n#        'left_neural_foraminal_narrowing_l3_l4',\n#        'left_neural_foraminal_narrowing_l4_l5',\n#        'left_neural_foraminal_narrowing_l5_s1',\n#        'left_subarticular_stenosis_l1_l2', 'left_subarticular_stenosis_l2_l3',\n#        'left_subarticular_stenosis_l3_l4', 'left_subarticular_stenosis_l4_l5',\n#        'left_subarticular_stenosis_l5_s1',\n#        'right_neural_foraminal_narrowing_l1_l2',\n#        'right_neural_foraminal_narrowing_l2_l3',\n#        'right_neural_foraminal_narrowing_l3_l4',\n#        'right_neural_foraminal_narrowing_l4_l5',\n#        'right_neural_foraminal_narrowing_l5_s1',\n#        'right_subarticular_stenosis_l1_l2',\n#        'right_subarticular_stenosis_l2_l3',\n#        'right_subarticular_stenosis_l3_l4',\n#        'right_subarticular_stenosis_l4_l5',\n#        'right_subarticular_stenosis_l5_s1', 'spinal_canal_stenosis_l1_l2',\n#        'spinal_canal_stenosis_l2_l3', 'spinal_canal_stenosis_l3_l4',\n#        'spinal_canal_stenosis_l4_l5', 'spinal_canal_stenosis_l5_s1'], range(25))}\n# id_to_label = {v:k for k,v in label_to_id.items()}\n# config.id2label = id_to_label\n# config.label2id = label_to_id\n\nmodel = RTDetrForObjectDetection.from_pretrained(\"/kaggle/input/spinal-stenosis-detection-demo/epoch9\", config = config,ignore_mismatched_sizes=True)\nsum(p.numel() for p in model.parameters() if p.requires_grad)\n\nimport pickle\nwith open('/kaggle/input/process-metadata/data.pkl', 'rb') as f:\n    data = pickle.load(f)\n#filter for only 'left_neural_foraminal_narrowing_l3_l4'\ndata_single_class = []\nfor pic in data:\n    if 2 in pic['labels']['class_labels']:\n        filt = pic['labels']['class_labels'] == 2\n        pic['labels']['class_labels'] = pic['labels']['class_labels'][filt]\n        pic['labels']['boxes'] = pic['labels']['boxes'][filt]\n        pic['classification'] = pic['classification'][filt]\n        pic['filename'] = pic['filename'].replace(\"val\",\"train\")\n        data_single_class.append(pic)\n\n#train, test split\ndata = data_single_class\nl = int(len(data)*.8)\ntrain = data[:l]\ntest = data[l:]","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2024-11-11T00:11:06.041243Z","iopub.execute_input":"2024-11-11T00:11:06.041647Z","iopub.status.idle":"2024-11-11T00:11:11.743443Z","shell.execute_reply.started":"2024-11-11T00:11:06.041605Z","shell.execute_reply":"2024-11-11T00:11:11.742199Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Let's take a look at the training data.\nIn each scan, an expert has made a diagnosis to the patient's L3-L4 veterbrate for spinal stenosis and annotated the location in the image that they used to make the diagnosis.\n\nIn the following images, the annotation is highlighted. It seems clear to me the location on the spine that the expert is looking for. Although less clear is how the diagnosis was made. But from a supervised leraning perspective, that could be easily done with a straight forward classification model after we localize the area of interest (the l3-l4 veterbrate). Obviously, it helps that the data is already labeled (even if it's just one label per series).","metadata":{}},{"cell_type":"code","source":"current_image = data[3] #change the indice to cycle through different images used to train the model\nimages = Image.open(current_image['filename']).convert(\"RGB\")\n\nplt.imshow(images, cmap=plt.cm.gray)\nplt.colorbar()\n\nplt.title(f'Normal Spine l3_l4')\nplt.xlabel('X-axis')\nplt.ylabel('Y-axis')\nx,y = current_image['labels']['boxes'][0]\nplt.plot(x,y, '.', markersize=50,fillstyle='none',label = f'localization')\ncurrent_image\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-11T01:53:08.948731Z","iopub.execute_input":"2024-11-11T01:53:08.949148Z","iopub.status.idle":"2024-11-11T01:53:09.436146Z","shell.execute_reply.started":"2024-11-11T01:53:08.949111Z","shell.execute_reply":"2024-11-11T01:53:09.435051Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Example of an MRI where the patient show severe spinal stenosis on L3_L4 veterbrate.","metadata":{}},{"cell_type":"code","source":"current_image = data[328]\nimages = Image.open(current_image['filename']).convert(\"RGB\")\n\nplt.imshow(images, cmap=plt.cm.gray)\nplt.colorbar()\nplt.title(f'Severe spinal stenosis l3_l4')\nplt.xlabel('X-axis')\nplt.ylabel('Y-axis')\nx,y = current_image['labels']['boxes'][0]\nplt.plot(x,y, '.', markersize=50,fillstyle='none',label = f'localization')\ncurrent_image","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-11T01:24:52.662089Z","iopub.execute_input":"2024-11-11T01:24:52.662543Z","iopub.status.idle":"2024-11-11T01:24:53.102205Z","shell.execute_reply.started":"2024-11-11T01:24:52.662502Z","shell.execute_reply":"2024-11-11T01:24:53.101006Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Let's look at an image that was annotated, but not included in the model's training.\nHow does the model perform on an image it has never seen. This image gives a clear view on the L3-l4 for spinal stenosis detection.","metadata":{"execution":{"iopub.status.busy":"2024-11-11T01:40:35.384129Z","iopub.execute_input":"2024-11-11T01:40:35.384546Z","iopub.status.idle":"2024-11-11T01:40:35.390155Z","shell.execute_reply.started":"2024-11-11T01:40:35.384510Z","shell.execute_reply":"2024-11-11T01:40:35.388823Z"}}},{"cell_type":"code","source":"t = test[69] #all images in the test dataset were not used to train the model.\nimage = Image.open(t['filename']).convert(\"RGB\")\ninputs = image_processor(images=image, return_tensors=\"pt\")\n\n\nwith torch.no_grad():\n    outputs = model(inputs['pixel_values'])\n\nresults = image_processor.post_process_object_detection(outputs, target_sizes=torch.tensor([image.size[::-1]]), threshold=0.4)\nfor result in results:\n    r = pd.DataFrame({'scores':result[\"scores\"].numpy(),'labels':result['labels'].numpy()})\n    #r = r.drop_duplicates('labels')\n    print('Top detection predictions')\n    for idx, row in r.iterrows():\n        score = row['scores']\n        label_id = row['labels']\n        box = result['boxes'][row.name]\n        score, label = score.item(), label_id.item()\n        box = [round(i, 2) for i in box.tolist()]\n        print(f\"label {idx} probability: {score:.2f}, box coordinates: {box}\")\n\nimage_data = image\nplt.imshow(image_data, cmap=plt.cm.gray)\nplt.colorbar()\nplt.title('Objective detection on validation set')\nplt.xlabel('X-axis')\nplt.ylabel('Y-axis')\nresult = results[0]\nfor idx, row in r.iloc[:100].iterrows():\n    score = row['scores']\n    label_id = row['labels']\n    box = result['boxes'][row.name]\n    score, label = score.item(), label_id.item()\n    box = [round(i, 2) for i in box.tolist()]\n    #print(f\"{model.config.id2label[label]}: {score:.2f} {box}\")\n\n    # Plot the DICOM image\n    # Calculate the corners of the square\n    #half_side = side_length / 2\n    enlarge_factor = 20\n    square_x = [box[0]-enlarge_factor, box[2]+enlarge_factor, box[2]+enlarge_factor, box[0]-enlarge_factor, box[0]-enlarge_factor]\n    square_y = [box[1]-enlarge_factor, box[1]-enlarge_factor, box[3]+enlarge_factor, box[3]+enlarge_factor, box[1]-enlarge_factor]\n\n    # Plot the square\n    \n    plt.plot(square_x, square_y, '-',label = f'prediction {idx}')\n    #plt.scatter([center_x], [center_y], color='blue')  # Mark the center point\n\n#actual\nfor label, xy in zip(t['labels']['class_labels'], t['labels']['boxes']):\n    plt.plot(xy[0], xy[1], '.',label = f'actual doctor anotation',\n             markersize=20,fillstyle='none' )    \n\nplt.legend(loc='upper left',bbox_to_anchor=(1.2, 1))\nplt.show()\nt","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-11T03:06:18.297958Z","iopub.execute_input":"2024-11-11T03:06:18.298387Z","iopub.status.idle":"2024-11-11T03:06:20.457150Z","shell.execute_reply.started":"2024-11-11T03:06:18.298349Z","shell.execute_reply":"2024-11-11T03:06:20.455864Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Now let's test our model on images that were not even annotated.","metadata":{"execution":{"iopub.status.busy":"2024-11-11T01:27:44.287103Z","iopub.execute_input":"2024-11-11T01:27:44.287532Z","iopub.status.idle":"2024-11-11T01:27:44.301696Z","shell.execute_reply.started":"2024-11-11T01:27:44.287496Z","shell.execute_reply":"2024-11-11T01:27:44.299281Z"}}},{"cell_type":"code","source":"meta = pd.read_csv(\"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_series_descriptions.csv\")\nmeta = meta[meta['series_description'] =='Sagittal T1'] #we are looking at ST1 view\nmeta","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-11T03:06:36.792373Z","iopub.execute_input":"2024-11-11T03:06:36.793252Z","iopub.status.idle":"2024-11-11T03:06:36.817806Z","shell.execute_reply.started":"2024-11-11T03:06:36.793206Z","shell.execute_reply":"2024-11-11T03:06:36.816740Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#let's look at the entire series of images from the MRI\nimport os\nimport pydicom\ndef read_dicom(dicom_path):\n    # Load DICOM file\n    dicom = pydicom.dcmread(dicom_path)\n    # Extract pixel array from DICOM file\n    pixel_array = dicom.pixel_array\n    # Normalize pixel values to the range [0, 255]\n    pixel_array = (pixel_array - np.min(pixel_array)) / (np.max(pixel_array) - np.min(pixel_array)) * 255\n    pixel_array = pixel_array.astype(np.uint8)\n    # Convert to a PIL image\n    image = Image.fromarray(pixel_array).convert(\"RGB\")\n    return image\n\n#selecting the entire series. The above prediction was one for the image on which the doctor made the annotation\n#we are now looking at the entire series from the MRI\ndirectory = f'/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/{t[\"study_id\"]}/{t[\"series_id\"]}' \nimages = os.listdir(directory)\nimages","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-11T03:06:50.291023Z","iopub.execute_input":"2024-11-11T03:06:50.292176Z","iopub.status.idle":"2024-11-11T03:06:50.314925Z","shell.execute_reply.started":"2024-11-11T03:06:50.292125Z","shell.execute_reply":"2024-11-11T03:06:50.313607Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Let's use the model on unlabel photos.\nThe expert only annotates one image per series. In practice, we want our model to analyze every image.\n\n\nThe model is used to segment each image in the series of MRI scans. In practice, we only select the images with highest probability.\n\nSome of the images without a great view, the model correctly does not detect any boxes at all. There are some false positives, however. Particularly around the l2-l3 region, likely due to how similar they look. For this particular problem, it's less problematic as spinal stenosis between the different veterbrates are highly correlated, but for some tasks that might be an issue.\n","metadata":{}},{"cell_type":"code","source":"for index in images:\n    print(index)\n    image = read_dicom(directory + '/' + index) \n    #saving the images in jpg before, distorted the images a bit. So we might get slightly numerical differences now that we are directly importing from dicom\n    inputs = image_processor(images=image, return_tensors=\"pt\")\n    \n    \n    with torch.no_grad():\n        outputs = model(inputs['pixel_values'])\n    \n    results = image_processor.post_process_object_detection(outputs, target_sizes=torch.tensor([image.size[::-1]]), threshold=0.4)\n    for result in results:\n        r = pd.DataFrame({'scores':result[\"scores\"].numpy(),'labels':result['labels'].numpy()})\n        #r = r.drop_duplicates('labels')\n        print('Top detection predictions')\n        for idx, row in r.iterrows():\n            score = row['scores']\n            label_id = row['labels']\n            box = result['boxes'][row.name]\n            score, label = score.item(), label_id.item()\n            box = [round(i, 2) for i in box.tolist()]\n            print(f\"label {idx} probability: {score:.2f}, box coordinates: {box}\")\n    \n    image_data = image\n    plt.imshow(image_data, cmap=plt.cm.gray)\n    plt.colorbar()\n    plt.title('Model to segment L3-L4 on the entire instance')\n    plt.xlabel('X-axis')\n    plt.ylabel('Y-axis')\n    result = results[0]\n    for idx, row in r.iloc[:100].iterrows():\n        score = row['scores']\n        label_id = row['labels']\n        box = result['boxes'][row.name]\n        score, label = score.item(), label_id.item()\n        box = [round(i, 2) for i in box.tolist()]\n        #print(f\"{model.config.id2label[label]}: {score:.2f} {box}\")\n    \n        # Plot the DICOM image\n        # Calculate the corners of the square\n        #half_side = side_length / 2\n        enlarge_factor = 20\n        square_x = [box[0]-enlarge_factor, box[2]+enlarge_factor, box[2]+enlarge_factor, box[0]-enlarge_factor, box[0]-enlarge_factor]\n        square_y = [box[1]-enlarge_factor, box[1]-enlarge_factor, box[3]+enlarge_factor, box[3]+enlarge_factor, box[1]-enlarge_factor]\n    \n        # Plot the square\n        \n        plt.plot(square_x, square_y, '-',label = f'prediction {idx}')\n        #plt.scatter([center_x], [center_y], color='blue')  # Mark the center point\n    \n    \n    plt.legend(loc='upper left',bbox_to_anchor=(1.2, 1))\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-11T03:06:55.075279Z","iopub.execute_input":"2024-11-11T03:06:55.075695Z","iopub.status.idle":"2024-11-11T03:07:26.562637Z","shell.execute_reply.started":"2024-11-11T03:06:55.075659Z","shell.execute_reply":"2024-11-11T03:07:26.561400Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}