BME470-570 Project Solved

Instructions

• In HW 1, we learned concepts in functional analysis for imaging through solving mathematical problems. In HW 2, we are learning about some of the computational tools for analysis of imaging systems such as singular-value decomposition and eigen analysis. In HW 3, we will learn about Fourier theory. Having learned these theoretical and computational tools, in this assignment, our objective will be to use these tools to analyze two imaging systems! Further, we will use these tools to conduct artificial-intelligence (AI)-based operations on the data obtained with these systems. We will conduct this analysis through a step-wise procedure.
• The assignment will be graded based on a presentation and a report.
• For each of the two imaging systems mentioned below, there are six sub-parts. The sixth sub-part will involve concepts from AI. However, this question is more challenging than the other questions due to multiple reasons such as this not being a course on AI, the AI-related question is more involved and longer and that question would require more computations. To ensure fairness while providing students with flexibility, we offer the following option to choose the questions a student wants to attempt:
  • –  The six sub-parts are divided into two categories, sub-parts (1)-(5) and subparts 6-7.
  • –  If you choose sub-parts (1)-(5), in your presentation, you would be randomly assigned to present one of these sub-parts. You would have to then present that sub-part indi- vidually. In your written report, you would have to provide solutions to all 5 sub-parts and for both the systems.
  • –  If you choose sub-part 6-7, you can collaborate and present that sub-part as a team of 3 people. In your written report, you would need to provide solution only to sub-part 6 and only for the system that you present.
  • –  You need to inform the instructor and the TAs by Nov. 12 if you want to work on sub-parts (1)-(5), or on sub-part (6)-(7). If you choose to work on sub-part (6)-(7), you would need to inform the composition of your team and if you have a preference for System 1 vs. System 2. You will be informed by Nov. 12 about the system that you will be required to present.
  • –  In case we do not hear from a student, we will assume that the student has chosen subparts (1)-(5) for their presentation and written report.
• The presentation and report components of this assignment are described below.
• Presentation:

1

• –  The objective of the presentation is to learn how to communicate mathematically com- plex ideas. This is a very important skill while presenting work on theoretical and computational imaging at venues such as conferences.
• –  The presentations will occur over the last three days of the class (Dec. 1, Dec. 6 and Dec. 8).
• –  In several of the questions, we ask that the student display the results for five objects. The requirement of five object functions is only for the presentation. Also, one of these objects should be the Shepp-Logan phantom. The other four object functions can be your choice.
• –  If the student chooses to attempt sub-parts (1)-(5), they will be assigned one sub-part within this category to present. The student will be assigned 8 minutes, of which 5- 6 minutes are assigned to presentation and the remaining time for a question-answer session. The questions will be asked by the students, TAs, and the instructor.
• –  For sub-part (6)-(7), the team of students will get a total of 20 minutes, with 15 minutes assigned for presentation and 5 minutes for questions. All members of the team need to present.
• –  In the presentation, the student should focus on outlining the concept they learned in the sub-part they are presenting. A guideline to organize the presentation is as follows:

  ∗ Introduce and very briefly go through how the problem was solved. ∗ Discuss the results.
  ∗ Describe the concept that was learned.

  More detailed guidelines have been provided in presentation format on Canvas.

• –  We have tried to ensure that sub-parts (1)-(5) have equal levels of difficulty. In case you choose this category of question to attempt, the students can mail their preferences on which of these sub-parts they would like to present to the instructor and teaching assistants by Nov. 17. To ensure that different students present different sub-parts, the sub-parts will be assigned to the students through a random-selection procedure that will account for the student preferences. The sub-parts will be assigned by around Nov. 22. Students who choose their own imaging system will need to mention this in their email and those students will be assigned one sub-part for that system.
• –  Students who have taken BME 470 and choose between subparts (1)-(5) can provide their preferences between sub-parts 1 − 4. Mention in your email that you have taken BME 470.
• –  The question-answer session at the end of each presentation is highly encouraged and will be part of the grading (both for the presenter and the rest of the students).

  • Report:

  • –  The report for this assignment is due Dec. 22, 2020 at 11:59 pm. There will be no

    extensions.

  • –  In the report, first answer all the sub-parts for one system before proceeding to the next

    system.

  • –  Start each sub-part on a new page.
  • –  In several of the questions below, we ask that the student display the results for five object functions. For the report, the results need to be displayed only for one object function: the Shepp-Logan phantom.

2

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 provid- ing 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 com- ments 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

3

system is given by

