3610

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 12, DECEMBER 2006

Voxel-Based 3-D Modeling of Individual Trees
for Estimating Leaf Area Density Using
High-Resolution Portable Scanning Lidar
Fumiki Hosoi and Kenji Omasa

Abstract—A method for accurate estimation of leaf area density(LAD) and the cumulative leaf area index (LAI) profiles of small
trees (Camellia sasanqua and Deutzia crenata) under different
conditions was demonstrated, which used precise voxel-based treemodels produced by high-resolution portable scanning lidar. In
this voxel-based canopy profiling (VCP) method, data for eachhorizontal layer of the canopy of each tree were collected fromsymmetrical azimuthal measurement points around the tree using
optimally inclined laser beams. The data were then converted intoa voxel-based three-dimensional model that reproduced the tree
precisely, including within the canopy. This precise voxel modelallowed the LAD and LAI of these trees, which have extremely
dense and nonrandomly distributed foliage, to be computed bydirect counting of the beam-contact frequency in each layer using a point-quadrat method. Corrections for leaf inclination andnonphotosynthetic tissues reduced the estimation error. A beamincident zenith angle near 57.5◦ offered a good correction for leaf
inclination without knowledge of the actual leaf inclination. Nonphotosynthetic tissues were removed by image-processing techniques. The best LAD estimations showed errors of 17% at theminimum horizontal layer thickness and of 0.7% at the maximum
thickness. The error of the best LAI estimations was also 0.7%.
Index Terms—Leaf area density (LAD), leaf area index (LAI),
portable scanning lidar, tree, voxel-based model.

I. I NTRODUCTIONTHE FOREST canopy has the important functional roles ofcycling material and energy through photosynthesis andtranspiration, maintaining forest microclimates, and providing
habitats for various species. Determining the vertical structureof the canopy is very important because the three-dimensional(3-D) composition of the canopy helps sustain these functionalroles. The vertical foliage structure is often represented by
the leaf area density (LAD) in each horizontal layer, where
LAD is defined as the total one-sided leaf area per unit layer
volume [1]. The leaf area index (LAI), which is defined as theleaf area per unit ground area, is calculated by integrating the
LAD data.
The point-quadrat method was first developed as an indirect
method of measuring LAD and LAI [2], [3]. In this method, aprobe with a sharp point is inserted into the canopy at a known

Manuscript received September 25, 2005; revised March 31, 2006.
The authors are with Graduate School of Agricultural and Life Sciences, The University of Tokyo, Tokyo 113-8657, Japan (e-mail: aomasa@
mail.ecc.u-tokyo.ac.jp).
Digital Object Identifier 10.1109/TGRS.2006.881743

inclination and azimuth angle, and the number of times thepoint contacts leaves or stems is counted. LAD and LAI canthen be estimated by calculating the contact frequency, which
is defined as the mean number of contacts per insertion. The
contact frequency is a useful quantity for estimating LAD and
LAI, because no assumptions on leaf spatial distribution, shape,or size are required [1]. However, because this method requires
many insertions of the probe into the actual canopy, it is very
laborious [4].
Another indirect method, the gap-fraction method, is widely
applied in field surveys, and it uses commercially available
tools such as a camera with a fish-eye lens or the LI-COR LAI2000 Plant Canopy Analyzer [5], [6]. This method allows automatic estimation of LAI without destruction of the plants andis less laborious. However, it requires the assumption that thefoliage distribution is random, which leads to underestimationof LAI when leaves are clumped or overlapping. Moreover, theestimate is accurate only for diffuse light filtered through the
canopy; direct sunlight introduces error.
Recently, light detection and ranging (lidar), a 3-D remotesensing technique, has been applied to forests. Lidar instruments mounted on aircraft—airborne lidar—can scan the forest
canopy of an entire region by emitting and receiving laser
pulses, and the results are not affected by varying light conditions [7]–[9]. Large-footprint airborne lidar has been used
to estimate the vertical canopy surface profile on the basisof the waveforms of returned pulses [10], [11]. Although thissystem has been proven useful on regional scale, the image
resolution and information on the inner part of the canopy
that it provides are limited. In contrast, commercially available
aircraft-mounted lidar instruments with fine spatial resolution
can provide precise profiles of canopy height [12]–[15]. This
type of lidar, however, receives basically only the first and last
returned pulses, with the result that little information on the
inner part of canopy is obtained [16].
A portable ground lidar system, however, can complement
the weak points of airborne lidar in forest applications. Its
fine spatial resolution and small beam size allow the inner
canopies of trees to be measured from the ground, making an
accurate estimation of LAD profiles possible. Portable lidar
instruments have been used previously for forest-canopy observation [17]–[19]. Radtke and Bolstad [18] and Parker et al. [19]used a nonscanning lidar (range finder) to emit vertical laser
pulses upward from the ground to estimate LAD by an optical point-quadrat method [20]. This method is a kind of

0196-2892/$20.00 © 2006 IEEEHOSOI AND OMASA: VOXEL-BASED 3-D MODELING OF INDIVIDUAL TREES FOR ESTIMATING LAD

3611

TABLE IPARAMETERS OF EACH OF THE TWO TREE SPECIMENS, C. sasanqua AND D. crenata

gap-fraction method rather than a standard point-quadrat
method, because LAD is estimated from the gap fractions
observed in each horizontal layer. The LAD is thus often
underestimated when the distribution of foliage is nonrandom[18], [21], so a calibration using measured LAI values isrequired to obtain absolute LAD values. In addition, when
the foliage density is high, vertical laser-beam penetration of thehorizontal layers is uneven because fewer laser beams reach theupper canopy. For example, Lovell et al. [16] reported obtaining
better measurements from the ground for the lower canopy
than the upper canopy with portable scanning lidar. Moreover,nonphotosynthetic tissues such as stem and branches are notdistinguishable from leaves, and it was also difficult to correct
for leaf-inclination angle by previous portable lidar methods.
Finally, measuring the horizontal differences in the foliage
is very labor intensive. Therefore, we developed a practicalmethod for accurate LAD estimation using a high-resolution
portable scanning lidar, which we refer to as the voxel-based
canopy profiling (VCP) method.
II. M ATERIALS AND M ETHODS

