## Abstract

Mathematical phantoms developed to synthesize realistic complex backgrounds such as those obtained when imaging biological tissue, play a key role in the quantitative assessment of image quality for medical and biomedical imaging. We present a modeling framework for the synthesis of realistic tissue samples. The technique is demonstrated using radiological breast tissue. The model employs a two-component image decomposition consisting of a slowly, spatially varying mean-background and a residual texture image. Each component is synthesized independently. The approach and results presented here constitute an important step towards developing methods for the quantitative assessment of image quality in medical and biomedical imaging, and more generally image science.

©1997 Optical Society of America

## 1. Introduction

The major goals of research in medical and biomedical imaging are to create better imaging systems, more accurate reconstructions, and to develop methods of image processing and analysis that utilize the most important information present in an image for accurate and timely diagnosis of disease. Realistic numerical models of human tissue and medical imaging systems are key components to achieving this goal. This paper specifically addresses an approach to the modeling of biological tissue, where the technique is demonstrated for radiological breast tissue, also referred to as mammographic tissue.

Biological tissues most often appear as highly complex backgrounds against which a physician may search for lesions and other specific features indicative of malignancy. In mammography, such features may be microcalcifications and architectural distortions. We demonstrate that a key component in the synthesis of medical and biomedical background images is to recognize that they can be thought of as a textured background with structures on both large and fine scales that may be synthesized independently.

The method for synthesis decomposes an image into a slowly, spatially varying
mean-background and a residual texture image.^{1,2} We shall refer to the slowly, spatially varying
mean-background as the mean background. We propose to model the mean background as a
stochastic process known as the lumpy background.^{3-5} The statistical properties of the lumpy background are well specified and
will be summarized in section 3.^{4,5} Thus, we shall also compare the statistical
properties of these two types of backgrounds. Some preliminary analysis of the first
and second order statistical properties of these two types of backgrounds will be
presented in section 6.

In the development of a mathematical phantom for complex backgrounds, an important question is whether the residual texture image can also be synthesized and statistically characterized. We shall demonstrate in this paper the synthesis of such texture backgrounds and also present results of a preliminary investigation of the first and second order statistical properties of the residual texture background and its synthesis.

The combination of two synthesis frameworks with known statistical properties, to represent the mean background and the residual texture background, respectively, leads to a digital tissue-phantom that may find important use in the development of methods for the quantitative assessment of image quality in medical and biomedical imaging. The overall synthesis approach may also find applications in various other areas of image science.

## 2. Synthesis of biological tissue via a two-component model

Some investigations into the statistics of texture backgrounds, specifically
mammographic backgrounds, have been conducted.^{6–9} It has been suggested that various classes of images, including mammograms,
have power spectra of the form 1/f^{α}. For mammograms, estimated
values of α in the range of 1.5 to 2 have been reported.^{8–9} A power-law spectrum exponent between 1.5 and 2 indicates that mammograms are
not fractals. A two-dimensional fractal would yield an exponent greater than 2. This
finding further suggests that such backgrounds cannot be synthesized using fractals.10 Therefore, while some investigations have demonstrated that
different mammographic tissue types can be classified according to their estimated
fractal dimension,7 some of these investigations have further demonstrated that
the addition of another parameter significantly yields improved classification.^{9,11} An additional complication with modeling
biological tissue as a fractal is that it is difficult to accurately estimate a
fractal dimension from digitized data.^{12}

The indication that the power spectra of mammograms and various natural images follow
some power law is significant and we shall relate our findings to these investigations.^{8–9} It has been demonstrated, however, that the power spectrum of a statistical
complex background is not a complete descriptor of the required background
statistics to predict human observer performance in various detection tasks:
specifically, two studies demonstrate that two sets of images with equal power
spectra, yet having Fourier spectra that differ in phase, yield different
detectability performance and thus require different predictive mathematical models.^{6,8} An ensemble of images with the same power spectrum
as that of another ensemble of images, but with a Fourier spectrum of random phase,
was obtained by filtering various realizations of white noise with the desired power
spectrum.

