{"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":"### Yet Another Notebook with Image Preprocessing :)\n\nThe basic idea comes from this:\n\n- We know dicom files has wide range of pixel values and while converting these images into model friendly image types such as \".jpg, .png\" we lose some amount of information; pixel values gets squashed between 0-255 range.\n\n- The problem above can be solved by using windowing ranges, there are many notebooks in previous medical imaging competitions using this technique. For different kind of tissues (soft, muscle, fat etc) we can adjust our window_center and window_with values to keep important signals in our images. Or you can use built-in method [David Roberts's notebook](https://www.kaggle.com/code/davidbroberts/mammography-apply-windowing/notebook).\n\n- I noticed that some of the default windowing ranges resulting undesirable images (at least for human eyes). So I decided to use manual windowing, but it seems like every machine_id has it's own windowing ranges registered into dicom files, altering these values gave better results on some machines. But we do not know which value is best for our models.\n\n- Then comes into last part of our preprocessing, converting dicoms resulting grayscale images with single channel, but with convolutional networks we can have more than one channels per image (usually 3 for pretrained models. I find myself asking **why don't we use this extra information to our benefit?** We can get different windowing ranges and stack them together into 3 channel images and feed them into our models...","metadata":{}},{"cell_type":"markdown","source":"## Dataset Link:\n\nYou can find the different sized versions of this approach created exactly as this notebook here:\n\n- 512x:\n    - https://www.kaggle.com/datasets/datafan07/multichannel-images-with-different-windows\n    \n    \n\n- 256x:\n    - https://www.kaggle.com/datasets/datafan07/multichannel-images-with-different-windows-256x","metadata":{}},{"cell_type":"code","source":"!pip install -qU python-gdcm pydicom pylibjpeg","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-12-03T16:34:09.596691Z","iopub.execute_input":"2022-12-03T16:34:09.597469Z","iopub.status.idle":"2022-12-03T16:34:23.863654Z","shell.execute_reply.started":"2022-12-03T16:34:09.597349Z","shell.execute_reply":"2022-12-03T16:34:23.862925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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 os\n\nimport glob\nimport cv2\nimport gdcm\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut, apply_modality_lut\n\nimport matplotlib.pyplot as plt","metadata":{"execution":{"iopub.status.busy":"2022-12-03T16:34:23.86532Z","iopub.execute_input":"2022-12-03T16:34:23.865684Z","iopub.status.idle":"2022-12-03T16:34:24.229799Z","shell.execute_reply.started":"2022-12-03T16:34:23.865657Z","shell.execute_reply":"2022-12-03T16:34:24.228701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# read training file\ntrain_df = pd.read_csv('/kaggle/input/rsna-breast-cancer-detection/train.csv')","metadata":{"execution":{"iopub.status.busy":"2022-12-03T16:34:24.231141Z","iopub.execute_input":"2022-12-03T16:34:24.231458Z","iopub.status.idle":"2022-12-03T16:34:24.317694Z","shell.execute_reply.started":"2022-12-03T16:34:24.231429Z","shell.execute_reply":"2022-12-03T16:34:24.316907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# create new column for dicom paths\ntrain_df['path'] ='/kaggle/input/rsna-breast-cancer-detection/train_images/'+ train_df['patient_id'].astype('str') +'/'+train_df['image_id'].astype('str') + '.dcm'\n","metadata":{"execution":{"iopub.status.busy":"2022-12-03T16:34:24.319228Z","iopub.execute_input":"2022-12-03T16:34:24.320196Z","iopub.status.idle":"2022-12-03T16:34:24.401806Z","shell.execute_reply.started":"2022-12-03T16:34:24.320171Z","shell.execute_reply":"2022-12-03T16:34:24.400813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# sample of metadata included in dicom file\npydicom.dcmread(train_df['path'][0])","metadata":{"execution":{"iopub.status.busy":"2022-12-03T16:34:24.402973Z","iopub.execute_input":"2022-12-03T16:34:24.403241Z","iopub.status.idle":"2022-12-03T16:34:24.479692Z","shell.execute_reply.started":"2022-12-03T16:34:24.403218Z","shell.execute_reply":"2022-12-03T16:34:24.478908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def window_image(img, window_center,window_width, intercept, slope, dcm, PhotometricInterpretation = \"MONOCHROME1\"):   \n    # copied/edited from https://www.kaggle.com/code/omission/eda-view-dicom-images-with-correct-windowing/notebook\n    \n    img = (img*slope +intercept) #for translation adjustments given in the dicom file. \n    img_min = window_center - window_width//2 #minimum HU level\n    img_max = window_center + window_width//2 #maximum HU level\n    img[img<img_min] = img_min #set img_min for all HU levels less than minimum HU level\n    img[img>img_max] = img_max #set img_max for all HU levels higher than maximum HU level\n    \n    if np.max(img) != 0:\n        img = (img - img_min) / (img_max - img_min)\n    img=(img * 255)\n    \n    if PhotometricInterpretation == \"MONOCHROME1\":\n        img = np.amax(img) - img\n        \n    else:\n        img = img - np.min(img)\n    \n    return img.astype('uint8')\n\ndef get_mean_of_dicom_field_as_int(x):\n    if type(x) == pydicom.multival.MultiValue: return int(np.mean(x))\n    else: return int(x)\n    \ndef get_windowing(data):\n    try:\n        dicom_fields = [data[('0028','1050')].value, #window center\n                        data[('0028','1051')].value, #window width\n                        data[('0028','1052')].value, #intercept\n                        data[('0028','1053')].value] #slope\n        return [get_mean_of_dicom_field_as_int(x) for x in dicom_fields]\n    except:\n        warnings.warn(\"Couldn't find dicom metadata for file, using default values.\")\n        return [1500.0,1000.0,0.0,1.0]","metadata":{"execution":{"iopub.status.busy":"2022-12-03T16:34:24.480703Z","iopub.execute_input":"2022-12-03T16:34:24.481092Z","iopub.status.idle":"2022-12-03T16:34:24.490381Z","shell.execute_reply.started":"2022-12-03T16:34:24.481068Z","shell.execute_reply":"2022-12-03T16:34:24.489299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def convert_dcm(df, resize=False, size=256, save=False, save_pth='/kaggle/working/', extension='png', w_range=(0.05, 0.10)):\n    dicom_raw = pydicom.dcmread(df['path'])\n    arr = dicom_raw.pixel_array\n    \n    window_center , window_width, intercept, slope = get_windowing(dicom_raw)\n    \n\n    X = window_image(arr, window_center, window_width, intercept, slope, dicom_raw, dicom_raw.PhotometricInterpretation)\n    Y = window_image(arr, window_center-window_center*w_range[0], window_width-window_width*w_range[1], intercept, slope,dicom_raw,  dicom_raw.PhotometricInterpretation)\n    Z = window_image(arr, window_center+window_center*w_range[0], window_width+window_width*w_range[1], intercept, slope,dicom_raw,  dicom_raw.PhotometricInterpretation)\n    \n    \n    im = np.concatenate([np.expand_dims(X, -1), np.expand_dims(Y, -1), np.expand_dims(Z, -1)], axis=-1)\n    \n    if resize:\n        im = cv2.resize(im, (size, size))    \n        \n    if save:\n        os.makedirs(save_pth + 'train_images/' +  df['patient_id'].astype('str') +'/', exist_ok=True)\n        path = save_pth + 'train_images/' +  df['patient_id'].astype('str') +'/' + df['image_id'].astype('str') +f'.{extension}'\n        \n \n        cv2.imwrite(path, im)\n    \n    return im","metadata":{"execution":{"iopub.status.busy":"2022-12-03T16:34:24.49153Z","iopub.execute_input":"2022-12-03T16:34:24.491836Z","iopub.status.idle":"2022-12-03T16:34:24.509112Z","shell.execute_reply.started":"2022-12-03T16:34:24.491809Z","shell.execute_reply":"2022-12-03T16:34:24.508416Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = train_df.sample(100, random_state=456).reset_index(drop=True)\n\nfig, axes = plt.subplots(5, 4,sharex=False, sharey=True, figsize=(32,32))\nfor i in range(5):\n    dcm = convert_dcm(df.iloc[i,:], resize=True, size=512, save=False, save_pth='/kaggle/working/', extension='png', w_range=(0.1, 0.25))\n    axes[i,0].imshow(dcm[:,:,0], cmap='gray')\n    axes[i,1].imshow(dcm[:,:,1], cmap='gray')\n    axes[i,2].imshow(dcm[:,:,-1], cmap='gray')\n    axes[i,3].imshow(dcm[:,:,:])\n    for j in range(4):\n        axes[i,j].axis('off')\n        axes[i,j].set_autoscale_on(False)\ncols=['Original Window', 'Lower Window', 'Upper Window', 'Stacked Windows']    \nfor ax, col in zip(axes[0], cols):\n    ax.set_title(col, fontsize=32)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-12-03T16:34:24.51023Z","iopub.execute_input":"2022-12-03T16:34:24.510667Z","iopub.status.idle":"2022-12-03T16:34:30.538629Z","shell.execute_reply.started":"2022-12-03T16:34:24.510642Z","shell.execute_reply":"2022-12-03T16:34:30.537618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## References\n\n* https://www.kaggle.com/code/dcstang/see-like-a-radiologist-with-systematic-windowing/notebook\n* https://www.kaggle.com/code/omission/eda-view-dicom-images-with-correct-windowing/notebook\n* https://www.kaggle.com/code/redwankarimsony/ct-scans-dicom-files-windowing-explained\n* https://www.kaggle.com/code/davidbroberts/mammography-apply-windowing\n* https://radiopaedia.org/articles/windowing-ct","metadata":{}},{"cell_type":"markdown","source":"## Ideas to Improve\n\n- I just take upper and lower bounds intuitively (some percent up and low ranges), you probably can get better hyperparameters than me easily.\n- You can increase number of channels with different windows.\n- You can apply different filters as new feature channels such as median blur etc.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}