HOME
PDF

Mathematical formulation of the HYPR algorithms

California State University, Fullerton. Summer 2008 page compiled on July 1, 2015 at 10:07pm

1 Introduction

This report is a summary of the HIghly constrained Back PRojection (HYPR) team work performed so far relating to the HYPR research project. We will describe the work done and results found.

The goals set for the HYPR project included formulating the HYPR algorithm and some of its variations (such as Wright-Huang HYPR (WH-HYPR), I-HYPR and HYPR-LR) in a mathematical framework which would allow the study and analyze of these algorithms in relation to other well known non-linear methods such as maximum a-posteriori (MAP) estimation and Maximum Likelihood Expectation Maximization (MLEM). These algorithms, like HYPR, use prior information on the object being reconstructed and they are extensively used in nuclear medicine where the data is intrinsically under sampled.

The initial period of this project, which this report reﬂects on, was spent becoming familiar with the HYPR algorithm and its connection to MLEM.  Towards this goal, the HYPR algorithm was formulated mathematically and schematic diagrams created which helped in its implementation. MATLAB simulation software was developed to enable more understanding of the algorithm and its behavior by running it on a number of test cases. Initial comparison between the original HYPR and the WH-HYPR made on a number of diﬀerent test conﬁguration which are described in detail in the simulation section below. The MLEM algorithm was implemented and compared the HYPR algorithm.

In addition, A mathematical connection between HYPR and Expectation Maximization (EM) is described and formulated.

2 Mathematical formulation of the HYPR algorithms

2.1 Original HYPR

2.1.1 Mathematical formulation

Please see the appendix for a complete description of the notation used in this section and throrught the rest of the report.

The mathematics of this algorithm will be presented by using the radon transform notation and not the matrix projection matrix notation.

The projection is obtained by applying radon transform on the image at some angle

When the original object image does not change with time one can drop the subscript from and just write

Next, the composite image is found from the ﬁltered back projection applied to all the as follows

Notice that the sum above is taken over and not over . Next, a projection is taken from at angle as follows

Then the unﬁltered back projection 2-D image is generated

And the unﬁltered back projection 2-D image is generated

Then the ratio of is summed and averaged over the time frame and the result multiplied by to generate a HYPR frame for the time frameHence for the time frame we obtain

2.2 Wright-Huang variation of HYPR

2.2.1 Mathematical formulation

This mathematics of this algorithm will be presented by using the radon transform notation and not the matrix projection matrix notation.

The projection is obtained by applying radon transform on the image at some angle

The composite image is found from the ﬁltered back projection applied to all the

Notice that the sum above is taken over and not over . Next a projection is taken from at angle as follows

Then the unﬁltered back projection 2-D image is generated

And the unﬁltered back projection 2-D image is generated

Now the set of and over one time frame are summed the their ratio multiplied by to obtain the HYPR frame

3 HYPR connection to Expectation Maximization

The following is a discussion of the Mathematics that connects the MLEM algorithm to HYPR.

According to O'Halloran's paper[1] the for Maximum-Likelihood Expectation-Maximization (MLEM) algorithm is mathematically equivalent to HYPR. The MLEM algorithm can be used in image reconstruction for medical purposes. Positron Emission Tomography (PET) and Single-Photon Emission Computed Tomography (SPECT) are two types of image reconstruction processes where the MLEM algorithm is used. The purpose here is to show that the MLEM algorithm will work for HYPR reconstructions.

The MLEM algorithm is a process that approximates the solution to

 (1)

In connection to HYPR, we can view as a forward projection matrix, as the original image being projected and the projection produced. The goal is to relate the above matrix based formulation to the radon transform based formulation seen above in the HYPR mathematical section, which is

 (2)

We formulate the ﬁrst iteration of the MLEM algorithm based on equation (1) and see how it can be translated into the HYPR process of image reconstruction. The ﬁrst step of MLEM is

 (3)

This can be rewritten in matrix form

If we replace by in (4), we obtain

 (4)

