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.¹²

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.¹³ 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.¹⁴ 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.¹⁵ 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.

 

Fig. 1. Mammography breast image decomposition: (a) The original sample. (b) The slowly, spatially varying mean-background. (c) The residual texture image.

Download Full Size | PDF

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.¹ 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.¹¹ 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).⁴ The first approach consisted of superimposing 2D Gaussian functions, referred to as Gaussian blobs, on a constant background of strength B₀. To keep the mathematics simple, we assumed 2D Gaussian blobs of constant amplitude b₀/π${{\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.⁴ A measure of lumpiness in the background is given by

where K̅/A_d is the mean number of Gaussian functions per mm².

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.¹⁷ 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.¹⁸ 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.¹⁷

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.²³ 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.¹⁸

 

Fig. 2. Illustration of the steerable pyramid transform used in the texture synthesis algorithm. The input image in the upper left corner would be either the texture sample or the white noise image. The output image in the upper right corner will be either the reconstruction of a decomposed image if only one input image is considered, or a synthesis image if two pyramid layers are combined as described in section 4. The left hand side of the pyramid is used for decomposing the two images and the right hand side of the pyramid is used for image reconstruction or synthesis.

Download Full Size | PDF

 

Fig. 3. Syntheses of a residual mammographic texture image: (a) a typical sample of a uniformly distributed white noise image used as a starting point for one synthesis; (b) original mammographic residual texture; (c) synthesis 1; (d) synthesis 2.

Download Full Size | PDF

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.¹⁹

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.¹³ 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.¹⁶ 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.¹⁵ 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.²⁶ An extension of Welch’s method to two dimensional noise images was detailed in Hanson.²⁷ 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:

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.²⁶

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⁶ counts²/(sec² 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.

 

Fig. 4. Average greylevel histograms over 18 images of an ensemble for five ensemble sets: (a) the original sample mammograms; (b) the mean backgrounds; (c) the lumpy backgrounds that best matched the mean backgrounds in visual appearance; (d) the residual texture images; and (e) the texture synthesis images.

Download Full Size | PDF

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.

 

Fig. 5. Average Power spectra over 18 images of an ensemble for five ensemble sets: (a) the original sample mammograms; (b) the mean backgrounds; (c) the lumpy backgrounds that best matched the mean backgrounds in visual appearance; (d) the residue texture images; and (e) the texture synthesis images.

Download Full Size | PDF

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.