Lorenzo Castiglia - 10578631
In this document I will address the work that I have done regarding the NL1 project, spanning across various image processing tasks. This project consists of two distinct sections, for both of which I will elaborate on the tasks to be performed, the theoretical aspects, the code implementation and the results obtained.
The first part of the project aimed at exploring some basic notions about image analysis and manipulation. The task was to manually implement the gaussian filter and the Sobel edge detector in the C language and apply them to 8 medical images saved in the bmp format.
My code is composed by two files: an header file named ‘bmp.h’ containing useful data structures and methods to load and save the bmp images and ‘project.c’ where the actual implementation is written.
The latter takes as input the path of the input image and the gaussian filter kernel size (which will be discussed in the next section) by means of the main function parameters. The output will be two new bmp images where in one the gaussian filter is applied (‘gauss_output.bmp’) and the other shows the output of the Sobel edge detector (‘sobel_output.bmp’).
1.2.1 Gaussian filter
[1] The gaussian filter is a 2D convolution operator which blurs an image, reducing noise and smoothing the edges. This filter is particularly useful as a preprocessing function for other filters types such as edge detection.
It uses a special kernel calculated using the Gaussian distribution formula, as shown in figure 1.1.
In my implementation the kernel is dynamically allocated, so that its size can be changed. The main function takes as first parameter the input image file path, the second parameter is optional and takes an integer which is used to set the kernel size. The default kernel size is 5, with which I obtained the most convincing results.
Some examples of the effects of this algorithm are shown in the figure 1.2 below.
1.2.2 Sobel filter
[2] The Sobel filter is a differential operator used as an edge detector for images. The algorithm emphasises regions of images with high spatial frequency, in other words how much the pixel intensity values change. Differently from the gaussian filter, this operator convolves using two different kernels, one for vertical edges and one for horizontal edges, as shown in figure 1.3. The two outputs obtained are combined to find the absolute magnitude of the gradient.
Some examples of the effects of this algorithm are shown in figures 1.4 and 1.5.
The output images however are particularly noisy and unpleasing. A way to solve this issue is to apply the gaussian filter in advance. Some examples follow.
In the results directory, 8 sub-directories are present, one for each input image. In the corresponding sub-directory there’s the original image and the various outputs generated by my code, including the combination of both filters as mentioned previously.
In addition, there’s the output of Matlab’s Sobel filter, which is implemented differently. In this case there’s a thresholding applied, forming a logical image where only bigger edges are displayed. This is another way to reduce noise in the output.
The aforementioned filters, together with many others, are often already implemented in high-level programming languages, such as Matlab and Python. However, using one line of code to execute them, although convenient, makes it difficult to actually understand how they work and in general how an image is manipulated.
I found that manually implementing these functions was a great way of getting started and introducing myself to the topic of image processing.
The second part of the project, building on the basic notions of the first part, consisted in the more advanced topic of image segmentation.
The starting point was the CT-ORG dataset [3], a collection of human bodies CT scans saved in the NIfTI1 file format, a medical data format for storing 3D scans. Once getting acquainted with this unusual file type and its visualisation techniques, the objective was to create a mask of the lungs portion of the CT scan.
To understand the methodology and the theory behind this task, I heavily relied on the research done when creating the CT-ORG dataset that I used. [4]
Note: to display every example image in this report I used 3D Slicer [5], a powerful and open source software package used in medical research.
Further information can be found in the official website: https://www.slicer.org
The following image shows an example slice from one scan of the CT-ORG dataset. The lungs are clearly visible as the two black air pockets inside the body.
I chose to write the code in Matlab because of the amount of documentation and utilities that already exist for this kind of subject. The script ‘processFile.m’ takes as input the path of the nifti image to segment, the outputs are two new nifti images: the lungs mask and the image of the lungs (original file combined with the mask).
In order to generate the mask, this script calls the function ‘lungsMask.m’, which performs the segmentation of the lungs.
Note: in order to work, the algorithm needs the ‘aimutil’ repository, a collection of scripts to handle medical images that can be found at the following link.
aimutil repository: https://github.com/bbrister/aimutil
2.2.1 Morphological detection and segmentation of the lungs
In order to perform this task, I followed the methodology presented in the aforementioned CT-ORG research [4], which I will briefly recap in this section.
These steps are based on the fact that the lungs are the two largest air pockets in the body.
- The air pockets are the voxels with values less than -150 Hounsfield Units, to extract them a thresholding is applied to the image.
- The air inside the hollow structure of the exam table where the patient is resting and the air outside the body has to be eliminated. In this step the exam table structure gets punctured with a prism shape to connect it to the air outside.
- The air outside the patient stretches until the border of the image, so in any axial slice these components of the mask are removed, which thanks to point 2, removes the air inside the table as well. Lungs are instead preserved.
- Other small air pockets in the body have to be removed as well. The image is opened using a spherical structuring element with diameter of 1cm.
- In the mask obtained, consider the two largest connected components, which are almost certainly the lungs.
- To undo the effect of erosion, a morphological reconstruction is performed
In figure 2.2 the computed mask of the lungs is shown, while in 2.3 the mask is applied to the original image to only show the lungs.
This algorithm is fairly robust in most cases present in the dataset and, as mentioned in the original article, it was a useful tool for generating a vast amount of segmented data to use for training a machine learning model.
However the results are not always comparable with manually segmented images, where a doctor can very precisely define the borders of every organ.
Most of the aspects of this second part of the project were completely new to me, starting from the NIfTI1 file format, the medical applications and the tools I used, therefore making it challenging yet particularly engaging.
In addition, I experienced firsthand the challenges of generating enough data to train a model and how some experts come up with a solution, which was something I had only discussed during theory lessons.
[1] Gaussian filter: https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm
[2] Sobel edge detector: https://homepages.inf.ed.ac.uk/rbf/HIPR2/sobel.htm
[3] CT-ORG dataset: https://wiki.cancerimagingarchive.net/display/Public/CT-ORG:+CT+volumes+with+multiple+organ+segmentations#61080890bcab02c187174a288dbcbf95d26179e8
[4] Rister, B., Yi, D., Shivakumar, K. et al. CT-ORG, a new dataset for multiple organ segmentation in computed tomography. Sci Data 7, 381 (2020). https://doi.org/10.1038/s41597-020-00715-8
[5] Fedorov A., Beichel R., Kalpathy-Cramer J., Finet J., Fillion-Robin J-C., Pujol S., Bauer C., Jennings D., Fennessy F.M., Sonka M., Buatti J., Aylward S.R., Miller J.V., Pieper S., Kikinis R. 3D Slicer as an Image Computing Platform for the Quantitative Imaging Network. Magnetic Resonance Imaging. 2012 Nov;30(9):1323-41. PMID: 22770690. PMCID: PMC3466397.
[6] Image segmentation introduction: https://www.analyticsvidhya.com/blog/2019/04/introduction-image-segmentation-techniques-python/
[7] Alternative Matlab tool for nifti visualisation and segmentation: https://it.mathworks.com/matlabcentral/fileexchange/74761-imtool3d_td?s_tid=srchtitle
[8] Nifti image manipulation tutorial: https://dartbrains.org/content/Introduction_to_Neuroimaging_Data.html