BME470-570-Project Solved

Figure 1: Demonstrating the procedure to display a sinogram

–    Collaboration on the assignment is fine, but any written material, including code, should be your own.

–    We recommend for the report to be typed.

–    Organize your solutions so that they are easy to follow. Some suggestions include providing titles for plots, plotting legible figures, and referring to the plots by figure numbers in your solution. For example, in questions that require plotting and/or providing comments on the plots, the template below could be used as a guideline:

∗ Provide theoretical details of the solution and conceptual understanding.

∗ Provide any implementation details. This is a great place to provide the code.

∗ Plot results. For each plot, provide title and caption.

∗ Provide comments on the results. In each comment, refer to the plot by Fig. number (e.g. Fig. 3a)

A general thought to keep in mind while organizing your solution is putting yourself in the shoes of the grader.

•    In some of the questions below, you are asked to display the data vectors as a sinogram. In this sinogram image, the pixels for each detector position are stacked in a vertical strip, with a grey level in each pixel corresponding to the data for that pixel, and the strips for successive angles are placed next to each other. An example is shown in Fig. 1.

•    Your entire submission should be a single file containing your answers to all the questions in PDF format. Do not submit multiple files or zip files.

•    All the code used for this assignment must be submitted. If you prefer to not include the code along with your solution, then include the code as an Appendix in your submitted file.

The key purpose of this project is to become familiar with analyzing imaging systems using the tools learned in class. We will consider two imaging systems in this project.