Fig. 1. Schematic diagram of the lidar-measurement technique for estimating
LAD and LAI using a voxel-based 3-D tree model. (a) Positions used for lidar
measurements around each tree. (b) Laser-beam inclination represented by thecentral zenith angle of the laser scan (θB ) and the zenith scan angle (∆θB ),(c) flow chart showing the steps of LAD and LAI computation, and(d) illustration of a voxel-based tree model.

A. Plant Materials and Their Characteristics

B. Measurement by Portable Scanning Lidar

We tested our method on two species of small tree: Camellia
sasanqua and Deutzia crenata (Table I). The actual LAD valuesof each horizontal layer of C. sasanqua and D. crenata were
measured by stratified clipping for layer thicknesses of 0.1 and0.05 m, respectively. For this measurement, all leaves in eachhorizontal layer of C. sasanqua and D. crenata were clippedand scanned into JPEG images using a commercially available
desktop scanner (FB636U, Canon, Inc.). The relationship between the number of pixels within the images and actual area(square meters) was determined by scanning a ruler together
with the leaves. The number of pixels of leaves within theimages was then converted into the actual leaf area of the leaves
using the relationship. The actual LAD in each horizontal layer
was obtained by dividing the actual area in each horizontal layer
by the horizontal layer thickness and the area of the projectionof the sample tree’s canopy onto the ground. These values
were then integrated with respect to height to obtain LAI. The
LAI and LAD of the two specimens were similar or somewhat
higher than those of forest stands. Therefore, the levels of LAI
and LAD seemed to be suitable for confirming the validity andthe robustness of our method.

A fine-resolution portable scanning lidar (TDS-130L 3-D
laser scanner, Pulstec Industrial Co., Ltd) was used to measureeach tree from several positions [Fig. 1(a)]. The range accuracy
was about 2 mm, and the scan resolution was about 1 mm. A
rotating mount run by a built-in stepping motor and a galvano
mirror within the lidar head facilitated horizontal and verticalscanning by the instrument.
Different laser-beam inclinations and lidar positions were
tested to obtain better information on the LAD profiles of eachhorizontal layer of the trees. Although the laser inclination,
which is defined by both the central zenith angle of the laser
scan θB and the zenith scan angle ∆θB [Fig. 1(b)], greatly
influences the accuracy of LAI and LAD estimation, at the
beginning of the experiment, the optimal laser inclination was
unknown. Therefore, the trees were scanned by using zenithscan angles from 8◦ to 14◦ at central zenith angles of 37.5◦ ,59.8◦ , 77.3◦ , and 180◦ for C. sasanqua and 16.2◦ , 55.7◦ , and180◦ for D. crenata (a central zenith angle of 180◦ means thatthe laser beam was emitted from directly above the top of thetree) from four positions [Fig. 1(a) and (b)]. The horizontal scanangle was set at about 5◦ to cover the width of the materials.

3612IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 12, DECEMBER 2006

Three-dimensional point cloud data for each central zenith
angle and each position were recorded. The tree element positions in the 3-D point cloud data were represented as points
in orthogonal coordinate systems. The measurements were firstcarried out for the intact, leaved trees, and then they wererepeated for only nonphotosynthetic tissues, after removal of
the leaves from the same trees.
C. Generation of the Voxel-Based 3-D Tree Model
The voxel-based 3-D tree model was generated in several
steps [Fig. 1(c), upper panel].
1) Registration Process: The complete data set for eachcentral zenith angle was composed of four-point cloud data, one
obtained from each of the four azimuthal positions. These four
data with their individual coordinate systems were registeredinto a single-point cloud data set with a common coordinate
system by using the iterative closest point algorithm [22].
2) Voxelization: A voxel-based tree model for each central
zenith angle was produced from the registered point cloud data
set. A voxel is defined as a volume element in a 3-D array. All
points within the registered data set were converted into voxel
coordinates by the following equations:X − Xmin(1)
i = Int∆iY − Yminj = Int(2)∆jZ − Zmin(3)
k = Int∆kwhere (i, j, k) are the voxel coordinates in the voxel array, Int
is a function to round off at one decimal place to the nearest
integer, (X, Y, Z) represent the point coordinates of the registered lidar data, (Xmin , Ymin , Zmin ) are the minimum values of(X, Y, Z), and (∆i, ∆j, ∆k) represent the voxel element size.
In this experiment, the voxel element size, which depends on
the range and scan resolution of the lidar, was 1 × 1 × 1 mm;
thus, the 3-D array was composed of 700 × 700 × 1600 voxels
for C. sasanqua and 600 × 600 × 600 voxels for D. crenata[see Fig. 1(d)]. Voxels corresponding to voxel coordinates
converted for all points within the registered data set wereassigned attribute 1. A voxel with attribute 1 represents a voxel,in which the laser beams were intercepted. Other attributes
for the remaining voxels, in which the laser beams were notintercepted, are assigned in next section.
D. LAD and LAI Computational Models
LAD and LAI were computed for four cases [Fig. 1(c), lower
panel], with or without corrections for leaf inclination andnonphotosynthetic tissues (stem and branches), of the voxelbased tree model.
1) Removal of Nonphotosynthetic Tissues: For each species,
two voxel-based tree models, one for the intact, leaved tree
and the other for only nonphotosynthetic tissues, were obtained
as described in Section II-C above, “Generation of the Voxel-