Knowing that complex backgrounds such as mammograms cannot be modeled either as a
fractal, or as filtered white noise to yield a given ensemble power spectrum, we
propose an alternative approach to modeling such backgrounds. The model is
established from knowledge of the anatomy of such tissues and their radiographic appearance.^{13} Radiographic contrast in mammography arises from differing attenuation
between tissues that comprise the breast. The breast is made essentially of a
mixture of fatty tissue, which appears dark on radiographs, connective and
epithelial tissues which produce bright radiographic appearances also referred to as
mammographic densities, and prominent ducts which yield cord-like structures or a
beaded appearance.^{11–13}

Overall, a set of mammographic radiographs from the same type of breast tissue may be described as a stochastic process with fairly large scale structures that account for mammographic densities on black backgrounds, and smaller scale structures that give the tissue the appearance of texture. The model we propose is thus based on the decomposition of such complex backgrounds into two components: a mean background (i.e. the slowly, spatially varying component defined earlier as the mean background) that accounts for large scale inhomogeneities in the background, and a texture image that characterizes the fluctuations of that image around the mean background.

A realization of the mean background is typically obtained by convolving a sample of
a mammographic sample image with a two-dimensional Gaussian kernel.^{1,2} An example of the sample image and the resulting
blurred image are shown in Figure 1a and 1b, respectively. The sample image is a 256 x 256 pixel
section extracted from a mammogram from the database of N. Karssemeijer of
University Hospital Nijmegen, The Netherlands.^{14} We propose to model the mean background as a lumpy background that we know to
be a wide-sense stationary stochastic process.^{4,15,16} We thus assume that the stochastic process
describing an ensemble of mean backgrounds extracted from several sub-images of a
set of mammograms is wide-sense stationary as well. This assumption will be fully
investigated in future work as it involves computation of an ensemble
autocorrelation function over a large number of images and the careful study of its properties.^{15} Such a validation is beyond the scope of this paper.

A sample of the residual texture image is obtained by subtracting the mean background from the original image. The residual texture image corresponding to the sample image shown in Figure 1a is shown in Figure 1c. We propose to synthesize the texture image using a four-layer pyramid framework described in section 4. Finally, we propose that various linearly weighted sums of the two model components, the lumpy background and the synthesized texture, yield mammograms with typical radiographic characteristics.

Current literature provides some support for this proposed model. First, Hunt and
Cannon demonstrated that natural scenes can be decomposed into intensity
fluctuations around a nonstationary ensemble mean where the mean is estimated by
blurring a typical ensemble member.^{1} Our model differs from that of Hunt
and Cannon as we hypothesize that the mean background arises from an underlying
wide-sense stationary process. Moreover, as it relates specifically to mammography,
several investigations by Byng and colleagues suggest that at least two parameters
are required to characterize mammograms: one parameter to describe the distribution
of breast tissue density as reflected by the brightness of the mammogram and another
parameter to characterize the texture.^{11} The mean background and texture components of the proposed decomposition are
reminiscent of the first and second parameters in Byng’s model,
respectively.

## 3. Synthesis of the slowly varying mean-background.

We propose to model the mean background as a stochastic process known as the lumpy
background which has been developed to specifically account for spatially varying
backgrounds in medical images as a result, for example, of anatomical structures.^{3–5} In the case of mammography, the mean background may account for the relative
amount of fat and densities in the breast tissue.

The lumpy background, detailed in Rolland (1990), was devised to be mathematically
tractable for computing signal to noise ratio predictions for various medical
imaging tasks.^{3,5,15} The lumpy background assumes wide-sense
stationarity, that is stationarity over the ensemble of images, where the
autocorrelation function is only a function of the shift variable **r**.
The second assumption is that the background autocorrelation function is a Gaussian
function. The power spectrum W(**ρ**) is then defined as the
Fourier transform of the autocorrelation function and is given by

where ρ is the 2D frequency variable in the Fourier domain conjugate to
**r**, r_{b} is the correlation length of the autocorrelation
function, and W(0) is the value of the power spectrum at zero frequency that we
refer to as the lumpiness measure.

