## SVD for imaging systems with discrete rotational symmetry |

Optics Express, Vol. 18, Issue 24, pp. 25306-25320 (2010)

http://dx.doi.org/10.1364/OE.18.025306

Acrobat PDF (2189 KB)

### Abstract

The singular value decomposition (SVD) of an imaging system is a computationally intensive calculation for tomographic imaging systems due to the large dimensionality of the system matrix. The computation often involves memory and storage requirements beyond those available to most end users. We have developed a method that reduces the dimension of the SVD problem towards the goal of making the calculation tractable for a standard desktop computer. In the presence of discrete rotational symmetry we show that the dimension of the SVD computation can be reduced by a factor equal to the number of collection angles for the tomographic system. In this paper we present the mathematical theory for our method, validate that our method produces the same results as standard SVD analysis, and finally apply our technique to the sensitivity matrix for a clinical CT system. The ability to compute the full singular value spectra and singular vectors could augment future work in system characterization, image-quality assessment and reconstruction techniques for tomographic imaging systems.

© 2010 Optical Society of America

## 1. Introduction

1. H. H. Barrett, J. N. Aarsvold, and T. J. Roney, “Null functions and eigenfunctions: tools for the analysis of imaging systems,” Prog. Clin. Biol. Res. **363**, 211–226 (1991). [PubMed]

3. A. K. Jorgensen and G. L. Zeng, “SVD-Based evaluation of multiplexing in multipinhole SPECT systems,” Int. J. Biomed. Imaging **2008**, 769195 (2008). [CrossRef]

6. G. Gullberg and G. Zeng, “A reconstruction algorithm using singular value decomposition of a discrete representation of the exponential radon transform using natural pixels,” IEEE Trans. Nucl. Sci. **41**(6), 2812–2819 (1994). [CrossRef]

7. S. Park and E. Clarkson, “Efficient estimation of ideal-observer performance in classification tasks involving high-dimensional complex backgrounds,” J. Opt. Soc. Am. A **26**(11), 59–71 (2009). [CrossRef]

8. S. Park, J. Witten, and K. Myers, “Singular vectors of a linear imaging system as efficient channels for the Bayesian ideal observer,” IEEE Trans. Med. Imaging **28**(5), 657–668 (2009). [CrossRef] [PubMed]

**H**that maps voxels in object space to detector pixels in image space. To minimize the error in this approximation, it is desirable to have a large number of voxels to sample object space causing the column dimension

*N*of the sensitivity matrix to be large. The row dimension of the sensitivity matrix is especially large in tomographic systems where

*M*is the product of the number detector elements

*P*and the number of collection angles

*J*. The consequent large dimensions of

**H**for tomographic medical imaging systems are problematic when using conventional SVD algorithms because they require adequate computer memory to compute and store

**HH**

^{†}or

**H**

^{†}

**H**.

## 2. The system operator and the symmetry operators

*f*(

**r**) of the three dimensional location vector

**r**, and

**R**

*is the matrix that describes a rotation by the angle*

_{j}*acting in object space is given by*

_{j}*Z*, the cyclic group of order

_{J}*J*. This group is represented in ℝ

^{3}by the set of matrices {

**I**,

**R**

_{1},...,

**R**

_{J−1}} which satisfy

**R**

_{j}**R**

*=*

_{k}**R**

_{j+k}and

**R**

_{1}by

_{1},..., 𝒯

_{J−1}}, which satisfy the same multiplication rules as the rotation matrices.

**m**given by

**m**= (

**p**,

*j*) = (

*p*,

_{x}*p*,

_{y}*j*). The vector

**p**identifies the location of a detector element in a given detector array, and the index

*j*identifies the detector array. These

*J*identical detector arrays are assumed to be spaced around the object at angles

*h*

**(**

_{p}**r**) is the detector sensitivity function for the detector located at position

**p**on the detector with

*j*= 0, then the detector sensitivity function for detector

**m**is given by

*h*

**(**

_{m}*r*) = 𝒯

_{j}h**(**

_{p}**r**). The mean data vector

**ḡ**has components given by which we will also write as

**ḡ**= ℋ

*f*. If the number of detector elements in one detector array is

*P*, then

**ḡ**is a concatenation of

*J P*-dimensional vectors

**ḡ**

*: Thus ℋ is an operator from some function space*

_{j}*U*to data space

*V*= ℝ

*. We represent this symbolically by ℋ :*

^{JP}*U*→