Based 3-D Tree Model.” The voxel-based tree model that excluded nonphotosynthetic tissues was obtained by subtracting
the nonphotosynthetic tissue model from the leaved tree model.
2) Computation of the Contact Frequency of the Laser
Beams in Each Horizontal Layer: For computation of LAD
and LAI, a plant region is set in a voxel array. The plant region
is defined as the region above an area covered by a projectionof the canopy on the horizontal plane [see Fig. 3(a)]. The
area covered by the projection of the canopy is produced in avoxel array by projecting all voxels with attribute 1 onto thehorizontal plane at k = 0. Voxels above the area are regarded
as voxels within the plant region, and they are used for LAD
and LAI computation.
LAD between height h and h + ∆H above the ground
LAD(h, ∆H) is given by
LAD(h, ∆H) =

1
∆H

mh+∆H



l(k)

(4)

k=mh

where l(k) is the LAI of the kth horizontal layer of the voxel
array within a plant region, ∆H is the horizontal layer thickness, and mh and mh+∆H are the voxel coordinates on thevertical axis equivalent to height h and h + ∆H in orthogonal
coordinates (h = ∆k × mh ). The LAI of the kth horizontal
layer l(k) is the product of the contact frequency N (k) of laser
beams in the kth layer and the coefficient α(θ), which corrects
for leaf inclination at laser incident zenith angle θl(k) = α(θ)N (k)= α(θ) ·nI (k)nI (k) + nP (k)(5)

where nI (k) is the number of laser beams intercepted by thekth layer, nP (k) is the number of laser beams passed throughthe kth layer, and nI (k) + nP (k) is the total number of incident
laser beams that reach the kth layer. α(θ) is expressed in termsof G(θ), which is the mean projection of a unit leaf area
on a plane perpendicular to the direction of the laser beam[1], [4]–[6]α(θ) =

cos θ
.
G(θ)

(6)

G(θ) is determined with the assumption that leaves are positioned symmetrically with respect to the azimuth, as follows(see Fig. 2):
1
G(θ) =
2π

2π π/2g(θL ) |cos(nB , nL )| dθL dϕL
0

0

π/2
=

g(θL )S(θ, θL )dθL
0

(7)

HOSOI AND OMASA: VOXEL-BASED 3-D MODELING OF INDIVIDUAL TREES FOR ESTIMATING LAD