The marked portion of the above equation can be viewed as the vector that is produced from unﬁltered back projection on the image produced by the ratio

Here the division is done element by element to produce the vector whose elements are the ratios of the respective elements of and

In HYPR the equation we want to tie to equation (3) above is as follows

 (5)

Where the represents an element by element multiplication, and the terms in (5) as deﬁned in the section of the HYPR mathematical derivation shown earlier. Hence for (5) and (4) to be equivalent,  We must have

Which represents the initial guess for the user image.  Therefore, by using for as an initial guess for the MLEM algorithm the above weighted term of the composite image , the MLEM algorithm will produce   which is a better approximation to the original image from that of the composite . And this is what the HYPR algorithm does. It uses the composite image to produce the HYPR image to approximate the original user image . Hence a one step of MLEM is equivalent to running HYPR for one time. Therefore, iterative HYPR algorithms can be seen as a multi step application of MLEM.

4 Software simulation and results

4.1 HYPR simulation

A software simulation written in MATLAB was designed and developed to enable more extensive HYPR testing of diﬀerent test conﬁgurations. The software is GUI based and all test results are saved in a tab-delimited plain text ﬁle to allow one further statistical analysis of the data generated by other software. The appendix contains a screen shot of the current version of the simulator (version 1.0).

4.1.1 Description of HYPR simulation and test results

This is a description of the diﬀerent tests performed. Both the original HYPR and the Wright-Huang HYPR (WH-HYPR) were run and results compared. In the following discussion, we use to mean number of projections in one time frame, and to mean the number of time frames. Hence the total number of projections is

The table below describes each test. In this table, a test with the letter 'a' represents the test being run using the original HYPR and a test with the letter 'b' represents the test being run using WH-HYPR. Each test was run under both the original HYPR and WH-HYPR.

The ﬁrst set of tests are designed to detect the eﬀect of Poisson noise on the accuracy of the HYPR algorithm as compared to the WH-HYPR. This was done for diﬀerent geometry of objects while keeping the number of projections per time frame and the Poisson noise parameter ﬁxed.

The second set of tests are designed to detect the eﬀect of Gaussian noise on the accuracy of the HYPR algorithm as compared to the WH-HYPR. This was done for diﬀerent geometry of objects while keeping the number of projections per time frame and the Gaussian distribution parameters (mean and variance) ﬁxed. The set of tests used here is smaller than the ﬁrst set due time limitation.

The third set of tests are designed to detect the eﬀect of increasing the number of projections on the accuracy of the HYPR and WH-HYPR. This was done under one ﬁxed conﬁguration and with the absence of noise.

The main measure of accuracy used was relative RMSE. This was calculated as follows: For each HYPR frame image generated, the set of user images which make up the time frame from which the HYPR frame was generated are averaged to obtain an average time frame image. Then the RMSE was obtained between these 2 images as follows: Assuming these are total pixels in each image, the error between each corresponding pixels is found as where is a pixel in the HYPR frame image and is the corresponding pixel in the averaged time frame image. This error is then squared. This was done for each pixel. The average of these square values is found, and the square root of the result is found. Hence

This quantity is normalized by dividing it by the mean intensity of the averaged time frame image found earlier. This gives a normalized RMSE value for each time frame generated. When there are more than one time frame generated then the average all these RMSE values are used to obtain one representative value of the RMSE for the test, and that is the value showed in the tables below for the purpose of comparing diﬀerent tests.

Other statistics are calculated to determine the algorithms accuracy. The relative error between the HYPR image and the averaged time frame is found using the standard formula for relative error. This measure however did not appear to be a good indicator for determining the accuracy of the HYPR image. Another statistical measure used is the histogram diﬀerence, which is found as follows. The histogram for each HYPR image frame and the corresponding histogram for the averaged time frame image are calculated and the diﬀerence between these histograms found. This measure appear to give a good indication of the performance of each test and correlated well with the RMSE measure used. These results are all written to the log ﬁle for further analysis, but are not currently taken into account in the following tests due to time limitation. Only the RMSE measure is currently used to determine the accuracy of the algorithm. The following sections describe each set of tests in more details.