(A)     System 1: 2-D tomographic system: Consider a simplified 2-D tomographic system. The field of view (FOV) is the unit disk in the plane (diameter = 1). The kernel for the imaging system is given by

 (

1 for

hm(x,y) =(1)

0,otherwise

In other words, each sensitivity function is equal to 1 in a horizontal strip of width 1/8, and zero outside this strip. The strips are non-overlapping and stack up to cover the whole unit disk. This corresponds to a 16 pixel detector oriented vertically with perfect collimation. To get the next 16 sensitivity functions we rotate these by π/8. We continue rotating by π/8 to generate a total of M = 16 × 8 sensitivity functions.

(B)     System 2: 2-D radial MRI: Consider a highly simplified 2D radial MRI system. The kernel for the imaging system is given by:

                                                                                       hm(r) = exp(−2πikm · r)                                                                    (2)

and the FOV is the disk centered at the origin with diameter 1. The first 16 km are (m−8,0) for m = 1,2,...,16. To get the next 16 km we rotate the first 16 by π/8 . We continue rotating by π/8 to generate a total of M = 16 × 8 sensitivity functions.

For each of the systems above, or for a system of your choice, perform the following:

1.    System modeling: Simulate the forward model for this system using the following expansion functions for representing the object:

(a)     16 × 16 pixels

(b)    32 × 32 pixels

(c)     Fourier-basis functions for the square of side-length 1 units for system A and system B.These would be of dimensions 32 × 32.

In the process, you will obtain the H matrices for each of the expansion functions. For 5 object functions in the unit disk and using these H matrices, display and compare the data vectors that you get from the three different discretizations.

2.    SVD analysis of DD operator: For each of the H matrices generated in part 1, perform the following:

(a)     Compute the SVD and plot the singular value spectra.

(b)    Plot some of the singular vectors in object space as images.

(c)     For each singular vector plotted above, plot the corresponding singular vectors in dataspace in the sinogram format.

(d)    For five object functions on the unit disk, use the SVD analysis to determine and thendisplay the measured and null components of the object.

(e)     Comment on all the results.

3.    Reconstruction with DD operator: For five object functions of your choice on the unit disk and each H matrix obtained in 1, generate the data vector and perform the following:

(a)     Use the SVD to perform a pseudoinverse reconstruction for each H matrix for five object functions of your choice.

(b)    Repeat with Poisson noise added to the data. To add Poisson noise in Matlab, you canuse the poissrnd command. Remember that the Poisson noise must be added to the data, not to the object.

(c)     You will find that adding noise leads to poor-quality reconstruction when all the SVDvalues are used. Explain why.

(d)    Reduce the number of SVD values by discarding the values with smaller magnitude andrepeat the reconstruction. This type of process is referred to as regularization.

4.    SVD analysis of CD operator: In this sub-part, we work directly with the original continuous-to-discrete (CD) operator H, the kernel of which is given by Eq. (1) and Eq. (2) for the two systems. Note that the operator cannot be discretized in the problems below.

(a)     Describe the procedure to compute the SVD. Hint: Use similar approach as problem 3 in Homework 2.

(b)    Plot the singular value spectra.

(c)     Plot the some of the singular vectors in object space as images. Note that you will needto discretize the singular vectors.

(d)    Plot some of the singular vectors in data space in the sinogram format.

(e)     For five object functions on the unit disk, use the SVD analysis to determine and thendisplay the measured and null components of the object. Again you will need to discretize the object functions.

5.    Reconstruction with CD operator: Here we will perform reconstruction with the original CD operator. Again note that the operator cannot be discretized in the problems below.

(a)     Use the SVD to perform a pseudoinverse reconstruction of the CD H operator for five object functions on the unit disk.

(b)    Repeat with Poisson noise added to the data. Remember that the Poisson noise mustbe added to the data, not to the object.

(c)     The results from part (b) above will likely be poor-quality reconstruction. Explain whythat is the case.

(d)    Implement the regularization-based strategy described above to address the issue ofpoor-quality reconstruction.

In the next question, we will explore the application of the concepts we learned in class in the emerging era of AI, and more specifically, deep learning. We will use the tool of convolutional neural network (CNN) in these questions, and some working knowledge of how to train and test CNNs will be assumed. There are multiple references to this though, and several opensource codes. The emphasis will not be on AI, but instead on the tools we learned in the class.

6.    Image segmentation:

Background: Images typically consist of multiple regions. For example, a medical image of a patient with a disease may consist of normal organs and abnormal organs. Similarly, a image of a road may consist of different cars, traffic signals, road signs and so on. Image segmentation is the process of partitioning an image into these different regions. This helps in analyzing the image to extract relevant information. For example, segmentation of a tumor in a medical image of a patient with cancer can quantify the volume of the tumor.

In this question, we will develop a supervised CNN-based approach to segment lesions from images. Supervised deep-learning-based approaches for image segmentation are broadly based on the idea that if a neural network is provided a population of images, and for each image, is provided the corresponding label that specifies the region to which each pixel in the image belongs, then the approach can learn how to segment the images. More specifically, the network is “trained” to classify each pixel and thus assign that pixel a label. By minimizing the distance between the true labels and the estimated labels (using a loss function), through an iterative process, the network becomes trained to segment a new image.

As we see, this is a multi-step process. In this question, we will go through each of these steps one by one.

(a) Population generation: The first step is to generate a population of objects. Note here that it is important that there should be variability in your population, else the network may just memorize how to segment the tumor in your training set, and would fail when given a new image with a new tumor.

The object will be assumed to consist of two parts, signal and background. The signal and background will be denoted by fs(r) and fb(r). If we denote the object as f(r), then we can write

                                                                                                   f(r) = fs(r) + fb(r)                                                                     (3)

We consider a simple circular model of the signal, fs(r). This is given by

                                                                                         fs(r) = Asrect(r)circ                                                          (4)

where As denotes the signal amplitude and where circ  denotes a function that is unity inside a circle centered at c = (cxs,cys) and of radius σs. Mathematically

 σs.

                                                                            circ                                                                                                           (5)

.

For the background fb(r) we will use a simplistic version of the lumpy background model.

This model denotes the objects as a collection of Gaussian lumps, and is given by

                                                       fb(r) = rect(                          (6)

where N denotes the total number of lumps, an and σn denote the amplitude and width of the lump, and (cx,n,cy,n) denote the center of the lump. We will consider that cx,n and cy,n) are both sampled from a uniform distribution. Consider that N and σn are fixed.

To generate an object, we follow the following process:

i.    Generate signal: Sample the signal center (cxs,cys) from a uniform distribution with range (−1/4,1/4). Sample σs from a uniform distribution with range (0,1/2). Fix As = 3. Insert these sampled values into Eq. (8). This leads to a continuous representation of the signal. Discretize this over 64 × 64 pixels. This yields the discretized signal.

ii.  Generate background: Fix the number of lumps N = 10 and for each lump, fix σn to σs and an = 1. Next, for each lump, sample cx,n and cy,n, each from a uniform distribution with range (−1/2,1/2). Insert these values into Eq. (10). This yields the continuous version of the background. Discretize this over 64 × 64 pixels. iii. Add the signal and the background. This yields the object over a 64 × 64 grid.

iv. Generate the segmentation labels: Label all the pixel values as 1, other than the pixels occupied by the signal, which will be labeled as 2. This yields the true segmentation.

Repeat the above three steps for J = 500 times will yield 500 object realizations. In each realization, in addition to saving the object, also save the true segmentations.

(b)    Image generation: Using the imaging system model for System A or System B, conduct the following:

i. Generating the system matrix: Simulate the forward model for the system using the

expansion function of the object as pixels with the following dimensions:

• 32 × 32 pixels

• 16 × 16 pixels

• 8 × 8 pixels ii. Computing the SVD: Perform the SVD of this system matrix to obtain the singular values and the singular vectors.

iii.    Obtaining the pseudo-inverse of the system operator: Use the SVD to obtain thepseudoinverse of the system operator.

iv.     Generating projection data: Apply the generated system matrix to the generatedobjects to compute 500 instances of the projection data.

v.       Reconstructing the image: Apply the pseudoinverse to the projection data to ob-tain 500 instances of the reconstructed images. Conduct this process for all three dimensions of the number of pixels mentioned in subpart 6(b)i.

We will split this dataset of 500 images into three parts, a training dataset, consisting of 300 images, a validation dataset consisting of 100 images, and a testing dataset, consisting of the remaining 100 images. For each of these reconstructed images, we have the true segmented object.

(c)     Training the CNN: Using any generic CNN-based segmentation code, train a CNN to segment the 300 images in the training dataset, using the 100 images in the validation test for ensuring that there has been no overfitting. One such code is available on the following website: link. The link directs you to a github page, where the file of interest is cnn model.py. You would be required to train one CNN for each of the three dimensions mentioned in subpart 6(b)i. For the sake of notation, we will denote the CNN trained with images of size K × K pixels by CNNK.

(d)    Studying the effect of pixel size: In this step, we will evaluate the performance of the three CNNs on the corresponding images in the testing set. To quantify segmentation performance, we typically use the figure of merit of dice similarity coefficient (DSC). In this question, we will also quantify performance on a simpler version of a clinical task, namely quantifying tumor area. Perform the following operations for all the three dimensions mentioned in subpart 6(b)i.

i.         Apply the trained CNNK to the 100 test images of dimension K × K pixels.

ii.       Compute the DSC between the true segmentation and the segmentation estimatedusing the CNN. Plot the variation in DSC as a function of the number of pixels. Comment on your results.

iii.     For each image, compute the true signal area (This would be ). Compare that to segmentation obtained with the CNN. Plot this as a function of the signal radius, σs. Show four plots on the figure, corresponding to the true area, and the three area values estimated at the three pixel values. Comment on your results.

7. Image denoising:

Background: Images are typically corrupted by noise. This noise could be due to various sources such as photon noise and electronic noise. The corruption by noise impacts image quality and thus, in the imaging community, there are several efforts on making the image look visually less noisy. This operation is referred to as image denoising. In this question, we will developed a supervised CNN-based approach for image denoising. These are based on the idea that if a neural network is provided a population of noisy images, and for image, is provided a less noisy version, then the algorithm can learn how to predict the less noisy images given a more noisy image.

A key question is how do we know if a denoising method is effective and better compared to other methods. To answer this question, we recall that in scientific and medical imaging, images are acquired for a certain task. An example of this could be that a medical image of a patient exhibiting symptoms of say cancer is acquired so that a radiologist could assess if the patient has a tumor. This is referred to as the detection task. Another task could be measurement of the tumor volume. This is referred to as a quantification task. Ideally, a denoising should be evaluated based on how well it performs in that task. In this question, we will also evaluate the performance of this denoising technique on the task of detecting a certain signal.

As we see, this is a multi-step process. In this question, we will go through each of these steps one by one.

(a) Population generation: The first step is to generate a population of objects. Note here that it is important that there should be variability in your population, else the network may just memorize how to denoise a specific set of noisy images, and would fail when given a new image.

As in the previous subpart, the object will be assumed to consist of two parts, signal and background. The signal and background will be denoted by fs(r) and fb(r). If we denote the object as f(r), then we can write

                                                                                                   f(r) = fs(r) + fb(r)                                                                     (7)

We consider a simple circular model of the signal, fs(r). This is given by

                                                                                         fs(r) = Asrect(r)circ                                                          (8)

where As denotes the signal amplitude and where circ  denotes a function that is unity inside a circle centered at c = (cxs,cys) and of radius σs. Mathematically

 σs.

                                                                            circ                                                                                                           (9)

.

For the background fb(r) we will use a simplistic version of the lumpy background model.

This model denotes the objects as a collection of Gaussian lumps, and is given by

                                                       fb(r) = rect(                       (10)

where N denotes the total number of lumps, an and σn denote the amplitude and width of the lump, and (cx,n,cy,n) denote the center of the lump. We will consider that cx,n and cy,n) are both sampled from a uniform distribution. Consider that N and σn are fixed.

To generate an object, we follow the following process:

i.    Generate signal: Let the signal be at the center of the object, i.e. if the coordinatesystem was designed such that the origin falls at the center, then the center would be 0,0). Fix σs = 1/4 and As = 3. Insert these sampled values into Eq. (8). This leads to a continuous representation of the signal. Discretize this over 64×64 pixels. This yields the discretized signal.

