Automated forest structure mapping from high resolution imagery based on directional semivariogram estimates

ELSEVIER

Automated Forest Structure Mapping from High Resolution Imagery Based on Directional Semivariogram Estimates Benott A. St-Onge’

and Frangois Cavauas f

A

that allows forest stands

new segmentation

identi$cation imagery

J

approach

on high spatial resolution

is presented.

Texture

information

rived by measuring

the range

monochrome

values in three

image

using a moving window. then

used to predict,

structure

parameters

growing

three regression

of the semivariogram different

on a per-pixel through

m) optical

was first de-

The semivariogram

oped for crown diameter, sure. A region

(~1

stand density, algorithm

directions ranges were

basis, three

regression

of

stand

equations develand crown clo-

was applied

to these

estimate images to identify the limits of

the forest stands. Calibration was made using artij&al

of the prediction

images created

by a geometri-

on the tree age, health,

tree height,

and other such

biodiversity.

in quality, sometimes in traditional

on arti$cial im-

as tone,

estimates within

associated

good

segmentation

and

per

Some errors were generated ing window

method

sometimes

stands because

stand

structure

estimates.

due to the fact that the mov-

overlapped

of the presence

forest vegetation or human

to real im-

sensor yielded

two difjerent

of areas covered

made structures.

forest by non-

The issue of

the moving window size and means to increase the precision of the

method

are

discussed.

OElsevier

Science

up-to-date

It is, however,

remote

sensing

large scale air photographs texture,

preters

(Leckie,

used to create and

pattern,

digital

(visual interpretation

context)

quirements

du Quebec

UniversitC

de Mont&al,

B Mont&al, Montr&l,

Address correspondence to Benoit A. St-Onge, Dept. de GBographie, Univ. du Quebec B Mont&al, C. P. 8888, Succursale Centre-Ville, Montreal, QuBbec, Canada H3C lP8. Receiwd 10 May 1996; revised 6 Novenlher 1996. REMOTE SENS. ENVIRON. 61:82-95 (1997) OElsevier Science Inc., 1997 655 Avenue of the Americas, New York, NY 10010

photointer-

1992). On the other hand, new methods forest maps from orbital remote

image

mapping,

of

due to the problems of retiring

processing

techniques

great capacity in reconnaissance forest

More-

based on visual criteria such

but seldom

sensing

have

shown

surveys, or small scale

meet

the information

re-

of large scale mapping. This situation will un-

doubtedly change with the advent of very high resolution orbital

sensors

promises

such as Earthwatch’s

to deliver I-m panchromatic

nm) as soon as 1997

(Tripp,

1995).

Indeed,

which

QuickBird,

imagery (450-900 However.

radically different

this new

image pro-

a pixel will no more repre-

sent the averaged signal reflected Universite

that large scale

difficult to find experts

with the replacement

cessing approaches.

“Department of Geography, Mont&al, QuCbec, Canada t Department of Geography, QuBbec, Canada

for

recent

digital information

recognized

increasingly

kind of data will demand

Inc., 1997

need

during

showing high error levels.

over, it is becoming

identijed

of the proposed

increased

official forest inventory maps are often outdated and vary

structure

with the METS-11 airborne

has

on such things as wood volume, habitat type, or potential

can be adequately

application

crown closure,

for each stand. The

information

years. Users now demand

ages and that average forest

stand are close to the actual values. Pre-

density,

characteristics

forest

boundaries

liminary

stands

forest inventory maps typically show the out-

species,

process.

ages acquired

that forest

Large-scale

line of forest stands and provide information

cal-optical

each delineated

It was found

equations

INTRODUCTION

by a number

of trees,

as is the case for Landsat Thematic Mapper-k images, but will only characterize other objects

a small part of a tree crown or of’

found in the forest. It follows that classical

image processing

techniques,

such as pixel-by-pixel

maxi-

mum likelihood classification, will most likely yield poor results. Thus, computer stand mapping will necessarily have to rely on processing

algorithms that analyze groups 0034.4257/9;/5 17.00 PI1 S0034-4257(96)00242-8

Automated Forest Structure Mapping

of contiguous pixels. This is even more critical for the structural aspects of the forest, such as density (number of trees per hectare) and crown closure (area covered by vertically projected tree crowns over total area, also called percent cover). These characteristics represent spatial averages of given quantities over a certain area and cannot be estimated on a per-pixel basis. Therefore, the problems associated with pixel-by-pixel classification techniques that have already been observed in the past for relatively high resolution imagery such as Landsat Thematic Mapper or SPOT HRV sensors (Moreno and Gandia, 1988; Cushnie, 1987) will certainly be much greater with the next generation of space images. Our general goals are to develop a forest mapping method able to 1) produce forest maps precise enough to derive wood volume estimates, 2) characterize the structure of forest ecosystems, 3) yield a map of forest structure assemblages useful for habitat and biodiversity evaluation, and 4) create a method for facilitating extensive inventories by generating automatic characterization of training sites for the classification of “conventional” space images (Lundsat-TM, SPOT-HRV, etc.). The specific objectives of the present study were to develop a set of algorithms designed to carry large-scale inventory of forest structure parameters, namely, tree size, tree density and percent cover, from very high resolution imagery, that is, images having a pixel size of approximately 1 m, and to assess their usefulness with simulated and real aerial images. To achieve automated forest inventory, the problem was divided in two subobjectives: the development of a stand delineation algorithm and the elaboration of a tree size, density, and percent cover estimation method. The objectives were pursued for specific forest conditions-monospecific softwood stands, flat terrain-and for trees having a shape, size, and density typical of the boreal forest. To achieve both automated stand delineation and forest characterization by image analysis, additional features, such as texture, must be taken into account along with the usual tone and color. Image texture has indeed been one of the most recommended criteria to characterize forest from aerial photographs (Spurr, 1960; Ktichler, 1967; Howard, 1970; UBC, 1971; Avery and Burkhart, 1983). The field of digital image processing has yielded numerous successful texture analysis algorithms, but, unfortunately, these were generally applied to images in which textures were visually very distinct, such as the ones found in the overused Brodatz photographic album (Brodatz, 1966). To our knowledge, it was never formally demonstrated that a texture measure based on second-order statistics, such as autocorrelation, the semivariogram, the cooccurrence matrix, etc., could differentiate two marginally distinct texture fields, a situation very often found on images of forested territory. The textural difference of two neighboring forest stands being indeed, in many situations, quite subtle. Nevertheless,