Two models of the lumpy backgrounds were presented in Rolland (1990).^{4} The
first approach consisted of superimposing 2D Gaussian functions, referred to as
Gaussian blobs, on a constant background of strength B_{0}. To keep the
mathematics simple, we assumed 2D Gaussian blobs of constant amplitude
b_{0}/π${{\mathrm{r}}_{\mathrm{b}}}^{2}$ and constant half-width
r_{b}. In this case, the lumpy background can be shown to be
mathematically specified as

where r_{j} is a random variable uniformly distributed over the background
area, and K is the number of Gaussian blobs in the background. Rolland further
showed that to yield a Gaussian autocorrelation function, the number of Gaussian
blobs K must also be a random variable with the mean of K equal to its variance. We
thus chose K to be Poisson distributed for this condition to be
satisfied.^{4} A measure of lumpiness in the background is given by

where K̅/A_{d} is the mean number of Gaussian functions per
mm^{2}.

By varying the background lumpiness, one can simulate various types of tissue inhomogeneities. We hypothesize that the width of the Gaussian blobs can be chosen to simulate existing correlation lengths of the tissue densities in various types of tissues.

## 4. Synthesis of the residual texture image

To yield a complete mathematical model of biological tissues, one may also account
for the residual texture shown in Figure 1c. Over the last year, we developed a four-layer
pyramid transform to synthesize texture backgrounds that can be best described by
looking at individual components: the pyramid transform, the image decomposition,
the histogram matching procedure, and finally the texture synthesis.^{17} The synthesis of the texture of biological tissues using the proposed
framework is presented here for the first time. We shall now describe each component
of the framework.

*Pyramid transform* The proposed algorithm for the synthesis of the
residual texture is based on a four-layer steerable pyramid transform. One layer of
the pyramid is depicted in Figure 2. Layers are connected by a factor-of-two decimation
of the image.^{18} Within each layer, the image is filtered by a set of bandpass filters and
followed by a set of orientation filters. The algorithm adopts a 4 (scales) x 4
(orientations: 0 degree, 45 degree, 90 degree and 135 degree) steerable filter bank.^{19–22} Details of the filters employed for the synthesis will be detailed elsewhere.^{17}

*Image Decomposition* The texture image is processed through the left
hand side of the pyramid transform shown in Figure 2. It is represented in Figure 2 as an input to the pyramid in the upper left corner.
In parallel, a realization of uniformly distributed white noise, referred thereafter
as white noise, is also processed by the same pyramid transform, that is, it is also
fed independently to the pyramid transform in the upper left corner. The role of the
white noise image is to provide a starting point for the synthesis.

*Histogram matching at multiple scales* After decomposition of a
texture sample and a realization image of white noise, the histograms of the subband
images (i.e. output images of the filters on the left hand side of the pyramid) of
the texture image and of the noise image are matched.^{23} Histogram matching is an image processing technique, specifically a point
operation, which modifies a candidate image so that its histogram matches that of a
model image.^{24–25}

*Image Synthesis* The histogram-matched noise subband-images obtained
at multiple scales are then recombined according to the right hand side of the
pyramid transform shown in Figure 2. This process yields a synthetic image such as that
shown in Figure 2c. If another realization of white noise is processed
instead, the synthesis yields another realization of the synthesized image. Such an
example is shown in Figure 2d. In our implementation, the set of filters used in
the decomposition and the reconstruction stages forms a quadrature-mirror filter bank.^{18}

In summary, the synthesis process yields a synthetic texture image that represents
one realization from a theoretical ensemble of “equal texture”
images (i.e. a textured image from an ensemble of images differs from another image
of the same ensemble pixel to pixel, yet their respective textures visually mimic
each other). For each new realization of the white noise image, a new realization of
the synthetic texture is obtained as desired. The development of the algorithm was
inspired from a previous approach to texture synthesis.^{19}

## 5. A proposed mathematical phantom

The synthesis of an ensemble of images M_{i}(x,y) according to the described
mathematical phantom can be established using an adaptive linear combination of
realizations from the two model components: a realization of a lumpy background
component denoted as L_{i}(x,y) and a realization of the synthesized texture
component denoted as T_{i}(x,y). The resulting synthesized image will then
be given by

