## Decomposition of images and objects into measurement and null components

Optics Express, Vol. 2, Issue 6, pp. 254-260 (1998)

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

Acrobat PDF (108 KB)

### Abstract

We derive and discuss an algorithm for separating measurement-and null-space components of an object with respect to a given imaging system. This algorithm requires a discrete-to-discrete approximation of a typically continuous-to-discrete imaging system, and problems associated with such an approximation are examined. Situations where knowledge of the measurement and null spaces is crucial for analyzing imaging systems are discussed. We conclude with two examples demonstrating the usefulness of this algorithm, even in the discrete-to-discrete approximation.

© Optical Society of America

## 1. Introduction

1. D.R. Gilland, B.M.W. Tsui, C.E. Metz, R.J. Jaszczak, and J.R. Perry, “An evaluation of maximum likelihood-EM reconstruction for SPECT by ROC analysis,” J. Nucl. Med. **33**, 451–457 (1992). [PubMed]

2. X-L. Xu, J-S. Liow, and S.C. Strother, “Iterative algebraic reconstruction algorithms for emission computed tomography: A unified framework and its application to positron emission tomography,” Med. Phys. , **20**, 1675–1684 (1993). [CrossRef] [PubMed]

3. S.C. Moore, S.P. Kijewski, S.P. Muller, and B.L. Holman, “SPECT Image Noise Power: Effects of Nonstationary Projection Noise and Attenuation Compensation,” J. Nucl. Med. **29**, 1704–1709, (1988). [PubMed]

4. H.H. Barrett, D.W. Wilson, and B.M.W. Tsui, “Noise properties of the EM algorithm: I. Theory,” Phys. Med. Biol. **39**, 833–846 (1994). [CrossRef] [PubMed]

## 2. Background and derivation

*, a continuous function of the position*

**f****r**, represents the radioisotope distribution within the patient, and the continuous-to-discrete imaging operator, 𝛨, has kernel

*h*

_{m}(

**r**) equal to the probability that a photon emitted at object position

**r**will be detected in detector pixel,

*m*. Elements of the discretely collected data vector,

**g**, are calculated as

*m*= 1,…,

*M*

*given*

**f****g**. However, it is clearly impossible to fully determine the infinite-dimensional object vector

*from an*

**f***M*-dimensional data vector

**g**. The usual method for circumventing this problem is to divide the continuous

*into a number of voxels approximately equal to the number of elements in*

**f****g**. We rewrite equation (1) as

**f**, the discrete approximation of

*has elements*

**f***n*= 1,…,

*N*, and

**H**, the matrix approximation of 𝛨, has elements

*v*

_{n}represents the

*n*

^{th}voxel in object space and the subscript on

**g**

_{a}indicates that this is only an approximation of the true projection process.

### 2.1 Null space, measurement space and singular-value decomposition

**H**would be invertible. However, it is generally found that even the approximation of the true imaging system is singular. This singularity results from the manner that tomographic systems sample the object space and the well-known fact that high frequencies are transferred very poorly in tomography. One way that the singularity of

**H**has been explored is with singular-value decomposition (SVD) theory. It can be shown that any

*m*×

*n*matrix,

**H**, can be decomposed as

*M*×

*M*

**U**and the

*N*×

*N*

**V**are nonsingular orthonormal matrices, and the

*M*×

*N*

**D**has zero elements only at the first

*p*locations on the diagonal.

**U**and

**V**are orthonormal, their column vectors,

**u**

_{i}and

**v**

_{i}, form an orthonormal basis for the projection space and discretized object space, respectively. Equation 3 can be rewritten, in terms of

**u**

_{i}and

**v**

_{i}, as

**v**

_{i}∙

**f**is the inner product of the

*i*

^{th}object-space eigenvector and the object. Thus, we see that the first

*p*

**v**

_{i}contribute to the projection operation and the rest do not. We can then therefore define an

*n*×

*n*measurement-space projector,

**M**, with

*p*rows equal to the first

*p*

**v**

_{i}, and

*N*-

*p*rows equal to the

**0**vector. The null-space projector,

**N**, is simply the identity matrix minus

**M**. The

**Mf**operation is then the projection of

**f**onto the measurement space, and

**Nf**is the projection of

**f**onto the null space.

**H**is

*d*

_{ii}for

*d*

_{ii}greater than zero, and zero otherwise. We note that the measurement- and null-space projectors can be written as

**M**=

**H**

^{+}

**H**and

**N**=

**I**-

**H**

^{+}

**H**.

**H**