83

computer texture analysis has shown interesting results when applied to forest mapping (Atkinson and Danson, 1988; Cohen et al., 1990; Nel et al., 1994). We have recently confirmed the great potential of textural analysis for the specific problem of stand structure mapping (StOnge and Cavayas, 1995) and delineating forest stands (St-Onge, 1994) that was foreseeable from previous studies (Yuan and Vlcek, 1989; Ryherd and Woodcock, 1990; Hay and Nieman, 1994; Kushwaha et al., 1994). The information contents of the texture of a wide variety of optical and radar images has been analyzed in previous research efforts (Laporte, 1983; Gougeon and Wong, 1987; Taberner et al., 1987; Edwards et al., 1988; Franklin and McDermid, 1994). Various methods were employed: simple local variance measures (Logan and Strahler, 1979), second-order spatial statistics such as the cooccurrence matrix and the semivariogram (Cohen et al., 1990; St-Onge, 1994), and neural network approach (Dreyer, 1993). Results in general show that texture measured at different scales reveals different types of information. For example, it is the mosaic of forest stands that tends to influence the texture of 30-m resolution Landsat Thematic Mapper images, thus giving an indication of the stand-to-stand variation, while it is the light and dark patterns created by individual tree crowns that are responsible for the texture appearing on high resolution aerial imagery (pixel size of 2 m or less), and that allow withinstand evaluation of tree size or density. It also appears that, in general, optical images contain more information than radar images and that second-order statistics offer more information extraction potential than simpler statistics. After a texture measure is applied, further processing is necessary to achieve stand delineation, a task for which different classical methods, such as edge detection, pixel-by-pixel classification, region growing, etc., can theoretically be employed. Because of natural forest continuums, the stand boundaries are often diffuse and ill-defined, even for an experienced photointerpreter. Thus, classical edge-detection algorithms applied on texture feature images derived from second-order statistic measurements are not a workable solution. The values of texture feature images could also be classified using conventional means such as maximum likelihood classification and the boundaries automatically extracted from the classified image. However, this method oversimplifies reality at the beginning of the process by imposing the discretization of the continuous values of tree size and density. Another problem comes from the fact that a pixel-by-pixel approach is used, where a “region-byregion” method would be better suited. Over the years, complex image segmentation problems such as stand delineation were partially overcome by using a combination of texture and region growing (Chen and Pavlidis, 1978; Enrich and Lai, 1979; PietikBnen and Rosenfeld, 1981; Laws, 1984; Chou et al., 1994; Spann and Wilson, 1985;

84

St-Onge and Caoayas

Beveridge et al., 1989). To our knowledge, the combination of texture and region growing methods have however rarely been applied to forest mapping problems. One notable exception is the work of Woodcock and Harward (1992) in which a centroid-linkage region growing algorithm (Tilton, 1983) was used to delineate forest stands from Thematic Mapper original bands and a derived texture channel. This approach was later tested using similar images (Shandley et al., 1996) or with multispectral SPOT and simulated 30-m resolution imagery (Ryherd and Woodcock, 1996). Region growing methods were also applied to textural channels derived from aerial or submeter resolution imagery (Gougeon and Wong, 1987; St-Onge et al., 1991). The meyhod reported here combines a texture estimation using multidirectional calculations of the semivariogram and a simple region growing algorithm to delineate and characterize forest stands. The two types of very high resolution imagery (pixel size of 0.5 m) that were used are described in the next section, followed by a detailed presentation of the algorithmns. An assessment of the precision of both the delineation and the stand structure estimations is then presented. The article is coneluded with a discussion of the method and its extension.

DATA The proposed segmentation approach was first developed on computer-simulated images of the forest and later tested on real digital aerial images. This choice was motivated by the impracticability of acquiring a real imageforest dataset where the structure parameters of interest varied one at the time and because of the high costs of in situ forest measurements. The 0.5-m monochrome images were generated from three-dimensional computer models of trees having a Lambertian reflectance as described in St-Onge and Cavayas (1995). Illumination was simulated by parallel rays and shadows were calculated from the 3-D model (tree shadows on the ground and tree shadows on trees). The trees were modeled by symmetrical functions, giving rise to shapes similar to softwood crown envelopes placed on sticks representing the boles. The data set was composed of series of 128 by 128 pixel individual images. Inside a series, either one of tree size, density, or percent cover, three parameters, was varied in a given range while all other variables were fixed. This geometrical-optical simulation procedure produces (Fig. la) realistic images having textural characteristics very close to those of similar natural images (StOnge and Cavayas, 1995). In order to assess the validity of the method in reallife situations, we used very low altitude (2400 ft) aerial imagery from the MEIS-II sensor. The MEIS-II pushbroom instrument was developed by the Canada Center for Remote Sensing (Till et al., 1987). With its ‘I-mrad this sensor can generate images of angular resolution,

pixel size as low as 30 cm. It has eight narrow spectral bands ranging from the blue to the infrared. In the summer of 1988, two flight lines over the Petawawa National Forest in Eastern Ontario, Canada produced 36 cm rcsol&ion images. A subscene of 1024 by 1024 pixels in the red band was taken from one of the MEIS-II scene. This spectral band behaves similarly to the opaque geometrical tree models used to generate forest images because of the strong absorption of red light by the chlorophyll activity. The subscene is composed mainly of softwood stands, but the presence of hardwood is not negligible. Three kinds of stands compose the subscene: natural stands, regular plantations, and irregular plantations.

ESTIMATING FOREST STRUCTURE USING SEMIVARIOGRAMS The Semivariogram

The texture measure used as the basis for the quent segmentation procedure is founded on the riogram. The semivariogram, hereafter called the gram” measures the average difference between separated by lag k, given in pixels or directly in meters:

subseserniva“variopixels grom1d

