Realistic anatomical images are useful for assessment and improvement of medical image quality. The use of synthesized images has the advantage of providing the user with a large number of independent samples, in a controlled environment. We propose a method to generate medical textures that are fully defined by a set of adjustable parameters and a random number generator, and which statistical properties are analytically tractable. This method, called the “clustered lumpy background”, is a generalization of the original lumpy background described by Rolland and Barrett (1992). A detailed application of the method in the case of mammography is presented. It is shown that the synthesized images are visually very similar and that their first and second order statistics can be considered as being equivalent.
© Optical Society of America
A main goal of research in medical imaging is to improve diagnostic accuracy. Whether the diagnostic decision is made by a physician visually inspecting the image, or a computer algorithm performing an image analysis, the accuracy is mainly limited by the image clutter consisting of the stochastic nature of photons (quantum noise) and the variability inherent in biological tissue. Performance for human visual detection/classification and image analysis algorithms such as texture analysis of mammography [1,2], osteoporosis , chest imaging [4,5], and ultrasound  are specific to the statistical properties of the images. It is for this reason that most observer studies [7–9] (ROC studies) and computed aided diagnosis  (CAD) algorithms development and testing are done using samples of real medical images pertinent to the imaging modality and anatomy.
However, one limitation of using real medical image backgrounds is there may be access to only a limited number of samples  and that the statistical properties of the images can be difficult to fully characterize. Within this context, developing a method to synthesize realistic looking textures might facilitate some of the evaluation/optimization of medical images for visual inspection and computer analysis. For example, for evaluation of CAD, the use of realistic looking synthetic textures might be an attractive option for a fast preliminary evaluation and optimization of algorithms prior to the testing in real medical images. In addition, many investigators are currently developing computer model observers that can predict human performance visually detecting simulated lesions in medical images [12–14]. These models could potentially be used for efficient evaluation/optimization of medical image quality. Tractable mathematical expressions to describe the statistics of the images highly simplify the calculation for these model observers . For this reason, many investigators have used computer-generated noise including white noise and filtered white noise [16,17]. However, these simulated backgrounds are visually very different from most images of medical or biological tissues containing anisotropic structures that are part of the anatomy and have to be taken into account in the diagnostic process. Others have applied the model observers to simulated lesions embedded in real x-ray angiographic images . However, for these cases calculation of model observer performance requires more complex sample driven methods and access to a sufficient number of sample images . The development of a method to create realistic looking textures with simple expressions characterizing the image statistics might facilitate development of these models for real practical use.
One attempt to estimate human and model performances with more realistic looking backgrounds than filtered white noise has been proposed by Rolland and Barrett . The images obtained, called “lumpy backgrounds”, consist of a random number of single structures called “blobs” located at random locations. The lumpy backgrounds have the advantage of being statistically tractable and stationary. Their main application for the method has been the simulation of variability in nuclear medicine images as part of an investigation of the resolution-noise trade-off for different aperture sizes in a pinhole imaging system. This approach has also proven useful in the investigation of image reconstruction algorithm .
In order to generate images that mimic complex biological images, Rolland and Strickland  proposed a method of texture synthesis based on the summation of a low frequency lumpy background and a pyramid framework algorithm with histogram equalization . Histograms of simulated and real images are identical, which is to be expected because of the histogram-matching step in the algorithm. However, the power spectra differ in a way that the real images have a monotonical power spectrum whereas the generated image have a power spectrum that contains ripples. An additional difficulty in the method of Rolland and Strickland is that no simple mathematical expressions are presented to characterize the statistical properties of the texture synthesis images and to relate them to the parameters of the simulation algorithm.
The present paper focuses on the development of a technique to generate textures that visually appear as real medical images but have simple tractable statistical properties. The method we describe is a generalization of the lumpy background proposed by Rolland and Barrett . We call this method the “clustered lumpy background” (CLB) because it clusterizes blobs within the image. It permits the generation of textures that visually resemble real images and have similar histograms and power spectra. The paper is organized in the following way. The theoretical section defines first and second order statistics, computation of the power spectra, and describes the CLB method. The results section presents an application of the method to mammography with image samples and statistical comparisons.
2.1 Statistics of an image
An image can be considered as a random variable described by its probability density function  (pdf). However, the complete pdf of a medical image is difficult if not impossible to compute. Image statistics are usually summerized by some of their first and second order statistics. The first order statistics of an image are defined relative to one point. It is the probability density of the brightness of that single point computed across all image realizations. Usually, it is assumed that the images are ergodic and stationary. Ergodicity means that the pdf computed across all images for a given point and the pdf computed across an infinite image at every point are equivalent . Stationarity means that the pdf does not depend on the image location. With these two assumptions, the first order statistics of an image can be derived directly from the probability density of the brightness of one image.
The density is usually presented as a histogram and can be quantified by some of its moments (or quantities derived from its moments) like mean, variance, skewness or kurtosis. By definition, these statistics only describe the distribution of the brightness in the image, without any reference to the spatial relation between different locations in the image.
The second order statistics of an image are defined in terms of the relationship between two points. They are fully determined by the joint probability density of the brightness at these two points. One of the most commonly used second order statistics is the covariance function or, if the image is considered stationary, the autocorrelation function. In this latter case, the Wiener-Kintchine theorem can be applied and the Fourier transform of the autocorrelation is the power spectrum.
Formally, stationarity requires an image of infinite extent. Although stationarity is never the case, this property can be used within the boundary of the image to estimate some statistical properties. The power spectrum of an image sample is estimated by the squared modulus of its Fourier transform. Because the power spectrum computed from a single sample image is noisy, the power spectrum is estimated as being the mean value of the sample power spectra computed from a larger group of images.
In the case of a digital image, the power spectrum is computed by using a discrete Fourier transform (DFT) that gives an estimate of the Fourier transform of the infinitely replicated image. For non-periodic images like real medical images, the replication effect acts as a sharp edge on each side of the image. This induces an artifact in the frequency domain in directions perpendicular to the edge. This can be avoided by windowing the image prior to computing its DFT. According to Papoulis , a window that results in a low bias spectrum is the following:
where r is the spatial coordinate, rw defines the size of the window, and the origin (r=0) is at the center of the image. The relationship between the power spectrum obtained when measured with a windowed image (Wwindow) and the actual “true” power spectrum (Wtrue) is given by :
where “~” is used to indicate the Fourier transform of the associated variable, and f is the spatial frequency.
2.2 Clustered lumpy background (CLB)
The CLB technique is a generalization of the lumpy background technique described by Rolland and Barrett . It permits generation of random backgrounds containing a texture that can be adjusted to have properties similar to natural structures like those observed in medical images. The CLB technique generates stationary images, with statistical properties that can be derived analytically.
The Rolland-Barrett lumpy background is generated by placing a random number (K) of delta functions at random points in the image. The number of points is selected using a Poisson probability distribution, whereas the locations (r k) are selected using a uniform probability distribution over the area of the image. Each delta function serves as the origin for a blob defined using a symmetrical function b(r). Formally, an image obtained by this method can be written as:
where r is the 2D spatial coordinate, and r k is the origin of the kth blob. In the original description of this method  this kind of background has been named type I lumpy background.
The CLB is generated by first producing a random uniform distribution of delta functions as in the type I lumpy background method. This serves as a “cluster map”. Each of these positions is used as an origin for a second process whereby a random number (Nk) of delta functions are randomly positioned. The number Nk in the kth cluster is a Poisson random variable and the positions r kn are selected using a pdf ϕ(r). The K delta functions are then removed from the image and replaced by K clusters of delta functions. Finally each delta function serves as the origin for a blob defined using a function b(r/a,R θ), where a is a random scaling factor, R θ is the rotation matrix that rotates the whole blob, and the angle θ is a random variable. A different angle and scale can be used for each of the M = Nk delta function locations. The blobs are then summed to create the image. Formally, the value of an image g at position r is then:
where the different variables are defined in Table 1. The blob asymmetry allows the local anisotropy observed in real tissues to be taken into account. The limited spread of the blobs within a cluster renders the biological middle-scale structure.
In the practical implementation, we did not consider any scale variation (akn=1). The rotation angle was selected using a uniform pdf (between 0 and 2π) and was considered constant within a given cluster (θkn=θk). Finally, the pdf of the blobs locations within a cluster was chosen Gaussian with standard deviation σϕ.
The power spectrum of such an image is given by (see Appendix for derivation):
where A is the area of a single image, Wb is the power spectrum component coming from the blob only and Ws is the power spectrum component coming from the interaction of the cluster shape and the blob. These two spectra are defined from the CLB parameters as follows:
For most reasonable and practical choices for ϕ(r), ϕ̃(f) will peak at f=0 with ϕ̃(0) =1. If this is the case, Ws(f) is then always smaller than Wb(f) and the influence of Ws is more important at low frequency. Eq. (5) shows that the shape of the image power spectrum is not affected by the mean number of clusters (K¯) but is influenced by the mean number of blobs per cluster (N¯). For a low value of N¯ , the shape of an individual blob dominates the image power spectrum whereas an increase of the value of N¯ modifies the shape of the image power spectrum towards Ws.
3.1 Mammographic power spectra
We measured the power spectra of a representative set of 32 mammograms free of disease according to the radiologist visual diagnosis. The images were obtained by a digital 510(K)-Bennett mammographic unit (Trex Medical Corporation, Waltham, MA-USA) (courtesy of University of California Los Angeles Medical Center). They have a pixel size of 0.04mm and are coded as 14 bits per pixel and the gray level is proportional to the x-ray exposure. The power spectra were measured on 1024×1024 pixel sample images (41×41mm2) with the window defined by Eq. (1). The mean image value was subtracted prior performing the DFT. These sample images were all located entirely within the breast excluding the chest wall. Fig. 1 shows the extreme values obtained for the measurement of the power spectra averaged over all orientations and displayed versus their radial spatial frequency. It is observed that the measured power spectra have a power law shape at low frequency with a plateau at high frequency. The slope of the power law is between -3.4 and -4.0. This is significantly higher in magnitude than our previous results , which were biased downward by the fact that no window was used.
3.2 Application of the CLB technique
As shown in the previous section, the simulation of a mammogram implies a power spectrum that has a power-law shape at low frequency and is white at high frequency. The constant part of the power spectrum can be simulated by a given quantity of white noise and the power-law part can be achieved with an asymmetrical exponential blob shown in Fig. 2(a) and defined in the following way:
where L(r) is the characteristic length of the exponential, and α and β are adjustable coefficients which effects are described below. The blob is first generated with characteristic length Lx and Ly in the x and y direction respectively. It is then rotated by the randomly selected angle θ having a uniform pdf between 0 and 2π. The blob function is not separable and the characteristic length L is equal to the “radius” of the ellipse having half-axes Lx and Ly (see Fig. 2(b)). For this particular application, we did not consider any scaling variation of the blobs. Therefore, the blobs contained in the image only differ by their orientations.
A value of θ lower than unity is necessary to give Wb a power law shape. The slope of the spectrum is only slightly affected by the actual value of β for β near 0.5. The ratio of Lx/Ly has only a small effect on the slope of Wb and can therefore be chosen in order to render the more or less fibrous appearance of the image. The value of α has a direct influence on the slope of the power spectrum: the larger the coefficient, the smaller the slope of the power spectrum. The parameters were chosen by trial and error using power spectrum and histogram analysis. Finally a visual inspection of some sample images was performed.
Fig. 3 shows some examples of simulated and typical real mammographic images displayed at a pixel size of 0.3mm. The chosen parameters are reported in Table 2 for a 128×128-pixel image (pixel size: 0.3mm). For a larger image, the number of clusters can be scaled proportionally to the image area. The simulated images shown were not selected from a larger group, but are simply the first four images created using the random number generator.
3.3 Statistical comparison
As explained in section 2.1, the first order statistics do not fully describe the pdf of an image but only the histogram of the brightness, without any reference to the spatial relation between different locations. The first order statistics parameters that characterize this pdf are presented in Table 3 for the real mammograms and CLBs. The first order statistics of these images were estimated from 286 non-overlapping areas taken from 32 real mammograms and 500 CLBs (image size: 128×128 pixels; pixel size: 0.3mm). For each individual image, the mean, the standard deviation, the skewness, and the kurtosis (with a definition such that a Gaussian distribution has a kurtosis equals to zero) were computed. Table 3 displays the means and standard deviations of these latter parameters. It can be seen that the values of both types of images are very similar.
The second order statistics are only analyzed in terms of the image power spectrum. Fig. 1 shows the CLB power spectrum compared to the extreme values measured on real mammograms.
The visual comparison of real mammograms and CLBs shows that the global characteristics of the mammograms are well reproduced. As the real mammograms, the CLBs have patterns appearing at different scales with more or less fibrous structure in-between. Nevertheless, with a sterner look, it is apparent that real mammograms have more variety and may contain in particular more fibrous structures than the simulated images.
The first order statistical analysis only shows minor differences between the two types of images. The mean values of the real mammograms are slightly higher than the CLBs’. The mean standard deviations of both types of images are very close and have the same amount of fluctuation. The two kinds of images are on average symmetric. The main difference lies in the fact that the real images skewness has a wider range of fluctuations. The negative kurtosis values indicate that both mammogram and CLB image brightness distributions are slightly less peaked than a Gaussian pdf. As with the skewness, the kurtosis has more variation in the real images than for the CLBs. Globally, this analysis shows that the CLB images have overall similar properties as the real mammograms. However, the real mammograms have more intrinsic variation from image to image. The CLB power spectrum is very similar to the real ones. It obeys a power law in the low frequency range and flattens at high frequency.
The main advantages of the CLB over previous methods proposed to synthesize textures are its simplicity and its statistical robustness. The images are stationary (within the boundary of the image) by definition and their statistical properties can be computed analytically. For example, this mathematical simplicity can greatly simplify the calculation of ensemble performance of model observers in signal detection tasks using closed form expressions [28,29] rather than having to use sample driven methods .
The main disadvantage of the CLB is the lack of a straightforward way to choose parameters for a given texture. Because the statistics as well as the visual appearance of the generated images are defined by all the CLB parameters, optimizing the simulation of a given texture may require some effort. However, the fact the statistics can be analytically computed speeds the process considerably since the formulas can be evaluated quickly and then checked against a sample of textures. In light of these issues, it should be noted that the parameters presented here for mammography might have another better solution.
A new method (CLB method) which builds on the Rolland-Barrett type I lumpy background is presented. It permits the generation of stationary random textures that are a reasonably good facsimile of realistic medical image backgrounds and have statistical properties that can be expressed analytically. A detailed application of the method has been presented for mammography. It is shown that the generated images have a close visual appearance to real mammographic backgrounds and that the first and second order statistics of both groups of images are very similar.
This method is of course not limited to mammography and could be generalized to textures in other applications by changing the parameters of the CLB algorithm, in particular the blob function or the cluster distribution.
5. Appendix: power spectrum of the CLB
In this appendix, we give a detailed power spectrum calculation for the CLB. A single realization of the CLB is presumed to have K clusters each of them containing Nk blobs. Each is Poisson distributed with mean K¯ and N¯ respectively. We restrict our calculation to the case where the blob orientation is constant within a cluster (θkn=θk), and where all the blobs have the same scaling factor (akn=1). The expectation value of the CLB defined by Eq. (4) is the following:
where A is the area of the image (considered sufficiently large so that we can neglect effects at the borders), is the integral of a single blob averaged over all angles, Θ is the domain of variation of the blobs centers inside a cluster, and ϕ(r) is the probability density of the blob centers inside a cluster. The expectation operation is expressed by brackets (“<⋯>”) with a subscript to indicate the variable over which the expectation is performed.
The covariance of the image is defined by the following relation:
For practical convenience, we will name the three terms of the right-hand side of Eq. (A.3) T1+T2+T3. The first term is:
where Rbθk (Rbθ (r - r’)= ∫A dr k b(r - r k ,R θ )b(r’-r k ,R θ)) is the autocorrelation function of bθk and is the autocorrelation of bθk averaged over all orientations θk. The second term of Eq. (A.3) can be written as:
where Sθ(r) (Sθ(r)= ∫Θdr’ϕ(r’)b(r - r’,R θ)) is the convolution of the cluster pdf and the blob function, is the autocorrelation function of Sθ(r) averaged across all orientations. The final step of Eq. (A.5) takes into account that the number of blobs in a cluster (N) is Poisson distributed. Finally, in the third term of Eq. (A.3) the fact that k is different than k’ implies that r kn and r k’n’ are independent as well as r k and r k’, and θk’ and θk’. This permits us to write T3 as:
where the last step has again been made possible by taking into account the Poisson distribution of the number of clusters (K). If we pool Eq. (A.4) to (A.6) into (A.3), and we take Eq. (A.1) into account, we can write the image covariance (Eq. (A.2)) as:
The Wiener-Kintchine theorem  says that the Fourier transform of the autocorrelation function is equal to the power spectrum of that function, which allows us to write the image power spectrum as stated in Eq. (5).
We are very grateful to Dr Carolyn Kimme-Smith from University of California Los Angeles who made the digital mammograms available for this study and Elias Pony Merhige for stimulating discussions. We would also like to thank an anonymous reviewer for the very constructive comments and suggestions, as well as suggesting the name “clustered lumpy background”. This work has been funded by the Swiss National Fund for Scientific Research and the National Institute of Health (NIH) ROI HL 53455.
References and links
1. A. Petrosian, H. P. Chan, M. A. Helvie, M. M. Goodsitt, and D. D. Adler, “Computer-aided diagnosis in mammography: classification of mass and normal tissue by texture analysis,” Phys. Med. Biol. 39, 2273–2288 (1994). [CrossRef] [PubMed]
2. H. P. Chan, W. Wei, M. A. Helvie, B. Sahiner, D.D. Adler, M. M. Goodsitt, and N. Petrick, “Computer-aided classification of mammographic masses and normal tissue: linear discriminant analysis in texture space,” Phys. Med. Biol. 40, 857–876 (1995). [CrossRef] [PubMed]
4. G. Revesz, H. L. Kundel, and M. A. Graber, “The influence of structured noise on the detection of radiologic abnormalities,” Am. J. Roentgenol. 9, 479–486 (1974).
5. J. Moshita, K. Doi, S. Katsuragawa, L. Monnier-Cholley, and H. MacMahon, “Computer aided diagnostic for interstitial infiltrates in chest radiographs: Optical-density dependence of texture measures,” Med. Phys. 22, 1515–1523 (1995). [CrossRef]
6. J. W. Allison, L.L. Barr, R. J. Massoth, G. P. Berg, B. H. Krasner, and B. S. Garra, “Understanding the process of quantitative ultrasonic tissue characterization,” RadioGraphics 14, 1099–1108 (1994). [PubMed]
7. T. Kobayashi, X. W. Xu, H. MacMahon, C. E. Metz, and K. Doi, “Effect of a computer-aided diagnosis scheme on radiologists’ performance in detection of lung nodules on radiographs,” Radiology 199, 843–848 (1996). [PubMed]
8. E. Samei, M. J. Flynn, and W. R. Eyler, “Simulation of subtle lung Nodules in projection chest radiography,” Radiology 202, 117–124 (1997). [PubMed]
9. A. J. Mendez, P. G. Tahoces, M. J. Lado, M. Souto, and J. J. Vidal, “Computer-aided diagnosis: automatic detection of malignant masses in digitized mammograms,” Med. Phys. 25, 957–964, (1998). [CrossRef] [PubMed]
10. W. Zhang, K. Doi, M. L. Giger, Y. Wu, R. M. Nishikawa, and R. A. Schmidt, “Computerized detection of clustered microcalcifications in digital mammograms using a shift-invariant artificial neural network,” Med. Phys. 21, 517–524 (1994). [CrossRef] [PubMed]
11. C. Kimme-Smith, M. McCombs, R. H. Gold, and L. W. Bassett, “Mammography fixed grid versus reciprocating grid: Evaluation using cadaveric breasts as test objects,” Med. Phys. 23, 141–147 (1996). [CrossRef] [PubMed]
12. R. F. Wagner and K. E. Weaver, “An assortment of image quality indices for radiographic film-screen combinations - can they be resolved?,” proceedings SPIE 35, 83–94 (1972). [CrossRef]
13. A. E. Burgess, “Statistically defined backgrounds: Performance of a modified nonprewhitening observer model,” J. Opt. Soc. Am. A 11, 1237–1242 (1994). [CrossRef]
14. M. P. Eckstein, C. K. Abbey, and J. S. Whiting, “Human versus model observers in anatomic backgrounds,” proceedings SPIE 3340, 16–26 (1998). [CrossRef]
17. K. J. Myers, H. H. Barrett, M. C. Borgstrom, D. D. Patton, and G. W. Seeley, “Effect of noise correlation on detectability of disk signals in medical imaging,” J. Opt. Soc. Am. A 2, 1752–1759 (1985). [CrossRef] [PubMed]
18. M. P. Eckstein and J. S. Whiting, “Visual signal detection in structured backgrounds. I. Effect of number of possible spatial locations and signal contrast,” J. Opt. Soc. Am. A 13, 1777–1787 (1996). [CrossRef]
22. E. P. Simoncelli, W. T. Freeman, E. H. Adelson, and D. J. Heeger, “Shiftable multi-scale transforms,” Trans. on Info. Theory, Special Issue on Wavelets 38, 587–607 (1992).
23. B. Picinbono, Random Signals and systems (Prentice Hall International, 1993), p.182.
24. A. Papoulis, Probability, random variables, and stochastic processes (McGraw-Hill, Inc, 1991), p.453.
25. A. Papoulis, Probability, random variables, and stochastic processes (McGraw-Hill, Inc, 1991), p.419.
26. J. P. Rolland, Factors influencing lesion detection in medical imaging (Ph.D. dissertation, University of Arizona, 1990).
27. F. O. Bochud, F. R. Verdun, C. Hessler, and J. F. Valley, “Detectability on radiological images: The effect of the anatomical noise,” proceedings SPIE 2436, 156–164 (1995). [CrossRef]
28. H. H. Barrett, J. L. Denny, R. F. Wagner, and K. J. Myers, “Objective assessment of image quality. II. Fisher information, Fourier crosstalk, and figures of merit for task performance,” J. Opt. Soc. Am. A 12, 834–852 (1995). [CrossRef]
29. A. E. Burgess, X. Li, and C. K. Abbey, “Visual signal detectability with two noise components: anomalous masking effects,” J. Opt. Soc. Am. A 14, 2420–2442 (1997). [CrossRef]
30. J. C. Dainty and R. Shaw, Image Science (Academic, London, 1974).