^{+}

**g**

_{a}is a least-squares solution to the problem stated by equation 3. However, a modeling error has been introduced by the discretization of equation 1, and a least-squares solution of the discretized problem does not imply a least-squares solution of the real continuous-to-discrete system. Although

**H**

^{+}is not the least-squares solution to the real imaging problem, it is very advantageous be able to calculate it when analyzing imaging systems or reconstruction algorithms. We shall present two examples demonstrating uses of

**H**

^{+}in section 3.

### 2.4 Derivation of the algorithm

**H**be stored in memory, and typical three-dimensional tomographic systems have

**H**matrices with a size of several gigabytes.

**f̂**of the discretized object is normally reconstructed by an iterative process (except in the case of the filtered-backprojection algorithm, which assumes a very simplified form of

**H**). Iterative algorithms allow the reconstruction process to be done “on the fly” so that the

**H**matrix does not need to be stored in memory. They also allow the possibility of stopping when a result of suitable accuracy has been achieved, thus saving the computation time required to reach convergence.

5. O.N. Strand, “Theory and methods related to the singular-function expansion and Landweber’s iteration for integral equations of the first kind,” SIAM J. Numer. Anal. , **11**, 798–825 (1974). [CrossRef]

**Nf**and

**Mf**for any given

**f**and

**H**. We shall present only an outline of the algorithm’s derivation here. (A full derivation is contained in [7

7. H.H. Barrett and D.W. Wilson, “An algorithm to calculate null and measurement spaces,” http://www.radiology.arizona.edu/~fastspec/nul_spc.pdf

**I**+

**X**]

^{-1}can be expanded in the Neuman series

**f**:

## 3. Examples

### 3.1 Example 1: FASTSPECT imaging resolution

9. H.H. Barrett, “The FASTSPECT imaging system,” http://www.radiology.arizona.edu/~fastspec.

10. I. Pang, “Wedge phantom for resolution measurement,” http://www.radiology.arizona.edu/~fastspec/wedge.htm.

**H**), 2) an improper model of the of

*h*

_{i}(

**r**) when estimating

**H**or 3) the reconstruction process. Since we physically measured

**H**, we were fairly confident that it accurately characterized a voxelized FASTSPECT system, but problems sometimes arise when we do the measurements. Additionally, scatter from within the wedge phantom was not included, and we had measured

**H**with larger than usual voxels. Also, as will demonstrate in example two, our reconstruction algorithm is sensitive to smoothing parameters.

**H**. Modeling errors, approximations due to voxelization, and reconstruction had now been removed from the problem. The only remaining effect was the measured response of the voxelized imaging system, so if the phantom’s measurement space showed quality similar to the reconstruction given in Fig. 1, we would assume that the reduced resolution was a function of the actual imaging system. It is observed in Figs. 2(b) and 2(c) that the measurement space of the digitized phantom looks much the same as the reconstructions, with most of the high-frequency information contained in the null space. From this we concluded that the poor resolution was an inherent system problem, and redesigned DCI FASTSPECT with smaller (two millimeter) pinholes. The measurement space of the digitized wedge phantom with respect to the redesigned system is given in Fig. 3(a). Fig. 3(b) shows reconstructed projections of the physical wedge phantom. It is seen that the system resolution is improved and that the null-space algorithm accurately predicts this improvement.

### 3.2 Example 2: Reconstruction using the null-space algorithm.

**Hf̂**, and the data collected by the imaging system,

**g**. This algorithm quickly arrives at an

**f̂**in good agreement with the projection data, but it in no way suppresses the null-space components of

**f̂**. Unless the reconstructed image is artificially smoothed, it will be of extremely poor visual quality. For this reason we perform some type of local smoothing of

**f̂**at each iteration of the reconstruction process. This diminishes the agreement between the estimate and the projection data, but increases the apparent visual quality of the image.

**Hf̂**to reach maximum agreement with

**g**. The resulting

**f̂**is given in Fig. 4. This estimate has very good agreement with the data collected by the imaging system but, clearly, little can be determined from it – including the phantom imaged in the experiment.

**Mf̂**(Fig. 5(a)) and

**Nf̂**(Fig. 5(b)) with respect to the

**H**used for reconstruction. It is seen that

**Mf̂**conveys more information about the object being imaged than

**f̂**alone, and that

**Nf̂**only served to confound the original estimate. We draw from this experiment the somewhat obvious conclusion that a reconstruction algorithm should not be allowed to proceed heedless of the null space.

## 4. Conclusions

## References and links