(1) It has the following characteristics and advantages: ;r) The measure relies on a limited set of assumptions in comparison of those of the autocorrelation function OI Fourier transform (Ramstein and Raffy, 1989); b) it is not tied to the choice of initial parameters such as relative position of the pixels and the matrix analysis method (contrast, entropy) needed to control cooccurrence algorithms; and, lastly, c) it yields texture coarseness values directly in terrain distance units. Moreover, the range of the variogram is rather iusensitive to variation of contrast from one image to the next such as those produced 1, sensor gain, emulsion speed or atmospheric attenuation (St-Onge and Cavayas, 1995). The variogram has beerr used in the past to estimate the optimal window sizes for texture or spectral measures in natural forests or plantations (Atkinson and Danson, 1988; Franklin and McDermid, 1994), to assess tree damage (Franklin et al., 1992) and to estimate stand structure parameters (Jupp et al., 1988; Woodcock et al., 1988). The typical variogram exhibits a rise that gradually slows to form a straight horizontal line (sill) due to the fact that increasing the distance between two pixels augments the chance of observing a great difference between their values. In many cases, the lag increase does not result in any more difference because the spatial autocorrelation between image values normally drops bcyond a certain distance (Fig. 1b). The height of the sill is normally proportional to the global image \.ariancc while the lag at which the sill is attained, called the

85

Automated Forest Strwtw~ Mapping

b) 1600

-

Theorical variogram at 45 degrees: range=335 cm, sill=1396 (DN)

1200 g

1000

z g ‘C z

800

Experimental degrees

variogram

at 45

Theoretical variogram at 90 degrees: range=316 cm, sill=1372 (DN)

600

Experimental degrees

400 -

200

200

400

600

800

at 90

Theoretical variogram at 135 degrees: rang-517 cm, sill=1245 (DN)

0 0

variogram

1000

Experimental degrees

distance (cm)

variogram

at 135

Figure 1. a) The three directions in which is measured the range of the variogram; b) experimental and adjusted theoretical variogram (spherical model) in three directions for the image shown in a).

“range,”

is generally

a very good

indicator

of texture

coarseness.

gated dark patches

on the ground. The variogram

value in the direction

The variogram can be calculated

range

will systematically

(At-

be higher than in the other directions.

and

sure increases,

the tree shadows will in most cases fall

Raffy, 1989) but gives a better representation when calculated in more than one direction (Cohen et al., 1990;

on other trees

and elongated

Woodcock

et al., 1988) because

this study, each directional

anisotropic

(Davis, 1981). The texture measure presented

kinson and Danson,

here was composed

1988;

Curran,

along transects

of the elongation

1988;

Ramstein

texture is in many cases

of three values given in ground me-

ters: the range of the variogram in the direction parallel to the sun rays at the moment of image acquisition, upa,.,

disappear,

yielding

first calculating experimental

the variogram,

variogram,”

retical variogram

patches

an isotropic

When

crown clo-

will progressively

texture

coarseness.

In

range value is obtained

by

usually referred

in one given direction.

model is then adjusted

to as “the A theo-

to the experi-

the range in the perpendicular direction, aBeT, and the range in the intermediate direction (45” from both parallel and perpendicular directions), urnd. These individually

mental variogram by least square fitting (Fig. lb). We used the widely accepted spherical model (see, e.g., Curran, 1988), a two-parameter function characterized by the range a and the sill c. The details of variogram esti-

characterize texture coarseness in different directions and, together, estimate texture anisotropy. This image characteristic is strongly influenced by crown closure (StOnge and Cavayas, 1995). Indeed, we can observe that, in open forest canopies, the shadows tend to form elon-

mation are presented in St-Onge and Cavayas (1995). We emphasize the fact that first-order statistics of the image are completely left out by the texture measure; in other words, the average reflectance of a forest stand or the global variance of the digital values corresponding to

86

St-Onge and Cavayas

the image of a given stand do not participate

in the texture Predicting

significantly

scores.

Stand Structure Parameters

Predicting stand structure parameters was done by using a set of three functions relating the directional range parameters to the structure parameters. These functions were derived from a multiple regression analysis having apor, ap,>,.,and a,& as independent variables and crown diameter CD, stand density SD, and crown closure CC, successively as dependent variables, The tapor-aper parameter was introduced to better capture coarseness anisotropy. The multiple regression analysis was run on an image dataset composed of 700 artificial images, as described in the next section, where CD ranged from 1 m to 8 m, SD from 50 trees/ha to 3000 trees/ha, and CC from 2% to 98%. Stepwise regression yielded the following equations: CD=O.62-0.50

a,,,+1.53

a,,,,,+0.33 lap,,r-a,,c,rl

(R”=0.93)

In SD=8.04-0.20 -1.45

a,,+0.02 ln a,&

In CC= -0.45-0.43

All R” were significant Mapping

Inlar,~,-al,p,l

(R’=O.86)

a,,,+0.27

+0.15 In &i
a,,,+0.15

(2)

a,,,,+0.09

(3) lnlal,O,-apt.~/

(R”=0.89)

dimensional image variogram y(ij,h,& is calculated for pixel (i,i), lag h, and direction B (determined in geographic azimuths: 0” = north). The polar coordinates constituted by the lag and the direction were converted to plane coordinates given by functions p(h,6) and v(hJ). expressing in pixels the x and y translation, respectively. For example, a lag of 2 pixels and a direction of 90” (horizontal) yields p(h,B)=2 and v(h,@=O, a lag of 42 ant1 a direction of 135”, gives ,u(h,@)=l and v(h,Q)=l, etc. The local variogram y(ij,,n(h,fl),u(h,tl)) for a moving window of size 2K-- 1 by 2L - 1 is given by-

where YLis the number of pixels pairs considered in the calculation. The local estimations of the directional variogram ranges given by the above equation are used to estimate the three stand structure parameters by applying Eqs. (2)-(4). The result was a set of three new images indicating respectively the estimate of local crown diameter in meters, local tree density in trees per hectare, and local crown closure in percentage (Fig. 2).

(4)

at F=O.OOl.

Local Stand Structure Parameters

Mapping the stand structure parameters was accomplished in two steps: 1) local estimation of al,a,, a,,,, and anridby calculating the variogram in a moving window and 2) translation of the variogram ranges into structure parameter using Eqs. (2)-(4) for each image pixel. Local estimation of the directional variogram ranges required that we define a window size for the local operator giving a good trade-off b etween an accurate stand structure estimation (obtainable from a large window) and good spatial precision of the resulting maps (for which a small window is preferable). It was also important to ascertain that the locally estimated range values reflected the contrast variations of one given stand at a time. Since few clues existed at the time, different window sizes were tested in order to assess the magnitude of the effect of window size on the quality of the results. This assessment was conveyed for the identification of the boundaries between the stands as well as the exactitude of stand structure estimation. Results associated with these different sizes are presented in a later subsection. Each range estimate was given by iteratively adjusting the spherical model to the locally calculated variograms, which makes for long computation runs. In order to minimize the length of calculations, the moving window was displaced by steps greater than 1 pixel. The local two-