Fig. 2. Schematic diagram of the projection of a unit leaf area on a plane
perpendicular to the direction of the laser beam. θ is the laser-beam incident
zenith angle, and θL is the leaf-inclination angle (the zenith angle of the normal
to the leaf surface), and ϕ and ϕL are the azimuth angles of the laser beamand the normal to the leaf surface, respectively. 
nB and nL are unit vectors
corresponding to the direction of the laser beam and the direction of the normal
nL )| is the absolute value of theto the leaf surface, respectively. | cos(nB , 
cosine of the angle between 
nB and nL , which is the projection of a unit leaf
area on a plane perpendicular to the direction of the laser beam 
nB : (A) θ ≤π/2 − θL and (B) θ > π/2 − θL .

where
1
S(θ, θL ) =
2π
1
=
2π

=

2π
|cos(nB , nL )| dϕL02π|cos θ cos θL +sin θ sin θL cos(ϕ − ϕL )| dϕL0cos θ cos θL ,

cos θ cos θL 1 + 2(tanπx−x)

for θ ≤ π2 − θL

x = cos−1 (cot θ cot θL ).

(9)

θL is the leaf-inclination angle (the zenith angle of the normal
to the leaf surface), and ϕ and ϕL are the azimuth angles ofthe laser beam and the normal to the leaf surface, respectively.
g(θL ) is the distribution function of the leaf-inclination angle;
it is independent of the azimuth angle of the normal to theleaf surface, assuming azimuth symmetry. | cos(nB , nL )| is theabsolute value of the cosine of the angle between two unit
vectors corresponding to the direction of the laser beam andthe direction of the normal to the leaf surface. This value is the
projection of a unit leaf area on a plane perpendicular to thedirection of the laser beam nB . S(θ, θL ) is the average of| cos(nB , nL )| with respect to the azimuth angle of the normal
to the leaf surface [2]. To use the actual measured distributionof leaf-inclination angles, (7) can also be expressed as follows:
G(θ) =Tq
Fig. 3. Schematic diagram of the laser beam tracing process in the voxel
array. (a) Laser beam tracing for determining voxel attributes (1, intercepted;
2, passed; and 3, not reached by any laser beam traces). The plant region isdefined as the region above an area covered by a projection of the canopyon the horizontal plane (between the broken lines). The bottom canopy layer
is the lowest horizontal layer that includes voxels corresponding to leaves.
(b) Example of attribute assignment to voxels within a horizontal layer.

then each class consists of a 5◦ interval. For example, q = 1,
q = 9, and q = 16 include the angles from 0◦ to 4◦ , 40◦ to 44◦ ,and 75◦ to 79◦ , respectively. g(q) is the distribution of the leafinclination-angle class q, which is a ratio of the leaf area belonging to class q to total leaf area; θL (q) is the midpoint angleof class q, which is the leaf-inclination angle used to represent
class q. From (4) and (5), the LAD of each horizontal layer isexpressed in terms of the product of the contact frequency sum
and the correction coefficient α(θ), on the assumption that α(θ)
is independent of the height h, for θ > π2 − θL
(8)g(q)S (θ, θL (q))

(10)

q=1

where q is the leaf-inclination-angle class and Tq is the total
number of leaf-inclination-angle classes. Thus, if there are
18 leaf-inclination-angle classes from 0◦ to 90◦ (Tq = 18),3613

LAD(h, ∆H) = α(θ) ·

1
∆H

mh+∆H



k=mh

nI (k).
nI (k) + nP (k)(11)

Voxels that intercepted the laser beams have already been assigned attribute 1 in the voxelization process (see Section II-C).
Therefore, nI (k) is obtained by counting the number of voxels
with attribute 1 in the kth layer of the voxel-based tree model.
To obtain nP (k), laser beams are traced in the tree model(see Fig. 3). The procedure starts from the bottom canopy
layer within a plant region, which is the lowest horizontal
layer, in which laser beams were intercepted by leaves (thevertical axis k is represented as mb at the bottom canopy
layer). First, all laser beams emitted from all lidar positions
are traced in the voxel array up to k = mb in accordancewith the actual lidar settings of the laser-beam angles. Next,
voxels that do not have attribute 1 in the bottom canopy
layer are considered. If one of the voxels is intersected by atleast one laser-beam trace, the voxel is assigned attribute 2.
A voxel with attribute 2 therefore represents a voxel throughwhich a laser beam (or laser beams) passes. If a voxel that does
not have attribute 1 in the bottom canopy layer is not intersectedby any laser-beam traces, the voxel is assigned attribute 3. A
voxel with attribute 3 represents a voxel that no laser beams cantouch. By these procedures, all voxels that do not have attribute
1 in the bottom canopy layer are assigned attribute 2 or 3. Third,

3614

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 12, DECEMBER 2006

traces of laser beams that intersect voxels with attribute 1 in thebottom canopy layer are not followed above that layer, because
beams that have been intercepted in the bottom layer cannot
reach any layers higher than the bottom canopy layer. Fourth,
the remaining traces, that is, those beams not intercepted in
the bottom canopy layer, are followed up to k = mb + 1 todetermine the voxel attributes of the layer at k = mb + 1. Inthat layer, the procedures used to determine the voxel attributesof the bottom canopy layer are repeated using the traces still
being followed. Similarly, traces that are not intercepted in thelayer at k = mb + 1 are followed up to k = mb + 2, and theprocedures used in the bottom canopy layer to determine voxel
attributes are also repeated in the layer at k = mb + 2. The
same procedures are repeated sequentially, layer by layer, up
to the layer that corresponds to the tree top until finally theattributes of all voxels up to canopy top have been determined.
Thus, nP (k) is obtained by counting the number of voxels withattribute 2 in the kth layer.
The cumulative LAI from height h to the height of the tree
top HT can also be used to characterize the vertical profile of thefoliage [18]. The cumulative LAI CLAI(h) is given by the sumof contact frequency for k = mh to k = mHT in (12), insteadof h to HT in height
CLAI(h) =

m
HT

k=mh

= α(θ) ·k=mhnI (k)nI (k) + nP (k)(12)

when h is 0 (ground level), CLAI(h) represents the LAI of theentire tree.
3) Correction for Leaf Inclination: Leaf inclination is corrected for by determining α(θ) as in (6). The distribution ofleaf inclination g(q) is required to determine α(θ) using (6),(8)–(10). g(q) can be obtained by measuring the actual leafinclination angle of each leaf directly, but this is a laborious
task. However, when the laser incident zenith angle θ is 57.5◦ ,α(θ) can be considered to be almost independent of leaf inclination: α(θ) ∼= 1.1 [1], [3]. In this case, leaf inclination is wellcorrected for without any actual leaf-inclination data by using
the following equation:
LAD(h, ∆H) ∼= 1.1 ×

1
∆H

mh+∆H



k=mh

nI (k).
nI (k) + nP (k)(13)

From (12), CLAI(h) for the particular case of θ = 57.5◦ isgiven by
CLAI(h) ∼= 1.1 ×

m
HT

k=mh

Error(case 1) = |MPE(case 1) − MPE(case 3)|+ MAPE(case 3)

nI (k).
nI (k) + nP (k)(14)

E. Analysis of Measured Results
1) Determination of the Optimal Laser Inclination: Fromthe voxel-based tree models excluding the nonphotosynthetic
tissues, LAD profiles for different laser inclinations were com-

(15)

Error(case 2) = |MPE(case 2) − MPE(case 4)|+ MAPE(case 4)

l(k)
m
HT


puted by (11), with α(θ) calculated by (6) and (8)–(10) fromthe leaf-inclination data obtained by actual measurement. The
optimal conditions of laser inclination for the two specimens
were determined by comparing the mean absolute percent
errors (MAPE) between the calculated LAD profile and thatobtained by stratified clipping.
2) Error Analysis of LAD and LAI Computations: The LAD
profiles computed using the optimal laser inclination werecompared among the four correction processes [cases 1 to 4,
Fig. 1(c)]. Cumulative LAI profiles were also estimated andcompared among the correction processes. The LAD andcumulative LAI profiles for cases 1 and 3, in which leaf
inclination was disregarded, were computed by (11) and (12)with α(θ) = 1. Those for cases 2 and 4 were computed withthe appropriate measurement-based values of α(θ). For optimallaser inclinations close to 57.5◦ , LAD profiles were estimatedby (13) and cumulative LAI profiles by (14).
To estimate the net errors of the LAD values for eachcorrection process, the following equations were defined:(16)

Error(case 3) = MAPE(case 3)(17)

Error(case 4) = MAPE(case 4)(18)