4.1.2 The ﬁrst set of tests

The appendix shows the output obtained from the above set of tests. We now present a summary of the results

Observation from running the ﬁrst set of tests  The original HYPR algorithm performed better in each test when noise is absent from projection. This occurs in either time varying or non-time varying conﬁguration. On the other hand, WH-HYPR performed better in each case when noise was present. This occurs in either time varying or non-time varying conﬁguration.

4.1.3 The second set of tests

These tests are a repeat of the ﬁrst set of tests, but with noise generated from normal distribution. Due to time limitation only test 2,6 and 10 described above are repeated since these 3 tests are good representative of the overall tests. The letter N is added to the test name to indicate the use of Normal distribution.

The appendix shows the output obtained from the above set of tests. Summary of the results is shown below. To clarify the nature of tests below a short description is given again below.

1.
Test 2N is a small ﬁxed disk, changes intensity linearly with time.
2.
Test 6N is one disk which moves vertically, oﬀ center.
3.
Test 10N is two disks separated from each others that move vertically across the image.

Observations from running the second set of tests  WH-HYPR performed better in all 3 cases. This correlated well with results found from the ﬁrst set of tests where it was observed that WH-HYPR performed better each time Poisson noise was added. In the above 3 tests, normal noise was added and it is observed that WH-HYPR performed better.

4.1.4 Third set of tests

As was mentioned earlier, the goal of these tests is to measure the relative accuracy of the algorithms on the same conﬁguration but with increasing number of projections per time frame. It is expected that the accuracy of each algorithm will improve, and we wish to obtain the measure of this improvement as a function of the number of projections per time frame.

For this purpose, the following test conﬁguration was used: small white disk moving vertically and oﬀ center, no noise added. One time frame was used and the following number of projections . These tests as named and respectively. The table below show the result of the tests.

Observations from running the third set of tests  As the number of projections per time frame increased, the accuracy of WH-HYPR improved. At high number of projections (over 512 per time frame) WH-HYPR bypassed original HYPR and became more accurate. It is not clear at this time if such high number of projections per time frame will conﬂict with other MRI requirements (sampling rate limitation or other issues), but the above shows that, even with the absence of noise, the WH-HYPR can become more accurate than the Original HYPR but at a cost of having large number of projections per time frame.

4.1.5 Conclusions drawn from HYPR test results

Original HYPR performed better than WH-HYPR when the number of projections is relatively low (below 256 per time frame) and when there was no noise present (noise added to projections taken from the original images). This occurred in all conﬁgurations (both objects moving in time or ﬁxed).

WH-HYPR performed better when noise is present (both Poisson and Normal noise) and for all number of projections and for all conﬁgurations.

In addition, WH-HYPR performed better when there was no noise added, but when the number of projections per time frame was increased.

These results seem to be a direct consequences of the fact that WH-HYPR sums the backprojection images over a time frame period before taking the ratio of these sums in order to obtain the mask image, while in the original HYPR the ratio for each backprojection images is ﬁrst found and the ratios added and averaged. More analysis will be needed to better understand this diﬀerence and to explain mathematically this observed diﬀerence between HYPR and WH-HYPR.

Since real MRI data is characterized by low SNR, this leads one to conclude that WH-HYPR should be the preferred choice between these 2 algorithms.

4.2 Expectation Maximization simulation

4.2.1 Description of simulation

The original HYPR algorithm was compared to 1 step of the MLEM algorithm. A time-invariant white disk with radius pixels centered in a by black image. diﬀerent projection angles were used (ordered using bit-reversed ordering), and the size of the window was set to 8 projections.

4.2.2 Results of simulation

The ﬁgures below are the actual images produced. The composite image, the HYPR reconstruction for the ﬁrst HYPR frame, and the corresponding MLEM image. The HYPR and the MLEM images are indistinguishable, although the mean absolute error is slightly higher for HYPR than for MLEM. More detailed comparisons of MLEM and HYPR are planned.