􏰈1 for 􏰀− 1 + m−1 􏰁 < y <= 􏰀− 1 + m 􏰁
hm(x,y) = 2 16 2 16 (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:

1. (a)  16 × 16 pixels
2. (b)  32 × 32 pixels
3. (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:

1. (a)  Compute the SVD and plot the singular value spectra.
2. (b)  Plot some of the singular vectors in object space as images.
3. (c)  For each singular vector plotted above, plot the corresponding singular vectors in data space in the sinogram format.
4. (d)  For five object functions on the unit disk, use the SVD analysis to determine and then display the measured and null components of the object.
5. (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.

4

2. (b)  Repeat with Poisson noise added to the data. To add Poisson noise in Matlab, you can use the poissrnd command. Remember that the Poisson noise must be added to the data, not to the object.
3. (c)  You will find that adding noise leads to poor-quality reconstruction when all the SVD values are used. Explain why.
4. (d)  Reduce the number of SVD values by discarding the values with smaller magnitude and repeat 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.

1. (a)  Describe the procedure to compute the SVD. Hint: Use similar approach as problem 3 in Homework 2.
2. (b)  Plot the singular value spectra.
3. (c)  Plot the some of the singular vectors in object space as images. Note that you will need to discretize the singular vectors.
4. (d)  Plot some of the singular vectors in data space in the sinogram format.
5. (e)  For five object functions on the unit disk, use the SVD analysis to determine and then display 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.

1. (a)  Use the SVD to perform a pseudoinverse reconstruction of the CD H operator for five object functions on the unit disk.
2. (b)  Repeat with Poisson noise added to the data. Remember that the Poisson noise must be added to the data, not to the object.
3. (c)  The results from part (b) above will likely be poor-quality reconstruction. Explain why that is the case.
4. (d)  Implement the regularization-based strategy described above to address the issue of poor-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 open- source 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

5

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

􏰄r − c􏰅
fs(r) = Asrect(r)circ σ (4)

s
where As denotes the signal amplitude and where circ􏰂r−c􏰃 denotes a function that is

σs
unity inside a circle centered at c = (cxs, cys) and of radius σs. Mathematically

􏰄r−c􏰅 􏰈1, |r−c|≤σs.

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

= 0, otherwise. (5) This model denotes the objects as a collection of Gaussian lumps, and is given by

s

􏰉N an 􏰆 (x−cx,n)2 (y−cy,n)2􏰇

(2πσn2)exp − 2σn2 − 2σn2 (6)

fb(r)=rect(r)

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.

n=1

6

2. 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.
3. Add the signal and the background. This yields the object over a 64 × 64 grid.
4. 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:

1. 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×8pixels

2. Computing the SVD: Perform the SVD of this system matrix to obtain the singular values and the singular vectors.
3. Obtaining the pseudo-inverse of the system operator: Use the SVD to obtain the pseudoinverse of the system operator.
4. Generating projection data: Apply the generated system matrix to the generated objects to compute 500 instances of the projection data.
5. 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. 7

2. Compute the DSC between the true segmentation and the segmentation estimated using the CNN. Plot the variation in DSC as a function of the number of pixels. Comment on your results.
3. For each image, compute the true signal area (This would be πσs2). 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

􏰄r − c􏰅
fs(r) = Asrect(r)circ σ (8)

s
where As denotes the signal amplitude and where circ􏰂r−c􏰃 denotes a function that is

σs
unity inside a circle centered at c = (cxs, cys) and of radius σs. Mathematically

􏰄r−c􏰅 􏰈1, |r−c|≤σs.

circ σ

s

= 0, otherwise. (9) 8

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(r)

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:

1. Generate signal: Let the signal be at the center of the object, i.e. if the coordinate system 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.
2. 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.
3. 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:

1. Generating the system matrix: Simulate the forward model for the system using the expansion function of the object as 32 × 32 pixels.
2. Computing the SVD: Perform the SVD of this system matrix to obtain the singular values and the singular vectors.
3. Obtaining the pseudo-inverse of the system operator: Use the SVD to obtain the pseudoinverse of the system operator.
4. Generating noiseless projection data: Apply the generated system matrix to the generated objects to compute 500 instances of the projection data. Let us refer to the generated noiseless projection data as g ̄.
5. Adding noise: In this step, we will be generating projection data at four different noise levels. For this purpose, we will be scaling the projection data generated in

the above step.

1. Low-noise version: Scale the generated projection so that the sum of g ̄ is 50,000.
2. Noisy version 1: Medium-noise version: Scale the generated projection so that

   the sum of g ̄ is 25,000.

3. Noisy version 2: High-noise version: Scale the generated projection so that the

   sum of g ̄ is 10,000.

4. Noisy version 3: Very high-noise version: Scale the generated projection so that

   the 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).

􏰉N an 􏰆 (x−cx,n)2 (y−cy,n)2􏰇

n=1

(2πσn2)exp − 2σn2 − 2σn2 (10)

9

vi. Reconstructing the image: Apply the pseudoinverse to the projection data at each noise 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.

1. Apply the trained CNNK to the 100 test images acquired at noise level K. Let us refer to the resultant images as “denoised” images.
2. Visually, indicate your perception of whether the “denoised” images are less noisy and 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
3. 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?

10