where MPE and MAPE are the mean percent error and themean absolute percent error, respectively, for each correction
process. Nonphotosynthetic tissues caused LAD to be overestimated, and leaf inclination caused LAD to be underestimated,so these errors tend to cancel each other out if they are notevaluated independently. Therefore, in cases 1 and 2, in which
nonphotosynthetic tissues were not removed, the error caused
by nonphotosynthetic tissues was separated for the net-error
estimation by the expressions |MPE(case 1) − MPE(case 3)|
or |MPE(case 2) − MPE(case 4)|, as in (15) and (16).
Another possible cause of error are the differences in thevertical positions of the boundaries between the horizontal
layers between the voxel-based tree model and the actual measurements, the magnitude of which depends on the horizontal
layer thickness ∆H. Therefore, (15)–(18) were used to analyze
the errors for different values of ∆H.
III. R ESULTS
The voxel-based 3-D model (that for C. sasanqua is shown in
Fig. 4) succeeded in reproducing each tree three-dimensionallyat the individual leaf scale, although some noise spots arevisible on the leaf edges [Fig. 4(c)].
To determine the optimal laser inclination, the LAD profilesof the voxel-based tree models that excluded the nonphotosynthetic tissues were compared among different laser inclinations.
To calculate LAD using (11), α(θ) was first computed with
(6) and (8)–(10) from the actual measured leaf-inclination
data. These procedures correspond to case 4 [Fig. 1(c)]. TheseHOSOI AND OMASA: VOXEL-BASED 3-D MODELING OF INDIVIDUAL TREES FOR ESTIMATING LAD

3615

Fig. 6. Comparison of LAD profiles for C. sasanqua between laser-derived
estimates (white squares) and stratified-clipping values (black squares) for eachof the four correction process cases [Fig. 1(c)]. MAPE—mean absolute percent
error. ∆H = 0.1 m.
Fig. 4. Example of a voxelized C. sasanqua measured at θB = 59.8◦ . Voxels
with attribute of 1 are shown by white dots. (a) Side view, (b) top view, and(c) close-up view.

Fig. 5. Comparisons of LAD profiles for C. sasanqua between laser-derived
estimates (white squares) obtained with different laser inclinations, θB , andstratified-clipping values (black squares). MAPE—mean absolute percent error.
∆H = 0.1 m. LAD was corrected for both leaf inclination and nonphotosynthetic tissues [case 4 in Fig. 1(c)].

Fig. 7. Comparison of the cumulative LAI profiles for C. sasanqua between
laser-derived estimates (white squares) and stratified-clipping values (blacksquares) for each of the four correction process cases [Fig. 1(c)]. ∆H = 0.1 m.

laser-derived estimates were then compared with stratifiedclipping values of the LAD profiles. The results for C. sasanqua at layer thickness ∆H = 0.1 m are shown in Fig. 5.
The error (MAPE) ranged from 17.2% to 55.3%. LAD was
generally underestimated in the upper layer of the canopy forθB = 37.5◦ , 59.8◦ , and 77.3◦ and in both the upper and lower
layers for θB = 180◦ . The best central zenith angle θB for
the laser scan was 59.8◦ for C. sasanqua and 55.7◦
for D. crenata. These angles are close to the particularangle θB = 57.5◦ , described in Section II-D3 above, “Correction for Leaf Inclination,” at which α(θ) is almost independentof leaf inclination.
The laser-derived estimates produced by the procedures ofcases 1–4 for C. sasanqua with θB = 59.8◦ were then comparedwith the stratified-clipping values of LAD profiles at ∆H =0.1 m (Fig. 6). Because θB was close to 57.5◦ , (13), in whichα(θ) = 1.1, was used to estimate LAD. The case 4 results fitted
the actual LAD profiles best because both correction processes,
removal of nonphotosynthetic tissues, and the correction forleaf inclination, were applied. The error (MAPE) was 17.4%,
almost the same as that for θB = 59.8◦ with α(θ) computedby (6) and (8)–(10) based on actual leaf inclination (Fig. 5).
The case 1 results were also relatively accurate even though
no correction processes were applied, because the error caused
by nonphotosynthetic tissues was compensated for due to leaf
inclination. The error in the case 2 results reflects nonphotosynthetic tissues, and that in the case 3 results is due to leaf
inclination. LAD was generally underestimated in the upper
part of the canopy and overestimated in the lower part. LAD
was also computed by (13) for D. crenata because its optimum
θB value was also near 57.5◦ . As with C. sasanqua, the best
result was obtained for case 4.

3616

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 12, DECEMBER 2006

Fig. 8. Relationship between the horizontal layer thickness and LAD errors ineach correction process for C. sasanqua: (a) No correction for leaf inclinationand (b) correction of leaf inclination. See Fig. 1(c) for descriptions of cases.
Errors were estimated by (15)–(18); εt , εl , and εd represent the errors caused
by nonphotosynthetic tissues, leaf inclination, and the layer boundary difference between the lidar measurement and stratified-clipping, respectively; γ andγ  represent the residual error in cases 3 and 4, respectively; and α(θ) is thecorrection coefficient for leaf inclination.

Fig. 9. Relationship between horizontal layer thickness and LAD errors ineach correction process for D. crenata. (a) No correction for leaf inclinationand (b) correction for leaf inclination. See Fig. 1(c) for descriptions of cases.
Parameters are as in Fig. 8.

