|
| 1 | +# Day 13 Medical-Image Segmentation |
| 2 | +Medical applications are among the most exciting use cases of image segmentation networks. |
| 3 | +In this exercise, you will study the publication |
| 4 | +["Towards Patient-Individual PI-RADS v2 Sector Map: |
| 5 | +CNN for Automatic Segmentation of Prostatic Zones |
| 6 | +from T2-Weighted MRI"](https://www.var.ovgu.de/pub/2019_Meyer_ISBI_Zone_Segmentation.pdf) |
| 7 | +by Meyer et al. |
| 8 | + |
| 9 | + |
| 10 | +1. To get started run |
| 11 | + |
| 12 | +```bash |
| 13 | +python ./data/download.py |
| 14 | +``` |
| 15 | + |
| 16 | +in your terminal. The script will download and prepare the medical scans and domain-expert |
| 17 | +annotations for you. |
| 18 | + |
| 19 | +Data loading and resampling work already. |
| 20 | + |
| 21 | +1. #### Find the bounding box roi as described below by finishing the `compute_roi` function. |
| 22 | +Once you have obtained the train and test data, you must create a preprocessing pipeline. |
| 23 | +Proceed to `src/util.py` and compute the so called region of interest. |
| 24 | +Meyer et al. define this region as: |
| 25 | + |
| 26 | +"The images were acquired by two different |
| 27 | +types of Siemens 3T MRI scanners (MAGNETOM Trio and Skyra) |
| 28 | +with a body coil. The ground truth segmentation of the prostate |
| 29 | +zones was created on the axial images with 3D Slicer [19] by a medical |
| 30 | +student and subsequently corrected by an expert urologist. All |
| 31 | +volumes were resampled to a spacing of 0.5 × 0.5 × 3 mm which |
| 32 | +corresponds to the highest in-plane resolution and maintains the rela- |
| 33 | +tion of in-plane to inter-plane resolution of the dataset. A bounding |
| 34 | +box ROI of the prostate was automatically extracted with help of |
| 35 | +sagittal and coronal T2w series: the ROI was defined as the intersecting |
| 36 | +volume of the three MR sequences." |
| 37 | + |
| 38 | +See wikipedia's [anatomical plane](https://en.wikipedia.org/wiki/Anatomical_plane) article for a description of the terminology. |
| 39 | +The plots below depict the situation for the 0004 scans: |
| 40 | + |
| 41 | + |
| 42 | + |
| 43 | + |
| 44 | + |
| 45 | +After computing the intersection of all tensors, we can consider i.e. slice 12 of the |
| 46 | +transversal scan: |
| 47 | + |
| 48 | + |
| 49 | + |
| 50 | +Your implementation needs to translate the array indices from local into global coordinate systems and back. |
| 51 | +In other words, we require a rotation and translation, or more formally |
| 52 | + |
| 53 | +$$ \mathbf{R}\mathbf{x} + \mathbf{o} = \mathbf{g} .$$ |
| 54 | + |
| 55 | +With a rotation matrix $\mathbf{R} \in \mathbb{R}^{3,3}$, the local coordinate vector $\mathbf{x \in \mathbb{R}^{3}}$, the offset $\mathbf{o} \in \mathbb{R}^{3}$, and the global coordinate line $\mathbf{g}$. |
| 56 | +Evaluate this transform for every coordinate box line. Use the `box_lines` function from the |
| 57 | +`util.py` module to generate a bounding box at the origin. All points in every line must be transformed using the above relationship. |
| 58 | + |
| 59 | +The region of interest is the overlap of all boxes in the global coordinate system. Use [np.amin](https://numpy.org/doc/stable/reference/generated/numpy.amin.html) and [np.amax](https://numpy.org/doc/stable/reference/generated/numpy.amax.html) to find roi-box points $\mathbf{r} \in \mathbb{R}^{3}$. |
| 60 | + |
| 61 | +To obtain array indices, transform all box points back into the local system. Or, more formally: |
| 62 | + |
| 63 | +$$ \mathbf{R}^{-1} \mathbf{r} - \mathbf{o} = \mathbf{x}_{\text{roi}} $$ |
| 64 | + |
| 65 | +With the inverse of the rotation matrix $\mathbf{R}^{-1}$ use [np.linalg.inv](https://numpy.org/doc/stable/reference/generated/numpy.linalg.inv.html) to compute it. $\mathbf{x}_{\text{roi}} \in \mathbb{R}^{3}$ is a point on the boundary of the local roi-box we seek. |
| 66 | +Transform all boundary points. |
| 67 | + |
| 68 | +Using the smallest and largest coordinate values of the roi box in |
| 69 | +local coordinates now allows array indexing. Following Meyer et al. we discard all but the axial `t2w` scans. |
| 70 | + |
| 71 | +Test your implementation by setting the if-condition wrapping the plotting utility in `compute_roi` to `True` and running vscode pytest `test_roi`. Remember to set it back to `False` afterwards. |
| 72 | + |
| 73 | +2. #### Implement the UNet. |
| 74 | + |
| 75 | +Navigate to the `train.py` module file in the `src` folder. |
| 76 | +Finish the `UNet3D` class, as discussed in the lecture. |
| 77 | +Use the [flax.linen.Conv](https://flax.readthedocs.io/en/latest/api_reference/flax.linen/_autosummary/flax.linen.Conv.html), [flax.linen.relu](https://flax.readthedocs.io/en/latest/api_reference/flax.linen/_autosummary/flax.linen.activation.relu.html), and [flax.linen.ConvTranspose](https://flax.readthedocs.io/en/latest/api_reference/flax.linen/_autosummary/flax.linen.ConvTranspose.html), to build your model. |
| 78 | + |
| 79 | +3. #### Implement the focal-loss |
| 80 | + |
| 81 | +Open the `util.py` module in `src` and implement the `softmax_focal_loss` function as discussed in the lecture: |
| 82 | + |
| 83 | +$$\mathcal{L}(\mathbf{o},\mathbf{I})=-\mathbf{I}\cdot(1-\sigma_s(\mathbf{o}))^\gamma\cdot\alpha\cdot\ln(\sigma_s(\mathbf{o})) $$ |
| 84 | + |
| 85 | +with output logits $\mathbf{o}$, the corresponding labels $\mathbf{I}$ and the softmax function $\sigma_s$. |
| 86 | + |
| 87 | +4. #### Run and test the training script. |
| 88 | + |
| 89 | +Execute the training script with by running `scripts/train.slurm` (locally or using `sbatch`). |
| 90 | + |
| 91 | +After training you can test your model by changing the `checkpoint_name` variable in `src/sample.py` to the desired model checkpoint and running `scripts/test.slurm`. |
| 92 | + |
| 93 | +#### Solution: |
| 94 | + |
| 95 | + |
| 96 | + |
| 97 | + |
| 98 | +5. #### (Optional) Implement mean Intersection-over-Union (mIoU) |
| 99 | + |
| 100 | +Open the `meanIoU.py` in `src` and implement the `compute_iou` function as discussed below. |
| 101 | +mIoU is the most common metric used for evaluating semantic segmentation tasks. It can be computed using the values from a confusion matrix as given below |
| 102 | + |
| 103 | +$$\text{mIoU} = \frac{1}{k} \sum_{c=0}^{k}\frac{TP_c}{TP_c+FP_c+FN_c}$$ |
| 104 | + |
| 105 | +where `k`, `TP`, `FP`, and `FN` are number of classes, True Positives, False Positives and False Negatives respectively. |
| 106 | +The mIoU value generally ranges from 0 to 1 with 0 means no intersection area between predicted segmentation map and ground truth map and 1 means its a perfect fit between these two. |
| 107 | +Generally any $\text{mIoU}>0.5$ is considered better. |
| 108 | + |
| 109 | +Run the script with |
| 110 | +```bash |
| 111 | +python -m src.meanIoU |
| 112 | +``` |
| 113 | + |
| 114 | +### Acknowledgments: |
| 115 | +We thank our course alumni Barbara Wichtmann, for bringing this problem to our attention. |
| 116 | +Without her feedback, this code would not exist. |
| 117 | + |
0 commit comments