- HOME

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

- 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.
- radon transform.
- 7.
- radon transform used at a projection angle .
- 8.
- When the projection angle is not constant but changes with time during the MRI acquiztion process, we add a subscript to indicate the angle at time instance .
- 9.
- radon transform used at an angle .
- 10.
- . radon transform applied to an image at angle . This results in a projection vector .
- 11.
- Forward projection matrix. The Matrix equivellent 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 .
- 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 equivelenet of applying the inverse radon transform in an unﬁltered mode to a projection (see item 12 above).
- 18.
- The psudo inverse of the forward projection matrix being multiplied by the projection vector . This is the matrix equivelenet 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.
- 25.
- Image ﬁdelity: " (inferred by the ability to discriminate between two images)" reference: The
relationship between image ﬁdelity and image quality by Silverstein, D.A.; Farrell, J.E

Sci-Tech Encyclopedia: Fidelity"The degree to which the output of a system accurately reproduces the essential characteristics of its input signal. Thus, high ﬁdelity in a sound system means that the reproduced sound is virtually indistinguishable from that picked up by the microphones in the recording or broadcasting studio. Similarly, a television system has a high ﬁdelity when the picture seen on the screen of a receiver corresponds in essential respects to that picked up by the television camera. Fidelity is achieved by designing each part of a system to have minimum distortion, so that the waveform of the signal is unchanged as it travels through the system. "

- 26.
- "image quality (inferred by the preference for one image over another)". Same reference as above
- 27.
- TE (Echo Time) "represents the time in milliseconds between the application of the 90 pulse and the peak of the echo signal in Spin Echo and Inversion Recovery pulse sequences." reference: http://www.fonar.com/glossary.htm
- 28.
- TR (Repetition Time) "the amount of time that exists between successive pulse sequences applied to the same slice." reference: http://www.fonar.com/glossary.htm

We will analyze the following HYPR algorithms:

- 1.
- Original HYPR
- 2.
- Wright-Hyper
- 3.
- I-HYPR

For each algorithm we show the schematic ﬂow chart and the mathematical description. At the end, we will run a simulation of each of the above 3 algorithm from the same set of images to describe the ﬁndings and compare the results.

The description of the simulation software and the results found are shown below in this report.

Before we discuss the Mathematical formulation, we need to better understand how to generate a set of projections from a well deﬁned image which we can express mathematically. For this purpose the following diagram shows the 3 possible cases in which projections can occur at.

We need to generate an analytical expression as a function of time of a simple object shape for the following cases.

In MRI, projection line refers to a line through the center of the k-space (a radial line). In the spatial domain one k-space projection line will be mapped to 2 projections (right and left projections). To clarify, the following diagram shows a disk at the center and the corresponding horizontal spatial projections. The projections are across the image and hence will generate a projection to the right and one to the left.

In CT scans however, the source is in one position and the detector is on another, hence we have forward projections only, not across the image as above. This is shown in the Kak and Slaney book

What all the above means, is that the Write-Huang paper, when they mention 8 projections, they mean 16 projections in the CT-sense. I.e. forward and backward projections. The simulation software written uses the CT projection convention. Hence use 16 projections to produce the same result as the Wright-Huag paper which uses 8 projections.

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

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

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

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

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

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

The following is the algorithm for original HYPR written in high level pseudo-code. For our implementation we will have the input to the HYPR algorithm as a set of original images . In other words, we will not start from the frequency domain, but from the spatial domain directly.

INPUT: N total number of projections

M number of time frames

set of 2D images {I}. The set size should be of size N*M

Theta_Initial,Theta_End -- initial and final angle (limited view)

OUTPUT: M number of HYPR 2D images

-- divide angles from Theta_Initial to Theta_Final into N*M equal divisions

Theta = linspace(Theta_Initial,Theta_Final,M*N)

FOR i from 1 to N LOOP

g(i) = radon(I(i),Theta(i)) -- generate forward projections

Filtered_Back_Projection(i) = iradon(g(i), Theta(i))

END LOOP

C = SUM( Filtered_Back_Projection ) -- composite image

Nr = N/M -- Number of projections per time frame

Starting_Projection = 1

FOR i from 1 to M LOOP

J(i) = Original_HYPR( C, g( Starting_Projection:Starting_Projection+Nr-1),

Theta( Starting_Projection:Starting_Projection+Nr-1) )

Starting_Projection = Starting_Projection + Nr

END LOOP

--

-- Original HYPR function

--

FUNCTION Original_HYPR( C, g, Theta) RETURNS HYPR_IMAGE

Nr = Length(g) --- Number of projections Per Frame

FOR i from 1 to Nr

gc(i) = radon(C, Theta(i))

PC(i) = iradon( gc(i), Theta(i), FILTER='NONE')

P(i) = iradon( g(i), Theta(i), FILTER='NONE')

Z(i) = P(i)/PC(i) -- element by element divison