IMAGE

SEGMENTATION

BY REGION

GROWING Region Growing Approach

Image segmentation refers to the division of the image plane in smaller units. Region growing consists in growing regions from starting points (sometimes homogeneous nuclei formed by pixel acretion based on a similarity criterion, by splitting the image recursively until all regions are homogeneous or by a combination of the two techniques). This segmentation approach is sometimes used concurrently with edge detection (see, e.g., Pavlidis and Liow, 1990). According to Haralick and Shapiro (1985), a segmented image should have the following properties: 1) Regions of a segmented image should be uniform and homogeneous with respect to the feature used; 2) adjacent regions should have significantly different values; 3) regions should be relatively simple and without many holes; and 4) boundaries of each region should be simple, not ragged, and be spatially accurate. In depth reviews of the most common methods to achieve these conditions can be found in Ballard and Brown (1982) and Haralick and Shapiro (1985). A recent survey of the application of image segmentation algorithms in remote sensing appears in Woodcock and Harward (1992). Many of the region growing algorithms presented in the literature dealt with only one spectral or textural channel at the time. Our situation is complicated by the

Automated Forest Structure Mapping

87

Figure 2. a) Example of a four stand image assemblage; b) resulting crown diameter value image for a): c) re. sulting density value image for a); d) resulting

fact that three

textural

channels

ously to analyze the forested considered:

simultane-

space. Two solutions were

1) Perform

the region

growing

separately

and create

the final map by in-

the channels tersecting

are needed

crown closure value image for a).

the three resulting

polygon layers; 2) perform

a single region growing using an integrated ture similarity

measure.

on each of

We suspected

forest struc-

that the former

approach might result in a map where a multiple of very small polygons would arise from the intersection individual layers and therefore proach. was performed

as fol-

lows: The first pixel (top left) of the overlaid forest structure images is considered, then compared

its multidimensional

value is

to that of each of its 8-connected

diate neighbors,

any neighboring

imme-

pixel having a good de-

gree of similarity is added to the first region and all rejected

pixels individually

form

starting

points

for new

regions. When new pixels are added to a growing region, their

values

are used to update

the average

structure

value of all the pixels forming that region. A region stops growing when no bordering when

all bordering

threshold.

pixels can be added, that is,

pixel can not

A compound

meet

similarity criterion

the

similarity

A was estab-

lished as follows:

+w&c-CC(i,j)l,

sure for a given region,

closure

would have a weight

vidually. This choice was motivated first two variables wood volume

10

that of the two first ones taken indiby the fact that the

have the strongest

influence

of a stand while crown closure

great part from the combination

on the results in

of crown diameter

and

stand density. Generalization

of the Regions

As will be showed in the next section, procedure

the segmentation

always gives rise to a certain

number

of very

small regions in addition to the main regions. Those gions appear mostly in the transition

zones between

refor-

est stand and are much too small to appear on large scale (typically 1:25,000) ple procedure

forest maps. We thus devised a sim-

to eliminate

place them by appropriate polygon

resulting

very small regions and to revalues. First, the area of each growing algorithm

is

computed.

Every region below a given area threshold

from the region

is

eliminated

by putting

all the pixels

inside

Lastly, all zero pixels of the image receive

it to zero. the value of

the closest nonzero pixels, having for effect that the main regions grow into the gaps left by removing This method

the smaller

is very similar to the “symmetric

fill operator”

of Pong et al. (1984).

APPLYING

THE TEXTURAL

(6)

CD, SD, and CC are, respectively,

values of crown diameter,

but that crown

segments.

A(i~)=oc~,lCD-CD(i~)l+os~ISD-SD(izj)l

where

A(i,j)

times less important

of the

opted for the second ap-

The region growing algorithm

and tree density would have equal weights in calculating

the average

stand density, and crown cloand CC(i,j)

CD(i,j),SD(ij)

respectively,

the structure

neighboring

pixels having image coordinates

parameter

are,

values of individual

i and j, and

ocD, usD, and occ are the relative weights. The decision rule is given by

APPROACH

SEGMENTATION

TO SIMULATED

IMAGES

Procedure The

segmentation

procedure

stand images generated ric-optical Each

technique

described

image was composed

by different

was first tested

by computer

on forest

using the geomet-

in the

second

section.

of four stands characterized

average tree height, average crown diameter,

accept pixel in a given segment

if

A(i@hreshold,

density,

reject pixel of a given segment

if

A(i,j)>threshold.

stands were of softwood type, having a shape that resem-

and

percent

cover.

bles that of most spruce The relative weights of forest structure

parameters

have

to be decided prior to the experimentation. In the present case, it has been decided that both crown diameter

diameter

trees.

Tree

composing

height

often used to characterize

in forest

stand (Halley

the

and crown

varied among each stand according

bull distribution, tribution

All trees

to a Wei-

tree size dis-

and Schreuder,

1977).

88

St-Onge

and Cavayas

Figure

3. Images in Figures 2b, 2c, and 2d respectively, smoothed by a 3 by 3 pixel mean filter.

The resolution of the 3-D stand model, as well as that of the resulting artificial images was of 0.5 m, the stand images having a size of 128 by 128 pixels and the composite forest images thus having a size of 256 by 256 pixels. Four different assemblages were generated. The first one comprised four very different stands and was devised to veritj, that the segmentation procedure would work in favorable conditions (situation a, Fig. 5). In the second assemblage, the four stands had the same density (800 trees/ha); only tree size (and inevitably percent cover) varied (situation b, Fig. 5). The reverse was tested with the third forest image were the average crown diameter was fixed at 4.1 m, but the density was increased by small steps of 200 trees/ha, from 600 to 1200 (situation c, Fig. 5). The last assemblage was devised to verify that size differences could also be detected when density and percent cover are quite small (situation d, Fig. !5). Figures 24 present an example of the application of our approach on artificial images. The local directional variogram was calculated using three different window sizes: 63 by 63 pixels, 47 by 47 pixels, and 31 by 31 pixels (Fig. 5). In each case, the horizontal and vertical displacement between two variogram calculations was of 21 were conpixels. The estimations of n,,,,,., (I,,~,, ant1 CL,,,,~~ verted to stand structure values using Eqs. (2)-(4), and the resulting images were filtered using an operator that averaged 3 by 3 local estimations (Fig. 3). The region growing algorithm was then applied to the three-parameter images followed by the generalization procedure (Fig. 4). Human intervention was needed at two stages: 1) to determine the threshold for A(ij) of Eq. (6) and 2) to determine the area threshold under which small regions are eliminated. The first threshold influences the size of the resulting regions; a low threshold value gives rise to