Fig. 7 shows cumulative LAI profiles of C. sasanqua forcases 1–4 at layer thickness ∆H = 0.1 m. These profiles show
the effects of the two correction processes more clearly, because
variations in each horizontal layer are muted, making the general trend more apparent. In case 2, in which the correction
for leaf inclination was applied but nonphotosynthetic tissues
retained, the cumulative LAI was consequently overestimatedat lower heights. In contrast, cumulative LAI in case 3 was underestimated because nonphotosynthetic tissues were removed
but the leaf-inclination correction was not applied. In case 4,
the match is excellent because both correction processes wereapplied. In case 1, the results are incidentally accurate because
the errors caused by the presence of nonphotosynthetic tissues
and by leaf inclination compensated for each other.
To clarify the independent effects of the two correction
processes, the LAD errors in cases 1–4 were estimated by(15)–(18) and then compared with the horizontal layer thickness ∆H for C. sasanqua (Fig. 8) and for D. crenata (Fig. 9). In
general, the larger the layer thickness became, the smaller the
LAD error was in all cases. At the maximum layer thickness,that is, ∆H equal to the tree height (C. sasanqua, 1.6 m;
D. crenata, 0.60 m), the tree was not divided into any horizontal layers; thus, the error caused by the difference in eachhorizontal layer boundary between the lidar measurement andstratified-clipping becomes zero.
For cases 1 and 3 of C. sasanqua [Fig. 8(a)], the error in
the layer boundary difference εd was 11% at ∆H = 0.1 m
compared with the error at ∆H = 1.6 m. Cases 2 and 4[Fig. 8(b)] were not suitable for estimating the net error ofdifferences in the layer boundaries because of the bias caused
by multiplication by the leaf-inclination correction coefficientα(θ). In case 3 [Fig. 8(a)], the error caused by leaf inclination
εl and a residual error γ remained at ∆H = 1.6 m. The error
caused by nonphotosynthetic tissues εt was also present incase 1 [Fig. 8(a)] at ∆H = 1.6 m. In contrast, the error at∆H = 1.6 m was less in case 2 [Fig. 8(b)] compared withthat of case 1 [Fig. 8(a)], owing to the leaf-inclination correction. Moreover, only a slight residual error γ  remained incase 4 [Fig. 8(b)] at ∆H = 1.6 m, because both corrections had
been applied. The error from the presence of nonphotosynthetic
tissues εt was estimated to be 10% by subtracting the case 3
error at ∆H = 1.6 m from the corresponding case 1 error. The
error due to leaf inclination εl was similarly estimated to be
9% by subtracting the case 3 error at ∆H = 1.6 m from thecorresponding case 4 error. The final residual error γ  in case 4
was 0.7%.
For D. crenata, in cases 1 and 3 [Fig. 9(a)], the error in
the layer boundary difference εd was 8% at ∆H = 0.05 m.
The error from nonphotosynthetic tissues εt was estimated to
be 3% by subtracting the case 3 error at ∆H = 0.60 m fromthe corresponding case 1 error. εt for D. crenata was smaller
than that for C. sasanqua, because the stem diameter of theformer was smaller than that of the latter (see Table I). The
leaf-inclination error εl was estimated to be 9% by subtracting
the case 3 error at ∆H = 0.60 m from the corresponding case 4
error. The final residual error γ  in case 4 was 5%.
To summarize, the best LAD estimations for C. sasanqua
were given in case 4, with an error of 17% at ∆H = 0.1 m
and of 0.7% at ∆H = 1.6 m. Similarly, for D. crenata, case 4
resulted in the best estimations of LAD with an error of 18%
at ∆H = 0.05 m and of 5% at ∆H = 0.60 m. The error in
the layer boundary difference εd might be improved by betteralignment of the boundaries used in lidar measurement withthose used for stratified clipping. The LAD error at maximum ∆H was the same as the overall LAI error for both
species. Therefore, the errors of 0.7% in C. sasanqua and of
5% in D. crenata were also the best results in the overall
LAI estimation.
IV. D ISCUSSION
A. Optimal Laser Inclination
The worst result was obtained for θB = 180◦ (Fig. 5). Although the lidar laser beams could illuminate the upper layerof the material at that inclination, for the most part, they could
not reach the inner layers because the upper layer obstructed
laser penetration. For θB = 37.5◦ , LAD in the middle andupper layers was underestimated because the laser beams were
obstructed by the lower layer. For θB = 77.3◦ and 59.8◦ , theinclination of the laser beams resulted in less obstruction tolaser penetration by the layers, so each layer was fully illuminated by the beams, thus yielding better results than the othertwo inclinations. Although the results for θB = 77.3◦ werecomparable to those for the optimal inclination of θB = 59.8◦ ,at θB = 77.3◦ , actual leaf-inclination data was required forthe leaf-inclination correction. In contrast, because θB = 59.8◦
is nearly equal to the particular angle of 57.5◦ [1], [3], leaf

HOSOI AND OMASA: VOXEL-BASED 3-D MODELING OF INDIVIDUAL TREES FOR ESTIMATING LAD

inclination can be well corrected for by setting α(θ) = 1.1[see (13)] without requiring any actual leaf-inclination data.
Thus, θB = 59.8◦ was not only optimal for beam penetration
but also made this correction easier.
B. Errors in LAD and LAI Estimation
LAD and cumulative LAI profiles for C. sasanqua in case 1
were apparently accurate, even though no correction process
was applied (Figs. 6 and 7), because the effects of nonphotosynthetic tissues and leaf inclination canceled each other out. The
error caused by nonphotosynthetic tissues was almost equal tothat from leaf inclination (both errors were about 10%). Theseerrors, however, should be separated for correct estimation of
LAD and LAI.
The accuracy of LAD estimation, particularly at smaller
values of ∆H, was affected by the differences in layer boundaries between stratified clipping and lidar measurement. For the
stratified clipping, the layer boundaries were measured using aruler, and the scale was different from that used for the lidar
measurement. The resulting slight differences in the vertical
positions of the boundaries between the voxel-based model andthe actual measurements led to an error that became larger at
smaller ∆H, when the layer boundary difference was not negligible relative to ∆H. This error could be decreased by betteralignment of the boundaries in the two measurement systems.
C. Advantages of the Voxel-Based Technique
Conventional measurement by portable lidar results in uneven laser-beam penetration into each horizontal layer because
penetration is less in the upper canopy [16]. In the present
VCP method, each horizontal layer was scanned fully andevenly by optimally inclined laser beams emitted from several
measuring points surrounding the specimen, thus solving the
problem of unevenness of laser-beam penetration. Moreover,in the VCP method, full laser illumination within the canopy
made it possible to count directly the contact frequency ineach layer. The contact frequency is a useful quantity for
estimating LAD and LAI, because no assumptions about leaf
spatial distribution, shape, or size are required [1]. Therefore,
underestimation caused by nonrandom foliage distribution, as
reported previously [18], [21], can be avoided.
The voxel-based representation of the trees obtained by the
VCP method has other advantages as well.
1) LAD and LAI computations are based not only on the
points at which laser beams are intercepted, but also
on those, through which the beams pass unobstructed.
Because the information about each point is represented
as a voxel with a particular attribution, voxel-based computation of LAD and LAI is easy.
2) The voxel-based 3-D array allows easy extraction of any
sites required for computation by retrieving the voxel
coordinates.
3) Overlapping points in the registered lidar data are represented as a single voxel; that is, the registered data
are merged through the voxelization process. Therefore,
overestimation of LAD caused by counting overlapping
points is avoided through voxelization.