where β ranges from 0 to 1. Such a combination will allow us to span a
wide range of tissue types with relative amounts of lumpy backgrounds and texture
backgrounds. We hypothesize that by choosing the parameters of the lumpy background
(i.e. correlation length and lumpiness values) and the types of texture to
synthesize, various tissue types as described by Wolfe, for example, can be synthesized.^{13} On a more theoretical basis, one can also study a wide range of combination
of such backgrounds by varying β and the parameters associated with each
component. Such a framework may naturally find application to other types of images
beside mammographic tissue.

## 6. Preliminary investigation of the model’s first- and second-order statistics

A useful mathematical phantom of biological tissue requires that task performance be
predicted accurately for some typical specific tasks.^{16} To address this issue, the ultimate test will be to conduct a set of
psychophysical studies using real images and mathematically simulated images. While
it has been shown that equalization of the first- and second-order statistics are
insufficient to predict detectability in complex backgrounds, we propose to conduct
a first investigation of such properties for the mammographic tissue, its proposed
two components, and their syntheses, to draw parallels in the equivalence of first-
and second-order statistics of the proposed synthesis model images compared to the
original images.

We shall take the greylevel histogram and the power spectrum estimates to represent
the first- and second-order statistical properties, respectively.^{15} The average histograms and power spectra for the five ensembles of images are
shown in Figure 4 and Figure 5, respectively. The procedure used to estimate a
power spectrum from an ensemble of images is the two dimensional extension of the
method due to Welch that we adapted to two-dimensional structured backgrounds.^{26} An extension of Welch’s method to two dimensional noise images was
detailed in Hanson.^{27} The extension to structured backgrounds involves computation of an ensemble
mean and subtraction of each image from the estimated ensemble-mean before
computation of the power spectrum. Rolland and Barrett used a similar method for
defining the power spectrum of the lumpy stochastic process, while the power
spectrum of the lumpy background was analytically, rather than numerically, computed.^{4–5} The procedure we adopted can be summarized in six steps:

- Given an ensemble of images, compute the ensemble mean
- Subtract from each image of the ensemble the ensemble mean to form a new ensemble set
- 3. Take the FFT of each image i from the ensemble to yield X
_{i}(m,n) - 4. Compute the normalized periodogram as |X
_{i}(m,n)|^{2}/(NxN) - 5. Compute the average periodogram over the ensemble of images
- 6. Plot the Log
_{10}of the power spectrum along the x-dimension.

In this preliminary investigation of the first and second-order statistical
properties, the power spectra were estimated using a limited number of images. All
ensembles of images had a total of eighteen images. For the mammogram power spectrum
estimation, the ensemble is formed of images from five mammograms where three to
four 256x256 overlapping sub-sample images were further extracted. Overlapping the
sub-sample images by half, in our case in both dimensions, was originally proposed
by Welch in order to reduce the variance of the power spectrum estimate.^{26}

The mean background was extracted from the mammographic images by convolving the
latter with a Gaussian kernel of standard deviation six pixels. For the lumpy
stochastic process, eighteen images were simulated of lumpiness value 10^{6}
counts^{2}/(sec^{2} pixels) (i.e. mean number of blobs was 400,
the strength of a blob b0 was 12800 counts/sec) and of correlation length rb equal
to 15 pixels. For the synthesized texture images, eighteen syntheses were obtained
using three of the residual texture images from the ensemble. The estimated power
spectra of mammographic backgrounds, its two proposed components, and its two model
components are shown in Figure 4.

Figure 4 demonstrates quasi-equivalent first order statistics for the two components of the mammogram images and their respective proposed models. The second order statistics represented as the power spectra are reported in Figure 5. While power spectra should be computed over a large number of images from a given image class for accurate estimations, the limited computations presented here may provide some indication of the similarity and differences experienced on a limited set of images for the respective ensembles and point to potential mismatch or agreement in the second order statistics for two ensemble sets.

