In the context of reconstruction of MRI images from K-space data sets, the following are two desirable properties which are difficult to achieve simultaneously: High spatial resolution and High temporal resolution.
The first requires longer acquisition time while the second requires less time be consumed acquiring each image. Hence the inherent conflict in achieving both simultaneously. One possible remedy is to undersample image acquisition (along a radial or other trajectories such as Cartesian) which results in speed up of data acquisition, hence improving the temporal resolution. Next, an appropriate image reconstruction method is applied to the acquired data which attempts to compensate for some of the effects of the image undersampling.
Due to undersampling, streaking artifacts will be present in the final image. These streaking artifacts become more visible the larger the undersampling. Radial undersampling is a preferred method of acquisition compared to using other trajectories such as Cartesian: "the aliasing artifacts from radial under-sampling usually appear as streaks, which are visually less distracting than the wrap-around artifacts obtained with Cartesian under-sampling."[?]
Mathematically, the problem of image reconstruction from undersampled K-space data is an inverse problem: "Mathematically, the reconstruction problem from sparse K-space samples is an ill-posed inverse problem with infinitely many solutions"[?]
Highly constrained backprojection reconstruction (HYPR) was introduced recently for the reconstruction of radially undersampled MRI images. HYPR is able to reconstruct these images with less visible artifacts while maintaining good SNR[?].
The method starts with the construction of one composite image made up from filtered backprojections of a large number of projections (each of these projections is the inverse Fourier transform of a corresponding radial lines from the K-space data set. This follows from the central slice theorem).
Since the composite image is made up of individual images collected over longer time period, it posses good spatial characteristics. In addition, its SNR is larger since a larger amount of images data is contained in it. However, the composite image temporal characteristics are poor since it combines images that were generated from varying time instances into a single image.
A weight image is then constructed from the ratio of a small number of unfiltered backprojections obtained from the original projections and from the composite image. Since the weight image is constructed from images that span a smaller time window than the case is with the composite image, the weight image posses good temporal characteristics. However, its spatial characteristic is poor due to the undersampling effect in the original data.
HYPR now generates a new image by multiplying the weight image with the composite image. This results in a HYPR image which combines the best characteristics found in the composite image and in the weight image, resulting in an image with good SNR, good temporal and spatial characteristics and with limited artifacts. The above process is repeated for the next HYPR image reconstruction until all the K-space data set is processed. The relationship between the weight image and the composite image and the resulting HYPR image is summarized in the following diagram.
Let \(I\) be a 2D image. Let \(R_{\theta \left ( t\right ) }\) be the forward projection matrix (implemented as Radon transformation in practice) which when applied to the image \(I\) at some projection angle \(\theta \) at time \(t\) will result in a 1D projection image \(P\left ( t\right ) \).
Since the image itself will change with time (blood flows, etc...) therefore we need to associated a time index \(t\) as well with the 2D image itself. Hence we write \(I\left ( t\right ) \) from now on, to indicate the 2D image at time \(t\). Notice that the image as a whole does not move (relative to fixed initerial frame of reference) but the image content can change as described above.
To avoid confusion in what follows, I will use the notation of \(\circledast \) to indicate a multiplication between 2 matrices element by element and will use \(\ast \) to indicate the normal matrix by matrix or matrix by vector multiplication. And will use \(/.\) to indicate division between vectors or matrices being done element by element.
We will analyze the following HYPR algorithms: Original HYPR, Wright-Hyper, I-HYPR and a new variation of I-HYPR based on on the Wright-HYPR. This new variation is different from I-HYPR since in I-HYPR the iteration is made on a composite image based on the Original HYPR construction while in this variation, the composite image is constructed is with the Write-HYPR algorithm. We will call this variation IW-HYPR.
For each algorithm we show the schematic flow chart and the mathematical description based on matrix formulation and an algorithm for the implementation. At the end, we will run a simulation of each of the above 4 algorithm from the same set of images and attempt to describe the finding and compare the results.
In the mathematical formula we derive an expression of the HYPR image as a function of the set of original images \(I\left ( t\right ) \) and the set of the forward projection matrices \(R_{\theta \left ( t\right ) }\).
Before we discuss the Mathematical formulation, we need to better undertstand 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 occure at.
We need to generate an analytical expression as a function of time of a simple object shape for each of the following 3 cases.
Given a set of projections \(P\left ( t\right ) \) over a number of time frames say \(M\) and assuming there are \(N\) projections generated per a time frame, the data set will consists of \(M\ast N\) projections in total. Within one time frame, a number of projections can be collected. These projections within each time frame will generate one HYPR frame using the HYPR algorithm. A time frame can have only one projection, but normally a time frame will have much larger number of projections.
For example, assuming there are \(M=10\) time frames, and \(N=20\) projections generated by time frame, then the total number of projections is \(200.\) In this case, we will obtain \(10\) HYPR frame images at the end.
HYPR starts by constructing a composite image from all the projections in the data set (\(200\) in the above example). This is done by computing the filtered backprojection of all the \(M\ast N\) projections into one image called \(C\).
Next, we process the projections from each time frame. For each time frame we generate a set of \(N\) projections from \(C\) by perform a radon tranform (forward projection) on \(C\) at the same angle corresponding to the projection being processed. This will generate set of projections called the \(P_{c}\) projections. Hence there will be \(N\) such projections per each time time.
Next the ratio of the each \(P\) projection over the corresponding \(P_{c}\) projection is found. This ratio is done pixel by pixel. Then each such ratio is multiplied by the composite image \(C\) generating a a new set of size \(M\ast N\) of weighted composite images. Now to generate a HYPR frame image, the average of these weighted compsite images is taken. The average is carried over each time frame at a time. Hence the first \(N\) weighted compsite images will generate one HYPR frame, and the next \(N\) weighted compsite images will generate the next HYPR frame and so on, resulting in \(M\) HYPR frames.
The following diagram illustrates the above where we used the original HYPR algorithm. Later on, we describe in detailes the different HYPR algorithms variations. In the following diagram we show 4 projections per time frame and a total of 3 time frames. This results in 3 HYPR frames.
The projection \(P\left ( t\right ) \) is obtained by applying forward projection on the image \(I\left ( t\right ) \), Hence we write\[ P\left ( t\right ) =R_{\theta \left ( t\right ) }\left [ I\left ( t\right ) \right ] \] Next, the set of \(\left \{ P\left ( t\right ) \right \} \) are combined and a filtered backprojection is applied to the result to generate a composite image \(C\)\[ C=\mathbf {FBP}\left [ {\displaystyle \sum \limits _{i=1}^{M\ast N}} R_{\theta \left ( t_{i}\right ) }\left [ I\left ( t_{i}\right ) \right ] \right ] \] Where \(\mathbf {FBP}\) is operator for the filtered backprojection.
The composite image \(C\) can be written as\begin {align*} C & ={\displaystyle \sum \limits _{i=1}^{N}} A_{i}^{+}\ast P_{i}\\ & ={\displaystyle \sum \limits _{i=1}^{N}} A_{i}^{+}\ast \left ( A_{i}\ast f_{i}\right ) \end {align*}
Where \(N\) in the number of projections. Now applying forward projection to \(C\) at angle \(\theta _{i}\) will generate a projection \(P_{c_{i}},\)Hence \[ P_{c_{i}}=A_{i}\ast C \] Let \(z_{i}\) be the ratio \(P_{i}/.P_{c_{i}}\) where this division is being carried out element by element between the two projections. Hence\begin {align*} z_{i} & =P_{i}/.P_{c_{i}}\\ & =\left ( A_{i}\ast f_{i}\right ) \ /.\ \left ( A_{i}\ast C\right ) \\ & =\left ( A_{i}^{T}\ast f_{i}\right ) \ /.\ \left ( A_{i}\ast \left ( {\displaystyle \sum \limits _{k=1}^{N}} A_{k}^{+}\ast \left ( A_{k}\ast f_{k}\right ) \right ) \right ) \end {align*}
Now apply unfiltered backprojection on the above projection ratio to obtain an unfiltered 2D image and then mutiply that with the composit image to obtain a HYPR frame \(F_{i}\)\begin {align*} F_{i} & =C\circledast \left ( A_{i}^{T}\ast z_{i}\right ) \\ & =\left ( {\displaystyle \sum \limits _{k=1}^{N}} A_{k}^{+}\ast \left ( A_{k}\ast f_{k}\right ) \right ) \circledast \left ( A_{i}^{T}\ast \left [ \left ( A_{i}^{T}\ast f_{i}\right ) \ /.\ \left ( A_{i}\ast \left ( {\displaystyle \sum \limits _{k=1}^{N}} A_{k}^{+}\ast \left ( A_{k}\ast f_{k}\right ) \right ) \right ) \right ] \right ) \end {align*}
Hence HYPR image is\begin {align} I & =\frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} F_{i}\nonumber \\ & =\frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} \left \{ \left ( {\displaystyle \sum \limits _{k=1}^{N}} A_{k}^{+}\ast \left ( A_{k}\ast f_{k}\right ) \right ) \circledast \left ( A_{i}^{T}\ast \left [ \left ( A_{i}^{T}\ast f_{i}\right ) \ /.\ \left ( A_{i}\ast \left ( {\displaystyle \sum \limits _{k=1}^{N}} A_{k}^{+}\ast \left ( A_{k}\ast f_{k}\right ) \right ) \right ) \right ] \right ) \right \} \tag {1} \end {align}
We see that the above expression for HYPR image depends only on \(f_{i}\) and \(A_{i}\).
Therefore, given a set of images \(f_{i}\) and a set of projection angles \(\theta _{i}\) we can compute the forward projection matrices \(A_{i}\) (analytically we can do this for simple shapes such as a disk rotating around a unit circle for example). And once the set of matrices \(A_{i}\) is computed, equation (1) could then be computed to obtain a HYPR image.
We can then generate the same HYPR images by performing backprojection (filtered an unfiltered) using the Fourier transform method. The algorithm for the backprojection (filtered) is known and given on page 62 of Kak and Slaney book which I will post a scan of. Or we could simply use the radon/iradon for the implementaion of \(A_{i}\),\(A_{i}^{T},A_{i}^{+}\), where \(A_{i}\) corresponds to applying radon on image \(f\) at angle \(\theta _{i},\) and \(A_{i}^{T}\) corresponds to applying iradon on \(g_{i}\)with filter ’none’ and \(A_{i}^{+}\)corresponds to applying iradon on projection \(g_{i}\) with a specified filter.
The projection \(g_{i}\) is obtained by applying forward projection on the image \(f\) at time \(i.\)Hence we write\[ g_{i}=A_{i}\ast f_{i}\] Where in the above the 2D image \(f_{i}\) needs to be first converted to a 1D vector as was done in the first assignment by folding the 2D image into a 1D column vector.
Let \(A_{i}^{+}\) be an psudoinverse of \(A_{i}\) which performs a filtered backprojection when applied on a projection vector \(g_{i}\) and let \(A_{i}^{T}\) be the transpose of the matrix \(A_{i}\) which performs an unfiltered backprojection on \(g_{i}\).
The composite image \(C\) can be written as\begin {align*} C & ={\displaystyle \sum \limits _{i=1}^{N}} A_{i}^{+}\ast g_{i}\\ & ={\displaystyle \sum \limits _{i=1}^{N}} A_{i}^{+}\ast \left ( A_{i}\ast f_{i}\right ) \end {align*}
where \(N\) in the number of projections. Now, each unfiltered backprojection is\begin {align*} P_{i} & =A_{i}^{T}\ast g_{i}\\ & =A_{i}^{T}\ast \left ( A_{i}\ast f_{i}\right ) \end {align*}
Now applying forward projection to \(C\) at angle \(\theta _{i}\) will generate a projection, call it \(r_{i}\). Hence \[ r_{i}=A_{i}\ast C \] Now applying an unfiltered backprojection on \(r_{i}\) will result in \(P_{c_{i}}\) hence\begin {align*} P_{c_{i}} & =A_{i}^{T}\ast r_{i}\\ & =A_{i}^{T}\ast \left ( A_{i}\ast C\right ) \end {align*}
Let \(z_{i}\) be the ratio \(P_{i}/.P_{c_{i}}\) where this division is being carried out element by element between the two images. Therefore\begin {align*} z_{i} & =P_{i}/.P_{c_{i}}\\ & =\left ( A_{i}^{T}\ast \left ( A_{i}\ast f_{i}\right ) \right ) \ /.\ \left ( A_{i}^{T}\ast \left ( A_{i}\ast C\right ) \right ) \\ & =\left ( A_{i}^{T}\ast \left ( A_{i}\ast f_{i}\right ) \right ) \ /.\ \left ( A_{i}^{T}\ast \left ( A_{i}\ast \left ( {\displaystyle \sum \limits _{k=1}^{N}} A_{k}^{+}\ast \left ( A_{k}\ast f_{k}\right ) \right ) \right ) \right ) \end {align*}
Hence one HYPR frame \(i\) is\begin {align*} F_{i} & =z_{i}\circledast C\\ & =\left ( A_{i}^{T}\ast \left ( A_{i}\ast f_{i}\right ) \right ) \ /.\ \left ( A_{i}^{T}\ast \left ( A_{i}\ast \left ( {\displaystyle \sum \limits _{k=1}^{N}} A_{k}^{+}\ast \left ( A_{k}\ast f_{k}\right ) \right ) \right ) \right ) \circledast \left ( {\displaystyle \sum \limits _{k=1}^{N}} A_{k}^{+}\ast \left ( A_{k}\ast f_{k}\right ) \right ) \end {align*}
Hence HYPR image is\begin {align} I & =\frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} F_{i}\nonumber \\ & =\frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} \left \{ \left ( A_{i}^{T}\ast \left ( A_{i}\ast f_{i}\right ) \right ) \ /.\ \left ( A_{i}^{T}\ast \left ( A_{i}\ast \left ( {\displaystyle \sum \limits _{k=1}^{N}} A_{k}^{+}\ast \left ( A_{k}\ast f_{k}\right ) \right ) \right ) \right ) \circledast \left ( {\displaystyle \sum \limits _{k=1}^{N}} A_{k}^{+}\ast \left ( A_{k}\ast f_{k}\right ) \right ) \right \} \tag {1} \end {align}
The difference between the original HYPR and Wright-HYPR can be seen in the following simplified diagram. We see than in the original HYPR, the ratio \(P/PC\) is performed on the 1-D projections, then the unfiltered backprojection is applied on the resulting 1-D set of images. In the Wright-HYPR algorithm, the unfiltered backprojection is first applied to the set of the 1-D projections and then the ratio is performed on the result 2-D set of images.
This is an iterative method where the composite image itself is updated and a new HYPR image determined with the hope of obtaining a better HYPR image (how to determine how many iterations? need to read more on this) as more iterations are performed. Let the number of iterations by \(M\) therefore this algorithm will generate \(C_{1},C_{2},\cdots ,C_{M}\) composite images.
This is the mathematical formulation for I-HYPR
The projection \(P_{i}\) is obtained by applying forward projection on the image \(f\) at time \(i.\)In the original HYPR, \(g_{i}\) is what is refered to \(P_{i}\) projection. Therefore\[ P_{i}=A_{i}\ast f_{i}\] The composite image \(C\) at iteration (1) can be written as\[ C_{1}={\displaystyle \sum \limits _{i=1}^{N}} A_{i}^{+}\ast P_{i}\] Where \(N\) in the number of projections. Now applying forward projection to \(C\) at angle \(\theta _{i}\) will generate a projection \(P_{c_{i}}\). Hence \[ P_{c(1)_{i}}=A_{i}\ast C_{1}\] Let \(z_{i}\) be the ratio \(P_{i}/.P_{c_{i}}\) where this division is being carried out element by element between the two projections. Hence\begin {align*} z_{i} & =P_{i}/.P_{c(1)_{i}}\\ & =P_{i}\ /.\ \left ( A_{i}\ast C_{1}\right ) \end {align*}
Now apply unfiltered backprojection on the above projection ratio to obtain an unfiltered 2D image and then mutiply that with the composit image to obtain a HYPR frame \(F_{i}\)\begin {align*} F_{i} & =C_{1}\circledast \left ( A_{i}^{T}\ast z_{i}\right ) \\ & =C_{1}\circledast \left ( A_{i}^{T}\ast \left [ P_{i}\ /.\ \left ( A_{i}\ast C_{1}\right ) \right ] \right ) \end {align*}
Hence HYPR image at interation (1) is\begin {align} I_{1} & =\frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} F_{i}\nonumber \\ & =\frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} C_{1}\circledast \left ( A_{i}^{T}\ast \left [ P_{i}\ /.\ \left ( A_{i}\ast C_{1}\right ) \right ] \right ) \tag {1} \end {align}
Now use \(I_{1}\) as the composite image for the next iteration, we obtain\[ C_{2}=I_{1}\] Hence\[ P_{c(2)_{i}}=A_{i}\ast C_{2}\] And \begin {align*} z_{i} & =P_{i}/.P_{c(2)_{i}}\\ & =P_{i}\ /.\ \left ( A_{i}\ast C_{2}\right ) \end {align*}
Now apply unfiltered backprojection on the above projection ratio to obtain an unfiltered 2D image and then mutiply that with the composit image to obtain a HYPR frame \(F_{i}\)\begin {align*} F_{i} & =C_{2}\circledast \left ( A_{i}^{T}\ast z_{i}\right ) \\ & =C_{2}\circledast \left ( A_{i}^{T}\ast \left [ P_{i}\ /.\ \left ( A_{i}\ast C_{2}\right ) \right ] \right ) \end {align*}
Hence HYPR image at interation (2) is\begin {align} I_{2} & =\frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} F_{i}\nonumber \\ & =\frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} C_{2}\circledast \left ( A_{i}^{T}\ast \left [ P_{i}\ /.\ \left ( A_{i}\ast C_{2}\right ) \right ] \right ) \tag {1} \end {align}
But \(C_{2}=I_{1}=\frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} C_{1}\circledast \left ( A_{i}^{T}\ast \left [ P_{i}\ /.\ \left ( A_{i}\ast C_{1}\right ) \right ] \right ) \) hence the above becomes\begin {align*} I_{2} & =\frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} \left ( \frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} C_{1}\circledast \left ( A_{i}^{T}\ast \left [ P_{i}\ /.\ \left ( A_{i}\ast C_{1}\right ) \right ] \right ) \right ) \circledast \\ & \left ( A_{i}^{T}\ast \left [ P_{i}\ /.\ \left ( A_{i}\ast \left ( \frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} C_{1}\circledast \left ( A_{i}^{T}\ast \left [ P_{i}\ /.\ \left ( A_{i}\ast C_{1}\right ) \right ] \right ) \right ) \right ) \right ] \right ) \end {align*}
But \(C_{1}={\displaystyle \sum \limits _{i=1}^{N}} A_{i}^{+}\ast P_{i}\), hence the above becomes\begin {align*} I_{2} & =\frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} \left ( \frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} \left ( {\displaystyle \sum \limits _{i=1}^{N}} A_{i}^{+}\ast P_{i}\right ) \circledast \left ( A_{i}^{T}\ast \left [ P_{i}\ /.\left ( A_{i}\ast \left ( {\displaystyle \sum \limits _{i=1}^{N}} A_{i}^{+}\ast P_{i}\right ) \right ) \right ] \right ) \right ) \circledast \\ & \left ( A_{i}^{T}\ast \left [ P_{i}\ /.\ \left ( A_{i}\ast \left ( \frac {1}{N}{\displaystyle \sum \limits _{i=1}^{N}} \left ( {\displaystyle \sum \limits _{i=1}^{N}} A_{i}^{+}\ast P_{i}\right ) \circledast \left ( A_{i}^{T}\ast \left [ P_{i}\ /.\ \left ( A_{i}\ast \left ( {\displaystyle \sum \limits _{i=1}^{N}} A_{i}^{+}\ast P_{i}\right ) \right ) \right ] \right ) \right ) \right ) \right ] \right ) \end {align*}
Repeate this processes by setting \(C_{3}=I_{2}\) and generate \(I_{3}\). Continue untill \(I_{M}\) where \(M\) is number of iterations or untill some other convergence criteria is achived.
Simplied version of the schematic diagram
Need to finish the mathematics of this version. (TODO)