END LOOP

mask = SUM(z)/Nr

RETURN( mask * C ) -- Element by element multiplication

END FUNCTION Original_HYPR

M number of time frames

set of 2D images {I}. The set size should be of size N*M

Theta_Initial,Theta_End -- initial and final angle (limited view)

OUTPUT: M number of HYPR 2D images

-- divide angles from Theta_Initial to Theta_Final into N*M equal divisions

Theta = linspace(Theta_Initial,Theta_Final,M*N)

FOR i from 1 to N LOOP

g(i) = radon(I(i),Theta(i)) -- generate forward projections

Filtered_Back_Projection(i) = iradon(g(i), Theta(i))

END LOOP

C = SUM( Filtered_Back_Projection ) -- composite image

Nr = N/M -- Number of projections per time frame

Starting_Projection = 1

FOR i from 1 to M LOOP

J(i) = Original_HYPR( C, g( Starting_Projection:Starting_Projection+Nr-1),

Theta( Starting_Projection:Starting_Projection+Nr-1) )

Starting_Projection = Starting_Projection + Nr

END LOOP

--

-- Original HYPR function

--

FUNCTION Original_HYPR( C, g, Theta) RETURNS HYPR_IMAGE

Nr = Length(g) --- Number of projections Per Frame

FOR i from 1 to Nr

gc(i) = radon(C, Theta(i))

PC(i) = iradon( gc(i), Theta(i), FILTER='NONE')

P(i) = iradon( g(i), Theta(i), FILTER='NONE')

Z(i) = P(i)/PC(i) -- element by element divison

END LOOP

mask = SUM(z)/Nr

RETURN( mask * C ) -- Element by element multiplication

END FUNCTION Original_HYPR

This mathematics of this algorithm will be presented by using the radon transform notation and not the matrix projection matrix notation. The conversion between the notation can be easily made by refering to the notation page at the end of this report.

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

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

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

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

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

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

The following is the algorithm for Wright HYPR written in high level pseudo-code. For our implementation we will have the input to the Wright-HYPR algorithm as a set of original images . In other words, we will not start from the frequency domain, but from the spatial domain directly.

INPUT: N total number of projections

M number of time frames

set of 2D images {I}. The set size should be of size N*M

Theta_Initial,Theta_End -- initial and final angle (limited view)

OUTPUT: M number of HYPR 2D images

-- divide angles from Theta_Initial to Theta_Final into N*M equal divisions

Theta = linspace(Theta_Initial,Theta_Final,M*N)

FOR i from 1 to N LOOP

g(i) = radon(I(i),Theta(i)) -- generate forward projections

Filtered_Back_Projection(i) = iradon(g(i), Theta(i))

END LOOP

C = SUM( Filtered_Back_Projection ) -- composite image

Nr = N/M -- Number of projections per time frame

Starting_Projection = 1

FOR i from 1 to M LOOP

J(i) = Wright_HYPR( C, g( Starting_Projection:Starting_Projection+Nr-1),

Theta( Starting_Projection:Starting_Projection+Nr-1) )

Starting_Projection = Starting_Projection + Nr

END LOOP

--

-- Wright HYPR function

--

FUNCTION Wright_HYPR( C, g, Theta) RETURNS HYPR_IMAGE

Nr = Length(g) --- Number of projections Per Frame

FOR i from 1 to Nr

gc(i) = radon(C, Theta(i))

PC(i) = iradon( gc(i), Theta(i), FILTER='NONE')

P(i) = iradon( g(i), Theta(i), FILTER='NONE')

END LOOP

mask = SUM(P) / SUM(PC) -- Element by element division

RETURN( mask * C ) -- Element by element multiplication

END FUNCTION Wright_HYPR

M number of time frames

set of 2D images {I}. The set size should be of size N*M

Theta_Initial,Theta_End -- initial and final angle (limited view)

OUTPUT: M number of HYPR 2D images

-- divide angles from Theta_Initial to Theta_Final into N*M equal divisions

Theta = linspace(Theta_Initial,Theta_Final,M*N)

FOR i from 1 to N LOOP

g(i) = radon(I(i),Theta(i)) -- generate forward projections

Filtered_Back_Projection(i) = iradon(g(i), Theta(i))

END LOOP

C = SUM( Filtered_Back_Projection ) -- composite image

Nr = N/M -- Number of projections per time frame

Starting_Projection = 1

FOR i from 1 to M LOOP

J(i) = Wright_HYPR( C, g( Starting_Projection:Starting_Projection+Nr-1),

Theta( Starting_Projection:Starting_Projection+Nr-1) )

Starting_Projection = Starting_Projection + Nr

END LOOP

--

-- Wright HYPR function

--

FUNCTION Wright_HYPR( C, g, Theta) RETURNS HYPR_IMAGE

Nr = Length(g) --- Number of projections Per Frame

