{"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>Determining MR Slice Orientation</h1>\n    - yet another MR notebook -\n    </div>\n    \n#### The goal of this notebook is to demonstrate how an MR slice is orientated relative to the patient's body.\n\n#### We'll use the DICOM tag ImageOrientationPatient to get positioning markers and draw them onto slices.\n\n#### - The position markers can be a single letter or a combination of up to three letters ..\n- **A** = Anterior\n- **P** = Posterior\n- **H** = Head\n- **F** = Feet\n- **L** = Left\n- **R** = Right\n\n#### - Position markers that are multiple letters indicate obliqueness (rotation) toward one or more axes. \n#### - I.E .. A letter marker of \"LH\" indicates the primary axis is *Left*, but the image is rotated in the *Head* direction.\n\n#### - More info can be found in the Image Plane Module of the DICOM Standard -> http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html","metadata":{}},{"cell_type":"markdown","source":"### - Define some functions","metadata":{}},{"cell_type":"code","source":"import matplotlib\nimport cv2\nimport numpy as np\nimport pydicom\nimport matplotlib.pyplot as plt","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:33.952614Z","iopub.execute_input":"2021-09-16T15:38:33.95306Z","iopub.status.idle":"2021-09-16T15:38:34.608596Z","shell.execute_reply.started":"2021-09-16T15:38:33.952943Z","shell.execute_reply":"2021-09-16T15:38:34.607144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def normalize_pixels(pixels):\n    pixels = pixels - np.min(pixels)\n    pixels = pixels / np.max(pixels)\n    pixels = (pixels * 255).astype(np.uint8)\n    return pixels","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:34.610356Z","iopub.execute_input":"2021-09-16T15:38:34.61079Z","iopub.status.idle":"2021-09-16T15:38:34.618278Z","shell.execute_reply.started":"2021-09-16T15:38:34.610742Z","shell.execute_reply":"2021-09-16T15:38:34.616907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_marker(iop, i):\n    \n    # Get the primary axis orientation markers. Get IOP vector values starting at 'i'\n    # Use the first three values for right, and last three for bottom.\n    Xx = 'R' if iop[i] < 0 else 'L'\n    Yx = 'A' if iop[i+1] < 0 else 'P'\n    Zx = 'F' if iop[i+2] < 0 else 'H'\n    \n    or1 = \"\"\n    absX = abs(iop[i])\n    absY = abs(iop[i+1])\n    absZ = abs(iop[i+2])\n\n    # Add oblique markers if not orthogonal\n    for i in range(0,3):\n        if absX > .0001 and absX > absY and absX > absZ:\n            or1 += Xx\n            absX = 0\n        elif absY > .0001 and absY > absX and absY > absZ:\n            or1 += Yx\n            absY=0\n        elif absZ > .0001 and absZ > absX and absZ > absY:\n            or1 += Zx\n            absZ = 0\n        else:\n            break\n            \n    return or1","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:34.620355Z","iopub.execute_input":"2021-09-16T15:38:34.620752Z","iopub.status.idle":"2021-09-16T15:38:34.632453Z","shell.execute_reply.started":"2021-09-16T15:38:34.620711Z","shell.execute_reply":"2021-09-16T15:38:34.631167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_markers(iop):\n    marker_right = get_marker(iop, 0)\n    marker_bottom = get_marker(iop, 3)\n    return marker_right, marker_bottom","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:34.634743Z","iopub.execute_input":"2021-09-16T15:38:34.635209Z","iopub.status.idle":"2021-09-16T15:38:34.645976Z","shell.execute_reply.started":"2021-09-16T15:38:34.635167Z","shell.execute_reply":"2021-09-16T15:38:34.644625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def draw_markers(pixels, marker_r, marker_b):\n    color = (255,255,255)\n    font = cv2.FONT_HERSHEY_SIMPLEX\n    img2_pixels = cv2.putText(pixels, marker_r,(int(pixels.shape[1] - 65),int(pixels.shape[0] / 2)), font, 1, color, 2, cv2.LINE_AA)\n    img2_pixels = cv2.putText(pixels, marker_b,(int(pixels.shape[1] / 2 - 20),int(pixels.shape[0] - 25)), font, 1, color, 2, cv2.LINE_AA)\n    return pixels","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:34.647634Z","iopub.execute_input":"2021-09-16T15:38:34.648175Z","iopub.status.idle":"2021-09-16T15:38:34.660198Z","shell.execute_reply.started":"2021-09-16T15:38:34.648117Z","shell.execute_reply":"2021-09-16T15:38:34.658662Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get the relavent tags and calculate the coordinates that specify the 'box' that one image projects onto another\ndef create_reference_line(src, dst):\n    dst_iop = dst.ImageOrientationPatient\n    dst_ipp = dst.ImagePositionPatient\n    src_iop = src.ImageOrientationPatient\n    src_ipp = src.ImagePositionPatient  \n\n    pos_x = []\n    pos_y = []\n    pos_z = []\n    row_pixel = []\n    col_pixel = []\n\n    row_len = int(dst.Rows * dst.PixelSpacing[1])\n    col_len = int(dst.Columns * dst.PixelSpacing[0])\n    \n    # Get the coordinates of the box\n    pos_x.append(src_ipp[0]) # Top Left Corner\n    pos_y.append(src_ipp[1])\n    pos_z.append(src_ipp[2])\n    pos_x.append(src_ipp[0] + src_iop[0] * row_len) # Top Right Corner\n    pos_y.append(src_ipp[1] + src_iop[1] * row_len)\n    pos_z.append(src_ipp[2] + src_iop[2] * row_len)\n    pos_x.append(src_ipp[0] + src_iop[0] * row_len + src_iop[3] * col_len) # Bottom Right Corner\n    pos_y.append(src_ipp[1] + src_iop[1] * row_len + src_iop[4] * col_len)\n    pos_z.append(src_ipp[2] + src_iop[2] * row_len + src_iop[5] * col_len)\n    pos_x.append(src_ipp[0] + src_iop[3] * col_len) # Bottom Left Corner\n    pos_y.append(src_ipp[1] + src_iop[4] * col_len)\n    pos_z.append(src_ipp[2] + src_iop[5] * col_len)\n    \n    for i in range (0,3):\n    \n        pos_x[i] -= dst_ipp[0]\n        pos_y[i] -= dst_ipp[1]\n        pos_z[i] -= dst_ipp[2]\n        \n        # Rotate the coordinates in 3D space\n        cp = int(dst_iop[0] * pos_x[i] + dst_iop[1] * pos_y[i] + dst_iop[2] * pos_z[i])\n        rp = int(dst_iop[3] * pos_x[i] + dst_iop[4] * pos_y[i] + dst_iop[5] * pos_z[i])\n\n        col_pixel.append(int(cp / float(dst.PixelSpacing[0])))\n        row_pixel.append(int(rp / float(dst.PixelSpacing[1])))\n    \n    return col_pixel, row_pixel","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:34.662321Z","iopub.execute_input":"2021-09-16T15:38:34.662805Z","iopub.status.idle":"2021-09-16T15:38:34.681472Z","shell.execute_reply.started":"2021-09-16T15:38:34.66275Z","shell.execute_reply":"2021-09-16T15:38:34.68012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### - Open two images from the same study, in different planes.","metadata":{}},{"cell_type":"code","source":"# Read in two images from the same study, different series.\nimg1 = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00021/FLAIR/Image-140.dcm')\nimg2 = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00021/T2w/Image-223.dcm')\n   \n# Crunch the pixels down to 8 bit so we can draw an 8 bit graphics on them\nimg1_pixels = normalize_pixels(img1.pixel_array)\nimg2_pixels = normalize_pixels(img2.pixel_array)","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:34.683079Z","iopub.execute_input":"2021-09-16T15:38:34.683581Z","iopub.status.idle":"2021-09-16T15:38:34.74411Z","shell.execute_reply.started":"2021-09-16T15:38:34.683526Z","shell.execute_reply":"2021-09-16T15:38:34.743134Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axes = plt.subplots(nrows=1, ncols=2,sharex=True, sharey=True, figsize=(15, 5))\nax = axes.ravel()\nax[0].set_title(f'{img1.ImageOrientationPatient} - FLAIR Axial')\nax[0].imshow(img1.pixel_array, cmap='gray');\nax[1].set_title(f'{img2.ImageOrientationPatient} - T2w Sagittal')\nax[1].imshow(img2.pixel_array, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:34.746401Z","iopub.execute_input":"2021-09-16T15:38:34.747144Z","iopub.status.idle":"2021-09-16T15:38:35.194774Z","shell.execute_reply.started":"2021-09-16T15:38:34.747102Z","shell.execute_reply":"2021-09-16T15:38:35.193928Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### - The slice on the left is Axial plane, the one on the right is Sagittal.\n#### - We can see from the ImageOrientationPatient tag values that these images are orthogonal to all three planes of the patient coordinate system.\n#### - The cosine values are all zero or one.\n#### - Now, we'll use the ImageOrientationTag to add position markers to the right and bottom of the image.","metadata":{}},{"cell_type":"code","source":"# Get markers and draw them on the pixels\nimg1_marker_r, img1_marker_b = get_markers(img1.ImageOrientationPatient)\nimg2_marker_r, img2_marker_b = get_markers(img2.ImageOrientationPatient)\n\nimg1_pixels = draw_markers(img1_pixels, img1_marker_r, img1_marker_b)\nimg2_pixels = draw_markers(img2_pixels, img2_marker_r, img2_marker_b)   \n    \n# Plot the results\nfig, axes = plt.subplots(nrows=1, ncols=2,sharex=True, sharey=True, figsize=(15, 5))\nax = axes.ravel()\nax[0].set_title(f'{img1.ImageOrientationPatient} - FLAIR Axial')\nax[0].imshow(img1_pixels, cmap='gray');\nax[1].set_title(f'{img2.ImageOrientationPatient} - T2w Sagittal')\nax[1].imshow(img2_pixels, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:35.196035Z","iopub.execute_input":"2021-09-16T15:38:35.196448Z","iopub.status.idle":"2021-09-16T15:38:35.602088Z","shell.execute_reply.started":"2021-09-16T15:38:35.196418Z","shell.execute_reply":"2021-09-16T15:38:35.601327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### - Now, we have Left and Posterior markers on the axial image, and the Feet/Posterior markers on the sagittal.\n#### - This tells us that we're viewing the axial slice as if we're standing at the patient's feet.\n#### - Let's plot out reference lines so we can visualize the orientation of the images onto the each other.","metadata":{}},{"cell_type":"code","source":"# Project reference line from img2 onto img1  \ncol_pixel, row_pixel = create_reference_line(img2, img1)\nfor i in range(0,2):\n    img1_pixels = cv2.line(img1_pixels, (col_pixel[i], row_pixel[i]), (col_pixel[i+1], row_pixel[i+1]), (255, 255, 255), 2)\n    \n# Project reference line from img1 onto img2\ncol_pixel, row_pixel = create_reference_line(img1, img2)\nfor i in range(0,2):\n    img2_pixels = cv2.line(img2_pixels, (col_pixel[i], row_pixel[i]), (col_pixel[i+1], row_pixel[i+1]), (255, 255, 255), 2)\n    \n# Plot the results\nfig, axes = plt.subplots(nrows=1, ncols=2,sharex=True, sharey=True, figsize=(15, 5))\nax = axes.ravel()\nax[0].set_title(f'{img1.ImageOrientationPatient} - FLAIR Axial')\nax[0].imshow(img1_pixels, cmap='gray');\nax[1].set_title(f'{img2.ImageOrientationPatient} - T2w Sagittal')\nax[1].imshow(img2_pixels, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:35.603354Z","iopub.execute_input":"2021-09-16T15:38:35.603759Z","iopub.status.idle":"2021-09-16T15:38:36.011749Z","shell.execute_reply.started":"2021-09-16T15:38:35.603729Z","shell.execute_reply":"2021-09-16T15:38:36.010642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### - Now we can see how the slices are spatially aligned to each other and in which orientation.\n#### - Let's try another study ..","metadata":{}},{"cell_type":"code","source":"# Read in two images from the same study, different series.\nimg1 = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00025/FLAIR/Image-350.dcm')\nimg2 = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00025/T1wCE/Image-55.dcm')\n   \n# Crunch the pixels down to 8 bit so we can draw an 8 bit graphics on them\nimg1_pixels = normalize_pixels(img1.pixel_array)\nimg2_pixels = normalize_pixels(img2.pixel_array)","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:36.013049Z","iopub.execute_input":"2021-09-16T15:38:36.013465Z","iopub.status.idle":"2021-09-16T15:38:36.058765Z","shell.execute_reply.started":"2021-09-16T15:38:36.013421Z","shell.execute_reply":"2021-09-16T15:38:36.0576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get markers and draw them on the pixels\nimg1_marker_r, img1_marker_b = get_markers(img1.ImageOrientationPatient)\nimg2_marker_r, img2_marker_b = get_markers(img2.ImageOrientationPatient)\n\nimg1_pixels = draw_markers(img1_pixels, img1_marker_r, img1_marker_b)\nimg2_pixels = draw_markers(img2_pixels, img2_marker_r, img2_marker_b)   \n\ncol_pixel, row_pixel = create_reference_line(img2, img1)\nfor i in range(0,2):\n    img1_pixels = cv2.line(img1_pixels, (col_pixel[i], row_pixel[i]), (col_pixel[i+1], row_pixel[i+1]), (255, 255, 255), 2)\n    \n# Project reference line from img1 onto img2\ncol_pixel, row_pixel = create_reference_line(img1, img2)\nfor i in range(0,2):\n    img2_pixels = cv2.line(img2_pixels, (col_pixel[i], row_pixel[i]), (col_pixel[i+1], row_pixel[i+1]), (255, 255, 255), 2)\n# Plot the results\nfig, axes = plt.subplots(nrows=1, ncols=2,sharex=True, sharey=True, figsize=(15, 5))\nax = axes.ravel()\nax[0].set_title(f'{img1.ImageOrientationPatient} - FLAIR Axial')\nax[0].imshow(img1_pixels, cmap='gray');\nax[1].set_title(f'{img2.ImageOrientationPatient} - T2w Sagittal')\nax[1].imshow(img2_pixels, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2021-09-16T15:38:36.060188Z","iopub.execute_input":"2021-09-16T15:38:36.060532Z","iopub.status.idle":"2021-09-16T15:38:36.482914Z","shell.execute_reply.started":"2021-09-16T15:38:36.060499Z","shell.execute_reply":"2021-09-16T15:38:36.481932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### - Now we have *Posterior* and *Left*/*Head* on the axial. Notice the reference lines aren't at right angles.\n#### - *Feet*/*Left*/*Posterior* and *Posterior*/*Right*/*Head* on the sagittal.\n#### - Notice the ImageOrientationPatient values are not whole numbers. This means the image is rotated along that axis.\n#### - The result can be visualized with the reference line .. The images are slightly obliqued.","metadata":{}},{"cell_type":"markdown","source":"### - Conclusion:\n\n#### - We can determine the laterality on slices as well as check their level of obliqueness by using the DICOM ImageOrientationPatient tag.\n#### - This can also allow us to relate coordinates between two slices more accurately.\n##### * Series need to have the same DICOM FrameOfReference tag to be spatially related. None of the studies in this competition have this tag, so they're not technically related. We'll assume they're close enough for now. \n\n#### Some of my other MR notebooks\n- Tumor Object Detection -> https://www.kaggle.com/davidbroberts/brain-tumor-object-detection\n- Determining DICOM image order -> https://www.kaggle.com/davidbroberts/determining-dicom-image-order\n- Determining MR image planes -> https://www.kaggle.com/davidbroberts/determining-mr-image-planes\n- Reference Lines on MR images -> https://www.kaggle.com/davidbroberts/mr-reference-lines\n- Manual VOI LUT on MR images -> https://www.kaggle.com/davidbroberts/manual-voi-lut-on-mr-images\n- Export DICOM Images by Plane -> https://www.kaggle.com/davidbroberts/export-dicom-series-by-plane/","metadata":{}}]}