{"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":"<div class='alert alert-info' style='text-align: center'><h1>Cropping Chest X-Rays</h1>\n    - yet another chest x-ray processing notebook -\n</div>\n\n### Simple CXR cropping technique\n\n- The idea here is to easily remove noise around the edges of a chest x-ray without using segmentation.\n- Since we know the edges of a CXR are generally very dark compared to the areas with tissue in them, we can remove pixel rows that are below a certain threshold and crop the image accordingly.\n- This is a rough idea and would require some tuning to work in production, but it could produce some very usable results.\n\n#### Let's grab an image, resize it and scale it to 8 bits (for simplicity), then attempt to crop down to the lung fields.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pydicom\nimport cv2\nimport matplotlib.pyplot as plt\nfrom skimage.transform import resize","metadata":{"execution":{"iopub.status.busy":"2021-07-02T20:46:13.272798Z","iopub.execute_input":"2021-07-02T20:46:13.273326Z","iopub.status.idle":"2021-07-02T20:46:14.872818Z","shell.execute_reply.started":"2021-07-02T20:46:13.273227Z","shell.execute_reply":"2021-07-02T20:46:14.871812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Grab a random DICOM file from the SIIM Covid-19 Detection set\nimg_file = \"../input/siim-covid19-detection/train/00792b5c8852/1f52bcb3143e/3fadf4b48db3.dcm\"\nimg = pydicom.dcmread(img_file)","metadata":{"execution":{"iopub.status.busy":"2021-07-02T20:49:30.569318Z","iopub.execute_input":"2021-07-02T20:49:30.56966Z","iopub.status.idle":"2021-07-02T20:49:30.582381Z","shell.execute_reply.started":"2021-07-02T20:49:30.569632Z","shell.execute_reply":"2021-07-02T20:49:30.581032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Resize the pixels\nw = int(img.pixel_array.shape[0] * .25)\nh = int(img.pixel_array.shape[1] * .25)\n\npx = img.pixel_array / 255\nimg = resize(px, (w, h), anti_aliasing=True).astype(float)\n\n# scale the pixels\nimg = (np.maximum(img,0) / img.max()) * 255.0\nimg = np.uint8(img)\n\n# TODO: Invert MONOCHROME1 images here.","metadata":{"execution":{"iopub.status.busy":"2021-07-02T20:49:32.611443Z","iopub.execute_input":"2021-07-02T20:49:32.6118Z","iopub.status.idle":"2021-07-02T20:49:32.878178Z","shell.execute_reply.started":"2021-07-02T20:49:32.611767Z","shell.execute_reply":"2021-07-02T20:49:32.877165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Display the original image\nplt.figure(figsize=(15,5))\nplt.imshow(img,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-07-02T20:49:35.505694Z","iopub.execute_input":"2021-07-02T20:49:35.506037Z","iopub.status.idle":"2021-07-02T20:49:35.867091Z","shell.execute_reply.started":"2021-07-02T20:49:35.506009Z","shell.execute_reply":"2021-07-02T20:49:35.866344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Make a binarized copy of the image\nthresh = 150\nimg_bin = cv2.threshold(img, thresh, 255, cv2.THRESH_BINARY)[1]\n\nplt.figure(figsize=(15,5))\nplt.imshow(img_bin,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-06-20T16:09:34.55568Z","iopub.execute_input":"2021-06-20T16:09:34.555986Z","iopub.status.idle":"2021-06-20T16:09:34.77066Z","shell.execute_reply.started":"2021-06-20T16:09:34.555957Z","shell.execute_reply":"2021-06-20T16:09:34.769717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Flip the image 90 degrees\nimg_bin = cv2.rotate(img_bin, cv2.cv2.ROTATE_90_CLOCKWISE)\nplt.figure(figsize=(15,5))\nplt.imshow(img_bin,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-06-20T16:09:34.771988Z","iopub.execute_input":"2021-06-20T16:09:34.77257Z","iopub.status.idle":"2021-06-20T16:09:34.987972Z","shell.execute_reply.started":"2021-06-20T16:09:34.772524Z","shell.execute_reply":"2021-06-20T16:09:34.986902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- With the image flipped, we can extract rows and calculate their mean. \n- Any row that is 'too dark' is probably not part of the patient and can be cropped out.","metadata":{}},{"cell_type":"code","source":"right = 0;\nleft = 0;\nline_thickness = 2\n\n# This is the value that specifies how bright a row is to consider it 'not the edge (too bright)'\nintensity_threshold = 190\n\n# Start at the bottom and work upwards checking the mean of pixels in every 10th row, this is the right side of the image\nfor i in range(img_bin.shape[0]-1,0,-10):\n    row_mean = img_bin[i].mean()\n    if row_mean > intensity_threshold:\n        right = i\n        \n        # Draw a line where we want to crop\n        cv2.line(img_bin, (0, i), (img_bin.shape[1], i), (0, img_bin.shape[1], 0), thickness=line_thickness)\n        break\n        \n# Start at the top and go down to find the left side\nfor i in range(0,img_bin.shape[0]-1,10):\n    row_mean = img_bin[i].mean()\n    if row_mean > intensity_threshold:\n        left = i\n        \n        # Draw a line where we want to crop\n        cv2.line(img_bin, (0, i), (img_bin.shape[1], i), (0, img_bin.shape[1], 0), thickness=line_thickness)\n        break","metadata":{"execution":{"iopub.status.busy":"2021-06-20T16:09:34.989408Z","iopub.execute_input":"2021-06-20T16:09:34.989912Z","iopub.status.idle":"2021-06-20T16:09:35.000067Z","shell.execute_reply.started":"2021-06-20T16:09:34.989866Z","shell.execute_reply":"2021-06-20T16:09:34.999109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Draw lines on the image where the mean intensity is > intensity_threshold\nplt.figure(figsize=(15,5))\nplt.imshow(img_bin,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-06-20T16:09:35.001093Z","iopub.execute_input":"2021-06-20T16:09:35.001371Z","iopub.status.idle":"2021-06-20T16:09:35.22012Z","shell.execute_reply.started":"2021-06-20T16:09:35.001345Z","shell.execute_reply":"2021-06-20T16:09:35.219098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- The two lines on the image represent where the row became much 'brighter'. This is where we'll crop.","metadata":{}},{"cell_type":"code","source":"# Rotate the image back to it's normal orientation\nimg_bin = cv2.rotate(img_bin, cv2.cv2.ROTATE_90_COUNTERCLOCKWISE)\n\nx1 = left\ny1 = 0\nx2 = right\ny2 = img_bin.shape[1]\n\n# Grab the region we identified from the binarized image\nimg_cropped = img_bin[y1:y2, x1:x2]\nplt.figure(figsize=(15,5))\nplt.imshow(img_cropped,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-06-20T16:09:35.221883Z","iopub.execute_input":"2021-06-20T16:09:35.222293Z","iopub.status.idle":"2021-06-20T16:09:35.415965Z","shell.execute_reply.started":"2021-06-20T16:09:35.22225Z","shell.execute_reply":"2021-06-20T16:09:35.415266Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Now we'll do the same process to find the best place to crop the top and bottom.","metadata":{}},{"cell_type":"code","source":"top = 0;\nbottom = 0;\n\n# Set some threshold values to specify what we consider edge vs patient\nbright_threshold = 240\ndark_threshold = 100\n\n# Start at the bottom and work upward\nfor i in range(img_cropped.shape[0]-1,0,-10):\n    row_mean = img_cropped[i].mean()\n    if row_mean < bright_threshold:\n        # Add 100 pixels of padding so we don't cut the costophrenic angles off\n        bottom = i + 100\n        cv2.line(img_cropped, (0, bottom), (img_cropped.shape[1], bottom), (0, img_cropped.shape[1], 0), thickness=line_thickness)\n        break\n        \n# Start at the top and go down\nfor i in range(0,img_cropped.shape[0]-1,10):\n    row_mean = img_cropped[i].mean()\n    if row_mean > dark_threshold:\n        top = i\n        cv2.line(img_cropped, (0, i), (img_cropped.shape[1], i), (0, img_cropped.shape[1], 0), thickness=line_thickness)\n        break\n        \nplt.figure(figsize=(15,5))\nplt.imshow(img_cropped,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-06-20T16:09:35.417041Z","iopub.execute_input":"2021-06-20T16:09:35.417478Z","iopub.status.idle":"2021-06-20T16:09:35.61652Z","shell.execute_reply.started":"2021-06-20T16:09:35.417433Z","shell.execute_reply":"2021-06-20T16:09:35.615543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x1 = 0\ny1 = top\nx2 = img_bin.shape[0]\ny2 = bottom\n\nimg_cropped = img_cropped[y1:y2, x1:x2]\nplt.figure(figsize=(15,5))\nplt.imshow(img_cropped,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-06-20T16:09:35.617756Z","iopub.execute_input":"2021-06-20T16:09:35.618031Z","iopub.status.idle":"2021-06-20T16:09:35.798933Z","shell.execute_reply.started":"2021-06-20T16:09:35.618005Z","shell.execute_reply":"2021-06-20T16:09:35.797884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Display the original image and the cropped section\nimg_cropped = img[top:bottom, left:right]\nplt.figure(figsize=(15,5))\n\nplt.subplot(1, 2, 1)\nplt.imshow(img,cmap=\"gray\");\n\nplt.subplot(1, 2, 2)\nplt.imshow(img_cropped,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-06-20T16:09:35.800684Z","iopub.execute_input":"2021-06-20T16:09:35.801108Z","iopub.status.idle":"2021-06-20T16:09:36.218751Z","shell.execute_reply.started":"2021-06-20T16:09:35.801073Z","shell.execute_reply":"2021-06-20T16:09:36.217631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- This isn't a perfect method, but hopefully someone can expand on it.","metadata":{}},{"cell_type":"markdown","source":"**Here are some other processing notebooks I made:**\n- Lung Segmentation Without CNN -> https://www.kaggle.com/davidbroberts/lung-segmentation-without-cnn\n- Applying filters to x-rays -> https://www.kaggle.com/davidbroberts/applying-filters-to-chest-x-rays\n- Rib supression on Chest X-Rays -> https://www.kaggle.com/davidbroberts/rib-suppression-poc\n- Manual DICOM VOI LUT -> https://www.kaggle.com/davidbroberts/manual-dicom-voi-lut\n- Apply Unsharp Mask to Chest X-Rays -> https://www.kaggle.com/davidbroberts/unsharp-masking-chest-x-rays\n- Bounding Boxes on Cropped Images -> https://www.kaggle.com/davidbroberts/bounding-boxes-on-cropped-images\n- Visualizing Chest X-Ray bit planes -> https://www.kaggle.com/davidbroberts/visualizing-chest-x-ray-bitplanes\n- DICOM full range pixels as CNN input -> https://www.kaggle.com/davidbroberts/dicom-full-range-pixels-as-cnn-input\n- Standardizing Chest X-Ray Dataset Exports -> https://www.kaggle.com/davidbroberts/standardizing-cxr-datasets","metadata":{}}]}