5 Future work

1.
Continue research into iterative HYPR and its connection to Expectation Maximization. Both analytically and through simulation.
2.
Work more on understanding the artifacts (modeling errors) and in the extreme, the pathological cases in which the HYPR algorithms will fail (worst-case scenario).
3.
Characterizing the noise ampliﬁcation and resolution of the HYPR algorithm through. Simulate HYPR algorithm with projection subjected to diﬀerent noise distributions and determine which variations of the algorithm are most accurate and under which  conditions.
4.
Investigate and implement a new Iterative HYPR variation proposed during work on this project which uses the Wright-Huang as its iterative step and compare to the current standard I-HYPR which uses the original HYPR and compare performance.
5.
Investigate possibility of a better measure to compare the accuracy of a HYPR image to a time frame image than was used in this report (RMSE), and if one is found, use the new measure for future testing.
6.
Obtained a mathematical description of the HYPR method based on the matrix formulation and not based on the radon transform. Apply for a simple geometrical shape which is time varying.

6 Appendix

6.1 nomenclature

1.
MLEM Maximum-Likelihood Expectation-Maximization
2.
PET Positron Emission Tomography
3.
SPECT Single-Photon Emission Computed Tomography
4.
A 2-D image. This represent the original user image at which the HYPR algorithm is applied to.
5.
When the original image content changes during the process, we add a subscript to indicate the image at time instance .
6.
7.
radon transform invoked at a projection angle .
8.
When the projection angle is not constant but changes with time during the MRI acquisition process, we add a subscript to indicate the angle at time instance .
9.
radon transform invoked at a projection angle .
10.
. radon transform applied to an image at angle . This results in a projection vector .
11.
Forward projection matrix. The Matrix equivalent to the radon transform .
12.
Estimate of an image .
13.
Multiply the forward projection matrix with an image estimate .
14.
Multiply the forward projection matrix with an image estimate to obtain a projection vector . Notice that for the inner dimensions of the matrix multiplication operation to be equal, this requires that the 2D image be linearized. In other words, the 2D image be written as a column vector.
15.
The inverse radon transform applied in unﬁltered mode to a projection which was taken at angle . This results in a 2D image.
16.
The inverse radon transform applied in ﬁltered mode to a projection which was taken at angle . This results in a 2D image.
17.
The transpose of the forward projection matrix multiplied by the projection vector . This is the matrix equivalent of applying the inverse radon transform in an unﬁltered mode to a projection (see item 12 above).
18.
The pseudo inverse of the forward projection matrix being multiplied by the projection vector . This is the matrix equivalent of applying the inverse radon transform in ﬁltered mode to a projection (see item 13 above).
19.
Composite image generated by summing all the ﬁltered back projections from projections of the original images . Hence
20.
The unﬁltered backprojection 2D image as a result of applying where is projection from user image taken at angle .
21.
The unﬁltered backprojection 2D image as a result of applying where is projection from the composite image taken at angle .
22.
Number of projections used to generate one HYPR frame image. This is the same as the number of projections per one time frame.
23.
The total number of projections used. This is the number of time frames multiplied by
24.
The HYPR frame image. A 2-D image generate at the end of the HYPR algorithm. There will be as many HYPR frame images as there are time frames.

References

[1]   Iterative projection reconstruction of time-resolved images using highly-constrained back-projection (HYPR) by Rafael L. O'Halloran, Zhifei Wen, James H. Holmes, Sean B. Fain

[2]   Highly Constrained Back projection for Time-Resolved MRI by C. A. Mistretta, O. Wieben, J. Velikina, W. Block, J. Perry, Y. Wu, K. Johnson, and Y. Wu

[3]   Principles of computerized Tomographic imaging by Kak and Staney

[4]   Highly Constrained Backprojection for Time-Resolved MRI by C. A. Mistretta, Wieben,z J, Velikina,W. Block,J. Perry,Y. Wu. K. ohnson and Y. Wu.