ii.  Generate background: Fix the number of lumps N = 10 and for each lump, fix σn to σs and an = 1. Next, for each lump, sample cx,n and cy,n, each from a uniform distribution with range (−1/2,1/2). Insert these values into Eq. (10). This yields the continuous version of the background. Discretize this over 64 × 64 pixels. iii. Add the signal and the background. This yields the object over a 64 × 64 grid.

Repeat the above three steps for J = 500 times will yield 500 object realizations.

(b)    Image generation: Using the imaging system model for System A or System B, conduct the following:

i.         Generating the system matrix: Simulate the forward model for the system using theexpansion function of the object as 32 × 32 pixels.

ii.       Computing the SVD: Perform the SVD of this system matrix to obtain the singularvalues and the singular vectors.

iii.     Obtaining the pseudo-inverse of the system operator: Use the SVD to obtain thepseudoinverse of the system operator.

iv.     Generating noiseless projection data: Apply the generated system matrix to thegenerated objects to compute 500 instances of the projection data. Let us refer to the generated noiseless projection data as g¯.

v.       Adding noise: In this step, we will be generating projection data at four differentnoise levels. For this purpose, we will be scaling the projection data generated in the above step.

A.      Low-noise version: Scale the generated projection so that the sum of g¯ is 50,000.