The power spectra are two dimensional and, while they appear to be fairly symmetrical, the x-dimension was reported here because we did not want to assume isotropy, as would have been assumed by taking a radial average. So the x-profiles of the power spectra are reported and compared. For the texture images, the spectra are both band-pass with similar shapes and an apparent peak around 0.05 cycles/pixel. This indicates that the texture synthesis process is capturing similar second-order statistics as those existing in the residual texture images extracted from the mammograms.

The power spectra for the lumpy background and mean background are different in shape, yet they are both low-pass. They mainly differ in that the power spectrum of the lumpy background is narrower than that of the mean-background. One complication in modeling the mean background as a lumpy process is to insure that all mammogram images belong to a “same class” of tissue density. While we attempted to select images according to such a criteria, future work will involve working with pre-classified images from experts. The correlation length and lumpiness values of the lumpy background were chosen to simulate backgrounds that best resemble the mean background of the mammograms. A more extensive investigation may lead to a better description of the mean background using the lumpy stochastic process. Furthermore, while it seems a reasonable assumption to state that the underlying stochastic process leading to the mean background of the mammograms is stationary, this property has not yet been fully investigated.15 Future work will include a more extensive study of the potential to model the mean background component with the lumpy process, as well as perhaps a more complex decomposition. The encouraging texture results suggests that the proposed framework has the potential to yield a useful mathematical phantom.

## 7. Conclusion

We presented a two-stage framework for the synthesis of complex texture backgrounds
and demonstrated its application to the synthesis of mammographic tissues. A tissue
sample was described as the sum of a slowly, spatially varying mean-background (i.e.
mean-background) and a residual texture image. We proposed to synthesize the
mean-background using a stochastic process known as the lumpy background that one of
the authors and colleagues established in previous investigations.^{3–5} Most importantly, we have shown in this paper that the residual texture image
can be synthesized using a four-layer pyramid decomposition framework we developed
over the last year.

Some preliminary investigation of the first and second order statistics of various ensemble of images were presented. The first order statistics compared were found to be equivalent and the power spectra of the residual texture and synthesized texture images were found to be similar. While a more extensive investigation is needed, the power spectrum of the mean background suggests that modeling the mean-background may require a more complex model than a simple lumpy background. From an investigation in progress not reported here, we postulate, however, that the lumpy background will likely play an important role in perhaps a more complex decomposition scheme.

Mathematical models of human tissue find key applications in the optimization and assessment of imaging systems, the generation of an ensemble of images for psychophysical studies for image quality assessment, and the extraction of economical statistics for description of detection performance in realistic complex backgrounds. Work in progress in our laboratory includes extraction of such economical statistical descriptions. Finally, realistic mathematical models of texture medical images may find useful application in the development of diagnostic teaching tools.

## Acknowledgments

We thank Marie Kijewski for pointing us to key references on power spectra computation, Kyle Myers for stimulating discussions about this work, Ha Evans for sharing some of her clinical expertise with the radiological appearance of mammographic tissues, and Nico Karssemeijer for providing images. We thank Liyun Yu for his assistance with programming for this project. Finally, we acknowledge Art Burgess for his insistence to include first and second order statistics computations which have significantly contributed to the objective validation of the proposed model. This work was supported in part by the University of Central Florida small grant program.

## References and links

**1. **B.R. Hunt and T. M. Cannon, “Nonstationary assumptions for Gaussian models of images,” IEEE Trans. on Sys., Man, and Cybern. , 876–882 (1976).

**2. **R.N. Strickland and H.I. Hahn, “Wavelet transforms for detecting microcalcifications in mammograms,” IEEE Trans. on Med. Imaging **15**, 218–229 (1996). [CrossRef]

**3. **K.J. Myers, J.P. Rolland, H.H. Barrett, and R.F. Wagner, “Aperture optimization for emission imaging: effect of a spatially varying background,” J. Opt. Soc. Am. A. **7**, 1279–1293 (1990). [CrossRef] [PubMed]

**4. **J.P. Rolland, “Factors influencing lesion detection in medical imaging,” Ph.D. Dissertation, University of Arizona, (1990).