a great number of very small regions while a high threshold value can sometimes result in the whole image becoming a single region. Some trials were necessq to empirically find a good region growing cutoff value. The degree of generalization of forest maps and its influence on the assessment of image processing methods for forest mapping is a complicated problem at best, for which no good solution exists. We have tried to find a threshold that was appropriate for the scale at which we apprehend the forest in this study. The cutoff value for the smoothing out of small regioL was somewhat easier to fiud, for it was clear from 3 visual analysis which regions corresponded to actual forest stands and which ones werra IXsidual (see, e.g., Fig. 4a). The results associated with this threshold can in fact be known in advance or (San I,(, fixed to agree with a given cartographic standard (r1.g.. no region can be smaller than x m”), which is not tile case for the first threshold. Results

The results of the segmentation are s1ww11in Figures 5a-5d and the corresponding stand parameter estimations for the 63 by 63 pixel window operators are pnhsented in Tables 14. Each of the four stands in all th(J forest images are ref&red to as 1, 2, %3,and 4, from left to right a1~1from top to bottom. The best segmentation results are obtained when a large window (63 by 63 pixels) is used to estimate the directional variogram. In general, the appropriateness of the limits given by the procedure decrease with local window size. The segmentatiol~ is better when the contrast between stands is at its ul;lx-imum, namely, in the first image (Fig. 5a). The rest&s for second image (Fig. 5b) show that, at least for high resolution simulated images, it is possible to find approximate limits between stands differing only in tree size us-

Automated

Forest

Structure

89

Mapping

Figure 5. Original four stand images, segmentation using a 63 by 63 pixel window, segmentation using a 47 by 47 pixel window, and segmentation using a 31 by 31 pixel window for situations a, b, c, and d, described in text.

ing texture

and region

growing.

some limits of the segmentation some encouraging

Figures

5c and 5d show

procedures

results the same time.

by 63 pixel window to measure were found between

while giving Using the 63

the variogram,

no limits

stands 1 and 2 in image 3. When

we also take into account

the results of figure Sd, it be-

comes

clear

that

precise

opened

forest

stands is difficult when only small differ-

ences in density are present.

segmentation

of

This observation

rated by the fact that the segmentation stands seems easier, as indicated

relatively

by the precise demarca-

stands 3 and 4 of Figure 5c. In low density

stands, trees tend to cluster if their distribution pens, the local density of trees is more subject

to varia-

tions than in higher density stands. This has a nonnegligible impact on the variance of the texture measure values for a given stand and, thus, on the success

of segmen-

tation. On the other hand, it is encouraging when the visual difference tion

Delineated

procedure

is able

by the Segmentation

to see that even

is very subtle, the segmentato

discriminate

Algorithm

Actual

Stand A Predicted

Stand B Actual

Stand B Predicted

Stand C Acfual

Stand C Predicted

5.1 100 20

4.7 123 15

4.1 500 54

4.5 712 66

2.0 500 14

2.1 434 15

Stund A

is close

to being random, which is the case here. When this hap-

is corrobo-

of high density

Table 1. Forest Structure Estimation for Stands Automatically with a 63 by 63 Pixel Window (Corresponding to Fig. 5a)

ml (m) SD (trees/ha) PC (%)

tion between

between

on Situation Stnnd D Actual -~ :3.0 2000 92

two

1 Stand D Predicted 3.0 1550 61

90

St-Onge

and Cavayas

Table 2. Forest Structure Estimation for Stands Automatically Delineated by the Segmentation Algorithm on Situation 2 with a 63 by 63 Pixel Window (Corresponding to Fig. 5b)

CD (m) SD (trees/ha) PC (%)

Stand A Actual

Stand A Predicted

Stand B Actual

Stand B Predicted

Stand C Actual

Stand C Predicted

Stand V

2.0 800 21

1.7 709 23

3.0 800 53

3.0 1153 57

4.1 800 73

3.8 920 63

5.1 800 90

Actual

Stand D Predicted ._.___~ -1.5 798 74

Table 3. Forest Structure Estimation for Stands Automatically Delineated by the Segmentation Algorithm on Situation 3 with a 63 by 63 Pixel Window (Corresponding to Fig. 5c) ..____ Stand A G B

Stand A S B

CD (m) SD (trees/ha) PC (%)

Stand C Actual

(Aoerage)

@vzr%y

Predicted

4.1

4.4

700

4.1

865 73

71

Stand C Predicted ___4.0

1000 87

stands differing only in density, when the density of trees is high. These results show that texture measures based on second-order statistics are capable in some cases of measuring small texture differences. The results of forest stand structure estimation are presented in Tables l-4 only for the best segmentation results (63 by 63 pixel window). Because the segmentation procedure did not recognize any frontier between stands 1 and 2 of Figure 5c, the true values for these stands were averaged in order to be compared with the predicted values which correspond to the space covered by both stands. The same was done with stands 3 and 4 of Figure 5d because of the poor results of the segmentation. Despite the fact that, for all other stands, the true and predicted outlines do not perfectly agree, estimates of forest stands are quite good. The predicted values for mean tree diameter are in general more accurate than for density and percent cover; The average error for tree diameter is indeed only 0.3 m. The amplitude of this error represents approximately 8% of the mean true diameter (3.6 m) calculated across all stands. The error levels are higher for density (130 trees/ha, or 18% of the mean density) and percent cover (10, or 19% of mean percent cover). This is probably due to the fact that, in general, variables related to the number of trees are more difficult to predict than variables related to tree size (StOnge, 1994). As indicated by formulas (2)-(4), the two

Stand D Actual -4.1

1024 73

I200 95

Stunrl II Predicted 3.6

latter structure variable estimates (SD and CC) tend to be less stable than those of crown diameter due to their nonlinear relationships with the semivariogram parameters. Nevertheless, the average performance is quite good, and the error is generally smaller than the class width found on most 1:25,000 maps. Canadian forest maps typically have a 2-m class interval in height (corresponding approximately, assuming the standard hardwood tree shape used in our study, to a 0.5-m crown diameter interval), 200 trees/ha class intervals in density, and 20% class intervals in percent cover.

APPLYING APPROACH

THE TEXTURAL SEGMENTATION TO MEIS-II IMAGES

The segmentation and structure parameter estimation were conducted on the subscenes using a 63 by 63 window moved by steps of 15 pixels to achieve greater spatial resolution. The direction of incident sun rays was identified from the shadows, and the directional variogram was calculated accordingly. The sun elevation was derived from the date and time of dav at which the image was acquired and the latitude. It was found to be approximately 40”, approaching the value of 45” used to generate the synthetic images on which are based the regression models (ey. 3-5). Figure 6 shows the resulting structure parameter images. The first segmentation trial

Table 4. Forest Structure Estimation for Stands Automatically Delineated hy the Segmentation Algorithm on Situation 4 with a 63 by 63 Pixel Window (Corresponding to Fig. 5d)

Cd (m) SD (trees/ha) PC (%)

..

I IJO 67

Stand A Predicted

Stand B Actual

Stand B Predicted

Stand C S I> (Amrage) Actual

Stand C G 1)

Stand A Actual 2.0 300 8

2.4 279 12

3.0 300 21

3.7 148 14

4.6 300 4:3

4.; :?I*54 AH

( Aoemgc)

Predicted

Automated Forest Structure Mapping

91

was conducted using a large threshold value, resulting in general partitioning of the image. Three regions were outlined (Fig. 7). The first region is a meadow, region 2 contains dense stands populated by small trees, and region 3 is composed of larger trees occurring in open stands. Table 5 shows the estimates for each region. The values indicated for region 1 are, of course, inexact but serve to show that the algorithm is not appropriate to discriminate forested from nonforested land, although the combination of the variogram and region growing appears to be a very good tool to delineate the forested land. Although the lack of ground-truth does not allow us to quantify the success of the algorithm, we can infer from a visual examination of Figure 7 that the segmentation is quite appropriate. The boundary between regions 2 and 3 follows very accurately the contour of an open plantation of large spruce trees. It also correctly follows the limit of a very dense plantation of small spruce trees. The size of the window is responsible for the offset between the boundary and the limit of the forest along the road. This problem could be alleviated by the use. of a digital map of roads and other man-made features in the context of GIS and remote sensing integration. A second segmentation was conducted using a narrower threshold for the similarity criterion [Eq. (6)] on the same image (Fig. 8). It yielded a finer partitioning of the image which is sometimes quite accurate. We can see, for example, in Figure 8 that two visibly distinct small plantations (zone 11 and zone 13) were correctly delineated. On the other hand, the limit of the forest was only poorly identified, and some obvious plantations, like the large Norway spruce plantation of zone 7 and part of zone 8, were correctly delineated only in their upper parts. It seems that the narrower threshold gives place to some degree of arbitrary segmentation that is difficult to relate to natural forest boundaries. DISCUSSION Forest mapping involves stand characterization and delineation, which can be viewed as two different but simultaneous operations when done manually. For the sake of simplicity, the task of translating these to a working algorithm was done by transforming the parallelism in a serial process: 1) local estimation of stand structure parameter, 2) delineation and generalization with a precision appropriate for a given scale, and 3) averaging of local structure values in order to generate a unique value for each stand and each structure parameter. This approach implies that an arbitrary area must be set for estimating stand density and tree size before the stand limits are known. We have shown that the greater the area, the more exact are the estimations. On the other hand, the risk of generating an estimation for an area encompassing parts of two potentially very different stands increases with the size of the area in which the initial variograms are cal-