B.      Noisy version 1: Medium-noise version: Scale the generated projection so thatthe sum of g¯ is 25,000.

C.      Noisy version 2: High-noise version: Scale the generated projection so that thesum of g¯ is 10,000.

D.     Noisy version 3: Very high-noise version: Scale the generated projection so thatthe sum of g¯ is 2,500.

. Next add Poisson noise to this data. In Matlab, this can be done by using the command poissrnd(gbar).

vi.     Reconstructing the image: Apply the pseudoinverse to the projection data at eachnoise level to obtain 500 instances of the reconstructed images. Conduct this process for all three dimensions of the number of pixels mentioned in subpart 7(b)v.

We will split this dataset of 500 images into three parts, a training dataset, consisting of 300 images, a validation dataset consisting of 100 images, and a testing dataset, consisting of the remaining 100 images.

Typically, in imaging, we never have access to the noiseless data. Thus, usually, a network is trained to estimate an image at the low noise levels, given the images at higher noise levels. Thus, in the above setup, for each of these reconstructed images, the low-noise version will be considered as the ground truth for the training process.

(c)     Training the CNN: Using any generic CNN-based denoising code, train a CNN to estimate the low-noise version given the (medium/high/very high) noisy image. You would assign 300 images in the training dataset and 100 images in the validation test for ensuring that there has been no overfitting. One such code is provided by us and can be used. You would be required to train one CNN for each of the three noise levels (medium/high/very high). For the sake of notation, we will denote the CNN trained for the images at “Noise level K” by CNNK.

(d)    Evaluating the denoising operation: In this step, we will evaluate the performance of the three CNNs on the corresponding images in the testing set. As mentioned above, the performance will be evaluated on a signal-detection task. Perform the following operations for all the three noise levels mentioned in Q. 7(b)v.

i.         Apply the trained CNNK to the 100 test images acquired at noise level K. Let us refer to the resultant images as “denoised” images.

ii.       Visually, indicate your perception of whether the “denoised” images are less noisyand look more similar to the low-noise image. Do you think the CNN is effective in suppressing noise? You can also study this using quantitative metrics such as root mean square error between the low-noise image and the true

iii.     We are providing you with code to evaluate performance on a signal-detection task.This code outputs the metric of area under the receiver operating characteristics curve (AUC). The values of the AUC lie between 0.5 and 1, and the higher the value, the more accurate the performance on the detection task. Compute and plot the AUC for all four noise levels mentioned in Q. 7(b)v. Does performance on the signal-detection task show a similar performance as the visual perception results?