- Email: [email protected]
Tracking of nonstationarities for texture fields
Signal Processing 14 (1988) 95-102 North-Holland
95
SHORT COMMUNICATION
TRACKING OF NONSTATIONARITIES FOR TEXTURE FIELDS* Todd REED and Harry WECHSLER Department of Electrical Engineering, University of Minnesota, Minneapolis, M N 55455, U.S.A. Received 19 November 1986 Revised 11 February 1987 and 1 July 1987
Abstract. We consider the problem of texture segmentation as one of tracking nonstationarities for 2-D signals. The underlying signal representation is the Wigner distribution. It is a co-joint representation of space/spatial-frequency and it enjoys both high locality and resolution. The Split-and-Merge algorithm tracks nonstationarities by outlining areas of homogeneity. Experiments were run on texture fields and the results are promising. We conclude by discussing issues related to our new approach and by suggesting directions for future research.
Zusammenfassung. Wir betrachten das Problem der Segmentation von Bildtexturen als Problem der Verfolgung von Nichtstationarit~iten in zweidimensionalen Signalen. Das Signal wird in Form der Wignerverteilung dargestellt. Diese ist eine Verbunddarsteilung von Ort und Ortsfrequenz; sic besitzt in beiden Bereichen ein hohes Aufl6sungsverm/bgen. Der Split-andMerge Algorithmus verfolgt Nichtstationarit,;iten, indem er homogene Gebiete kennzeichnet. Die experimentelle Untersuchung verschiedener Texturfelder ergibt vielversprechende Resultate. AbschlieBend diskutieren wit einige Gesichtspunkte in Verbindung mit unseren neuen Verfahren sowie m6gliche Fragestellungen fiir kiJnftige Untersuchungen. R6sum6. Nous consid6rons le probl~me de la segmentation de texture comme celui de la poursuite des non stationnarit6s pour des signaux bidimensionnels. La repr6sentation des signaux sous-jacente est la distribution de Wigner. C'est une repr6sentation conjointe espace/fr6quences spatiales et elle est pr6cise en position et en r6solution. L'algorithme de 'split-andmerge' suit ies non stationnadt6s en d61imitant les r6gions uniformes. Des exp6riences ont 6t6 effectu6es sur des images de texture et les r6sultats sont prometteurs. Nous concluons en discutant des cons6quences li6es ~ notre nouvelle approche et en sugg6rant des directions de recherche futures.
Keywords. Wigner distribution, texture segmentation, nonstationary signal.
1. Introduction Texture is often described as being generated by some basic elements referred to as primitives. These elements are repeated at positions and orientations determined by some placement rules. While this definition may be useful when describing artificially, man-made textures, such as checkerboards and tilings, it is quite difficult to * This work was supported in part by DARPA and monitored through the MIT Lincoln Laboratory under Contracts Nos. CX-7401 and CX-7916.
apply it to natural textures, which exhibit random appearance. In fact, textures may be found in a continuous spectrum ranging from purely deterministic (regular grids, etc.) to purely stochastic (clouds, etc.). Texture analysis can be considered at three basic levels of increasing difficulty: (a) classification/discrimination, (b) description, and (c) segmentation. Classification is the determination of to which of a number of known types of texture a given sample belongs. Discrimination involves the task of deciding if two samples come
0165-1684/88/$3.50 © 1988, Elsevier Science Publishers B.V. (North-Holland)
96
T. Reed, H. Wechsler / Texture segmentation
from the same texture. Both tasks involve the formulation of feature-based models, the parameters of which may be used to classify/discriminate texture samples. This is the most straightforward area of texture analysis, and the one in which the majority of activity has concentrated [34]. Description is more difficult to accomplish, as it is easy to find textures that are perceptually distinctive, but whose differences are difficult to quantify. Finally, texture segmentation, the most difficult level of texture analysis, is the partitioning of an image based on texture. The detection of boundaries can be thought of as tracking the transitions between approximately stationary regions of a nonstationary 2-D signal. It should be pointed out that image segmentation based on the picture value function is a specific case of texture segmentation. Specifically, nonuniform gray levels, like salt-andpepper, can still belong to the same area and as such no segmentation should occur. A formal definition for segmentation [29] is given below. 1.1. Definition. Let X denote the grid of sample points for a given image. Let Y be a subset of X containing at least two points. Then, a uniformity predicate P ( Y ) is one which assigns the value true or false to Y depending only on properties o f the picture value function f ( i , j ) as defined over Y. Furthermore, P has the property that if Z is a nonempty subset of Y, then P ( Y ) = T-> P ( Z ) = T.
on the way the 'seeds' are planted and the way the regions are grown and tracked. Complexity issues related to nondecidability of a predicate being uniform and the number of possible segmentations for a given uniform predicate are further discussed by Gurari and Wechsler [16]. There are few methods attempting texture segmentation [11, 12, 32]. Segmentation methods may range from fully supervised to fully unsupervised in nature. Most research to date deals with the supervised case where a priori modelling about expected textures is provided. Bayesian probability and the maximum likelihood approach are then used as is the case in remote sensing. Intermediate methods assume no prior knowledge about the number of textures in the image, or their identity. It is this level of texture segmentation that we wish to address in this paper. The only available attempt to achieve segmentation at this level is that by Pentland [32] using the fractal concept. However, Peleg [31] disputes the idea that textures can be modelled as fractals. Fully unsupervised texture segmentation (UTS), in which absolutely no information is known a priori (e.g., about potential boundary locations), and which involves no human intervention (e.g., through selection of thresholds) remains essentially unexplored, and it is the goal of this paper to start research toward solving it.
2. Process modelling 1.2. Definition. A segmentation of the grid X for a uniformity predicate P is a partition of X into disjoint nonempty subsets {Xi}i"_-~ such that: (a) [._]X , = X, (b) Vl<~i<~n, X~ is connected and P ( X i ) = T, (c) P is false on the union of any number of adjacent members of the induced partition. It should be noted that Definition 1.2 does not imply uniqueness of the segmentation or that the number o f regions, n, is as small as possible, unless some minimizing criteria, like in simulated annealing, is further imposed. The segmentation depends Signal Processing
The seminal work of Marr [24] considered computational vision (CV) as an informationprocessing task and it defines the three levels at which any (machine) vision system carrying out a visual task must be understood. There is a first, basic computational theory level which specifies the goal of computation, why it is appropriate, and what the strategy is by which it can be carried out. Next levels deal with representation (REP) and algorithm (ALG) (i.e., how is the computational theory implemented in terms of input/output representations and transformational procedures),
T. Reed, H. Wechsler / Texture segmentation
and hardware realization. Most of the state-of-theart research in CV comes from a synergetic approach which draws heavily from both psychophysics and neurophysiology. The segmentation task (Definition 1.2) suggests choosing a representation which can implicitly encode nonstationarities (i.e., signal changes) and an algorithm which can keep track of such changes and make a record of them. The choice of an adequate representation is tightly coupled to the uniformity predicate P. In texture segmentation, and especially in UTS, the task is further complicated by the cell unit problem [37]. The image resolution is intimately linked to the representation being evaluated, since computation is valid only over a specific range of resolution windows. There seems to be a consensus regarding the attributes required for adequate representations, even that the agreement does not extend to how the attributes can be attained. Specifically, one should consider the following arguments: (a) Multiresolution. This idea comes from the belief in the availability of 4-5 channels for human processing [40]. It has been implemented through the use of the pyramid concept and it has been used by Witkin [41] for signal analysis. Such a concept as multiresolution also naturally leads to spectral analysis. (b) Locality. There is enough evidence to believe that spectral analysis plays a major role in (human) vision [33, 36], dissent notwithstanding. However, empirical studies [10, 38] showed that frequency analysis compared poorly against other methods. The reason behind the poor performance is the lack of locality in analysis. The work done by Gagalowicz [14], and more recently by Caelli [5] and Julesz [21] showed conclusively the relevance of local analysis for texture synthesis/ analysis. In retrospect, one can now understand why statistics (like the 2nd order), equivalent to power spectral analysis, performed better. The reason is that the statistics were estimated on a local basis rather than on a whole image. The short-term periodogram is one attempt to rectify the lack of locality for spectral analysis. The Wigner distribu-
97
tion (WD), to be discussed in the next section, is one further step towards improved locality and resolution [20, 22]. There are two types of algorithms which can be used on representations (fitting the above specifications) for tracking signal changes and to yield an appropriate segmentation. Split-andMerge (S&M) [6, 17] and relaxation are the two generic l~ypes of algorithms. The use of S&M is detailed in the next section while reasons for using relaxation are discussed in the concluding section.
3. Representation and algorithm We consider texture as a 2-D nonstationary random process, i.e., a process with space varying statistics or, alternatively, space-varying frequency contents. The analysis of nonstationary processes is not new. For the 1-D case (e.g., speech analysis), the short-time periodogram [1] has been used for this purpose. However, it has recently been shown [26] that the estimation of time-varying frequency characteristics (again for I-D) can be better achieved via the conjoint time-frequency representation referred to as the Wigner-Ville spectrum, or Wigner distribution (WD) [7, 8, 9, 35, 39]. The use of the WD for 2-D image processing was first suggested by Jacobson and Wechsler [18, 19]. The formal WD definition is given below:
Wf(x,u)= f f~o f(x+ lot)f*(x-½ot) × exp{-jau} dot. The WD can track nonstationarities as seen in the next example. Assume the function
f(t) = A exp{jate/2}. This function is often referred to as a chirp. Its instantaneous frequency increases linearly with time. The corresponding WD is given by
Wr(t, o,) = Ial=2~(o, -
at).
Thus, the WD of a chirp is concentrated about the instantaneous frequency for any time, as desired. Vol. 14, No. 1, January 1988
98
T. Reed, H. Wechsler / Texture segmentation I
Fig. 1. The composite image of Example 1.
Fig. 3. The lowpass filtered image, as input to the PWD.
T h e n e e d to o p e r a t e o n d a t a a n d / o r a n a l y t i c a l c o m p l e x i t y r e q u i r e s e s t i m a t i n g the W D . W e refer to these n u m e r i c a l l y c a l c u l a t e d r e p r e s e n t a t i o n s as pseudo-Wigner distributions (PWD). A PWD estimator used successfully by Martin and F l a n d r i n [26, 27] to t r a c k 1-D signal n o n s t a t i o n arities, was e x t e n d e d b y the p r e s e n t a u t h o r s to 2-D Signal Processing
Fig. 2. The 64 x 64 pixel image as input to the anti-aliasingfilter.
Fig. 4. The segmented 32 x 32 image with PWD edge effect removed, as found using Method 1. a n d is s h o w n b e l o w : P W ( m , n, p, q) =4
N2--1
NI-I
E
~,
k=-- N2+ 1 I=-Nt+1 M2--1 MI--I
x
~.
~
r=--M2+l s=--Ml+l
IhN,,N2(k,/)l
2
g~,,M~(r, s)
T. Reed, H. Wechsler / Texture segmentation
99
Fig. 5. The composite image of Example 2 (256 x 256).
Fig. 6. The 64 x 64 pixel image as input to the anti-aliasmg tilter.
Fig. 7. The image after lowpass filtering.
Fig. 8. The segmented 32 x32 center of the image as found using Method 2.
xf(m+r+k, n+s+l) xf*(m+r-k,n+s-l) f .[2~rkp 27rlq\) x exp/-Jk
P + O )J'
where p = 0 , + l . . . . , + ( N 2 - 1), q = 0 , + 1 , . . . , +(N1-1), P=2N2-1,
Q=2N~-I, Vol. 14, No. 1, January 1988
100
T. Reed, H. Wechsler / Texture segmentation
m and n are integers, and n and g are window functions. There are a number of issues related to the use of the WD and they are considered next. (a) Aliasing. As seen from the PWD definition, the product of the data and its complex conjugate is calculated. The spectrum of this product is twice as wide as the spectrum of the original data, yet its sampling period is the same. If the original data is sampled at the Nyquist rate, aliasing will result. The method we use is to bandlimit the data so that the sampling rate is twice the Nyquist rate. (We could blur the image by defocussing the camera and achieve the same result.) For future work we are considering the use o f the Hilbert transform [28], implicit spline interpolation [2] a n d / o r the CSF (contrast sensitivity function). (b) Window functions. The PWD is calculated using normed rectangular windows. Improvement in performance, like avoiding edge artifacts, can be expected from using Gaussian a n d / o r Kaiser windows. (c) Computational complexity. The PWD definition is expensive in terms of computer implementation. Furthermore, the data storage requirements place an additional burden. Data compression thus seems to be a real need once the PWD has been calculated. KLT/eigenvalue analysis [23] was suggested for WD. Other methods like the maximum spectral component(s) [41] a n d / o r a fractal approximation will be considered as well. The fractal approximation will enjoy some invariant characteristics, like those related to perspective. In future applications, the highly parallel nature of the PWD calculation may be exploited, so that processing time may be reduced by nearly a factor o f N for N parallel processors. A systolic processor for computing the PWD for onedimensional signals is discussed in [13]. The PWD representation is stored in a quadtree structure. A Split-and-Merge algorithm [30] is then used for segmentation. The uniformity predicates P consider homogeneity in terms of similarity of the local spectral signatures as defined by PWD. The thresholds used are a function of the tree level Signal Processing
occupied by the candidate regions, with the threshold decreasing as one moves toward the top of the tree. The varying threshold accounts for increased averaging over larger candidate regions. If one considers the difference between spectral information of region pairs at a given tree level as a random variable of the Normal distribution, then the threshold for that level falls within the 80-90% range of the tail distribution.
4. Experimental results We implemented a system for texture segmentation which incorporates the ideas of the previous section concerning representation and algorithm. Specifically, there are three basic steps: (i) Preprocessing ( anti-aliasing ). (ii) Local spectra are derived from the WD. The first alternative, Method 1, considers the average local spectra for the uniformity predicate P. Regions are merged when the total difference between the average frequency content, averaged over each of the candidate regions and summed over the frequencies in the discrete spectrum, falls below a threshold. The second approach, Method 2, relates P to the variance of those harmonics which contain most of the energy. Regions are merged when the variance of the frequency content over the merged regions falls below a threshold for all of the main harmonics. (iii) Split-and-Merge segmentation algorithm, where regions are split or merged based on homogeneity with respect to P. The threshold is a function of the quadtree level. In Example 1, we consider an image composed of textures from [4]. Fig. 1 shows a 128 x 128 pixel sample of each texture. Counter-clockwise from the upper right, the textures are of paper, paper, grass, and fur (D57, D57, D9, and D93 in [4]). Fig. 2 is a 64 x 64 version of the image, at the same resolution. This is the image input to the antialiasing filter. Fig. 3 is the input to the PWD, filtered to prevent aliasing. Fig. 4 is the segmented 32 x 32 image, as found by using Method 1. This
T. Reed, H. Wechsler / Texture segmentation
smaller image results from removing the edges of the 64 x 64 image prior to merging, in order to eliminate edge effects due to the PWD. Example 2 involves another composite of textures from [4]. Fig. 5 is a 256 x 256 image illustrating the textures used. Counter-clockwise from the upper right, the textures are of wood, paper, and fur (D71, D57, and D93 in [4]). Fig. 6 is a 64 x 64 version of the image, of the same resolution as the previous figure. This image is lowpass filtered, with the result shown in Fig. 7. Finally, Fig. 8 is the segmented (32 x 32) center of the image, as found using Method 2. Although perfect segmentation was not achieved, the results obtained are still impressive. We allowed for any number of full gray scale textures and no a priori knowledge a n d / o r modelling for the textures which might be present in the image. Due to computational reasons we were limited to small image sizes which severely restrict the amount of detail that can be examined and used for fine discrimination.
5. Conclusions
In this paper we have considered the problem of texture segmentation. The main component of the approach suggested here is that of analyzing textural information by using local spectral representations based on the Wigner distribution. The experiments reported showed that our method holds promise. The most obvious areas of future work are the application of our method to larger and more complex images, data compression, and possibly to UTS. There are several modifications to the segmentation algorithms that are o f immediate interest as well. The use of the Hilbert transformation to calculate the 'analytic image' as a means to reduce aliasing, rather than lowpass filters, will retain detail that is now lost. Another area of contemplated change is that of using relaxation as suggested by Caelli [5] and Grossberg and Mingolla
101
[15] instead of Split-and-Merge. (The Split-andMerge processes all levels of the multisolution pyramid level in a constant fashion and averaging (i.e., loss of detail) is the means for moving up the tree (i.e., to merge).) Finally, we plan to consider the issues of invariance and data fusion. The results reported here should be made invariant to linear transformations and perspective distortion. Data fusion is concerned with making the segmentation process dependent on several sensory representations in retinotopic correspondence.
References [1] J. Allen and L. Rabiner, "A unified approach to short-time Fourier analysis and synthesis", Proc. IEEE, Vol. 65, No. 11, November 1977, pp. 1558-1564. [2] G. Boudreaux-Bartels and T.W. Parks, "Reducing aliasing in the Wigner distribution using implicit spline interpolation", Proc. IEEE Internat. Conf. on Acoustics, Speech, and SignalProcessing, Boston, MA, April 1983, pp. 1438-1441. [3] G.F. Boudreaux-Bartels and T.W. Parks, "Time-varying filtering and signal estimation using Wigner distribution synthesis techniques", IEEE Trans. Acoust., Speech, Signal Process., Vol. ASSP-34, No. 3, June 1986, pp. 442-451. [4] P. Brodatz, Textures--A Photographic Album for Artists and Designers, Dover, New York, 1966. [5] T. Caelli, "'Three processing characteristics of visual texture segmentation", Spatial Vision, Vol. 1, No. 1, 1985, pp. 19-30. [6] P. Chen and T. Pavlidis, "Segmentation by texture using a co-occurrence matrix and a split-and-merge algorithm", Proc. 4th Internat. Joint Conf. on Pattern Recognition, Kyoto, Japan, November 7-10, 1978, pp. 565-569. [7] T.A.C.M. Claasen and W.F.G. Mecklenbrauker, "The Wigner distribution--A tool for time-frequency signal analysis, Part I: Continuous-time signals", Philips J. Res., Vol. 35, No. 3, 1980, pp. 217-250. [8] T.A.C.M. Claasen and W.F.G. Mecklenbrauker, "The Wigner distribution--A tool for time-frequency signal analysis, Part II: Discrete-time signals", Philips J. Res., Vol. 35, Nos. 4/5, 1980, pp. 276-300. [9] T.A.C.M. Claasen and W.F.G. Mecklenbrauker, "The Wigner distribution--A tool for time-frequency signal analysis, Part III: Relations with other time-frequency signal transformations", Philips J. Res., Vol. 35, No. 6, 1980, pp. 372-389. [10] R. Connors and C. Harlow, "A theoretical comparison of texture algorithms", [EEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-2, No. 3, May 1980, pp. 204-222. Vol. 14, No. 1, January 1988
102
T. Reed, H. Wechsler / Texture segmentation
[11] H. Derin and W. Cole, "Segmentation of textured images using Gibbs random fields", Computer Vision, Graphics & Image Process., Vol. 35, 1986, pp. 72-98. [12] H. Derin, H. Elliot, R. Cristi and D. Geman, "Bayes smoothing algorithms for segmentation of binary images modeled by Markov random fields", IEEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-6, No. 6, November 1984, pp. 707-720. [13] T. Durrani, R. Chapman and T. Willey, "Systolic processor for computing the Wigner distribution", Electronic Lett., Vol. 19, No. 13, June 1983, pp. 476-477. [14] A. Gagalowicz, "A new method for texture fields synthesis: Some applications to the study of human vision", IEEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-3, No. 5, September 1981, pp. 520-533. [15] S. Grossberg and E. Mingolla, "Neutral dynamics of perceptual grouping: Textures, boundaries, and emergent segmentations", Perception & Psychophysics, Vol. 38, No. 2, 1985, pp. 141-171. [16] E. Gurari and H. Wechsler, "On the difficulties involved in the segmentation of pictures", IEEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-4, No. 3, May 1982, pp. 304-306. [17] S. Horowitz and T. Pavlidis, "Picture segmentation by a tree traversal algorithm", Z ACM, Vol. 23, No. 2, April 1976, pp. 368-388. [18] L. Jacobson and H. Wechsler, "The Wigner distribution and its usefulness for 2-D image processing", in: Proc. 6th lnternat. Joint Conf. on Pattern Recognition, Munich, Fed. Rep. Germany, October 1982. [19] L. Jacobson and H. Wechsler, "A theory for invariant object recognition in the frontoparallel plane", IEEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-6, No. 3, May 1984, pp. 325-331. [20] L. Jacobson and H. Wechsler, "Joint spatial/spatialfrequency representations for image processing", Proc. SPIE / Cambridge Internat. Conf. on Intelligent Robots and Computer Vision, Cambridge, MA, September 1985, pp. 11-18. [21] B. Julesz and J. Bergen, "Textons, the fundamental elements in preattentive vision and perception of textures", Bell Syst. Tech. J., Vol. 62, No. 6, July-August 1983, pp. 1619-1645. [22] J.J. Kulikowski, S. Marcelja and P.O. Bishop, "Theory of spatial position and spatial frequency relations in the respective fields of simple cells in the visual cortex", Biol. Cyber., Vol. 43, 1982, pp. 187-198. [23] N.M. Marinovic and G. Eichmann, "Feature extraction and pattern classification in space-spatial frequency domain", Proc. SPIE / Cambridge Internat. Conf. on Intelligent Robots and Computer Vision, Cambridge, MA, September 1985, pp. 19-26. [24] D. Marr, Vision, Freeman, New York, 1982.
Signal Processing
[25] W. Martin, "Measuring the degree of non-stationarity by using the Wigner-Ville spectrum", Proc. Internat. Conf. on Acoustics, Speech, and Signal Processing, San Diego, CA, March 19-21, 1984, pp. 41B.3.1-4. [26] W. Martin and P. Flandrin, "Analysis of non-stationary processes: Short-time periodograms versus a pseudoWigner estimator", in: H. Schussler, editor, EUSIPCO-83, North-Holland, Amsterdam, 1983, pp. 455-458. [27] W. Martin and P. Flandrin, "Detection of changes of signal structure by using the Wigner-Ville spectrum", Signal Processing, Vol. 8, No. 2, April 1985, pp. 215-233. [28] W. Martin and P. Flandrin, "Wigner-Ville spectral analysis of nonstationary processes", IEEE Trans. Acoust., Speech, Signal Process., Vol. 33, No. 6, December 1985, pp. 1461-1470. [29] T. Pavlidis, Structural Pattern Recognition, Springer, New York, 1977, Chap. 4, pp. 68-70. [30] T. Pavlidis, Algorithms for Graphics and Image Processing, Computer Science Press, Rockville, MD, 1982, Chap. 6, pp. 113-116. [31] S. Peleg, J. Naor, R. Harley and D. Avnir, Multiple Resolution Texture Analysis and Classification, Tech. Rept. TR1306, Dept. of Computer Science, Univ. of Maryland, College Park, MD, 1983. [32] A. Pentland, "Fractal-based description of natural scenes", IEEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-6, No. 6, November 1984, pp. 661-674. [33] D. Pollen and S. Ronner, "Visual cortical neurons as localized spatial frequency filters", IEEE Trans. Syst. Man & Cybernet., Vol. SMC-13, No. 5, September/October 1983, pp. 907-915. [34] L. Van Gool, P. Dewaele and A. Oosterlinck, "Texture analysis anno 1983", Computer Vision, Graphics & Image Process., Vol. 29, 1985, pp. 336-357. [35] J. Ville, "Theorie et applications de la notion de signal analytiqne", Cables et Trans., Vol. 2, No. 1, 1948, pp. 61-74. [36] A. Watson, H. Barlow and J. Robson, "What does the eye see best?" Nature, Vol. 302, 1983, pp. 419-422. [37] H. Wechsler, "Texture analysis--A survey", Signal Processing, Vol. 2, No. 3, July 1980, pp. 271-282. [38] J.S. Weszka, C.R. Dyer and A. Rosenfeld, "A comparison study of texture measures for terrain classification", IEEE Trans. Syst. Man & Cybernet., Vol. SMC-6, No. 4, 1976, pp. 269-286. [39] E. Wigner, "On the quantum correction for thermodynamic equilibrium", Phys. Rev., Vol. 40, June 1932, pp. 749-759. [40] H.R. Wilson and J.R. Bergen, "A four mechanism model for threshold spatial vision, Vision Research, Vol. 19, 1979, pp. 19-32. [41] A.P. Witkin, "Recovering surface shape and orientation from texture", Artificial Intelligence (Special Volume on Computer Vision), Vol. 17, Nos. 1-3, August 1981, pp. 17-45.
95
SHORT COMMUNICATION
TRACKING OF NONSTATIONARITIES FOR TEXTURE FIELDS* Todd REED and Harry WECHSLER Department of Electrical Engineering, University of Minnesota, Minneapolis, M N 55455, U.S.A. Received 19 November 1986 Revised 11 February 1987 and 1 July 1987
Abstract. We consider the problem of texture segmentation as one of tracking nonstationarities for 2-D signals. The underlying signal representation is the Wigner distribution. It is a co-joint representation of space/spatial-frequency and it enjoys both high locality and resolution. The Split-and-Merge algorithm tracks nonstationarities by outlining areas of homogeneity. Experiments were run on texture fields and the results are promising. We conclude by discussing issues related to our new approach and by suggesting directions for future research.
Zusammenfassung. Wir betrachten das Problem der Segmentation von Bildtexturen als Problem der Verfolgung von Nichtstationarit~iten in zweidimensionalen Signalen. Das Signal wird in Form der Wignerverteilung dargestellt. Diese ist eine Verbunddarsteilung von Ort und Ortsfrequenz; sic besitzt in beiden Bereichen ein hohes Aufl6sungsverm/bgen. Der Split-andMerge Algorithmus verfolgt Nichtstationarit,;iten, indem er homogene Gebiete kennzeichnet. Die experimentelle Untersuchung verschiedener Texturfelder ergibt vielversprechende Resultate. AbschlieBend diskutieren wit einige Gesichtspunkte in Verbindung mit unseren neuen Verfahren sowie m6gliche Fragestellungen fiir kiJnftige Untersuchungen. R6sum6. Nous consid6rons le probl~me de la segmentation de texture comme celui de la poursuite des non stationnarit6s pour des signaux bidimensionnels. La repr6sentation des signaux sous-jacente est la distribution de Wigner. C'est une repr6sentation conjointe espace/fr6quences spatiales et elle est pr6cise en position et en r6solution. L'algorithme de 'split-andmerge' suit ies non stationnadt6s en d61imitant les r6gions uniformes. Des exp6riences ont 6t6 effectu6es sur des images de texture et les r6sultats sont prometteurs. Nous concluons en discutant des cons6quences li6es ~ notre nouvelle approche et en sugg6rant des directions de recherche futures.
Keywords. Wigner distribution, texture segmentation, nonstationary signal.
1. Introduction Texture is often described as being generated by some basic elements referred to as primitives. These elements are repeated at positions and orientations determined by some placement rules. While this definition may be useful when describing artificially, man-made textures, such as checkerboards and tilings, it is quite difficult to * This work was supported in part by DARPA and monitored through the MIT Lincoln Laboratory under Contracts Nos. CX-7401 and CX-7916.
apply it to natural textures, which exhibit random appearance. In fact, textures may be found in a continuous spectrum ranging from purely deterministic (regular grids, etc.) to purely stochastic (clouds, etc.). Texture analysis can be considered at three basic levels of increasing difficulty: (a) classification/discrimination, (b) description, and (c) segmentation. Classification is the determination of to which of a number of known types of texture a given sample belongs. Discrimination involves the task of deciding if two samples come
0165-1684/88/$3.50 © 1988, Elsevier Science Publishers B.V. (North-Holland)
96
T. Reed, H. Wechsler / Texture segmentation
from the same texture. Both tasks involve the formulation of feature-based models, the parameters of which may be used to classify/discriminate texture samples. This is the most straightforward area of texture analysis, and the one in which the majority of activity has concentrated [34]. Description is more difficult to accomplish, as it is easy to find textures that are perceptually distinctive, but whose differences are difficult to quantify. Finally, texture segmentation, the most difficult level of texture analysis, is the partitioning of an image based on texture. The detection of boundaries can be thought of as tracking the transitions between approximately stationary regions of a nonstationary 2-D signal. It should be pointed out that image segmentation based on the picture value function is a specific case of texture segmentation. Specifically, nonuniform gray levels, like salt-andpepper, can still belong to the same area and as such no segmentation should occur. A formal definition for segmentation [29] is given below. 1.1. Definition. Let X denote the grid of sample points for a given image. Let Y be a subset of X containing at least two points. Then, a uniformity predicate P ( Y ) is one which assigns the value true or false to Y depending only on properties o f the picture value function f ( i , j ) as defined over Y. Furthermore, P has the property that if Z is a nonempty subset of Y, then P ( Y ) = T-> P ( Z ) = T.
on the way the 'seeds' are planted and the way the regions are grown and tracked. Complexity issues related to nondecidability of a predicate being uniform and the number of possible segmentations for a given uniform predicate are further discussed by Gurari and Wechsler [16]. There are few methods attempting texture segmentation [11, 12, 32]. Segmentation methods may range from fully supervised to fully unsupervised in nature. Most research to date deals with the supervised case where a priori modelling about expected textures is provided. Bayesian probability and the maximum likelihood approach are then used as is the case in remote sensing. Intermediate methods assume no prior knowledge about the number of textures in the image, or their identity. It is this level of texture segmentation that we wish to address in this paper. The only available attempt to achieve segmentation at this level is that by Pentland [32] using the fractal concept. However, Peleg [31] disputes the idea that textures can be modelled as fractals. Fully unsupervised texture segmentation (UTS), in which absolutely no information is known a priori (e.g., about potential boundary locations), and which involves no human intervention (e.g., through selection of thresholds) remains essentially unexplored, and it is the goal of this paper to start research toward solving it.
2. Process modelling 1.2. Definition. A segmentation of the grid X for a uniformity predicate P is a partition of X into disjoint nonempty subsets {Xi}i"_-~ such that: (a) [._]X , = X, (b) Vl<~i<~n, X~ is connected and P ( X i ) = T, (c) P is false on the union of any number of adjacent members of the induced partition. It should be noted that Definition 1.2 does not imply uniqueness of the segmentation or that the number o f regions, n, is as small as possible, unless some minimizing criteria, like in simulated annealing, is further imposed. The segmentation depends Signal Processing
The seminal work of Marr [24] considered computational vision (CV) as an informationprocessing task and it defines the three levels at which any (machine) vision system carrying out a visual task must be understood. There is a first, basic computational theory level which specifies the goal of computation, why it is appropriate, and what the strategy is by which it can be carried out. Next levels deal with representation (REP) and algorithm (ALG) (i.e., how is the computational theory implemented in terms of input/output representations and transformational procedures),
T. Reed, H. Wechsler / Texture segmentation
and hardware realization. Most of the state-of-theart research in CV comes from a synergetic approach which draws heavily from both psychophysics and neurophysiology. The segmentation task (Definition 1.2) suggests choosing a representation which can implicitly encode nonstationarities (i.e., signal changes) and an algorithm which can keep track of such changes and make a record of them. The choice of an adequate representation is tightly coupled to the uniformity predicate P. In texture segmentation, and especially in UTS, the task is further complicated by the cell unit problem [37]. The image resolution is intimately linked to the representation being evaluated, since computation is valid only over a specific range of resolution windows. There seems to be a consensus regarding the attributes required for adequate representations, even that the agreement does not extend to how the attributes can be attained. Specifically, one should consider the following arguments: (a) Multiresolution. This idea comes from the belief in the availability of 4-5 channels for human processing [40]. It has been implemented through the use of the pyramid concept and it has been used by Witkin [41] for signal analysis. Such a concept as multiresolution also naturally leads to spectral analysis. (b) Locality. There is enough evidence to believe that spectral analysis plays a major role in (human) vision [33, 36], dissent notwithstanding. However, empirical studies [10, 38] showed that frequency analysis compared poorly against other methods. The reason behind the poor performance is the lack of locality in analysis. The work done by Gagalowicz [14], and more recently by Caelli [5] and Julesz [21] showed conclusively the relevance of local analysis for texture synthesis/ analysis. In retrospect, one can now understand why statistics (like the 2nd order), equivalent to power spectral analysis, performed better. The reason is that the statistics were estimated on a local basis rather than on a whole image. The short-term periodogram is one attempt to rectify the lack of locality for spectral analysis. The Wigner distribu-
97
tion (WD), to be discussed in the next section, is one further step towards improved locality and resolution [20, 22]. There are two types of algorithms which can be used on representations (fitting the above specifications) for tracking signal changes and to yield an appropriate segmentation. Split-andMerge (S&M) [6, 17] and relaxation are the two generic l~ypes of algorithms. The use of S&M is detailed in the next section while reasons for using relaxation are discussed in the concluding section.
3. Representation and algorithm We consider texture as a 2-D nonstationary random process, i.e., a process with space varying statistics or, alternatively, space-varying frequency contents. The analysis of nonstationary processes is not new. For the 1-D case (e.g., speech analysis), the short-time periodogram [1] has been used for this purpose. However, it has recently been shown [26] that the estimation of time-varying frequency characteristics (again for I-D) can be better achieved via the conjoint time-frequency representation referred to as the Wigner-Ville spectrum, or Wigner distribution (WD) [7, 8, 9, 35, 39]. The use of the WD for 2-D image processing was first suggested by Jacobson and Wechsler [18, 19]. The formal WD definition is given below:
Wf(x,u)= f f~o f(x+ lot)f*(x-½ot) × exp{-jau} dot. The WD can track nonstationarities as seen in the next example. Assume the function
f(t) = A exp{jate/2}. This function is often referred to as a chirp. Its instantaneous frequency increases linearly with time. The corresponding WD is given by
Wr(t, o,) = Ial=2~(o, -
at).
Thus, the WD of a chirp is concentrated about the instantaneous frequency for any time, as desired. Vol. 14, No. 1, January 1988
98
T. Reed, H. Wechsler / Texture segmentation I
Fig. 1. The composite image of Example 1.
Fig. 3. The lowpass filtered image, as input to the PWD.
T h e n e e d to o p e r a t e o n d a t a a n d / o r a n a l y t i c a l c o m p l e x i t y r e q u i r e s e s t i m a t i n g the W D . W e refer to these n u m e r i c a l l y c a l c u l a t e d r e p r e s e n t a t i o n s as pseudo-Wigner distributions (PWD). A PWD estimator used successfully by Martin and F l a n d r i n [26, 27] to t r a c k 1-D signal n o n s t a t i o n arities, was e x t e n d e d b y the p r e s e n t a u t h o r s to 2-D Signal Processing
Fig. 2. The 64 x 64 pixel image as input to the anti-aliasingfilter.
Fig. 4. The segmented 32 x 32 image with PWD edge effect removed, as found using Method 1. a n d is s h o w n b e l o w : P W ( m , n, p, q) =4
N2--1
NI-I
E
~,
k=-- N2+ 1 I=-Nt+1 M2--1 MI--I
x
~.
~
r=--M2+l s=--Ml+l
IhN,,N2(k,/)l
2
g~,,M~(r, s)
T. Reed, H. Wechsler / Texture segmentation
99
Fig. 5. The composite image of Example 2 (256 x 256).
Fig. 6. The 64 x 64 pixel image as input to the anti-aliasmg tilter.
Fig. 7. The image after lowpass filtering.
Fig. 8. The segmented 32 x32 center of the image as found using Method 2.
xf(m+r+k, n+s+l) xf*(m+r-k,n+s-l) f .[2~rkp 27rlq\) x exp/-Jk
P + O )J'
where p = 0 , + l . . . . , + ( N 2 - 1), q = 0 , + 1 , . . . , +(N1-1), P=2N2-1,
Q=2N~-I, Vol. 14, No. 1, January 1988
100
T. Reed, H. Wechsler / Texture segmentation
m and n are integers, and n and g are window functions. There are a number of issues related to the use of the WD and they are considered next. (a) Aliasing. As seen from the PWD definition, the product of the data and its complex conjugate is calculated. The spectrum of this product is twice as wide as the spectrum of the original data, yet its sampling period is the same. If the original data is sampled at the Nyquist rate, aliasing will result. The method we use is to bandlimit the data so that the sampling rate is twice the Nyquist rate. (We could blur the image by defocussing the camera and achieve the same result.) For future work we are considering the use o f the Hilbert transform [28], implicit spline interpolation [2] a n d / o r the CSF (contrast sensitivity function). (b) Window functions. The PWD is calculated using normed rectangular windows. Improvement in performance, like avoiding edge artifacts, can be expected from using Gaussian a n d / o r Kaiser windows. (c) Computational complexity. The PWD definition is expensive in terms of computer implementation. Furthermore, the data storage requirements place an additional burden. Data compression thus seems to be a real need once the PWD has been calculated. KLT/eigenvalue analysis [23] was suggested for WD. Other methods like the maximum spectral component(s) [41] a n d / o r a fractal approximation will be considered as well. The fractal approximation will enjoy some invariant characteristics, like those related to perspective. In future applications, the highly parallel nature of the PWD calculation may be exploited, so that processing time may be reduced by nearly a factor o f N for N parallel processors. A systolic processor for computing the PWD for onedimensional signals is discussed in [13]. The PWD representation is stored in a quadtree structure. A Split-and-Merge algorithm [30] is then used for segmentation. The uniformity predicates P consider homogeneity in terms of similarity of the local spectral signatures as defined by PWD. The thresholds used are a function of the tree level Signal Processing
occupied by the candidate regions, with the threshold decreasing as one moves toward the top of the tree. The varying threshold accounts for increased averaging over larger candidate regions. If one considers the difference between spectral information of region pairs at a given tree level as a random variable of the Normal distribution, then the threshold for that level falls within the 80-90% range of the tail distribution.
4. Experimental results We implemented a system for texture segmentation which incorporates the ideas of the previous section concerning representation and algorithm. Specifically, there are three basic steps: (i) Preprocessing ( anti-aliasing ). (ii) Local spectra are derived from the WD. The first alternative, Method 1, considers the average local spectra for the uniformity predicate P. Regions are merged when the total difference between the average frequency content, averaged over each of the candidate regions and summed over the frequencies in the discrete spectrum, falls below a threshold. The second approach, Method 2, relates P to the variance of those harmonics which contain most of the energy. Regions are merged when the variance of the frequency content over the merged regions falls below a threshold for all of the main harmonics. (iii) Split-and-Merge segmentation algorithm, where regions are split or merged based on homogeneity with respect to P. The threshold is a function of the quadtree level. In Example 1, we consider an image composed of textures from [4]. Fig. 1 shows a 128 x 128 pixel sample of each texture. Counter-clockwise from the upper right, the textures are of paper, paper, grass, and fur (D57, D57, D9, and D93 in [4]). Fig. 2 is a 64 x 64 version of the image, at the same resolution. This is the image input to the antialiasing filter. Fig. 3 is the input to the PWD, filtered to prevent aliasing. Fig. 4 is the segmented 32 x 32 image, as found by using Method 1. This
T. Reed, H. Wechsler / Texture segmentation
smaller image results from removing the edges of the 64 x 64 image prior to merging, in order to eliminate edge effects due to the PWD. Example 2 involves another composite of textures from [4]. Fig. 5 is a 256 x 256 image illustrating the textures used. Counter-clockwise from the upper right, the textures are of wood, paper, and fur (D71, D57, and D93 in [4]). Fig. 6 is a 64 x 64 version of the image, of the same resolution as the previous figure. This image is lowpass filtered, with the result shown in Fig. 7. Finally, Fig. 8 is the segmented (32 x 32) center of the image, as found using Method 2. Although perfect segmentation was not achieved, the results obtained are still impressive. We allowed for any number of full gray scale textures and no a priori knowledge a n d / o r modelling for the textures which might be present in the image. Due to computational reasons we were limited to small image sizes which severely restrict the amount of detail that can be examined and used for fine discrimination.
5. Conclusions
In this paper we have considered the problem of texture segmentation. The main component of the approach suggested here is that of analyzing textural information by using local spectral representations based on the Wigner distribution. The experiments reported showed that our method holds promise. The most obvious areas of future work are the application of our method to larger and more complex images, data compression, and possibly to UTS. There are several modifications to the segmentation algorithms that are o f immediate interest as well. The use of the Hilbert transformation to calculate the 'analytic image' as a means to reduce aliasing, rather than lowpass filters, will retain detail that is now lost. Another area of contemplated change is that of using relaxation as suggested by Caelli [5] and Grossberg and Mingolla
101
[15] instead of Split-and-Merge. (The Split-andMerge processes all levels of the multisolution pyramid level in a constant fashion and averaging (i.e., loss of detail) is the means for moving up the tree (i.e., to merge).) Finally, we plan to consider the issues of invariance and data fusion. The results reported here should be made invariant to linear transformations and perspective distortion. Data fusion is concerned with making the segmentation process dependent on several sensory representations in retinotopic correspondence.
References [1] J. Allen and L. Rabiner, "A unified approach to short-time Fourier analysis and synthesis", Proc. IEEE, Vol. 65, No. 11, November 1977, pp. 1558-1564. [2] G. Boudreaux-Bartels and T.W. Parks, "Reducing aliasing in the Wigner distribution using implicit spline interpolation", Proc. IEEE Internat. Conf. on Acoustics, Speech, and SignalProcessing, Boston, MA, April 1983, pp. 1438-1441. [3] G.F. Boudreaux-Bartels and T.W. Parks, "Time-varying filtering and signal estimation using Wigner distribution synthesis techniques", IEEE Trans. Acoust., Speech, Signal Process., Vol. ASSP-34, No. 3, June 1986, pp. 442-451. [4] P. Brodatz, Textures--A Photographic Album for Artists and Designers, Dover, New York, 1966. [5] T. Caelli, "'Three processing characteristics of visual texture segmentation", Spatial Vision, Vol. 1, No. 1, 1985, pp. 19-30. [6] P. Chen and T. Pavlidis, "Segmentation by texture using a co-occurrence matrix and a split-and-merge algorithm", Proc. 4th Internat. Joint Conf. on Pattern Recognition, Kyoto, Japan, November 7-10, 1978, pp. 565-569. [7] T.A.C.M. Claasen and W.F.G. Mecklenbrauker, "The Wigner distribution--A tool for time-frequency signal analysis, Part I: Continuous-time signals", Philips J. Res., Vol. 35, No. 3, 1980, pp. 217-250. [8] T.A.C.M. Claasen and W.F.G. Mecklenbrauker, "The Wigner distribution--A tool for time-frequency signal analysis, Part II: Discrete-time signals", Philips J. Res., Vol. 35, Nos. 4/5, 1980, pp. 276-300. [9] T.A.C.M. Claasen and W.F.G. Mecklenbrauker, "The Wigner distribution--A tool for time-frequency signal analysis, Part III: Relations with other time-frequency signal transformations", Philips J. Res., Vol. 35, No. 6, 1980, pp. 372-389. [10] R. Connors and C. Harlow, "A theoretical comparison of texture algorithms", [EEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-2, No. 3, May 1980, pp. 204-222. Vol. 14, No. 1, January 1988
102
T. Reed, H. Wechsler / Texture segmentation
[11] H. Derin and W. Cole, "Segmentation of textured images using Gibbs random fields", Computer Vision, Graphics & Image Process., Vol. 35, 1986, pp. 72-98. [12] H. Derin, H. Elliot, R. Cristi and D. Geman, "Bayes smoothing algorithms for segmentation of binary images modeled by Markov random fields", IEEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-6, No. 6, November 1984, pp. 707-720. [13] T. Durrani, R. Chapman and T. Willey, "Systolic processor for computing the Wigner distribution", Electronic Lett., Vol. 19, No. 13, June 1983, pp. 476-477. [14] A. Gagalowicz, "A new method for texture fields synthesis: Some applications to the study of human vision", IEEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-3, No. 5, September 1981, pp. 520-533. [15] S. Grossberg and E. Mingolla, "Neutral dynamics of perceptual grouping: Textures, boundaries, and emergent segmentations", Perception & Psychophysics, Vol. 38, No. 2, 1985, pp. 141-171. [16] E. Gurari and H. Wechsler, "On the difficulties involved in the segmentation of pictures", IEEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-4, No. 3, May 1982, pp. 304-306. [17] S. Horowitz and T. Pavlidis, "Picture segmentation by a tree traversal algorithm", Z ACM, Vol. 23, No. 2, April 1976, pp. 368-388. [18] L. Jacobson and H. Wechsler, "The Wigner distribution and its usefulness for 2-D image processing", in: Proc. 6th lnternat. Joint Conf. on Pattern Recognition, Munich, Fed. Rep. Germany, October 1982. [19] L. Jacobson and H. Wechsler, "A theory for invariant object recognition in the frontoparallel plane", IEEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-6, No. 3, May 1984, pp. 325-331. [20] L. Jacobson and H. Wechsler, "Joint spatial/spatialfrequency representations for image processing", Proc. SPIE / Cambridge Internat. Conf. on Intelligent Robots and Computer Vision, Cambridge, MA, September 1985, pp. 11-18. [21] B. Julesz and J. Bergen, "Textons, the fundamental elements in preattentive vision and perception of textures", Bell Syst. Tech. J., Vol. 62, No. 6, July-August 1983, pp. 1619-1645. [22] J.J. Kulikowski, S. Marcelja and P.O. Bishop, "Theory of spatial position and spatial frequency relations in the respective fields of simple cells in the visual cortex", Biol. Cyber., Vol. 43, 1982, pp. 187-198. [23] N.M. Marinovic and G. Eichmann, "Feature extraction and pattern classification in space-spatial frequency domain", Proc. SPIE / Cambridge Internat. Conf. on Intelligent Robots and Computer Vision, Cambridge, MA, September 1985, pp. 19-26. [24] D. Marr, Vision, Freeman, New York, 1982.
Signal Processing
[25] W. Martin, "Measuring the degree of non-stationarity by using the Wigner-Ville spectrum", Proc. Internat. Conf. on Acoustics, Speech, and Signal Processing, San Diego, CA, March 19-21, 1984, pp. 41B.3.1-4. [26] W. Martin and P. Flandrin, "Analysis of non-stationary processes: Short-time periodograms versus a pseudoWigner estimator", in: H. Schussler, editor, EUSIPCO-83, North-Holland, Amsterdam, 1983, pp. 455-458. [27] W. Martin and P. Flandrin, "Detection of changes of signal structure by using the Wigner-Ville spectrum", Signal Processing, Vol. 8, No. 2, April 1985, pp. 215-233. [28] W. Martin and P. Flandrin, "Wigner-Ville spectral analysis of nonstationary processes", IEEE Trans. Acoust., Speech, Signal Process., Vol. 33, No. 6, December 1985, pp. 1461-1470. [29] T. Pavlidis, Structural Pattern Recognition, Springer, New York, 1977, Chap. 4, pp. 68-70. [30] T. Pavlidis, Algorithms for Graphics and Image Processing, Computer Science Press, Rockville, MD, 1982, Chap. 6, pp. 113-116. [31] S. Peleg, J. Naor, R. Harley and D. Avnir, Multiple Resolution Texture Analysis and Classification, Tech. Rept. TR1306, Dept. of Computer Science, Univ. of Maryland, College Park, MD, 1983. [32] A. Pentland, "Fractal-based description of natural scenes", IEEE Trans. Pattern Analysis & Machine Intelligence, Vol. PAMI-6, No. 6, November 1984, pp. 661-674. [33] D. Pollen and S. Ronner, "Visual cortical neurons as localized spatial frequency filters", IEEE Trans. Syst. Man & Cybernet., Vol. SMC-13, No. 5, September/October 1983, pp. 907-915. [34] L. Van Gool, P. Dewaele and A. Oosterlinck, "Texture analysis anno 1983", Computer Vision, Graphics & Image Process., Vol. 29, 1985, pp. 336-357. [35] J. Ville, "Theorie et applications de la notion de signal analytiqne", Cables et Trans., Vol. 2, No. 1, 1948, pp. 61-74. [36] A. Watson, H. Barlow and J. Robson, "What does the eye see best?" Nature, Vol. 302, 1983, pp. 419-422. [37] H. Wechsler, "Texture analysis--A survey", Signal Processing, Vol. 2, No. 3, July 1980, pp. 271-282. [38] J.S. Weszka, C.R. Dyer and A. Rosenfeld, "A comparison study of texture measures for terrain classification", IEEE Trans. Syst. Man & Cybernet., Vol. SMC-6, No. 4, 1976, pp. 269-286. [39] E. Wigner, "On the quantum correction for thermodynamic equilibrium", Phys. Rev., Vol. 40, June 1932, pp. 749-759. [40] H.R. Wilson and J.R. Bergen, "A four mechanism model for threshold spatial vision, Vision Research, Vol. 19, 1979, pp. 19-32. [41] A.P. Witkin, "Recovering surface shape and orientation from texture", Artificial Intelligence (Special Volume on Computer Vision), Vol. 17, Nos. 1-3, August 1981, pp. 17-45.