Figure 6. a) Crown diameter image; b) density image, and c) crown closure image for the MEIS-II image shown in Figure 7.

92

St-Onge and Chaps

Figure tation.

culated. drops,

7. 1024 by 1024 pixel red band MEIS-II

In that case, the stability of structure yielding

a high within-stand

make segmentation method

texture

composed

This will assure a certain

of the texture primitive

of small

level of repetition

within a single unit of measure

with images

posed of small or medium-sized process

results

For this reason, the algorithm we pro-

posed works better relaxation

that will

will give rise to better

when stands have a uniform

of the variogram.

estimation

variance

quite difficult. This explains that the

we developed

primitives.

image with a wsolution

of dense

stands com-

trees, We believe that a

involving the iteration

of steps

and 3 could bring better results by progressively

1, 2,

adapting

the size, and possibly the shape, of moving window inside which the local variogram The discrepancy

observed

of the measure

cases, predictable.

The tree inodels used for image syu-

fer trees, but the MEIS-II tion of hardwood Synthetic

trees

the presence

It is known, for example, that

the overall image structure effect

For the sake ot

of gaps was not taken into ac-

count in the geometrical-optical

process presented

in this

study but should eventually be included in further devrlopments of the technique. The fact that other objects, such as roads, clearings,

to real images was, in many

of a stand aud will

on the variogram.

the presence

the natural image also induce

from the

in three). in their in-

of canopy gaps due to local tree falls will

have a certain simplicity,

one

and real forest stands also differed

ternal spatial arrangement. influence