3617

4) Voxelized lidar data can be treated as a 3-D image, sovalues can be computed by image-processing techniques.
For example, nonphotosynthetic tissues were easily removed by subtracting the voxelized nonphotosynthetic
tissues from the voxel-based tree model using image
processing.
D. Potential for Field Use and Other Applications
The significance of this VCP method is that LAD and LAIof trees can be estimated accurately without destruction oftrees. This significance is greatly enhanced when the method
is extended in field studies to larger trees. Optimally inclined
laser beams with a zenith angle near 57.5◦ contributed to thebetter estimation of LAD in this paper. This particular angle
also offered better estimation in previous field studies that used
the gap-fraction method because of the insensitivity of thatangle to leaf inclination [1], [23], [24]. Thus, this particularangle becomes a criterion for optimal beam inclination.
In the present technique, it is necessary to remove only
nonphotosynthetic tissues from lidar data of leaved trees. Although the leaves were removed manually in this experiment,
the natural defoliation of deciduous trees would permit easy
data collection of only nonphotosynthetic tissues. Developinga method for removing nonphotosynthetic tissues using difference in reflectance of each part, i.e., leaf, branch, and trunk, of
evergreen trees is a future challenge.
This precise voxel-based tree model is expected to be usefulfor obtaining other types of information about a tree. For
example, Sinoquet and Rivet [25] computed branching pattern,
fruit and shoot distribution, and plant topology using a precise
3-D model representing all vegetative components, which they
constructed with an electromagnetic digitizer. Practical use of
their method is limited, however, because a very large number of elements must be measured manually point by point.
However, the present voxel-based model can easily provide
estimates of various tree parameters by a 3-D image analysis.
The measurement range of our lidar system (less than 10 m)
is insufficient for field application to large trees. However, acurrently commercially available lidar system has a measurement range of up to several hundred meters, and this system
could be used to test the validity of our method for large-tree
measurements.
V. C ONCLUSION
The VCP method for accurate estimation of LAD and cumulative LAI profiles, which used precise voxel-based treemodels produced by high-resolution portable scanning lidar,
was demonstrated. Sufficient amounts of evenly distributed dataon each horizontal layer throughout the canopy were collectedby using optimally inclined laser beams and by taking measurements from symmetrical azimuthal measurement points.
The data were converted into a voxel-based 3-D model that
reproduced each tree precisely, including within the canopy.
This precise voxel model made it possible to compute LAD and
LAI by directly counting the contact frequency in each layer ofthe experimental trees.

3618IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 12, DECEMBER 2006