FOR i from 1 to Nr

gc(i) = radon(C, Theta(i))

PC(i) = iradon( gc(i), Theta(i), FILTER='NONE')

P(i) = iradon( g(i), Theta(i), FILTER='NONE')

END LOOP

mask = SUM(P) / SUM(PC) -- Element by element division

RETURN( mask * C ) -- Element by element multiplication

END FUNCTION Wright_HYPR

This table shows the diﬀerence between the original HYPR and Wright-HYPR

This table shows the diﬀerence when we express the HYPR frame as function of original image and the operator only

A Matlab simulation software was written to implement the HYPR algorithms as described above. A User interface developed to make it easier for the user to run diﬀerent simulations using diﬀerent parameters. The appendix describes the user interface in more details.

The following are the initial results of the simulation. relative RMSE (root mean square error relative to the user image frame) was computed for the HYPR frame generated and the corresponding original image over the same time frame used to generate the HYPR frame (The average of the original image over the time frame was used if the time frame contained more than one snap shot of the user image). One time frame was used, and hence one HYPR frame was generated. The number of projections was increased, and the RMSE was calculated each time. The following table summarizes the ﬁnding.

Initial results shows that the original HYPR algorithm results in smaller RMSE than the Wright-Huang-HYPR algorithm using the Wright-Huang disk (WH disk), as the user image sampled for all the simulation. The WH disk is a disk of 50 pixels in diameter whose density is increased linearly from 1 to 128 units in time. It is also found that for later in time generated HYPR frames (i.e. HYPR frames generated from user images whose intensities were large), the RMSE is also larger than earlier HYPR images. This shows the larger the intensity of the original image, the less approximate is the HYPR image to the original image. This needs to be investigated more to see the reason why the HYPR algorithm is sensitive to intensity.

In addition to measuring the RMSE between each HYPR frame and corresponding image frame, the mean of each frame is calculated.

Another result found is that the Iterative HYPR algorithm using the original HYPR algorithm produces the least RMSE with 5 iterations. This shows that I-HYPR does improve the ﬁnal quality of the HYPR images. Convergence seems to be fast. After 4 iterations, the RMSE error did not change too much (please see table below). Most of the image improvement (as measured by lower RMSE) occurred in the ﬁrst 3 iterations.

The following table lists the RMSE and means found when running the original and Wright-Huang HYPR and with running I-HYPR for 5 iterations.

16 Frames were generated for each simulation. All simulations used 16 projections (or 8 per MRI deﬁnition) per a time frame, over full view (0-180 degrees) and used 16 time frames.

The following table shows the convergence of the I-HYPR. The RMSE is measured as more iterations are performed.

16 Frames were generated for each simulation. All simulations used 16 projections (or 8 per MRI deﬁnition) per a time frame, over full view (0-180 degrees) and used 16 time frames. Notice that converges seem to occur after 4 steps. (not counting the initial zero step to create the ﬁrst composite image).

This diagram shows the result of running the Original HYPR algorithm on the WH disk. Notice that the composite image intensity (red line below) is averaged over all time frames. This image shows the last HYPR frame and it is compared with the average of the user images over the same time frame. Notice that the HYPR frame is close in intensity to the corresponding user image, but some variations on the circle edges show up.

This is the same result as above, but for the Wright-Huang HYPR algorithm.

This is the same result for I-HYPR with 5 iterations

Notice the above spatial proﬁle. The HYPR image (black line) is closer the true line than the earlier results using original and WH HYPR. However, we notice an interesting artifact here, a Gibbs like phenomena around the edges of the disk shows up where the HYPR converges is not coming to exact at those places. This needs more investigation.

- 1.
- Implement iterative HYPR using the Wright-Huang version of HYPR as its iterative step. Compare results.
- 2.
- Run all the above simulation for user images with small object close to each others (small disks) and change intensity to simulate time changes. Compare quality of HYPR images reconstructed with those when the images are further apart.
- 3.
- Run all the above simulation for user image which moves in time and changes intensity as well. Use a DISK that moves in diﬀerent conﬁgurations.
- 4.
- Investigate the Gibb-like phenomena observed around the edges of the disk when using I-HYPR.
- 5.
- Investigate why original HYPR produces a better RMSE than Wright-Huang, while the temporal proﬁle of the last HYPR frame has more oﬀset from the true proﬁle as compared to the Wrigth-Huang proﬁle. This seems counter intituive. It is possible I am doing something wrong in the simulation?

The following is a screen shot of the Matlab simulation software.

- 1.
- 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
- 2.
- Iterative projection reconstruction of time-resolved images using HYPR by O'Halloran et.all
- 3.
- Time-Resolved MR Angiography With Limited Projections by Yuexi Huang1,and Graham A. Wright
- 4.
- GE medical PPT dated 6/6/2008
- 5.
- Book principles of computerized Tomographic imaging by Kak and Staney