Pattern Recognition Letters 19 Ž1998. 735–739
Self-similar texture characterization using a Fourier-domain maximum likelihood estimation method C.-Y. Wen a
a,)
, R. Acharya
b
Department of Computer Science and Information Management, Hung-Kuang Institute of Technology, Sha-Lu, 34 Chung-Chie Rd., Taichung, Taiwan b Department of Electrical and Computer Engineering, State UniÕersity of New York at Buffalo, Buffalo, NY, USA Received 8 September 1997; revised 13 March 1998
Abstract A Maximum Likelihood Estimator ŽMLE. has been applied to estimating the Hurst parameter H on a self-similar texture image. Much of the work done so far has concentrated on the spatial domain. In this paper, we propose an approximate MLE method for estimating H in the Fourier domain. This method saves computational time and can be applied to estimating the parameter H directly from the Fourier-domain raw data collected by the Magnetic Resonance Imaging ŽMRI. scanner. We use synthetic fractal datasets and a human tibia image to study the performance of our method. q 1998 Elsevier Science B.V. All rights reserved. Keywords: Maximum likelihood estimator; Fourier domain; Fractals; Textures
1. Introduction Fractals have been successfully used in the characterization of many natural textures. Fractional Brownian Motion ŽFBM. has been used to model self-similar textures ŽPeitgen and Saupe, 1988.. While using the fractal model, the most important procedure is to measure the Hurst parameter H, which is directly related to the fractal dimension. A Maximum Likelihood Estimator ŽMLE. has been used to estimate the Hurst parameter H on a selfsimilar texture image ŽLundahl et al., 1986; Benhamou et al., 1994., which is in the spatial domain. )
Corresponding author. E-mail:
[email protected].
However, for some applications, such as Magnetic Resonance Imaging ŽMRI., its collected data are the Fourier transform of the tomographic image ŽAcharya et al., 1995.. We hope to directly work with the collected data before obtaining the tomographic image. In this work, we propose an approximate Fourier-domain MLE method and apply it to estimating the parameter H. In this Fourier-domain MLE method, we use the discrete-time Fourier Transform to whiten discrete Fractional Gaussian Noise ŽDFGN. Žpractically, we use the discrete Fourier transform ŽDFT. to approximate., so that we can avoid the extensive computation to obtain the inverse matrix and determinant of the covariance matrix compared to the spatial-domain MLE method. Similar work was done by Wornell and Oppenheim Ž1992.. They
0167-8655r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 8 6 5 5 Ž 9 8 . 0 0 0 5 1 - 8
C.-Y. Wen, R. Acharyar Pattern Recognition Letters 19 (1998) 735–739
736
exploited the Wavelet transform as a whitening filter for 1rf processes and derived an approximate MLE, using wavelet coefficients. We use synthetic fractal datasets and compare both spatial-domain and Fourier-domain MLE methods. The results show that we can obtain similar results with both methods. We also applied both methods on images obtained by scanning a human tibia in the MRI scanner.
Fractional Brownian Motion ŽFBM. is given as ŽMandelbrot and Ness, 1968; Voss, 1986. BH Ž 0 . s b 0 ,
H . ½
0
G Ž H q 1r2
y`
Ž t ys.
Hy 12
H0 Ž t y s .
q
t
Hy 1r2
Hy 1r2
d BŽ s.
5
d BŽ s. ,
Ž 1.
where 0 - H - 1, B Ž t . is an ordinary Brownian motion and b 0 is an arbitrary real number. For H s 0.5, FBM is an ordinary Brownian motion. From Eq. Ž1., we can see the difference signal B H Ž t 2 . y BH Ž t 1 . s
1
G Ž Hq
. ½H
t2
1 2
y`
Ž t2 y s .
Hy 1r2
Hy` Ž t y s .
y
t1
1
DFGN has zero mean and its autocorrelation can be written as rw kx s
s2
d BŽ s.
Hy 1r2
< k q 1 < 2 H y 2 P < k < 2 H q < k y 1 < 2 H , Ž 4.
2
3. Maximum Likelihood Estimator (MLE)
`
SŽ v . s
y Ž ys .
Ž 3.
DFGN is a Gaussian stationary process with zero mean and autocorrelation r w k x s r wyk x in Eq. Ž4.. Its Power Spectral Density ŽPSD. was obtained by Papoulis Ž1991. and Porat Ž1994.,
BH Ž t . y BH Ž 0 . 1
x w n x s B w n x y B w n y 1x .
where s 2 is the variance of x w n x.
2. Fractional Brownian Motion (FBM)
s
Gaussian Noise ŽFGN.. Discrete FBM can be defined as Bw n x and discrete FGN ŽDFGN. as
d BŽ s.
Ž 2.
is a strict-sense stationary process and has a Gaussian distribution with zero mean and variance VH < t 2 y t 1 < 2 H , where VH is a constant. FBM is a nonstationary process and is difficult to analyze. However, its increments are a strict-sense stationary process which has been termed Fractional
Ž 5.
ksy`
Because r w k x s r wyk x, SŽ v . is real, symmetric in v , and nonnegative. The discrete-time Fourier transform of DFGN, X Ž v ., is a Gaussian nonstationary white noise with zero mean and autocovariance E X Ž u. X Ž Õ .
)
s 2 p S Ž u. d Ž u y Õ .
Ž 6.
where yp - u,Õ - p. Practically, Eq. Ž5. and ŽEq. Ž6.. can be just approximated because of the finite data length. Let us define x N w n x as the windowed time series of DFGN Ž x w n x. xN w nx s
5
Ý r w k x eyj k v .
½wx
x n , 0,
n s 0,1, . . . , N y 1, other,
Ž 7.
where N is the total data length Žduration of the window.. We have x w n x s lim x N w n x . N™`
Ž 8.
The discrete-time Fourier transform of the windowed time series of x N w n x is given by XN Ž u . s
Ny1
Ý
ns0
x N w n x eyj u n ,
Ž 9.
C.-Y. Wen, R. Acharyar Pattern Recognition Letters 19 (1998) 735–739
where yp - u - p. The mean of X N Ž u. can be written as E XN Ž u . s
Ny1
Ý
E x N w n x eyj u n .
Ž 10 .
If u s Õ, let l s m y n and Eq. Ž11. can be written as
s
Ý Ý
S Ž u . s lim
N™`
rN w m y n x eyjŽ u myÕ n. ,
Ž 11 .
where yp - u,Õ - p, and
½
rw kx, 0,
1 N
where a is integer.
,
/
rN w l x eyju l .
Ž 14 .
Ž 15 .
If u / Õ, then Ew X N Ž u. X N Ž Õ . ) x / 0. However, the values for u / Õ are much less than the values for u s Õ, as N is large. Let us define
Ž 16 .
and
)
)
N
E < XN Ž u . < 2 .
Ž 12 .
s E X N Ž u q 2 ap . X N Ž Õ q 2 ap .
R N w i , j x s rN w < i y j < x
k s 0," 1, . . . ," Ž N y 1 . , other,
where r w k x is given in Eq. Ž4.. Because rN w k x s rN wyk x, Ew X N Ž u. X N Ž Õ . ) x is real and E XN Ž u . XN Ž Õ .
ž
1y
We can see that, when N ™ `, Eq. Ž5. can be obtained by
)
ms0 ns0
rN w k x s
Ý
lsyNq1
We can see that, when N ™ `, X N Ž u. also has zero mean. The autocovariance of X N Ž u. can be written as
Ny1 Ny1
Ny1
E < XN Ž u . < 2 s N P
ns0
E XN Ž u . XN Ž Õ .
737
Ž 13 .
L N w p,q x s E X N
ž / ž / 2p p N
XN
2pq N
)
,
Ž 17 .
where i, j, p,q s 0,1, . . . , N y 1. Without losing generality, we let s 2 s 1 and take H s 0.4 for example to compare Eqs. Ž16. and Ž17. with plots Žsee Fig. 1.. For R N w i, j x ŽFig. 1Ža.., the difference between diagonal and non-diagonal values is small Žy0.2 ; 1.0., while the one for L N w p,q x ŽFig. 1Žb.. is large Ž0.0 ;
Fig. 1. Covariance matrix, H s 0.4, s 2 s 1, N s 100. Ža. R N w i, j x; Žb. L N w p,q x.
C.-Y. Wen, R. Acharyar Pattern Recognition Letters 19 (1998) 735–739
738
Table 1 Results for applying both spatial-domain and Fourier-domain MLE methods to the generated datasets True H
Spatial
Fourier
Mean H
STD
Mean H
STD
0.2 0.3 0.4 0.5 0.6 0.7 0.8
0.193 0.291 0.396 0.492 0.591 0.693 0.788
0.05 0.06 0.07 0.07 0.08 0.08 0.08
0.192 0.291 0.390 0.490 0.587 0.685 0.778
0.05 0.06 0.08 0.08 0.09 0.10 0.10
140.0.. The difference for L N w p,q x between diagonal and non-diagonal values begins larger when N increases. We can approximate L N w p,q x as a diagonal matrix when N is large enough. Let us go back to derive the Fourier-domain MLE estimator. Because of linearity, we can assume that the DFT
of DFGN still has a Gaussian distribution when N is large, and its PDF can be written as pŽ X ; H . 1 1 ) s exp y X NT Ly1 N XN , Nr2 1r2 2 Ž 2 p . < LN <
½
5
Ž 18 .
where L N is the autocovariance matrix given in Eq. Ž17.. Let us define 1 LXN s 2 L N Ž 19 . s and after some manipulation, we obtain log p Ž X ; H . s y
2 N
log 2 p y 1
N 2
log
ž
) X NT LXy1 N XN
N
/
log < LXN < . Ž 20 . 2 2 has to be inverted and its determinant has to be y
LXN
N
y
Fig. 2. An MRI image of the human tibia. We choose 34 scan lines Žblack. to estimate their H values with both spatial-domain and Fourier-domain MLE methods.
C.-Y. Wen, R. Acharyar Pattern Recognition Letters 19 (1998) 735–739
computed. When N is large, we can approximate LXN as a non-zero diagonal matrix. We can obtain its inverse and determinant directly. While a spatial-domain MLE needs extensive computation to obtain the inverse and determinant of the covariance matrix ŽOŽ N 2 .., a Fourier-domain MLE saves more computational time ŽOŽ N ...
4. Experimental results 4.1. Mathematical simulation When we compute the PSD of DFGN ŽEq. Ž5.., at least 100 terms must be included in the sum ŽLundahl et al., 1986.. So, we chose the data length to be 100. One hundred datasets Ždata length N s 100. for true H s 0.2, 0.3, . . . ,0.8 were generated by the Cholesky decomposition method ŽLundahl et al., 1986.. The Fourier data was generated by transformation of the Cholesky decomposition results from the generation algorithm. We used those datasets to compare both spatial-domain and Fourier-domain MLE methods. The experimental results are presented in Table 1. We can see that, when H ( 0.5, both methods almost obtain similar results; however, when H ) 0.5, the effect of approximation error increases. 4.2. Imaging tibia in the MRI scanner An isolated tibia with varying degree of trabecular structure was imaged in the MRI scanner. The trabecular structure exhibits itself as a rich texture in the MRI images. We applied both MLE methods to estimating the H value in textural regions of the human tibia Žsee Fig. 2.. For the spatial-domain test, we used directly the MRI images Ž20 spatial slices.. For the Fourier-domain test, we used transformation of the individual slices. The estimated H mean values we obtained were very similar: 0.247825 with the spatial-domain MLE and 0.243243 with the Fourier-domain MLE.
739
5. Conclusion FBM has been used to describe many medical images, such as human calcaneus radiographs. A MLE has been applied to estimate a fractal parameter H on those images. However, for some applications, such as MRI, its collected data are the Fourier transform of the tomographic image. In this paper, we propose a Fourier-domain MLE method which can be applied to those Fourier-domain raw data instead of images. From computer simulations, we find that the Fourier-domain MLE can obtain performance similar to the spatial-domain MLE when H ( 0.5. The Fourier-domain MLE can save computational time without extensive computation in order to obtain the inverse and determinant of a matrix. This is not so in a spatial-domain MLE. We have applied our algorithms on a real human tibia image scanned in the MRI scanner. References Acharya, R., Wasserman, R., Stevens, J., Hinojosa, C., 1995. Biomedical imaging modalities: A tutorial. Computerized Medical Imaging and Graphics 19 Ž1., 3–25. Benhamou, C.L., Lespessailles, E., Jacquet, G., Harba, R., Jennane, R., Loussot, T., Tourliere, D., Ohley, W., 1994. Fractal organization of trabecular bone images on calcaneus radiographs. J. Bone and Mineral Res. 9 Ž12., 1909–1918. Lundahl, T., Ohley, W.J., Kay, S.M., Siffert, R., 1986. Fractional Brownian motion: a maximum likelihood estimator and its application to image texture. IEEE Trans. Medical Imaging 5 Ž3., 152–161. Mandelbrot, B.B., Ness, H.W.V., 1968. Fractional Brownian motion, fractional noises and applications. SIAM Review 10, 422–436. Papoulis, A., 1991. Probability, Random Variables, and Stochastic Process. McGraw-Hill, New York. Peitgen, H., Saupe, D., Eds., 1988. The Science of Fractal Images. Springer, New York. Porat, B., 1994. Digital Processing of Random Signals: Theory and Methods. Prentice-Hall, Englewood Cliffs, NJ. Voss, R.F., 1986. Characterization and measurement of random fractals. Physica Scripta 13, 27–32. Wornell, G.W., Oppenheim, A.V., 1992. Estimation of fractal signals from noisy measurement using wavelets. IEEE Trans. Signal Processing 40 Ž3., 611–623.