image shows a good propor(approximately

is calculated. the results asso-

segmeu-

thesis were designed to imitate the crown shape of coni-

between

ciated with artificial images and those stemming application

of’ :36 cm with general

the segmentation

and meadows,

also appear on

some amount

of error in

process. The use of properly stnictured

digital topographical

maps stored in a geographical

infor-

Fore,st Structure Mapping

Automated

Figure tation.

8.

1024 by 1024 pixel red band MEIS-II

image with a resolution

system could help focus the variogram analysis on

mation

forest

stand and reduce

the “contamination”

by other

types of objects. method used to smooth out poly-

gons too small to be considered sible for some imprecision Indeed,

in the delineation

actual stands. The method

gaps originating

from the removal

Table 5. Predicted

we used to fill the of those

in that it allocates

one or the other neighboring This simplification

of reality could be replaced

by a reanalysis of raw structure founded

on distance

image pixels to

stands solely on a distance values by a decision

rule

but also “spectral” criteria.

of stands.

those polygons arise mostly near the boundaries

between

Three

as stands is also respon-

of 36 cm with fine segmen-

gons is very arbitrary criterion.

The generalization

93

small poly-

Values for CD, SD, and PC for the

Forest

CONCLUSION The variogram calculated

Region

Predicted SD (t/ha)

Predicted PC

1

4.2

635

68%

2

3.3

3

7.2

331 lS2

57% 45%

was used to

forest images and

regression

generated

models built from computer

images

were used to convert the variogram range values to stand structure

Predicted CD (m)

in three directions

measure the texture of high resolution

parameters.

When applied to artificial images, a

region growing algorithm conducted yielded a good delineation

on the range images

of the stands while giving rel-

atively accurate

estimates

sity and percent

cover. It was demonstrated

cases (dense

stands),

of crown diameter,

the method

stand den-

that, in some

used was able to cor-

94

St-Onge and Cavayas

rectly identify the boundary between two neighboring stands showing only minor visual and structural differences. The results are generally better when percent cover is high. When applied on real aerial images of the forest, the method yielded less precise results, in part due to the presence of other objects (roads, clearings) and a difference virtual trees (softwood only) and the real trees appearing on the MEIS-II images (hardwood and softwood). The results of combined forest characterization and delineation are nevertheless very good. The average difference between estimates and true forest value becomes negligible if one classifies the values. Our characterization and segmentation methods are based solely on monospectral values. It has the advantage of being impervious to linear modifications of contrasts such as those usually induced by atmospheric variations (St-Onge and Cavayas, 1995). Complete forest maps will, nevertheless, necessitate the use of multispectral information to facilitate the identification of forest species. Perhaps the most interesting aspect of our study is that the method we developed offers the possibility of fine tuning the relative importance of a certain number of criteria in the segmentation process. Although external factors such as terrain slope and viewing conditions could impede on the success of the algorithm, it is reasonable to conclude that, with the support of GIS data, automated detailed mapping of the forest could be attained in a foreseeable future, and we believe that the method we presented here shows the feasibility of automated forest stand mapping based on image processing. It clears the way for the elaboration of more precise and more complete methods that will exploit the potential of the next generation of l-m resolution space sensors. This new data combined with recent techniques such as the ones we described will greatly facilitate commercial forest inventories and provide a very fine stand characterization useful to habitat and biodiversity studies. The authors wish to thank the Canada Centre for Remote Sensing for the permission to use the MEIS II images. Our appreciation also goes to the Fetawawa National Forest Institute staff, particularly to Dr. Don Leckie for providing access ta useful ground data and to Dr. Frango& Gougeon for his assistance in transferring the MEIS II images.

REFERENCES Atkinson, P. M., and Damon, F. M. (1988), Spatial resolution for remote sensing of forest plantations, In Proceedings of IGARSS ‘88 Symposium, Edinburgh, Scotland, 13-16 September, ESA SP-284, IEEE 88CH2497-6, pp. 221-223. Avery, E. A., and H. E. Burkhart (1983) Forest measurements, McGraw Hill, New York. Ballard, D. H., and Brown, C. M. (1982), Computer Vision, Prentice-Hall, Englewood Cliffs, NJ.

Beveridge, J. R., Griffith, J., Kohler, R. R., Hanson, A. R., and Riseman, E. M. (1989), Segmenting images using localized histograms and region merging. Int. J. Comput. Vision 2:311-347. Brodatz, P. (1966), Textures: A Photographic Album fur Artists and Designers, Dover, Toronto. Chen, P. C., and Pavhdis, T. (1978), Segmentation by texture using a co-occurrence matrix and a split-and-merge algorithm, Technical Report 237, Princeton University Princeton, NJ. Chou, J., Weger, R. C., Ligtenberg, J. M., Kuo, K. S., Welch, R. M., and Breeden, P. (1994), Segmentation of polar scenes using multi-spectral texture measures and morphological filtering. Int. J. Remote Sens. 15:1019-1036. Cohen, W. B., Spies, T. A., and Bradshaw, G. A. (1990), Semvariograms of digital imagery for analysis of conifer canopy structure. Remote Sens. Environ. 34:167-178. Curran, P. J. (1988) The semi-variogram in remote sensing: an introduction. Remote Sens. Environ. 24:493507. Cushnie, J. L. (1987), The interactive effect of spatial resolution and degree of internal variability within land-cover types on classification accuracies. Int. J. Rem&c Sens. 8: 15-29. Davis, L. S. (1981), Polarograms: a new tool for image texture analysis. Pattern Recognition 13:219-223. Dreyer, P. (1993), Classification of land cover using optimized neural nets on SPOT data. t’hotogramm. Eng. Remote Sens. 59:617-621. Edwards, G. R., Landry, R., and Thompson, K. P. B. (1988), Texture analysis of forest regeneration sites in high resolution SAR imagery. In International Geoscience and Remote Sensing Symposium, pp. 1355-1360. Enrich, R. W., and Lai, P. F. (1979), Texture region growing based on a structural model of texture, Virginia Polytechnic Institute and State University, Department of Computer Science, Blacksburg, VA. Franklin, S. E., Bowers, W., Hudack, J., and McDemud, G. J. (1992), Estimating structural damage in balsam fir stands using semivariance. In Proceedings of the 15th Canadian Symposium on Remote Sensing, Vancouver, British Columbia, June, pp. 96-99. Franklin, S. E., and McDermid, G. J. (1994), Empirical relations between digital SPOT HRV and CASI spectral response and lodgepole pine (Pinus contorta) forest stand parameters. Int. J. Remote Sens. 14:2331-2348. Gougeon, F. A., and Wong, A. K. C. (1987), Towards a knowledge-based spectral and textural signature recognition. In Proceedivqs of the 1Ith Canadiun Symposium on Remote Sensing, Waterloo, Ontario, 22-25 June, pp. 565-578. Hafley, W. I., and Schreuder, H. T. (1977). Statistical distributions for fitting diameter and height data in even aged stands. Can. J. For. Res. 7:481489. Harahck, R. M., and Shapiro, L. G. (1985), image segmentation techniques. Comput. Vi/is. Graph. Image Process. 29:100-132. Hay, G. J ,, and Niemann, K. 0. (1994), Visualizing 3-D texture: a three-dimensional structural approach to model forest texture. Con. J. Remote Sens. 20:90-101. Howard, J. A. (1970), Aerial Photo-Ecology, Eisevier, Anisterdam. Jupp, 1). L. B., Strahler, A. H., and Woodcock, C. E. (1988),

