Image fidelity: " (inferred by the ability to discriminate between two images)"
reference: The relationship between image fidelity 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 fidelity 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 fidelity 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. "
We will analyze the following HYPR algorithms:
For each algorithm we show the schematic flow 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 findings 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 defined 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 \(R\) notation and not the matrix projection matrix \(H\) notation.
The projection \(s_{t}\) is obtained by applying radon transform \(R\) on the image \(I_{t}\) at some angle \(\phi _{t}\)\[ s_{t}=R_{\phi _{t}}\left [ I_{t}\right ] \]
When the original object image does not change with time then we can drop the subscript \(t\) from \(I_{t}\) and just write \(s_{t}=R_{\phi _{t}}\left [ I\right ] \)
The composite image \(C\) is found from the filtered back projection applied to all the \(s_{t}\)\[ C={\displaystyle \sum \limits _{i=1}^{N}} R_{\phi _{t_{i}}}^{f}\left [ s_{t_{i}}\right ] \] Notice that the sum above is taken over \(N\) and not over \(N\). Next a projection \(s_{c}\) is taken from \(C\) at angle \(\phi \) as follows\[ s_{c_{t}}=R_{\phi _{t}}\left [ C\right ] \] The the unfiltered back projection 2-D image \(P_{t}\) is generated\[ P_{t}=R_{\phi _{t}}^{u}\left [ s_{t}\right ] \] And the unfiltered back projection 2-D image \(P_{c_{t}}\) is found\[ P_{c_{t}}=R_{\phi _{t}}^{u}\left [ s_{c_{t}}\right ] \] Then the ratio of \(\frac {P_{t}}{P_{c_{t}}}\) is summed and averaged over the time frame and multiplied by \(C\) to generate a HYPR frame \(J\) for the time frame\(.\)Hence for the \(k^{th}\) time frame we obtain\begin {align*} J_{k} & =C\ \left ( \frac {1}{N_{p}}{\displaystyle \sum \limits _{i=1}^{N_{p}}} \frac {P_{t_{i}}}{P_{c_{t_{i}}}}\right ) \\ & =\frac {1}{N_{p}}\left ( {\displaystyle \sum \limits _{i=1}^{N}} R_{\phi _{t_{i}}}^{f}\left [ s_{t_{i}}\right ] \right ) \ {\displaystyle \sum \limits _{j=1}^{N_{p}}} \frac {R_{\phi _{t_{j}}}^{u}\left [ s_{t_{j}}\right ] }{R_{\phi _{t_{j}}}^{u}\left [ s_{c_{t_{j}}}\right ] } \end {align*}
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 \(I_{i}\). 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
This mathematics of this algorithm will be presented by using the radon transform \(R\) notation and not the matrix projection matrix \(H\) notation. The conversion between the notation can be easily made by referring to the notation page at the end of this report.
The projection \(s_{t}\) is obtained by applying radon transform \(R\) on the image \(I_{t}\) at some angle \(\phi _{t}\)\[ s_{t}=R_{\phi _{t}}\left [ I_{t}\right ] \]
When the original object image does not change with time then we can drop the subscript \(t\) from \(I_{t}\) and just write \(s_{t}=R_{\phi _{t}}\left [ I\right ] \)
The composite image \(C\) is found from the filtered back projection applied to all the \(s_{t}\)\[ C={\displaystyle \sum \limits _{i=1}^{N}} R_{\phi _{t_{i}}}^{f}\left [ s_{t_{i}}\right ] \] Notice that the sum above is taken over \(N\) and not over \(N\). Next a projection \(s_{c}\) is taken from \(C\) at angle \(\phi \) as follows\[ s_{c_{t}}=R_{\phi _{t}}\left [ C\right ] \] The the unfiltered back projection 2-D image \(P_{t}\) is generated\[ P_{t}=R_{\phi _{t}}^{u}\left [ s_{t}\right ] \] And the unfiltered back projection 2-D image \(P_{c_{t}}\) is found\[ P_{c_{t}}=R_{\phi _{t}}^{u}\left [ s_{c_{t}}\right ] \]
Now the set of \(P_{t}\) and \(P_{c_{t}}\) over one time frame are summed the their ratio multiplied by \(C\) to obtain the \(k^{th}\) HYPR frame \begin {align*} J_{k} & =C\ \frac {{\displaystyle \sum \limits _{i=1}^{N_{p}}} P_{t_{i}}}{{\displaystyle \sum \limits _{i=1}^{N_{p}}} P_{c_{t_{i}}}}\\ & =C\ \frac {{\displaystyle \sum \limits _{i=1}^{N_{pr}}} R_{\phi _{t}}^{u}\left [ s_{t}\right ] }{{\displaystyle \sum \limits _{i=1}^{N_{pr}}} R_{\phi _{t}}^{u}\left [ s_{c_{t}}\right ] } \end {align*}
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 \(I_{i}\). 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
This table shows the difference between the original HYPR and Wright-HYPR
| HYPR frame \(J\) in original HYPR | \(\frac {1}{N_{pr}}C\ {\displaystyle \sum \limits _{i=1}^{N_{pr}}} \frac {P_{i}}{P_{c_{i}}}\) |
| HYPR frame \(J\) in Wright HYPR | \(C\ \frac {{\displaystyle \sum \limits _{i=1}^{N_{pr}}} P_{i}}{{\displaystyle \sum \limits _{i=1}^{N_{pr}}} P_{c_{i}}}\) |
This table shows the difference when we express the HYPR frame as function of original image \(I\) and the operator \(R\) only
| HYPR frame \(J\) in original HYPR | \(\frac {1}{N_{pr}}\left ( {\displaystyle \sum \limits _{i=1}^{N}} R_{\phi _{i}}^{+}\left [ R_{\phi _{i}}\left [ I_{i}\right ] \right ] \right ) \ {\displaystyle \sum \limits _{j=1}^{N_{pr}}} \frac {R_{\phi _{j}}^{\ast }\left [ R_{\phi _{j}}\left [ I_{j}\right ] \right ] }{R_{\phi _{j}}^{\ast }\left [ R_{\phi _{j}}\left [ {\displaystyle \sum \limits _{\mu =1}^{N_{pr}}} R_{\phi _{\mu }}^{+}\left [ R_{\phi _{\mu }}\left [ I_{\mu }\right ] \right ] \right ] \right ] }\) |
| HYPR frame \(J\) in Wright HYPR | \(=\left ( {\displaystyle \sum \limits _{i=1}^{N}} R_{\phi _{i}}^{+}\left [ g_{i}\right ] \right ) \ \frac {{\displaystyle \sum \limits _{j=1}^{N_{pr}}} R_{\phi _{j}}^{\ast }\left [ R_{\phi _{j}}\left [ I_{j}\right ] \right ] }{{\displaystyle \sum \limits _{j=1}^{N_{pr}}} R_{\phi _{j}}^{\ast }\left [ R_{\phi _{j}}\left [ {\displaystyle \sum \limits _{\mu =1}^{N}} R_{\phi _{\mu }}^{+}\left [ g_{\mu }\right ] \right ] \right ] }\) |
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 different simulations using different 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 finding.
| Number of projections | Original HYPR | Wright-Huang HYPR |
| \(16\) | \(2.172299\) | \(1.897816\) |
| \(128\) | \(0.884654\) | \(0.932832\) |
| \(256\) | \(0.779807\) | \(0.817493\) |
| \(512\) | \(0.743475\) | \(0.767085\) |
| \(700\) | \(1.241370\) | \(1.167796\) |
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 final 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 first 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 definition) 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 definition) 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 first 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 profile. 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.
The following is a screen shot of the Matlab simulation software.