**5. **J.P. Rolland and H.H. Barrett, “Effect of random background inhomogeneity on observer detection performance,” J. Opt. Soc. Am. A. **9**, 649–658 (1992). [CrossRef] [PubMed]

**6. **M.G.A. Thomson and D.H. Foster, “Role of second- and third-order statistics in the discriminability of natural images,” J. Opt. Soc. Am. A. **14**(9), 2081–2090 (1997). [CrossRef]

**7. **C. Caldwell and M. Yaffe, “Fractal analysis of mammographic parenchemal pattern,” Phys. Med. Biol. **35**, 235–247 (1990). [CrossRef] [PubMed]

**8. **F.O. Bochud, F. R. Verdun, C. Hessler, and J.F. Valley, “Detectability on radiological images: the influence of anatomical noise,” Proc. SPIE **2436**, 156–165 (1995). [CrossRef]

**9. **B. Zheng, Y.H. Chang, and D. Gur, “Adpative computer-aided diagnosis scheme of digitized mammograms,” Acad. Radiol. **3** (10), 806–814 (1996). [CrossRef] [PubMed]

**10. **M.F. Barnsley. *Fractals Everywhere*. (Academic Press, San Diego, CA, 1988)

**11. **J.W. Byng, MJ. Yaffe, G.A. Lockwood, L.E. Little, D.L. Tritchler, and N.F. Boyd, “Automated analysis of mammographic densities and breast carcinoma risk,” Cancer **80**(1), 66–74 (1997). [CrossRef] [PubMed]

**12. **B. Dubuc, C.R. Carmes, C. Tricot, and S.W. Zucker, “The variation method: a technique to estimate the fractal dimension of surfaces,” Proc. SPIE **845**, 241–248 (1987). [CrossRef]

**13. **J.N. Wolfe, “Breast patterns as an index of risk for developing breast cancer,” Am. J.Roentgenol. **126**, 1130–1139 (1976).

**14. ** The Nijmegen database is available by anonymous FTP from ftp://figment.csee.usf.edu /pub/mammograms/nijmegen-images

**15. **A. Papoulis. *Probablity, Random Variables, and Stochastic Processes*. (Mc Graw-Hill, NY, 1991).

**16. **H.H. Barrett, J. Yao, J.P. Rolland, and K.J Myers, “Model observers for assessment of image quality,” Proc. Natl. Acad. Sci. USA **90**, 9758–9765 (1993). [CrossRef] [PubMed]

**17. **J.P. Rolland and L. Yu, “A four-layer pyramid framework for statistical texture synthesis,” (in preparation).

**18. **P.P. Vaidyanathan. *Multivariate systems and filter banks*. (Prentice Hall, NJ, 1993).

**19. **D.J. Heeger and J.R. Bergen, “Pyramid-based texture analysis/synthesis,” Compt. Graph.. , 229–238 (1995).

**20. **E.P. Simoncelli and E.H. Adelson, “Subband transforms”. *In Subbands Image Coding*, (Kluwer Academic Publishers, J.W. Woods, eds., MA 1991).

**21. **E.P. Simoncelli, W.T. Freeman, E.H. Adelson, and D.J. Heeger, “Shiftable multi-scale transforms,”. IEEE Trans. on Info. Theory , Special Issue on Wavelets **38**, 587607 (1992).

**22. **P. Perona, “Deformable kernels for early vision,” IEEE Trans. Pattern Analysis and Machine Intelligence , **17**(5), 448–499 (1995). [CrossRef]

**23. **J. W. Woods. *Subband Image Coding*. (Kluwer Academic Publishers, MA,1991).

**24. **W.K. Pratt, *Digital Image Processing*. (John Wiley & Sons, NY, 1991).

**25. **K.R. Castleman, *Digital Image Processing*. (Prentice Hall, NJ,1996).

**26. **P.D. Welch, “The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms,” IEEE Trans. Audio Electroacoust. **15**, 70–73 (1967). [CrossRef]

**27. **K.M. Hanson, “Detectability in computed tomographic images,” Med. Phys. **6**(5), 441–451 (1979). [CrossRef] [PubMed]