The ability to automatically identify areas of homogeneous texture present within a greyscale image is an important feature of image processing algorithms. This article describes the R package LS2Wstat which employs a recent wavelet-based test of stationarity for locally stationary random fields to assess such spatial homogeneity. By embedding this test within a quadtree image segmentation procedure we are also able to identify texture regions within an image.
This paper provides an introduction to the LS2Wstat package (Taylor and Nunes 2014), developed to implement recent statistical methodology for the analysis of (greyscale) textured images. Texture analysis is a branch of image processing concerned with studying the variation in an image surface; this variation describes the physical properties of an object of interest. The key applications in this field, namely discrimination, classification and segmentation, are often dependent on assumptions relating to the second-order structure (variance properties) of an image. In particular many techniques commonly assume that images possess the property of spatial stationarity (Gonzalez and Woods 2001). However, for images arising in practice this assumption is often not realistic, i.e. typically the second-order structure of an image varies across location. It is thus important to test this assumption of stationarity before performing further image analysis. See Figure 1 for examples of textured images. For a comprehensive introduction to texture analysis, see (Bishop and Nasrabadi 2006) or (Petrou and Sevilla 2006).
Recently, (Taylor et al. {in press}) proposed a test of spatial stationarity founded on the locally stationary two-dimensional wavelet (LS2W) modelling approach of (Eckley et al. 2010). The LS2W modelling approach provides a location-based decomposition of the spectral structure of an image. The \(Bootstat_{LS2W}\) test proposed by (Taylor et al. {in press}) uses a statistic based on an estimate of the process variance within a hypothesis testing framework, employing bootstrap resampling under the null hypothesis assumption of stationarity to assess its significance.
Given a test of spatial stationarity for random fields, it is natural to consider how this might be usefully applied within a problem such as texture segmentation. The ability to determine non-stationarity and the presence of localised textured regions within images is important in a broad range of scientific and industrial applications, including product evaluation or quality control purposes. Possible areas of use for the methods described in this article include identifying uneven wear in fabrics (Taylor et al. {in press}; Chan and Pang 2000) and defect detection on the surface of industrial components (Wiltschi et al. 2000; Bradley and Wong 2001) or natural products (Funck et al. 2003; Pölzleitner 2003). For a review of texture segmentation, see (Pal and Pal 1993).
Readily available implementations for stationarity assessment have, up until now, been restricted to the time series setting; examples of such software in this context include the R add-on packages urca (Pfaff 2008; Pfaff and Stigler 2013), CADFtest (Lupi 2009) and locits (Nason 2013b,a).
Below we describe the package LS2Wstat which implements the
spatial test of stationarity proposed by (Taylor et al. {in press}). The package has
been developed in R and makes use of several functions within the
LS2W package (Eckley and Nason 2011, 2013). The article
is structured as follows. We begin by describing details of simulation
of LS2W and other processes. An overview of the \(Bootstat_{LS2W}\) test
of stationarity is then given, focussing in particular on the function
TOS2D
. We then illustrate the application of the test on both
simulated and real texture images. Finally, the article concludes by
describing how the algorithm might be embedded within a quadtree image
splitting procedure to identify regions of texture homogeneity within a
multi-textured image.
Before describing the implementation of the work proposed in
(Taylor et al. {in press}), we first explain how to generate candidate textured
images using the LS2Wstat package. Several different spatially
stationary and non-stationary random fields can be generated with the
simTexture
function. See the package help file for full details of the
processes available.
To demonstrate the LS2Wstat implementation, throughout this
article we consider a realisation of a white noise process with a
subregion of random Normal deviates in its centre with a standard
deviation of 1.6. This simulated texture type is called NS5
, and is
one of several textures which can be simulated from the package. In
particular, we consider an image of dimension \(512\times 512\) with a
central region that is a quarter of the image, i.e. a dimension size of
\(128\times 128\). This can be obtained as follows:
> library("LS2Wstat")
> set.seed(1)
> X <- simTexture(n = 512, K = 1, imtype = "NS5", sd = 1.6, prop = 0.25)[[1]]
> image(plotmtx(X), col = grey(255:0/256))
The simTexture
function returns a list of length K
with each list
entry being a matrix representing an image of dimension n
\(\times\) n
with the chosen spectral structure. In this case, since K = 1
, a list
of length 1 is returned. The simulated image X
is shown in
Figure 2. Note in particular that visually, one can
discern that the image consists of two subtly different texture types.
Firstly, the centre of the image has one particular form of second order
structure. The second texture structure can be seen in the remainder of
the image. Throughout the rest of this article we shall apply the
approach of (Taylor et al. {in press}) to this image.
We now briefly introduce the LS2W random field model of (Eckley et al. 2010) together with some associated notation, before describing the implementation of the test of stationarity proposed in (Taylor et al. {in press}). For an accessible introduction to wavelets, please see (Prasad and Iyengar 1997), (Vidakovic 1999) or (Nason 2008).
The LS2W process model is defined by
\[\label{eq:ls2wproc} X_{\mathbf{r}} = \sum_{l} \sum_{j=1}^{\infty}\sum_{\mathbf{u}} w^{l}_{j,\mathbf{u}}\psi^{l}_{j,\mathbf{u}}(\mathbf{r})\xi^{l}_{j,\mathbf{u}} \, , \tag{1}\]
for directions \(l=h, v \mbox{~or~} d\) and spatial locations \(\mathbf{r}\), where \(\{\xi^l_{j, \mathbf{u}}\}\) is a zero-mean random orthonormal increment sequence; \(\{\psi^l_{j,\mathbf{u}}\}\) is a set of discrete nondecimated wavelets and \(\{w^l_{j,\mathbf{u}} \}\) is a collection of amplitudes, constrained to vary slowly over locations of an image (Eckley et al. 2010). In the above definition, we assume the image is of dyadic dimension, i.e. we have \(\mathbf{r}=(r,s)\) with \(r, s\in\{1,\dots,2^J\}\) and where \(J\) is the coarsest observed scale.
(Eckley et al. 2010) also define the local wavelet spectrum (LWS) associated with an LS2W process. The LWS for a given location \(\mathbf{z}=\left(\frac{r}{2^J},\frac{s}{2^J}\right)\in (0,1)^2\), at scale \(j\) in direction \(l\) is \(S^l_j(\mathbf{z})\approx w^l_j(\mathbf{u}/\mathbf{R})^2\). The LWS provides a decomposition of the process variance at (rescaled) locations \(\mathbf{z}\), directions \(l\), and wavelet scales \(j\). In practice the LWS is usually unknown and so needs to be estimated (see Eckley et al. 2010 for details). Spectral estimation using the LS2W model is implemented in R in the add-on package LS2W (Eckley and Nason 2013). The LS2Wstat routines described below thus have a dependence on some functions from the LS2W package.
Next we turn to describe the implementation of a test of stationarity
within the LS2Wstat package. We focus on describing the
\(Bootstat_{LS2W}\) approach implemented in the LS2Wstat package,
referring the interested reader to (Taylor et al. {in press}) for details of other
tests which might be employed. Throughout this section let us assume
that we have some image \(X_{\mathbf{r}}\) (as in
Figure 2), whose second-order structure we wish to test
for spatial stationarity. We assume that \(X\) is an LS2W process with
associated unknown spectrum, \(S_{j}^{\ell}\) for \(j=1,\ldots,J\) and
\(\ell=v\), \(h\) or \(d\). Since the model in (1) assumes the
process has zero mean, if necessary the image can be detrended. This can
be done in R, for example, by using the core stats package
function medpolish
, which implements Tukey’s median polish technique
(Tukey 1977).
Under the null hypothesis of stationarity, the wavelet spectrum will be constant across location for each scale and direction. Motivated by this fact (Taylor et al. {in press}) formulate a hypothesis test for the stationarity of the image \(X_{\mathbf{r}}\) with \[\begin{aligned} H_0 : & \ S_{j}^{\ell}(\mathbf{z}) \mbox{ is constant across $\mathbf{z}$ for all $j$ and $\ell$}, \\ H_A : & \ S_{j}^{\ell}(\mathbf{z}) \mbox{ is not constant across $\mathbf{z}$ for some $j$ or $\ell$}. \end{aligned}\] Hence, a test statistic for the hypothesis should measure how much the wavelet spectrum for an observed image differs from constancy. (Taylor et al. {in press}) suggest using the average scale-direction spectral variance as a test statistic to measure the degree of non-stationary within an image, where the variance is taken over pixel locations, that is:
\[\label{eq:tos} T\left\{ \hat{S}_j^{\ell}(\mathbf{z})\right\} = \frac{1}{3J}\sum_{\ell} \sum_{j=1}^{J} \mbox{var}_{\mathbf{u}}\left( \hat{S}_{j,\mathbf{u}}^{\ell}\right). \tag{2}\]
In practice this statistic is computed based on an (unbiased) estimate
of the local wavelet spectrum, produced by the LS2W function
cddews
(see the documentation in LS2W for details on optional
arguments to this function). For the (square) image X
, the test
statistic is calculated using the function avespecvar
as follows:
> TSvalue <- avespecvar(cddews(X, smooth = FALSE))
> TSvalue
1] 0.2044345 [
Since the spectrum characterises the second-order structure of the observed random field (and hence its stationarity properties), (Taylor et al. {in press}) suggest determining the p-value of the hypothesis test by performing parametric bootstrapping. This corresponds to sampling LS2W processes assuming stationarity under the null hypothesis, and comparing the observed test statistic to that of the simulated LS2W processes under stationarity. For pseudo-code of this algorithm, please see Algorithm 1.
This bootstrap algorithm is performed with the LS2Wstat function
TOS2D
. The function has arguments:
The image you want to analyse.
A binary value indicating whether the image should be detrended
before applying the bootstrap test. If set to TRUE
, the image is
detrended using Tukey’s median polish method.
The number of bootstrap simulations to carry out. This is the value \(B\) in the pseudocode given above. By default this takes the value 100.
This specifies the test statistic function to be used within the
testing procedure to measure non-stationarity. The test statistic
should be based on the local wavelet spectrum and by default is the
function avespecvar
representing the statistic (2).
A binary value indicating whether informative messages should be printed.
Any optional arguments to be passed to the LS2W function
cddews
. See the documentation for the cddews
function for more
details.
Note that TOS2D
uses the LS2W process simulation function
LS2Wsim
from the LS2W R package to simulate bootstrap
realizations under the null hypothesis. The output of TOS2D
is a list
object of class "TOS2D"
, which describes pertinent quantities
associated with the bootstrap test of stationarity. The object contains
the following components:
The name of the image tested for stationarity.
A vector of length nsamples + 1
containing each of the test
statistics calculated in the bootstrap test. The first element of
the vector is the value of the test statistic calculated for the
original image itself.
The statistic used in the test.
The bootstrap p-value associated with the test.
In particular, the object returns the measure of spectral constancy in
the entry statistic
, together with the p-value associated with the
stationarity test (in the p.value
component).
An example of the function call is
> Xbstest <- TOS2D(X, nsamples = 100)
Note that the p-value returned within the "TOS2D"
object is computed
using the utility function getpval
, which returns the parametric
bootstrap p-value for the test from the bootstrap test statistics
provided by counting those test statistic values less than
\(T^{\mathit{obs}}\) (see Davison et al. 1999 for more details). In other
words, the p.value
component is obtained by the following call:
> pval <- getpval(Xbstest$samples)
0.204
Observed bootstrap is -value is 0.00990099 p
This p-value can then be used to assess the stationarity of a textured image region.
Information on the "TOS2D"
class object can be obtained using the
print
or summary
S3 methods for this class. For example, using the
summary method, one would obtain
> summary(Xbstest)
2D bootstrap test of stationarity
object of class TOS2D----------------------------------
summary=======
: X
data: 0.204
Observed test statistic-value: 0.01 bootstrap p
Alternatively, the print
method for the "TOS2D"
class prints more
information about the Xbstest
object. Note that the function
internally calls the summary
method for "TOS2D"
objects:
2D bootstrap test of stationarity
object of class TOS2D----------------------------------
summary=======
: X
data: 0.204
Observed test statistic-value: 0.01
bootstrap p
: 100
Number of bootstrap realizations: avespecvar spectral statistic used
To demonstrate the test of stationarity further, we now provide some
other textured image examples. Firstly, we consider a Haar wavelet
random field with a diagonal texture, an example of a LS2W process as
described in (Eckley et al. 2010). The realisation of the process (shown
in Figure 3) is simulated using the simTexture
function
with the command:
> Haarimage <- simTexture(512, K = 1, imtype = "S5")[[1]]
The test of stationarity of (Taylor et al. {in press}) performed on the image
Haarimage
with the function TOS2D
reveals that the image is
spatially stationary as expected, with a high p-value associated to the
test.
> Haarimtest <- TOS2D(Haarimage, smooth = FALSE, nsamples = 100)
> summary(Haarimtest)
2D bootstrap test of stationarity
object of class TOS2D----------------------------------
summary=======
: Haarimage
data: 0.631
Observed test statistic-value: 0.673
bootstrap p
: 100
Number of bootstrap realizations: avespecvar spectral statistic used
As another example of a textured image, we construct an image montage
using two of the textures shown in Figure 1 from the
package LS2W. The montage, montage1
, is shown in
Figure 4.
Note that since this image may not have zero mean as assumed by the LS2W
model (1), we detrend the montage first using the
medpolish
function in the stats package.
> data(textures)
> montage1 <- cbind(A[1:512, 1:256], B[, 1:256])
> montage1zm <- medpolish(montage1)$residuals
The TOS2D
test indicates that the texture montage is non-stationary:
> montage1zmtest <- TOS2D(montage1zm, smooth = FALSE, nsamples = 100)
> summary(montage1zmtest)
2D bootstrap test of stationarity
object of class TOS2D----------------------------------
summary=======
: montage1zm
data: 0
Observed test statistic-value: 0.01
bootstrap p
: 100
Number of bootstrap realizations: avespecvar spectral statistic used
In this section we describe embedding a test of stationarity into a quadtree algorithm to identify regions of spatial homogeneity within a textured image. This segmentation approach is similar in spirit to, e.g., (Spann and Wilson 1985) or (Pal and Pal 1987) which use homogeneity measures within a quadtree structure. We first give details of the quadtree implementation, and subsequently describe functions to aid illustration of quadtree decompositions.
Each image is further split if deemed non-stationary, which is
determined by a test of stationarity such as TOS2D
. After the first
subdivison of an image, each sub-image is of size \(n/2 \times n/2\). The
sizes of the regions halve in size at each progressive division but
increase in number. The R function in LS2Wstat which creates the
quadtree structure described in Algorithm 2 is imageQT
. The
function has inputs:
The image to be decomposed with the quadtree algorithm.
A function for assessing regions of spatial homogeneity, for example
TOS2D
.
The testing size of sub-images below which we should not apply the
function test
.
The significance level of the \(Bootstat_{LS2W}\) test, with which to assess spatial stationarity of textured regions.
Any other optional arguments to test
.
As an illustration of using the imageQT
function, consider the code
below to decompose the (non-stationary) input image X
. We use the
function TOS2D
to assess the regions of spatial homogeneity although
the imageQT
function allows other functions to be supplied.
> QTdecX <- imageQT(X, test = TOS2D, nsamples = 100)
The output of the imageQT
function is a list object of class
"imageQT"
with components:
The index representation of the non-stationary images in the quadtree decomposition.
The results of the stationarity testing (from the test
argument)
during the quadtree decomposition. The results giving FALSE
correspond to those non-stationary sub-images contained in the
indl
component and the results giving TRUE
correspond to the
stationary sub-images, i.e. those contained in the indS
component.
The index representation of the stationary images in the quadtree decomposition.
This particular way of splitting an image has a convenient indexing representation to identify the position of subregions within an image. If a (sub)image is subdivided into quadrants, we assign it a base 4 label as follows: 0 – top-left quadrant; 1 – bottom-left quadrant; 2 – top-right quadrant; 3 – bottom-right quadrant. By continuing in this manner, we can assign an index to each tested subregion, with the number of digits in the index indicating how many times its parent images have been subdivided from the “root” of the tree (the original image). This indexing is illustrated for the quadtree decomposition given in the example in Figure 5.
Examining the quadtree decomposition of the image X
using the print
S3 method for the "imageQT"
class, we have
> print(QTdecX)
2D quadtree decomposition
object of class imageQT-------------------------
summary=======
: X
data-stationary sub-images:
Indices of non"0" "1" "2" "3" "03" "12" "21" "30"
-images:
Indices of stationary sub"00" "01" "02" "10" "11" "13" "20" "22" "23" "31" "32" "33" "030" "031" "032" "033"
"120" "121" "122" "123" "210" "211" "212" "213" "300" "301" "302" "303"
: 64 minimum testing region
The resl
component gives the results of the test of stationarity for
all sub-images tested during the quadtree procedure, reporting FALSE
for the non-stationary sub-images and TRUE
for the stationary ones:
> QTdecX$resl
1] FALSE
[2]]
[[1] FALSE FALSE FALSE FALSE
[3]]
[[1] TRUE TRUE TRUE FALSE TRUE TRUE FALSE TRUE TRUE FALSE TRUE TRUE
[13] FALSE TRUE TRUE TRUE
[4]]
[[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[16] TRUE [
By performing the quadtree algorithm given in Algorithm 2, it is possible to decompose images into regions indicating regional stationarity. Note that if a texture discrimination procedure is used to classify the output from the stationarity quadtree algorithm, the image segmentation method can be seen as a split and merge technique.
Suppose we have performed the quadtree decomposition. The
LS2Wstat package includes an S3 plot
method for "imageQT"
objects to plot the output for the "imageQT"
class and optionally a
classification of those textured regions. If the classification output
is plotted (class = TRUE
), each textured region is uniquely coloured
according to its texture group. The function has arguments:
A quadtree decomposition object, such as output from imageQT
.
Vector of class labels associated to the subimages produced by the quadtree decomposition.
A value for any unclassified values in a quadtree decomposition.
A Boolean value indicating whether to plot the classification of the quadtree subimages.
A Boolean value indicating whether to plot the quadtree decomposition.
We now illustrate the use of this function with the example given in
Figure 2. Suppose the textured regions identified by the
quadtree algorithm in the QTdecX
object have been classified according
to some texture discrimination procedure. For the purposes of this
example, we suppose that the 28 regions of stationarity in QTdecX
(see
Figure 5) have been classified as coming from two groups
according to the labels
> texclass <- c(rep(1, times = 15), rep(c(2, 1, 1), times = 4), 1)
> texclass
1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 2 1 1 2 1 1 1 [
Using the output from the quadtree technique (QTdecX
) and the texture
classification vector texclass
, we can use the quadtree plotting
function for "imageQT"
objects as follows:
> plot(QTdecX, texclass, class = TRUE, QT = TRUE)
> plot(QTdecX, texclass, class = TRUE, QT = FALSE)
The quadtree decomposition from this example is shown in Figure [fig:CT1]a; the same decomposition is shown together with the texture classification in Figure 6b.
We also consider an image montage using the textures from the package
LS2W. The montage Y
is shown in Figure 7. Prior
to performing the quadtree decomposition, we detrend the image.
> data(textures)
> Y <- cbind(A[1:512, 1:256], rbind(B[1:256, 1:256], C[1:256, 1:256]))
> Yzm <- medpolish(Y)$residuals
Similarly to above, we can now perform a quadtree decomposition of the
image Y
:
> QTdecYzm <- imageQT(Yzm, test = TOS2D, nsamples = 100)
> print(QTdecYzm)
2D quadtree decomposition
object of class imageQT-------------------------
summary=======
: Yzm
data-stationary sub-images:
Indices of non
-images:
Indices of stationary sub"0" "1" "2" "3"
: 64 minimum testing region
The function imageQT
initially assesses that the image is indeed
non-stationary, and then proceeds to analyse sub-images of the montage.
The algorithm stops the quadtree decomposition after the first
decomposition level, since it judges all quadrants of the image to be
stationary, described by the indices "0", "1", "2", and "3".
In this article we have described the LS2Wstat package, which implements some recent methodology for image stationarity testing (Taylor et al. {in press}). Our algorithm is most useful as a test of homogeneity in textures which are visually difficult to assess. We have also extended its potential use by embedding it within a quadtree implementation, allowing assessment of potentially multi-textured images. The implementation is demonstrated using simulated and real textures throughout the paper.
We thank Aimée Gott for suggestions on an early version of the package. We would also like to thank Heather Turner, two anonymous referees and the Editor for helpful comments which have resulted in an improved manuscript and package.
LS2Wstat, LS2W, urca, CADFtest, locits
Econometrics, Finance, TimeSeries
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Nunes, et al., "A Multiscale Test of Spatial Stationarity for Textured Images in R", The R Journal, 2014
BibTeX citation
@article{RJ-2014-002, author = {Nunes, Matthew A. and Taylor, Sarah L. and Eckley, Idris A.}, title = {A Multiscale Test of Spatial Stationarity for Textured Images in R}, journal = {The R Journal}, year = {2014}, note = {https://rjournal.github.io/}, volume = {6}, issue = {1}, issn = {2073-4859}, pages = {20-30} }