In addition, corrections were applied for leaf inclinationand the presence of nonphotosynthetic tissues, thus reducing
the estimation error. A beam incident zenith angle near 57.5◦
made it possible to apply a good correction for leaf inclination
without measuring the actual leaf inclination. By using a voxelized model, the nonphotosynthetic tissues were successfully
removed by image-processing techniques. Therefore, it waspossible to obtain good estimations of LAD and LAI with only
small errors. This paper showed the fundamental method. The
significance of the method is greatly enhanced when the method
is applied in field studies to large trees, but further development
is needed before it is ready for field application. We believethat additional work will enhance the validity of the method forvarious field applications.
R EFERENCES[1] M. Weiss, F. Baret, G. J. Smith, I. Jonckheere, and P. Coppin, “Review ofmethods for in situ leaf area index (LAI) determination Part II. Estimationof LAI, errors and sampling,” Agric. For. Meteorol., vol. 121, no. 1/2,
pp. 37–53, Jan. 2004.
[2] J. Warren-Wilson, “Inclined point quadrats,” New Phytol., vol. 59, no. 1,
pp. 1–8, 1960.
[3] ——, “Estimation of foliage denseness and foliage angle by inclined point
quadrats,” Aust. J. Bot., vol. 11, no. 1, pp. 95–105, 1963.
[4] I. Jonckheere, S. Fleck, K. Nackaerts, B. Muys, P. Coppin, M. Weiss, and
F. Baret, “Review of methods for in situ leaf area index determination
Part I. Theories, sensors and hemispherical photography,” Agric. For.
Meteorol., vol. 121, no. 1/2, pp. 19–35, 2004.
[5] J. M. Norman and G. S. Campbell, “Canopy structure,” in Plant Physiological Ecology: Field Methods and Instrumentation, R. W. Pearcy,
J. Ehleringer, H. A. Mooney, and P. W. Rundel, Eds. London, U.K.:
Chapman & Hall, 1989, pp. 301–326.
[6] J. M. Welles and J. M. Norman, “Instrument for indirect measurement ofcanopy architecture,” Agron. J., vol. 83, no. 5, pp. 818–825, 1991.
[7] M. Nilsson, “Estimation of tree heights and stand volume using an airborne lidar system,” Remote Sens. Environ., vol. 56, no. 1, pp. 1–7,
Apr. 1996.
[8] E. Næsset, “Determination of mean tree height of forest stands using
airborne laser scanner data,” ISPRS J. Photogramm., vol. 52, no. 2,
pp. 49–56, Apr. 1997.
[9] J. E. Means, S. A. Acker, D. J. Harding, J. B. Blair, M. A. Lefsky,
W. B. Cohen, M. E. Harmon, and W. A. McKee, “Use of large-footprint
scanning airborne lidar to estimate forest stand characteristics in the Western Cascades of Oregon,” Remote Sens. Environ., vol. 67, no. 3, pp. 298–
308, Mar. 1999.
[10] M. A. Lefsky, W. B. Cohen, G. G. Parker, and D. J. Harding, “Lidar
remote sensing for ecosystem studies,” Bioscience, vol. 52, no. 1, pp. 19–
30, Jan. 2002.
[11] D. J. Harding, M. A. Lefsky, G. G. Parker, and J. B. Blair, “Laser altimeter canopy height profiles: Methods and validation for closed-canopy,
broadleaf forests,” Remote Sens. Environ., vol. 76, no. 3, pp. 283–297,
Jun. 2001.
[12] K. Omasa, Y. Akiyama, Y. Ishigami, and K. Yoshimi, “3-D remote sensingof woody canopy heights using a scanning helicopter-borne lidar systemwith high spatial resolution,” J. Remote Sens. Soc. Jpn., vol. 20, no. 4,
pp. 394–406, 2000.
[13] K. Omasa, G. Y. Qiu, K. Watanuki, K. Yoshimi, and Y. Akiyama, “Accurate estimation of forest carbon stocks by 3-D remote sensing of individual
trees,” Environ. Sci. Technol., vol. 37, no. 6, pp. 1198–1201, Mar. 2003.
[14] T. Brandtberg, T. A. Warner, R. E. Landenberger, and J. B. McGraw,“Detection and analysis of individual leaf-off tree crowns in small footprint, high sampling density lidar data from the eastern deciduous forest
in North America,” Remote Sens. Environ., vol. 85, no. 3, pp. 290–303,
2003.

[15] M. Maltamo, K. Eerikäinen, J. Pitkänen, J. Hyyppä, and M. Vehmas,“Estimation of timber volume and stem density based on scanning laser
altimetry and expected tree size distribution functions,” Remote Sens.
Environ., vol. 90, no. 3, pp. 319–330, Apr. 2004.
[16] J. L. Lovell, D. L. B. Jupp, D. S. Culvenor, and N. C. Coops,“Using airborne and ground-based ranging lidar to measure canopy
structure in Australian forests,” Can. J. Remote Sens., vol. 29, no. 5,
pp. 607–622, 2003.
[17] K. Omasa, Y. Urano, H. Oguma, and Y. Fujinuma, “Mapping of tree
position of Larix leptolepis woods and estimation of diameter at breast
height (DBH) and biomass of the trees using range data measured by aportable scanning lidar,” J. Remote Sens. Soc. Jpn., vol. 22, no. 5, pp. 550–
557, 2002.
[18] P. J. Radtke and P. V. Bolstad, “Laser point-quadrat sampling for estimating foliage-height profiles in broad-leaved forests,” Can. J. For. Res.,
vol. 31, no. 3, pp. 410–418, 2001.
[19] G. G. Parker, D. J. Harding, and M. L. Berger, “A portable LIDAR system
for rapid determination of forest canopy structure,” J. Appl. Ecol., vol. 41,
no. 4, pp. 755–767, Aug. 2004.
[20] R. H. MacArthur and H. S. Horn, “Foliage profile by vertical measurements,” Ecology, vol. 50, no. 5, pp. 802–804, 1969.
[21] J. D. Aber, “Foliage-height profiles and succession in Northern hardwood
forests,” Ecology, vol. 60, no. 1, pp. 18–23, 1979.
[22] P. J. Besl and N. D. McKay, “A method for registration of 3-D shapes,”
IEEE Trans. Pattern Anal. Mach. Intell., vol. 14, no. 2, pp. 239–256,
Feb. 1992.
[23] R. Bonhomme, C. V. Grancher, and P. Chartier, “Use of hemisphericalphotographs for determining leaf area index of young crops,” Photosynthetica, vol. 8, no. 3, pp. 299–301, 1974.
[24] S. G. Leblanc and J. M. Chen, “A practical scheme for correcting multiple
scattering effects on optical LAI measurements,” Agric. For. Meteorol.,vol. 110, no. 2, pp. 125–139, Dec. 2001.
[25] H. Sinoquet and P. Rivet, “Measurement and visualization of the architecture of an adult tree based on a three-dimensional digitising device,”
Trees, vol. 11, no. 5, pp. 265–270, 1997.

Fumiki Hosoi received the M.E. degree from University of Tokyo, Tokyo, Japan, in 1995, where hehas been working toward the Ph.D. degree in the
Department of Agricultural and Life Sciences since2005.
In 1995, he joined Opto-Technology Laboratory,
Furukawa Electric Company, Ltd., Ichihara, Chiba,
Japan.

Kenji Omasa received the M.E. degree from Ehime
University, Ehime, Japan, in 1975 and the Ph.D.
degree in engineering from the University of Tokyo,
Tokyo, Japan.
He is currently a Professor with the Department ofAgricultural and Life Sciences, University of Tokyo.
His research interests include imaging of cell levelto plants, remote sensing, modeling of ecosystems,
analysis of global change effects on ecosystems andinformation engineering for biological and environmental systems.