*V*. The decomposition of

**ḡ**given above corresponds to the decomposition of

*V*into an orthogonal direct sum of subspaces

*V*, all of which are

_{j}*P*-dimensional: A different decomposition of

*V*into subspaces arises from symmetry considerations and will be described below.

**ḡ**

*can be thought of as the result of an operator ℋ*

_{j}*acting on the object function: We will consider ℋ*

_{j}*to be an operator that takes an object function and gives a vector in the subspace*

_{j}*V*of the full data space In other words, the result of applying this operator is the vector Physically, the operator ℋ

_{j}*would be obtained by covering all of the apertures except the*

_{j}*j*one. We will assume that the range of ℋ

^{th}*is all of*

_{j}*V*. We may then write the action of the system operator ℋ on an object function as with ℋ

_{j}*:*

_{j}*U*→

*V*and

*range*(ℋ

*) =*

_{j}*V*.

_{j}**S**

*. The matrix*

_{j}**S**

_{1}is described by the equation and the matrices

**S**

*satisfy*

_{j}*=*

_{j}**S**

*ℋ*

_{j}_{0}𝒯

_{–j}. This equation says that we can get the projection on the

*j*detector array by rotating the object by an angle of

^{th}*detector array, and then moving this*

^{th}*P*-dimensional data vector back to the correct location in the

*JP*-dimensional vector

**ḡ**.

**S**

*ℋ = ℋ𝒯*

_{k}*. In mathematical terms we say the operator ℋ intertwines the group representation {ℐ, 𝒯*

_{k}_{1},...,𝒯

_{J}_{–1}} with the representation {

**I**,

**S**

_{1},...,

**S**

_{J}_{–1}}. Intuitively the symmetry condition tells us that rotating the object

*J*. Before we show this, we will introduce the inner products that we will use for object space and data space.

## 3. Inner products and adjoint operators

*S*is a support region for all object functions and is invariant under the rotations in the symmetry group. For example,

*S*could be a circular cylinder whose axis is the axis of the rotations in the group. In data space, we have the standard inner product We will have occasion to make use of complex object functions and data vectors, which is why the complex conjugate operation appears in these inner products.

^{†}, by the equation (ℋ

*f*,

**g**)

*= (*

_{V}*f*, ℋ

^{†}

**g**)

*. This operator is a map from data space to object space: ℋ*

_{U}^{†}:

*V*→

*U*. In terms of the detector sensitivity functions, this adjoint operator is given by This operator is also called the backprojection operator for the imaging system described by ℋ. We can use the the system operator given above to write Notice that we used the unitary property of the symmetry operators here. The adjoint operator

*V*. Symbolically, we write ℋ :

_{j}*V*→

*U*and

*nullspace*

## 4. Projection operators

*operators and data vectors into components that are eigenvectors of the*

_{j}**S**

*matrices. This process is essentially a Fourier decomposition with respect to the group of rotations. To this end define operators 𝒬*

_{j}*and*

_{k}**P**

*, for*

_{k}*k*= 0, 1,...,

*J*– 1 via the equations [9] These operators satisfy the idempotent property: 𝒬

*𝒬*

_{k}*= 𝒬*

_{k}*and*

_{k}**P**

_{k}**P**

*=*

_{k}**P**

*. They are also Hermitian:*

_{k}*Ũ*is the range of 𝒬

_{j}*, and*

_{j}*Ṽ*is the range of

_{j}**P**

*, then we have the corresponding orthogonal decompositions of object space and data space This suggests that we define the operators ℋ̃*

_{j}*by ℋ̃*

_{k}*= ℋ𝒬*

_{k}*=*

_{k}**P**

*ℋ, where the second equality here follows from the symmetry property of ℋ. These operators are maps from object space to data space with the following ranges and null spaces:*

_{k}*range*(ℋ̃

*) =*

_{k}*Ṽ*and

_{k}*nullspace*

*maps*

_{j}*Ũ*to

_{k}*Ṽ*and vanishes on

_{k}*Ũ*for

_{j}*j*≠

*k*.

*maps data space to object space and has the following properties:*

_{k}*nullspace*

*range*

*JP*×

*JP*matrix ℋℋ

^{†}has block diagonal form as well: Each block

*Ṽ*to

_{j}*Ṽ*and vanishes on

_{j}*Ṽ*for

_{k}*k*≠