1. | D.R. Gilland, B.M.W. Tsui, C.E. Metz, R.J. Jaszczak, and J.R. Perry, “An evaluation of maximum likelihood-EM reconstruction for SPECT by ROC analysis,” J. Nucl. Med. |

2. | X-L. Xu, J-S. Liow, and S.C. Strother, “Iterative algebraic reconstruction algorithms for emission computed tomography: A unified framework and its application to positron emission tomography,” Med. Phys. , |

3. | S.C. Moore, S.P. Kijewski, S.P. Muller, and B.L. Holman, “SPECT Image Noise Power: Effects of Nonstationary Projection Noise and Attenuation Compensation,” J. Nucl. Med. |

4. | H.H. Barrett, D.W. Wilson, and B.M.W. Tsui, “Noise properties of the EM algorithm: I. Theory,” Phys. Med. Biol. |

5. | O.N. Strand, “Theory and methods related to the singular-function expansion and Landweber’s iteration for integral equations of the first kind,” SIAM J. Numer. Anal. , |

6. | K.M. Hanson, “Bayesian and Related Methods in Image Reconstruction from Incomplete Data,” in |

7. | H.H. Barrett and D.W. Wilson, “An algorithm to calculate null and measurement spaces,” http://www.radiology.arizona.edu/~fastspec/nul_spc.pdf |

8. | P. Gilbert, “Iterative Methods for the Three-Dimensional Reconstruction of an Object from Projections,” J. Theor. Biol. |

9. | H.H. Barrett, “The FASTSPECT imaging system,” http://www.radiology.arizona.edu/~fastspec. |

10. | I. Pang, “Wedge phantom for resolution measurement,” http://www.radiology.arizona.edu/~fastspec/wedge.htm. |

**OCIS Codes**

(100.6950) Image processing : Tomographic image processing

(110.6960) Imaging systems : Tomography

**ToC Category:**

Focus Issue: Tomographic image reconstruction

**History**

Original Manuscript: December 9, 1997

Published: March 16, 1998

**Citation**

Donald Wilson and Harrison Barrett, "Decomposition of images and objects into measurement and null components," Opt. Express **2**, 254-260 (1998)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-2-6-254

Sort: Journal | Reset

### References

- D.R. Gilland, B.M.W. Tsui, C.E. Metz and R.J. Jaszczak, J.R. Perry, "An evaluation of maximum likelihood-EM reconstruction for SPECT by ROC analysis," J. Nucl. Med. 33, 451-457 (1992). [PubMed]
- X-L. Xu, J S. Liow and S.C. Strother, "Iterative algebraic reconstruction algorithms for emission computed tomography: A unified framework and its application to positron emission tomography," Med. Phys., 20, 1675- 1684 (1993). [CrossRef] [PubMed]
- S.C. Moore, S.P. Kijewski, S.P. Muller and B.L. Holman, "SPECT Image Noise Power: Effects of Nonstationary Projection Noise and Attenuation Compensation," J. Nucl. Med. 29, 1704-1709, (1988). 4. H.H. Barrett, D.W. Wilson and B.M.W. Tsui, "Noise properties of the EM algorithm: I. Theory," Phys. Med. Biol. 39, 833-846 (1994). [PubMed]
- H.H. Barrett, D.W. Wilson and B.M.W. Tsui, "Noise properties of the EM algorithm: I. Theory," Phys. Med. Biol. 39, 833-846 (1994). [CrossRef] [PubMed]
- O.N. Strand, "Theory and methods related to the singular-function expansion and Landwebers iteration for integral equations of the first kind," SIAM J. Numer. Anal., 11, 798-825 (1974). [CrossRef]
- K.M. Hanson, "Bayesian and Related Methods in Image Reconstruction from Incomplete Data," in Image Recovery, Theory and Application, Henry Stark, ed. (Academic, Orlando, Fla. 1987)
- H.H. Barrett and D.W. Wilson, "An algorithm to calculate null and measurement spaces," http://www.radiology.arizona.edu/~fastspec/nul_spc.pdf
- P. Gilbert, "Iterative Methods for the Three-Dimensional Reconstruction of an Object from Projections," J. Theor. Biol. 36, 105-117 (1972). [CrossRef] [PubMed]
- H.H. Barrett, "The FASTSPECT imaging system," http://www.radiology.arizona.edu/~fastspec.
- I. Pang, "Wedge phantom for resolution measurement," http://www.radiology.arizona.edu/~fastspec/wedge.htm.

## 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.

OSA is a member of CrossRef.