{"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>Drawing Reference Lines On MR Images</h1>\n    - yet another MR notebook -\n    </div>\n    \n#### The goal of this notebook is to demonstrate how images within a study are spatially related to one another.\n\n- We'll use the ImageOrientationPatient and ImagePositionPatient tags to calculate and draw a reference line.\n- This is a simple demo of the technique. There are some scaling and spatial orientation steps that should be added for completeness, but this will do in the context of this comp.","metadata":{}},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pydicom\nimport matplotlib.pyplot as plt\nimport cv2","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:25:25.946416Z","iopub.execute_input":"2021-09-01T19:25:25.947018Z","iopub.status.idle":"2021-09-01T19:25:26.590478Z","shell.execute_reply.started":"2021-09-01T19:25:25.946975Z","shell.execute_reply":"2021-09-01T19:25:26.589496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Define a function that determines the reference line position from the position and orientation tags of a 'source' and 'destination' image.\n\n- Since we're projecting a 2D image into 3D space, we'll actually make a box that represents all four sides of the source image.\n- If the two images are orthogonal to each other in at least two planes, the sides of the 2D image will overlap and we'll see a single line drawn on the destination image.\n- It's important to understand the difference in the patient's coordinates (3D) and the image's coordinates (2D). \n- Here is a detailed explanation of how the ImagePositionPatient and ImageOrientationPatient tags work to orient the image. -> http://dicomiseasy.blogspot.com/2013/06/getting-oriented-using-image-plane.html","metadata":{}},{"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-01T19:25:26.592059Z","iopub.execute_input":"2021-09-01T19:25:26.592571Z","iopub.status.idle":"2021-09-01T19:25:26.607254Z","shell.execute_reply.started":"2021-09-01T19:25:26.592518Z","shell.execute_reply":"2021-09-01T19:25:26.606418Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Now we'll read in a couple images from the same study, but different series. \n- We want series in different planes, as it doesn't make sense to display reference lines on coplanar images.\n- Technically, we should verify the FrameOfReferenceUID tags match, which don't exist in this datatset, so we'll pretend they do.","metadata":{}},{"cell_type":"code","source":"# Read in two images from the same study, different series.\nsrc_img = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00006/T1w/Image-15.dcm')\ndst_img = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00006/T2w/Image-200.dcm')\n\n# Print these two tags so we can get used to seeing them. Nerds will be able to determine orientation by just looking at the IOP tags (after some time).\nprint(\"Src IOP: \", src_img.ImageOrientationPatient)\nprint(\"Dest IOP: \", dst_img.ImageOrientationPatient)\n\n# Crunch the destination image's pixels down to 8 bit so we can draw an 8 bit line on it\npixels = dst_img.pixel_array\npixels = pixels - np.min(pixels)\npixels = pixels / np.max(pixels)\npixels = (pixels * 255).astype(np.uint8)\n\n# Get the line (box) coords\ncol_pixel, row_pixel = create_reference_line(src_img, dst_img)\n\n# Draw lines on the destination image\nfor i in range(0,2):\n    dst_img = cv2.line(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=(14, 7))\nax = axes.ravel()\n\nax[0].set_title('Source Image (T1w - Axial)')\nax[0].imshow(src_img.pixel_array, cmap='gray');\nax[1].set_title('Destination Image (T2w - Sagittal)')\nax[1].imshow(dst_img, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:25:26.609297Z","iopub.execute_input":"2021-09-01T19:25:26.609822Z","iopub.status.idle":"2021-09-01T19:25:27.175668Z","shell.execute_reply.started":"2021-09-01T19:25:26.609769Z","shell.execute_reply":"2021-09-01T19:25:27.174553Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### We can see from the printout above the images that the source image IOP cosines are all either 0 or 1. The destination cosines are not whole numbers.\n- This means the source image **is** at right angles relative to the patient's body.\n- The destination image **is not** positioned at right angles to the patient's body.\n\n#### Now, we'll plot the same two images, but reverse their order so we can see this.","metadata":{}},{"cell_type":"code","source":"# Read in two images from the same study (technically, we should verify the StudyInstanceUID and SeriesInstanceUIDs match, but we'll assume they do here)\nsrc_img = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00006/T2w/Image-200.dcm')\ndst_img = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00006/T1w/Image-15.dcm')\n\nprint(\"Src IOP: \", src_img.ImageOrientationPatient)\nprint(\"Dest IOP: \", dst_img.ImagePositionPatient)\n\n# Crunch the destination image's pixels down to 8 bit so we can draw an 8 bit line on it\npixels = dst_img.pixel_array\npixels = pixels - np.min(pixels)\npixels = pixels / np.max(pixels)\npixels = (pixels * 255).astype(np.uint8)\n\n# Get the line (box) coords\ncol_pixel, row_pixel = create_reference_line(src_img, dst_img)\n\n# Draw lines on the destination image\nfor i in range(0,2):\n    dst_img = cv2.line(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=(14, 7))\nax = axes.ravel()\n\nax[0].set_title('Source Image (T2w - Sagittal)')\nax[0].imshow(src_img.pixel_array, cmap='gray');\nax[1].set_title('Destination Image (T1w - Axial)')\nax[1].imshow(dst_img, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:25:27.177328Z","iopub.execute_input":"2021-09-01T19:25:27.177661Z","iopub.status.idle":"2021-09-01T19:25:27.630876Z","shell.execute_reply.started":"2021-09-01T19:25:27.177627Z","shell.execute_reply":"2021-09-01T19:25:27.629616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### We can see when we project the T2w sagittal series onto the T1w axial series, it's tilted slightly.\n- This is because the technologist adjusted the reconstruction of this series to fit the patient who wasn't laying perfectly straight inside the magnet.\n\n#### Let's add another series from this study","metadata":{}},{"cell_type":"code","source":"# Read in two images \nsrc_img = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00006/FLAIR/Image-455.dcm')\ndst_img = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00006/T2w/Image-200.dcm')\n\nprint(\"Src IOP: \", src_img.ImageOrientationPatient)\nprint(\"Dest IOP: \", dst_img.ImagePositionPatient)\n\n# Crunch the destination image's pixels down to 8 bit so we can draw an 8 bit line on it\npixels = dst_img.pixel_array\npixels = pixels - np.min(pixels)\npixels = pixels / np.max(pixels)\npixels = (pixels * 255).astype(np.uint8)\n\n# Get the line (box) coords\ncol_pixel, row_pixel = create_reference_line(src_img, dst_img)\n\n# Draw lines on the destination image\nfor i in range(0,2):\n    dst_img = cv2.line(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=(14, 7))\nax = axes.ravel()\n\nax[0].set_title('Source Image (FLAIR - Axial)')\nax[0].imshow(src_img.pixel_array, cmap='gray');\nax[1].set_title('Destination Image (T2w - Sagittal)')\nax[1].imshow(dst_img, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:25:27.632213Z","iopub.execute_input":"2021-09-01T19:25:27.632513Z","iopub.status.idle":"2021-09-01T19:25:28.085611Z","shell.execute_reply.started":"2021-09-01T19:25:27.632481Z","shell.execute_reply":"2021-09-01T19:25:28.084379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### This one is slightly tilted as well. \n- It's difficult to tell by the line, but you can see from the IOP numbers, that neither of these series are orthogonal to the patient.\n\n#### Let's go back to the T2w (sagittal) and T1w (coronal) combo and plot a bunch of reference lines.","metadata":{}},{"cell_type":"code","source":"# Load the destination image\ndst = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00006/T1w/Image-11.dcm')\npixels = dst.pixel_array\npixels = pixels - np.min(pixels)\npixels = pixels / np.max(pixels)\npixels = (pixels * 255).astype(np.uint8)\n\n# Load a bunch of images from the source series and plot their lines onto the dest image\nfor i in range(100, 300, 10):\n    src = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00006/T2w/Image-' + str(i) + '.dcm')\n\n    # Get the line (box) coords\n    col_pixel, row_pixel = create_reference_line(src, dst)\n\n    # Draw lines on the destination image\n    for i in range(0,2):\n        pixels = cv2.line(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(figsize=(14, 7))\nplt.imshow(pixels, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2021-09-01T19:25:28.087158Z","iopub.execute_input":"2021-09-01T19:25:28.087472Z","iopub.status.idle":"2021-09-01T19:25:28.577051Z","shell.execute_reply.started":"2021-09-01T19:25:28.087439Z","shell.execute_reply":"2021-09-01T19:25:28.576018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### This shows the orientation of every tenth image of the T2w sagittals on the T1w axial, starting at image 100 and ending at image 300.\n- You could loop through the entire series, remove empty images etc.","metadata":{}},{"cell_type":"markdown","source":"#### Some of my other MR notebooks\n- Tumor Object Detection -> https://www.kaggle.com/davidbroberts/brain-tumor-object-detection\n- Determining MR Slice Orientation -> https://www.kaggle.com/davidbroberts/determining-mr-slice-orientation\n- Determining MR image planes -> https://www.kaggle.com/davidbroberts/determining-mr-image-planes\n- Determining DICOM image order -> https://www.kaggle.com/davidbroberts/determining-dicom-image-order\n- Manual VOI LUT on MR images -> https://www.kaggle.com/davidbroberts/manual-voi-lut-on-mr-images\n- Standardizing MR Images -> https://www.kaggle.com/davidbroberts/standardizing-mr-images\n- Export DICOM Images by Plane -> https://www.kaggle.com/davidbroberts/export-dicom-series-by-plane","metadata":{}}]}