*j*. This fact, together with the orthogonality of the decompositions of object and data space, implies that the SVD problem for the system operator reduces to finding the eigenvalues and eigenvectors for the

*J*matrices

*P*-dimensional space to itself and therefore the corresponding eigenvalue problems reduce to finding the eigenvalues for

*P*×

*P*matrices. The original

*JP*×

*JP*eigenvalue problem for the SVD of ℋ has now been reduced to a set of

*JP*×

*P*eigenvalue problems, which is more computationally tractable, especially for large values of

*J*.

## 5. Reduction of dimension for SVD

**P**

*ℋ𝒬*

_{k}*ℋ*

_{k}^{†}=

**P**

*ℋ ℋ*

_{k}^{†}, where we have used the symmetry property of ℋ and the orthogonal projection properties of

**P**

*. We may write this last form out as a triple sum: Now we let*

_{k}*q*=

*p*–

*l*and re-index the sum over

*l*to get Next we let

*r*=

*j*+

*p*–

*q*and re-index the sum over

*j*to arrive at By separating the exponential factor, we can separate the sums. The end result is

*Ṽ*, so we will write such a vector as

_{k}**P**

_{k}**g**in the eigenvalue equation:

*Ṽ*may be reproduced from its values on the first detector, we may assume that

_{k}**g**has the following form Now we look at the eigenvalue equation

*J*

^{2}is replaced with

*J*because

**P**