Automated Forest Structure Mapping

Auto-correlation and regularization in digital images: I. Basic theory. IEEE Trans. Geosci. Remote Sens. 26:463473. Ktichler, A. W. (1967), Vegetation Mapping, Ronald Press, New York. Kushwaha, S. P. S., Kuntz, S., and Oesten, G. (1994), Apphcations of image texture in forest classification. Int. J. Remote Sens. 15:2273-2284. Laporte, J. M. (1983), Texture analysis on SPOT simulations. in 17th Symposium on Remote Sensing of the Environment, pp. 124331252. Laws, K. I. (1984), Goal-directed textured-image segmentation, in Computer Vision Research and Its Applications to Automated Cartography, Technical Report, Defense Advanced Research Projects Agency, Arlington, VA. Leckie, D. G. (1992), Application of airborne multispectral scanning to forest inventory mapping. In Proceedings of the International Forum on Airborne M&spectral Scanning for Forestry and Mapping (with Emphasis on MEIS), Petawawa National Forestry Institute, Information Report PI-X-113, pp. 86-93. Logan, T. L., and Strahler, A. H. (1979), Use of a standard deviation based texture channel for Landsat Classification of forest strata. In Machine Processing of Remotely Sensed Data Symposium, p. 395. Moreno, J. F., and Gandia, S. (1988), The influence of the high spatial resolution of SPOT HRV data on the separability of spectral classes. In International Geoscience and Remote Sensing Symposium, p. 597. Nel, E. M., Wessman, C. A., and Veblen, T. T. (1994), Digital and visual analysis of Thematic Mapper imagery for differentiating old growth from younger spruce-fir stands. Remote Sens. Environ. 48:291301. Pavlidis, T., and Liow, Y. T. (1990), Integrating region growing and edge detection. IEEE Trans. Pattern Anal. Mach. lntell. 12~225-233. Pietikainen, M., and Rosenfeld, A. (1981), Image segmentation by texture using pyramid node linking. IEEE Trans. Syst. Man Cybemeti. 11:822-825. Pong, T.-C., Shapiro, L. G., Watson, L. T., and Harahck, R. M. (1984), Experiments in segmentation using a facet model region grower. Comput. Vis. Graph. Image Process. 25:1-23. Ramstein, G., and Raffy, M. (1989) Analysis of the structure of radiometric remotely-sensed images. Int. 1. Remote Sens. 10:1.049-1073. Ryherd, S. L., and Woodcock, C. E. (1990), The use of texture in image segmentation for the definition of forest stand boundaries. In Proceedings of the 23rd International Symposium on Remote Sensing of Environment, Bangkok, Thailand, 18-25 April, Vol. 2, pp. 12091213. Ryherd, S., and Woodcock, C. (1996), Combining spectral and texture data in the segmentation of remotely sensed images. Photogramm. Eng. Remote Sens. 62:181-194.

95

Shandley, J., Franklin, J., and White, T. (1996) Testing the Woodcock-Harward image segmentation algorithm in an area of southern California chapparral and woodland vegetation. lnt. J. Remote Sens. 5:983-1004. Spann, M., and Wilson, R. (1985), A quad-tree approach to image segmentation which combines statistical and spatial information. Pattern Recognition 18:257-269. Spurr, S. H. (1960), Forest Inventory, Ronald Press, New York. St-Onge, B. A. (1994), L’apport de la texture des images numeriques de haute resolution a la cartographic forestibre automatisee, Ph.D. thesis, Department of Geography, University of Montreal, 251 pp. St-Onge, B. A., and Cavayas, F. (1995) Estimating forest stand structure from high resolution image7 using the directional variogram. lnt. J Remote Sens. 16:1999-2021. St-Onge, B. A., Cavayas, F., and Teillet, P. M. (1991) Etude de la signature spatiale des couverts forestiers par modelisation geometrique-optique. In Proceedings of the 5th lnternational Colloqium on Physical Measurements and Signatures in Remote Sensing Courchevel, France, 1418 January, ESA SP-319, pp. 671-674. Tabemer, M. J., Philpot, W. D., and Philipson, W. R. (1987), Spectral and spatial characterization of orchards in NewYork State using Thematic Mapper imagery, in Proceedings of the 21st Symposium on Remote Sensing of Environment, Ann Arbor, MI, 26-30 October, pp. 1061-1072. Till, S. M., Neville, R. A., Leckie, D. G., and Strome, W. M. (1987) Advanced airborne electrooptical imager. In Proceedings of the 21st International Symposium on Remote Sensing of Environment, Ann Arbor, MI, 26-30 October, pp. 4147. Tilton, J. C. (1983) Multiresolution spatially constrained clustering of remotely sensed data on the massively parallel processor, Proceedings, Tenth International Symposium on Machine Processing of Remotely Sensed Data, 12-14 June, Purdue University, W. Lafayette, Indiana, 297-304. Tripp, B. (1995), EarlyBird satellite expected to sharpen focus on commercial remote sensing industry. Earth Ohs. Mag. (10):4648. UBC (1971), Forestry Handbook for British Columbia, The Forest Club, UniversitC de Colombie-Britannique, Vancouver. Woodcock, C., and Harward, V. J. (1992), Nested-hierarchical scene models and image segmentation. lnt. J. Remote Sens. 13:31673187. Woodcock, C. E., Strahler, A. H., and Jupp, D. L. B. (1988), The use of variograms in remote sensing: l-Scene models and simulated images. Remote Sens. Environ. 25:323348. Yuan, X., and Vlcek, J. (1989), Cl assification of coniferous forest species using spectral and textural analysis of MEIS II data. In Proceedings of the ASPRS-ACSM Annual Meeting, Baltimore, MD, 2-7 April, Vol. 3, pp. 128-139.