*has a 1/*

_{k}*J*factor. These last two vectors are equal iff

*P*×

*P*linear system since the operator maps

*V*

_{0}to

*V*

_{0}.

**g**

_{0}=

**Zg**so that

**Z**is a

*P*×

*M*matrix. With

**g**given as above, we also have

**g**=

**Z**

^{†}

**g**

_{0}. Our eigenvalue equation is now

**Z**, we get

*, we can write this as*

_{k}*J*(

**Z**ℋ

_{0}𝒬

*)(*

_{k}**Z**ℋ

_{0}𝒬

*)*

_{k}^{†}

**g**

_{0}=

*λ*

**g**

_{0}. Therefore, the eigenvalue problem for each

*k*is equivalent to the SVD of the operator

**Z**ℋ

_{0}𝒬

*. Finally, note that*

_{k}**Z**ℋ

_{0}𝒬

*=*

_{k}**Z**ℋ𝒬

*=*

_{k}**ZP**

*ℋ so that we could also do the SVD for these last two operators.*

_{k}*by*

_{k}*P*×

*P*matrix. The eigenvalues are real and non-negative. To get the singular vectors in the full data space, we use:

*f*in object space are then given by:

_{kl}*. This operator rotates the function through each of the*

_{k}*J*angles, multiplies each rotated function by a phase factor determined by

*k*, and then sums. The symmetry properties of the resulting singular functions are determined by

*k*.

## 6. Examples

**H**, for the simulated x-ray CT system was computed using a discrete-to-discrete forward model to map from object space

**U**to image space

**V**. The object space was discretized into

*N*= 256 × 256 voxels with an imposed circular field of view. The image space was measured with a

*P*= 32 detector element camera collecting at

*J*= 90 equally spaced angles over 360 degrees of rotation about the center of the field of view. As discussed in Section 5, the reduced dimension SVD therefore led to

*J*SVD computations of the

*P*×

*N*matrices

*JP*×

*N*matrix

**H**. In both cases, the algorithm used to compute singular data for either the reduced dimension or the standard dimension matrix was the SVD function in the MATLAB software package [10

10. E. Anderson, Z. Bai, and C. Bischof, *LAPACK Users’ Guide* (Society for Industrial Mathematics, 1999). [CrossRef]

**H**on this machine required 15.86 minutes using the reduced dimension SVD algorithm and 473.55 minutes using the standard dimension SVD algorithm. Our reduced dimension algorithm was therefore approximately 30 times faster than the standard dimension algorithm at computing the full SVD of the simulated sensitivity matrix. The speed of our algorithm could be further improved by trivially distributing the SVD computation of the

**A**matrices onto the cores of a multi-core machine such as a Mac Pro computer.

_{k}*u*can be expressed mathematically in terms of the inner product and the Kronecker delta as [2] Shown in Table 1 and Table 2 are the magnitudes of the inner products for selected singular vectors that were computed using the reduced dimension SVD algorithm. The magnitudes of the inner products of singular vectors with the same index are unity while the magnitudes of the inner products of singular vectors with differing indices are close to zero. The error in the inner products is due to the estimations we made when simulating

_{n}**H**for the proof-of-concept system. In our forward model, we used subvoxels and monte carlo techniques when quantifying the values of

**H**at each collection angle. Consequently, our forward model does not strictly exhibit the property of discrete rotational symmetry discussed in Section 2. To further assess how this modeling error manifested itself in the SVD data, we compared singular vectors produced by the reduced dimension and standard dimension techniques. For this analysis we chose singlets

*u*

_{1},

*u*

_{6}and

*u*

_{15}. We subtracted the singlet output by the reduced dimension SVD from the singlet output by the standard SVD resulting in an image of the error for the singlet. The results of this subtraction are shown in Fig. 3. The structure of these images shows that that there are contributions from additional singular vectors within each singlet that are a consequence of not preserving discrete rotational symmetry in our forward model. Note that the magnitude of the inner product between two singular vectors of differing indices that came from a reduced dimension SVD computation of the same

**A**matrix (

_{k}*e.g. u*

_{1}and

*u*

_{6}) are in fact zero to within the numerical precision of the computer.

14. S. Steckmann, M. Knaup, and M. Kachelrieß, “High performance cone-beam spiral backprojection with voxel-specific weighting,” Phys. Med. Biol. **54**(12), 3691–3708 (2009). [CrossRef] [PubMed]

*R*= 600 mm and a center of rotation to center of detector array distance of

_{F}*R*= 450 mm. The slice thickness for each of the 12 detector column was 0.5 mm and each detector row had 672 channels of 0.9 mm width. This detector geometry was modeled as a rectangular array with

_{D}*P*= 8064 detector elements. The detector array collected data at J = 180 equally spaced angles over 360 degrees. Object space was centered at the center of rotation and discretized using

*N*= 512 × 512 × 5 square voxels of dimension 0.6741 mm. This choice of voxel parameters ensured that a ray traveling from the source through the edge of object space was observed by the detector array. A circular field of view was also imposed within each 512 × 512 slice through object space to preserve rotational symmetry. An illustration of the geometry for the modeled system is shown in Fig. 4.

*J*SVD computations of

*P*×

*N*matrices. The singular data for a given

*P*×

*N*matrix was computed using a power methods algorithm [3

3. A. K. Jorgensen and G. L. Zeng, “SVD-Based evaluation of multiplexing in multipinhole SPECT systems,” Int. J. Biomed. Imaging **2008**, 769195 (2008). [CrossRef]

## 7. Conclusion

## Acknowledgments

## References and links

1. | H. H. Barrett, J. N. Aarsvold, and T. J. Roney, “Null functions and eigenfunctions: tools for the analysis of imaging systems,” Prog. Clin. Biol. Res. |

2. | H. Barrett and K. Myers, |

3. | A. K. Jorgensen and G. L. Zeng, “SVD-Based evaluation of multiplexing in multipinhole SPECT systems,” Int. J. Biomed. Imaging |

4. | Y. Hsieh, G. Zeng, and G. Gullberg, “Projection space image reconstruction using strip functions to calculate pixels more natural for modeling the geometric response of the SPECT collimator,” IEEE Trans. Med. Imaging |

5. | G. Zeng and G. Gullberg, “An SVD study of truncated transmission data in SPECT,” IEEE Trans. Nucl. Sci. |

6. | G. Gullberg and G. Zeng, “A reconstruction algorithm using singular value decomposition of a discrete representation of the exponential radon transform using natural pixels,” IEEE Trans. Nucl. Sci. |

7. | S. Park and E. Clarkson, “Efficient estimation of ideal-observer performance in classification tasks involving high-dimensional complex backgrounds,” J. Opt. Soc. Am. A |

8. | S. Park, J. Witten, and K. Myers, “Singular vectors of a linear imaging system as efficient channels for the Bayesian ideal observer,” IEEE Trans. Med. Imaging |

9. | M. Hamermesh, |

10. | E. Anderson, Z. Bai, and C. Bischof, |

11. | J. Aarsvold, “Multiple-pinhole transaxial tomography: a model and analysis,” Ph. D. Dissertation (University of Arizona, 1993). |

12. | J. Aarsvold and H. Barrett, “Symmetries of single-slice multiple-pinhole tomographs,” in |

13. | P. Varatharajah, B. Tankersley, and J. Aarsvold, “Discrete models and singular-value decompositions of single-slice imagers with orthogonal detectors,” in |

14. | S. Steckmann, M. Knaup, and M. Kachelrieß, “High performance cone-beam spiral backprojection with voxel-specific weighting,” Phys. Med. Biol. |

15. | E. Isaacson and H. Keller, |

**OCIS Codes**

(110.2960) Imaging systems : Image analysis

(110.3000) Imaging systems : Image quality assessment

(110.6955) Imaging systems : Tomographic imaging

**ToC Category:**

Imaging Systems

**History**

Original Manuscript: August 13, 2010

Revised Manuscript: November 11, 2010

Manuscript Accepted: November 11, 2010

Published: November 19, 2010

**Citation**

Eric Clarkson, Robin Palit, and Matthew A. Kupinski, "SVD for imaging systems with discrete rotational symmetry," Opt. Express **18**, 25306-25320 (2010)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-24-25306

Sort: Year | Journal | Reset

### References

- H. H. Barrett, J. N. Aarsvold, and T. J. Roney, "Null functions and eigenfunctions: tools for the analysis of imaging systems," Prog. Clin. Biol. Res. 363, 211-226 (1991). [PubMed]
- H. Barrett, and K. Myers, Foundations of Image Science (John Wiley and Sons, 2004).
- A. K. Jorgensen, and G. L. Zeng, "SVD-based evaluation of multiplexing in multipinhole SPECT systems," Int. J. Biomed. Imaging 2008, 769195 (2008). [CrossRef]
- Y. L. Hsieh, G. L. Zeng, and G. T. Gullberg, "Projection space image reconstruction using strip functions to calculate pixels more "natural" for modeling the geometric response of the SPECT collimator," IEEE Trans. Med. Imaging 17(1), 24-44 (1998). [CrossRef] [PubMed]
- G. Zeng, and G. Gullberg, "An SVD study of truncated transmission data in SPECT," IEEE Trans. Nucl. Sci. 44(1), 107-111 (1997). [CrossRef]
- G. Gullberg, and G. Zeng, "A reconstruction algorithm using singular value decomposition of a discrete representation of the exponential radon transform using natural pixels," IEEE Trans. Nucl. Sci. 41(6), 2812-2819 (1994). [CrossRef]
- S. Park, and E. Clarkson, "Efficient estimation of ideal-observer performance in classification tasks involving high-dimensional complex backgrounds," J. Opt. Soc. Am. A 26(11), 59-71 (2009). [CrossRef]
- S. Park, J. M. Witten, and K. J. Myers, "Singular vectors of a linear imaging system as efficient channels for the bayesian ideal observer," IEEE Trans. Med. Imaging 28(5), 657-668 (2009). [CrossRef] [PubMed]
- M. Hamermesh, Group Theory and its Application to Physical Problems (Dover Publications, 1989).
- E. Anderson, Z. Bai, and C. Bischof, LAPACK Users’ Guide (Society for Industrial Mathematics, 1999). [CrossRef]
- J. Aarsvold, "Multiple-pinhole transaxial tomography: a model and analysis," Ph. D. Dissertation (University of Arizona, 1993).
- J. Aarsvold, and H. Barrett, "Symmetries of single-slice multiple-pinhole tomographs," in Conference Record of the 1996 IEEE NSS/MIC (IEEE, 1997), vol. 3, pp. 1673-1677. [CrossRef]
- P. Varatharajah, B. Tankersley, and J. Aarsvold, "Discrete models and singular-value decompositions of single-slice imagers with orthogonal detectors," in Conference Record of the 1998 IEEE NSS/MIC (IEEE, 1999), vol. 2, pp. 1184-1188.
- S. Steckmann, M. Knaup, and M. Kachelrieß, "High performance cone-beam spiral backprojection with voxel-specific weighting," Phys. Med. Biol. 54(12), 3691-3708 (2009). [CrossRef] [PubMed]
- E. Isaacson, and H. Keller, Analysis of Numerical Methods (Dover Publications, 1994).

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

### Supplementary Material

» Media 1: MOV (40 KB)

» Media 2: MOV (55 KB)

» Media 3: MOV (57 KB)

» Media 4: MOV (66 KB)

» Media 5: MOV (53 KB)

» Media 6: MOV (56 KB)

» Media 7: MOV (71 KB)

» Media 8: MOV (59 KB)

» Media 9: MOV (70 KB)

» Media 10: MOV (59 KB)

« Previous Article | Next Article »

OSA is a member of CrossRef.