Context-sensitive Filtering of Terrain Data based on Multi Scale Analysis

Steffen Goebbels and Regina Pohle-Fr¨ohlich

Faculty of Electrical Engineering and Computer Science, Niederrhein University of Applied Sciences,

Reinarzstr. 49, 47805 Krefeld, Germany

Keywords:

Surface Triangulation, Polygon Count Reduction, Wavelet Transform.

Abstract:

This paper addresses triangle count reduction for regular terrain triangulations. Supported by additional infor-

mation available from cadastre data, non-linear low pass ﬁltering based on wavelet coefﬁcients is applied. The

aim is to generate a simpliﬁed view of a city model’s terrain so that still landmarks can be easily recognized.

Height data, that are not critical for visual impressions, are compressed with lower resolution than height data

belonging to landmarks. Results are presented for the German city of Krefeld.

1 INTRODUCTION

Airborne laser scanning of terrain results in a huge

amount of vertices that is not immediately suitable

for drawing a surface triangulation. In this paper an

approach is discussed where additional information

from register data about ground usage is used for data

reduction. From that information we derive precision

that is necessary to draw certain areas. Goal is that the

terrain still should be recognizable while the number

of triangles is minimized.

Wavelet transforms are an often used method for

point cloud compression, ﬁltering and extraction of

bare earth, cf. (Vu and Tokunaga, 2004; Wei and Bar-

tels, 2006; Xu et al., 2007; Hu et al., 2014). Besides

polygon merging, polygon elimination, vertex elim-

ination, edge collapsing, and different retiling tech-

niques, multi-resolution wavelet analysis is one of

the main methods for mesh reduction, cf. (Kalber-

matten, 2010, pp. 45–77), (Mocanu et al., 2011;

W¨unsche, 1998; Wiman and Yuchu, 2009). For ex-

ample (Olanda et al., 2014; Bjørke and Nilsen, 2003;

Bruun and Nilsen, 2003) give an overview of wavelet

techniques for simpliﬁcation of digital terrain mo-

dels. These techniques do not use additional infor-

mation besides height data. Studies were carried

out with biorthogonal average interpolating wavelets.

The property of conserving average height values at

different scales yields good approximationsof terrain.

As wavelet transform allows for local frequency

ﬁltering, it is predestined for local, context-sensitive

data reduction. However, knowledge of local proper-

ties of ground is needed. For example, good visible

sharp edges can be detected from height data using

curvature coefﬁcients. Such edges should be main-

tained and not eliminated by ﬁltering. In the example

of a city model we do not have to detect such patterns

in height data as we can use additional information

from a geographic information system (GIS)

1

. These

cadastre data describe ground usage in terms of poly-

gons that stand for streets, railway lines, bridges, etc.

This approach to local ﬁltering based on digital maps

has been proposed before, for example in (Mahler,

2001, p. 58) in the context of coast lines. New is a set

of heuristic rules to determine the local level of ﬁlter-

ing using a simple ﬁltering algorithm that differenti-

ates between low frequencies describing terrain on a

coarse scale and high frequencies needed for details.

Whereas for example in (Bjørke and Nilsen, 2003)

one biorthogonal wavelet transform is used to ﬁlter

the whole frequencyspectrum, we apply two different

transforms to minimize triangle count:

Local areas like embankments, shallow hills and

heaps that afford a higher quality of rendering will

be ﬁltered using the simplest average interpolating

wavelet that is the discontinuous Haar wavelet. We

apply non-linear ﬁltering using thresholds for wavelet

coefﬁcients. In our scenario the work with continuous

orthogonal or biorthogonal wavelets is not necessary.

In contrast, local support of the Haar wavelet helps to

be very location speciﬁc.

Haar wavelet transform leads to artifacts like visi-

ble steps if only low frequencies remain after ﬁltering.

Therefore, areas which do not require detailed ren-

1

Case study is based on geographical data, which are

protected by copyright: Geobasisdaten der Kommunen und

des Landes NRW

c

Geobasis NRW 2014.

106

Goebbels S. and Pohle-Fröhlich R..

Context-sensitive Filtering of Terrain Data based on Multi Scale Analysis.

DOI: 10.5220/0005252201060113

In Proceedings of the 10th International Conference on Computer Graphics Theory and Applications (GRAPP-2015), pages 106-113

ISBN: 978-989-758-087-1

Copyright

c

2015 SCITEPRESS (Science and Technology Publications, Lda.)

dering will be low pass ﬁltered via continuous linear

spline approximation (that also can be interpreted as

a biorthogonalwavelet transform based on B-splines).

To be more speciﬁc, we do this spline based low pass

ﬁltering for the whole scenario and not only for ar-

eas with little detailed information. Reason is, that

local areas with relevant details are not only detected

heuristically on the basis of context information but

also based on deviations between the surface of spline

approximation and unﬁltered heights. For these ar-

eas we triangulate with heights obtained from ﬁltering

with the Haar wavelet. Otherwise, heights for trian-

gulation are taken from spline approximation. Spline

approximation gives us a degree of freedom to fur-

ther decrease triangle count while preserving a more

detailed triangulation of areas that are displayed in

higher quality. We modify the precision of spline in-

terpolation to ﬁnd a minimum number of triangles.

2 TRIANGULATION BASED ON

FILTERED HEIGHTS

Our algorithm consists of two steps. In the ﬁrst step,

height data are ﬁltered with spline approximation and

wavelet transform as described in subsequent sec-

tions. In a second step, mesh simpliﬁcation can be

performed more efﬁcient due to the adjusted heights.

For the second step a variety of well established algo-

rithms like edge collapsing can be used, cf. (Garland

and Heckbert, 1997; Silva and Gomes, 2004). Also,

the simpliﬁcation can be made dependent on distance

to the viewer by using geometry clipmaps as intro-

duced in (Losasso and Hoppe, 2004). However, to test

our approach we decided to use a very basic algorithm

that is described in this section. The algorithm works

with an equidistant mesh where 400 square kilometers

are represented by (5 ·2

9

+ 1)

2

grid points. Triangu-

lation is done for each square kilometer separately on

the basis of squares that are split up into two triangles

using a diagonal. We divide a square recursively into

four smaller squares of equal size in case the surfaces

of triangles deviate more than ε from previously ﬁl-

tered heights. A predeﬁned accuracy level ε = 0.2 m

is a good compromise ﬁtting with the accuracy of un-

derlying laser scanning data, but we also show num-

bers for ε = 0.05 and ε = 0.1 m. Because mesh con-

nectivity changes, larger ε might lead to visible gaps

between surface triangles.

To avoid unsightly edges when displaying diag-

onal structures, we use following rule to choose the

diagonal for splitting up a square into two triangles.

We compute the four absolute differences between the

arithmetic mean of all four heights of square vertices

Figure 1: Triangulation based on an equidistant mesh.

and each of the four heights. If P is a vertex with a

maximum difference, then we select the diagonal for

which P is not end point. Further down, this rule is

described with formula (1). Because no orientation of

the diagonal is a-priori favored over the other, we do

not get a type-1 triangulation for which this would be

the case. A typical outcome is shown in Figure 1.

3 LINEAR SPLINE

APPROXIMATION

Multi scale analysis on the basis of piecewise lin-

ear, continuous functions with local support seems to

be adequate in connection with drawing triangulated

surfaces. Therefore, in this section we reduce num-

ber of triangles with the means of a one-dimensional

equidistant L

2

best linear continuous spline approx-

imation on different scales. The result is a coarse

model of terrain. Our approach to handle visible de-

tailed structures is described in the next section. Pure

edge collapsing methods like (Garland and Heckbert,

1997; Silva and Gomes, 2004) simplify the mesh in

planar zones. Using spline approximation we signif-

icantly simplify in non planar zones as well. If such

zones are not important for visual impression, we will

not increase their level of detail later.

Using continuous, piecewise linear hat function

Λ(x) :=

1+ x : −1 ≤ x < 0

1−x : 0 ≤ x ≤ 1

0 : x 6∈ [− 1,1]

we deﬁne scales V

i

, V

i

⊂V

i+1

,

V

i

:=

(

m·2

i

∑

k=0

c

k

Λ

i,k

(x) : c

0

,. .. ,c

m·2

i

∈ R

)

,

with Λ

i,k

(x) := m·2

i

Λ(m·2

i

·x−k)

[0,1]

where “|

[0,1]

”

denotes the restriction of a function to the real interval

[0,1] and m is an odd natural number – in our case

m = 5 because of (5·2

9

+ 1)

2

grid points.

Each scale has a B-spline basis of translated and

scaled hat functions. They do not constitute an or-

thogonal set of functions but can be seen in context of

biorthogonal wavelets. For the sake of simplicity, we

directly compute best L

2

approximations and do not

use terms of biorthogonal wavelets.

Context-sensitiveFilteringofTerrainDatabasedonMultiScaleAnalysis

107

We simplify a function

f

i+1

=

m·2

i+1

∑

k=0

c

i+1,k

Λ

i+1,k

(x) ∈V

i+1

to a coarser function f

i

=

∑

m·2

i

k=0

c

i,k

Λ

i,k

(x) ∈ V

i

by

computing the best approximation with regard to L

2

inner product < f,g>:=

R

1

0

f(x) ·g(x) dx. It is well

known, that coefﬁcients c

i,k

form the unique solution

of normal equations

<Λ

i,0

,Λ

i,0

> . .. <Λ

i,0

,Λ

i,m·2

i

>

.

.

.

.

.

.

<Λ

i,m·2

i

,Λ

i,0

> ... <Λ

i,m·2

i

,Λ

i,m·2

i

>

·

c

i,0

.

.

.

c

i,m·2

i

=

< f

i+1

,Λ

i,0

>

.

.

.

< f

i+1

,Λ

i,m·2

i

>

.

By considering the representation of f

i+1

as a lin-

ear combination of functions Λ

i+1,k

, and by solving

elementary integrals, we get

1

6

2 1 0 0 ... 0 0 0

1 4 1 0 ... 0 0 0

0 1 4 1 ... 0 0 0

.

.

.

0 0 0 0 ... 1 4 1

0 0 0 0 ... 0 1 2

·

c

i,0

.

.

.

c

i,m·2

i

=

5

12

1

2

1

12

0 0 0 0 ... 0

1

12

1

2

5

6

1

2

1

12

0 0 ... 0

0 0

1

12

1

2

5

6

1

2

1

12

... 0

.

.

.

0 0 0 0 0 0 0 ...

1

12

0 0 0 0 0 0 0 ...

5

12

·

c

i+1,0

.

.

.

c

i+1,m·2

i+1

where the tridiagonal matrix on the left side belongs

to R

(m·2

i

+1)×(m·2

i

+1)

, and the matrix on the right side

is element of R

(m·2

i

+1)×(m·2

i+1

+1)

. Thus, we can start

on scale n with a representation

f

n

(x) =

m·2

n

∑

k=0

c

n,k

Λ

n,k

(x)

where c

n,k

= a

k

/(m·2

n

) for m·2

n

+1 givenheight val-

ues (a

0

,a

1

,. .. ,a

m·2

n

). By solving normal equations

we successively get coefﬁcients c

i,k

that can be inter-

preted as heights m ·2

i

c

i,k

on the coarser grit of scale

i, 0 ≤i < n.

Tridiagonal systems of linear equations can be

solved with Thomas-algorithm, cf. (Mooney and

Swift, 1999, p. 235), so that we get a pro-

cedure APPROX(c

i+1,0

,. .. ,c

i+1,m·2

i+1

) that returns

(c

i,0

,c

i,1

,. .. ,c

i,m·2

i

). Algorithm 1 shows how we ap-

ply this one-dimensional best approximation to all

rows and then to columns of a matrix. This is similar

to the simpliﬁed non-standard wavelet transform that

will be discussed in Section 4. In order to ensure that

the result consists of heights and not of spline coefﬁ-

cients we have to normalize with factors 1/(m·2

i+1

)

and m·2

i

, respectively, for example

B[k, 0:m·2

i

]:=m·2

i

· APPROX(

1

m·2

i+1

·A

i+1

[k, 0:m·2

i+1

]).

As best approximation is a linear calculation, factors

can be shortened to

1

2

, cf. Algorithm 1.

Algorithm 1: Height data reduction. Height matrix A

i+1

∈

R

(m·2

i+1

+1)×(m·2

i+1

+1)

is mapped to A

i

∈ R

(m·2

i

+1)×(m·2

i

+1)

using one-dimensional best approximations of rows and

columns.

procedure HEIGHTDATAREDUCTION(A

i+1

)

for k = 0 : m·2

i+1

do

B[k, 0:m·2

i

]:=

1

2

·APPROX(A

i+1

[k, 0:m·2

i+1

])

for k = 0 : m·2

i

do

A

i

[0:m·2

i

,k] :=

1

2

·APPROX(B[0:m·2

i+1

,k])

return A

i

Figure 2: Reduction of triangles with three steps of linear

spline approximation (triangles are transparent).

Computing of spline approximation consists of

n−i calls of Algorithm 1 so that we get a (m·2

i

)×(m·

2

i

)-matrix A

i

of ﬁltered heights. Figure 2 shows the

result of three calls in our example. The number n−i

of calls will be chosen with respect to the result of

detail analysis such that the overall number of trian-

gles becomes small. From A

i

we derive approximate

heights h

i

( j,k) for all (m ·2

n

+ 1) ×(m ·2

n

+ 1) grid

points ( j,k) through interpolation. The given point

( j,k) lies within the coarser net’s square with ver-

tices (x

0

,y

0

), (x

0

+ 1,y

0

), (x

0

,y

0

+ 1), (x

0

+ 1,y

0

+ 1)

where x

0

and y

0

are down rounded integers deﬁned

by x

0

:= ⌊j/2

n−i

⌋ if j < m ·2

n

and y

0

:= ⌊ k/2

n−i

⌋

if k < m · 2

n

. For end points with i = m ·2

n

set

x

0

:= m·2

i

−1, and for k = m·2

n

let y

0

:= m·2

i

−1.

In order to compute an approximate height with re-

spect to scale i, we have to triangulate the square by

choosing a diagonal as described in Section 2. Let

h := (A

i

[x

0

,y

0

] + A

i

[x

0

+ 1, y

0

] +

+A

i

[x

0

,y

0

+ 1] + A

i

[x

0

+ 1, y

0

+ 1])/4

be the average of heights. If

max{|A

i

[x

0

,y

0

] −h|,|A

i

[x

0

+ 1, y

0

+ 1] −h|} (1)

≤ max{|A

i

[x

0

+ 1, y

0

] −

h|,|A

i

[x

0

,y

0

+ 1] −h|}

then the diagonal connects (x

0

,y

0

) and (x

0

+1, y

0

+1).

Otherwise (x

0

,y

0

+ 1) and (x

0

+ 1, y

0

) are connected.

GRAPP2015-InternationalConferenceonComputerGraphicsTheoryandApplications

108

The surface of this triangulation at point (i, j) de-

ﬁnes height h

i

( j,k). As a result, we get a coarse ap-

proximation of terrain that lacks detailed structures.

Please also note that even in grid points the heights

might be changed signiﬁcantly by spline approxima-

tions (see ﬁgure 2). Compared to original data, these

changes might lead to increased detail in some areas

after the next processing steps (see lower left corner

of pictures in Figure 3).

4 DISCRETE WAVELET

TRANSFORM

We now deal with areas like heaps, bridge ramps

and embankments that have to be shown in more de-

tail. Despite higher accuracy is required, these areas

also can be smoothened using a low pass ﬁlter based

on discrete wavelet transform (DWT) for the Haar

wavelet. In one dimension Algorithm 2 shows how to

transform a vector (a

0

,. .. ,a

m·2

n

−1

) (not dealing with

last element a

m·2

n

). In each step arithmetic means b

j

of pairs of values a

2j

and a

2j+1

are computed. They

represent data on a coarser scale. By adding wavelet

coefﬁcients b

j+N/2

to the mean value b

j

, a

2j

is recon-

structed. Subtracting b

j+N/2

leads to a

2j+1

. Wavelet

coefﬁcients can be interpreted as amplitudes of local

frequencies. Haar scaling and wavelet functions form

an orthonormal basis if correctly normed, cf. (Stroll-

nitz et al., 1995). If wavelet coefﬁcients are seen as

scalar factors regarding this wavelet basis, then one

has to divide by

√

2 instead of 2 in Algorithm 2 and

gets coefﬁcients 2

−i/2

c

k

as a result of Algorithm 2.

Algorithm 2: One-dimensional DWT for size m·2

n

.

procedure STEP(~a)

N :=SIZE(~a)

for j = 0 : N/2−1 do

b

j

:=

a

2j

+a

2j+1

2

, b

j+N/2

:=

a

2j

−a

2j+1

2

return

~

b

procedure DWT(~a, m, n)

if n ≤ 0 then return~a

~

b := STEP(~a)

~c[0:m·2

n−1

−1] :=DWT(

~

b[0:m·2

n−1

−1], m, n−1)

~c[m·2

n−1

: m·2

n

−1] :=

~

b[m·2

n−1

: m·2

n

−1]

return~c

The one-dimensional transform often is extended

to two dimensions by iteration. Another approach,

that covers iteration as well but is not discussed here,

is the use of multi dimensional wavelets. In the lit-

erature two iterative methods are prominent: stan-

dard and non-standard transform, cf. (Strollnitz et al.,

1995). The standard transform is a tensor product

method that applies the one-dimensional transform to

all rows of a matrix A :=

a

j,k

resulting in matrix B:

B = [DWT((a

j,0

,a

j,1

,. .. ,a

j,m·2

n

−1

))]

j=0,...,m·2

n

−1

.

In a second step each column of B gets transformed

into a column DWT((b

0,k

,b

1,k

,. .. ,b

m·2

n

−1,k

)), 0 ≤

k ≤m·2

n

−1 of resulting matrix C. However, for our

purposeof reducing height differencesit appears to be

a disadvantage of the standard transform that wavelet

coefﬁcients themselves get transformed during col-

umn transforms. For image manipulation instead of

the standard wavelet transform often a non-standard

variant is used that is slightly faster (but of same linear

order). The non-standard transform iterates between

row and column transformation steps and excludes

previously computed wavelet coefﬁcients in follow-

on steps, see Algorithm 3.

Algorithm 3: Two-dimensional non-standard DWT.

procedure NST DWT(A)

m·2

n

×m·2

n

:= SIZE(A)

if n ≤ 0 then return A

for j = 0 : m·2

n

−1 do

B[ j,0:m·2

n

−1] := STEP(A[ j, 0:m·2

n

−1])

for j = 0 : m·2

n

−1 do

⊲ Use 0 : m·2

n−1

−1 for simpliﬁed transform.

C[0:m·2

n

−1, j] := STEP(B[0:m·2

n

−1, j])

D[0 : m·2

n−1

−1,0 : m·2

n−1

−1] :=

NST

DWT(C[0:m·2

n−1

−1, 0:m·2

n−1

−1]))

D[m·2

n−1

: m·2

n

−1,0 : m·2

n

−1] :=

C[m·2

n−1

: m·2

n

−1,0 : m·2

n

−1]

D[0 : m·2

n−1

−1,2

n−1

: m·2

n

−1] :=

C[0 : m·2

n−1

−1,m·2

n−1

: m·2

n

−1]

return D

Even with the non-standard transform, there is

one column transformation step applied to previously

computed coefﬁcients of rows. Therefore, we choose

to use a simpliﬁed version of non-standard transform

that has been described in (Kopp and Purgathofer,

1998). The only change to Algorithm 3 is that the

second loop stops at m·2

n−1

−1 instead of m·2

n

−1.

The number of column transforms then only is half

the number of row transforms. The advantage is that

height differences immediately become apparent in

terms of coefﬁcients. The associated inverse compu-

tation for a single component a

j,k

reads:

a

j,k

= c

⌊j2

−n

⌋,⌊k2

−n

⌋

(2)

+

n−1

∑

i=0

[c

m·2

i

+⌊j2

i−n

⌋,⌊k2

i−n

⌋

ψ( j2

i−n

−⌊j2

i−n

⌋)

+c

⌊j2

i+1−n

⌋,m·2

i

+⌊k2

i−n

⌋

ψ(k2

i−n

−⌊k2

i−n

⌋)]

Context-sensitiveFilteringofTerrainDatabasedonMultiScaleAnalysis

109

where ψ is the Haar wavelet function

Ψ(x) :=

1, 0 ≤x <

1

2

,

−1,

1

2

≤x < 1,

0, else

that determines whether wavelet coefﬁcients have to

be added or subtracted.

Non-linear low pass ﬁltering consists of cut-

ting off higher summands depending on coefﬁcients

c

m·2

i

+⌊j2

i−n

⌋,⌊k2

i−n

⌋

and c

⌊j2

i+1−n

⌋,m·2

i

+⌊k2

i−n

⌋

during re-

construction of height data with equation (2). To

this end we need context-sensitivesequences of thres-

holds. For a point ( j,k) let (t

0

,t

1

,. .. ,t

n−1

) be such

a sequence of non-negative numbers. We investigate

two variants of non-linear ﬁltering:

Variant 1: Let

˜a

j,k

:= c

⌊j2

−n

⌋,⌊k2

−n

⌋

(3)

+

n−1

∑

i=0

f

j,k

c

m·2

i

+⌊j2

i−n

⌋,⌊k2

i−n

⌋

,i

·ψ( j2

i−n

−⌊j2

i−n

⌋)

+

n−1

∑

i=0

f

j,k

c

⌊j2

i+1−n

⌋,m·2

i

+⌊k2

i−n

⌋

,i

·ψ(k2

i−n

−⌊k2

i−n

⌋)

be a component of the ﬁltered matrix where depend-

ing on point ( j,k) and frequency i

f

j,k

(x,i) :=

x : |x| ≥ t

i

0 : |x| < t

i

.

This variant allows to eliminate intermediate frequen-

cies while higher frequencies might remain. Because

high frequencies have a higher impact on number of

triangles, we compare this approach with

Variant 2: Tails of reconstruction sum in equation

(2) (i.e. high frequencies) get cut off. Let M

1

be the

set of indices i ∈{0,...,n −1} with

c

⌊j2

i+1−n

⌋,m·2

i

+⌊k2

i−n

⌋

≥t

i

and M

2

be the subset of {0,... ,n −1} where

c

m·2

i

+⌊j2

i−n

⌋,⌊k2

i−n

⌋

≥t

i

.

Then we cut off by using maximum values m

1

:=

max(M

1

∪{−1}) and m

2

:= max(M

2

∪{−1}):

˜a

j,k

:= c

⌊j2

−n

⌋,⌊k2

−n

⌋

(4)

+

max{m

1

−1,m

2

}

∑

i=0

c

m·2

i

+⌊j2

i−n

⌋,⌊k2

i−n

⌋

·ψ( j2

i−n

−⌊j2

i−n

⌋)

+

max{m

1

,m

2

}

∑

i=0

c

⌊j2

i+1−n

⌋,m·2

i

+⌊k2

i−n

⌋

·ψ(k2

i−n

−⌊k2

i−n

⌋).

If we alternate between reconstruction steps for

columns and rows beginning with highest frequency

Figure 3: First row: unﬁltered height data of a motorway,

last row: ﬁlter (4) applied with ε = 0.2 m and δ = 4.0 m

where height deviations are marked with colors in left pic-

ture: The bigger the difference, the darker the color, red

marks points where unﬁltered heights are larger than ﬁltered

heights, blue indicates the opposite.

n−1 going down to 0 then the formula starts adding

summands when encountering a ﬁrst coefﬁcient that

exceeds its threshold.

In contrast to JPEG 2000 we do not quantize

wavelet coefﬁcients. For triangle count reduction

quantization does not appear to work as well as omis-

sion of coefﬁcients. All non-zero coefﬁcients of

Haar wavelets deﬁne height alterations that might

lead to additional triangles. But quantization does re-

duce storage space and transmission bandwidth, see

(Olanda et al., 2014).

5 HEURISTICS BASED LOW

PASS FILTERING

Our aim is to preserve salient structures like ramps,

surroundings of bridges, hills, and shore lines. Poly-

gons covering areas of these structures are available

from GIS data. Whereas ﬁelds can be drawn with re-

duced precision, a more detailed drawing is needed

for streets and railway lines. This section deals with

heuristics to implement context-sensitiveﬁltering that

takes care of such differences.

If a difference of unﬁltered height and the result

of spline approximation at a point ( j,k) is less than

0.3 m, then this difference is not disturbing. The

threshold of 0.3 m is a parameter different than the

accuracy ε that is used for mesh simpliﬁcation on the

basis of previously ﬁltered heights. If the difference

exceeds 0.3 m, in certain circumstances we compute

a more speciﬁc height using low pass wavelet trans-

form. In this case, depending on terrain usage and

height of (j,k), we divide coordinates ( j, k) into three

main classes:

1. Class 1: A point ( j,k) is near to a point P, iff P

is located within a square covering 9

2

grid points

GRAPP2015-InternationalConferenceonComputerGraphicsTheoryandApplications

110

Figure 4: Triangulation based on unﬁltered data (ﬁrst row)

and based on non-linear ﬁltereddata (second row) with ﬁlter

(4) and ε = 0.2 m, δ = 4.0 m. Last row shows differences

in the same manner as in Figure 3.

with ( j, k) at its center. The size of this neigh-

borhood was determined experimentally. Smaller

squares lead to poorer quality whereas larger

neighborhoods increase number of triangles un-

necessarily. In our context, high precision is re-

quired for those points ( j,k) that are near to a

point P which fulﬁlls one of following conditions.

• P is next to a shore line of a river or lake.

• P is part of a railway embankment. If ( j,k) is

an interior point of the embankment then even

extra precision is needed to avoid non smooth

railway tracks.

• P lies within the border polygon of a bridge.

Surroundings of bridges need to be accurate

so that bridges ﬁt into terrain. In our exam-

ple, even the ﬁnest triangle size is too coarse

to avoid gaps between bridges and their ramps.

This problem is ﬁxed by extending bridges.

• P lies within the polygon of a street that inter-

sects at least one polygon of a bridge. Such

streets can be part of ramps.

2. Class 2: Terrain of our city model is ﬂat apart

from some heaps and local hills. Therefore, zones

outside class 1 with unﬁltered heights that dif-

fer more than a parameter δ from the heights of

spline approximation should be drawn with inter-

mediate precision. Hence heap contours remain

clearly visible. If a forest grows on a hill then even

less precision might be needed as trees hide the

ground. In our example, hilly terrain covered by

trees is limited, so that special treatment of forests

does not improve triangle count signiﬁcantly. It is

not performed.

1 2 3 4 5 6

0.5

1

1.5

2

2.5

x 10

6

Figure 5: Triangle count for different numbers of spline it-

erations and accuracy levels ε, see Table 3. Broken line:

ε = 0.2 m, solid line: ε = 0.05 m (δ = 4.0 m, ﬁlter (4)).

3. Class 3: All remaining terrain can be drawn with

poor precision, i. e., heights are taken directly

from spline approximation.

We apply the simpliﬁed non-standard wavelet

transform for the grid with m·2

n

= 5·2

9

points. Low

pass ﬁltering for points of classes 1 and 2 is done us-

ing speciﬁc sets of threshold values (t

0

,t

1

,. .. ,t

n−1

).

Height values resulting from low pass ﬁltering re-

place heights computed by spline approximation only

if their distance to spline heights exceeds 0.3 m (less

can’t be seen).

For δ = 4 m good looking results were obtained

with threshold values from Table 1. We determined

the values experimentally. Instead, one could deﬁne

a distance measure between exact and ﬁltered heights

(for example an L

2

distance) together with separate

local error bounds for areas of class 1 and class 2 as

constraints. Then separately for both classes, varia-

tion of thresholds can yield a set of values for which

the number of triangles becomes minimal.

Our data set contains an error where a hole ap-

pears in the middle of the harbor - probably an artifact

from a ship. Generally, heights of grid points within

polygons deﬁning areas of water need to be constant,

because water surfaces should be ﬂat. Linear spline

approximation and wavelet ﬁltering do not guarantee

this outcome. This is achieved by assigning an av-

erage value to all grid points of water polygons in a

post-processing step.

6 RESULTS

We count triangles but we do not present a global

measure for quality of ﬁltered terrain. Because we al-

low to over-simplifyareas that are uncritical for visual

impressions, measures like scene-wide L

2

metrics that

involve only height differences between ﬁltered and

unﬁltered heights are not sufﬁcient. A comparison of

contour lines as done in (Bjørke and Nilsen, 2003)

Context-sensitiveFilteringofTerrainDatabasedonMultiScaleAnalysis

111

Table 1: Thresholds.

coefﬁcient 0 1 2 3 4 5 6 7 8

class 1 0 0.1 0.1 0.2 0.2 0.2 0.3 0.3 0.4

class 1, inner points of railway embankments 0 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

class 2 0 0.1 0.2 0.2 0.3 0.4 0.5 0.5 0.5

class 2, forest 0 0.1 0.2 0.2 0.3 20 20 20 20

Table 2: Triangle count for δ = 4 m and three iterations of spline approximation.

triangulation triangles test data triangles spline approx. triangles, detail triangles. detail

accuracy level ε (no ﬁltering) with unﬁltered detail ﬁltered with (3) ﬁltered with (4)

0.2 m 2901422 1221938 962012 940214

0.1 m 4564400 1405856 1046510 1035626

0.05 m 6144422 1434602 1067888 1058636

Table 3: Inﬂuence of spline approximation on triangle count for δ = 4 m and wavelet ﬁlter (4), see Figure5.

accuracy level ε = 0.05 m accuracy level ε = 0.2 m

6144422 triangles without compression 2901422 triangles without compression

spline triangles spline triangles triangles spline triangles

iterations approximation ﬁltered detail approximation ﬁltered detail

1 1820300 2081120 1201376 1450412

2 544388 1158722 420140 984206

3 170210 1058636 135488 940214

4 48218 1165586 44372 1049498

5 12422 1444874 12308 1312340

6 3164 1838468 3140 1690142

might be a better tool but it lacks the ability to distin-

guish between “important” and “unimportant”.

Table 2 shows numbers of triangles in our case

study. Last two columns correspond to results of ﬁl-

tering as described. Compression rates up to a fac-

tor of 5.8 can be reached. We have added the mid-

dle column in order to demonstrate positive effects of

low pass Haar wavelet ﬁltering of detailed areas. It

shows triangle count when using exact heights instead

of heights resulting from wavelet ﬁltering of points

belonging to class 1 and 2. Heights of all other points

still are taken from spline approximation. The com-

parison of the two ﬁlter variants (3) and (4) shows that

(4) produces slightly better results. That is no surprise

as highest non-zero coefﬁcients determine size of tri-

angles. Non-zero coefﬁcients of frequencies in the

middle do not have this effect.

Finally, we have to determine a good value for the

number of iterations of spline approximation. This

number in turn depends on the rendering of detailed

areas. Too many iterations lead to larger areas where

height deviations exceed δ. The number of triangles

saved through spline approximation then is outnum-

bered by the amount of additional triangles needed to

show details. On the other hand, too few iterations

generate a too detailed spline approximation and too

many triangles. Analysis of detailed areas then be-

comes obsolete. In our setting, for δ = 4 m the low-

est triangle count is obtained with three iterations, see

Table 3 and Figure 5. The smaller δ the larger are

areas of class 2 and the less iterations of spline ap-

proximation can be used for simpliﬁcation. How δ

inﬂuences the number of spline approximations and

triangle count is shown in Table 4 and Figure 6.

With respect to different landscapes, the outcome

of our approach to ﬁltering is shown in Figures 3 and

4.

1 2 3 4 5 6

0.5

1

1.5

2

2.5

3

3.5

x 10

6

Figure 6: Triangle count depending on iterations of spline

approximation for different height parameters δ (ﬁlter (4),

ε = 0.2 m), see Table 4. Thick solid line: δ = 1 m, dashed

line: δ = 2 m, thin solid line: δ = 4 m, dash-dot line: δ = 6

m.

GRAPP2015-InternationalConferenceonComputerGraphicsTheoryandApplications

112

Table 4: Triangle count depending on iterations of spline approximation for different height parameters δ (ﬁlter (4), ε = 0.2

m), see Figure 6.

iterations → 1 2 3 4 5 6

δ = 1.0 1529450 1274276 1440734 1813544 2460890 3038612

δ = 2.0 1458482 1044968 1075772 1275776 1676342 2203256

δ = 4.0 1450412 984206 940214 1049498 1312340 1690142

δ = 6.0 1450058 978740 917552 995828 1190306 1489028

7 CONCLUSIONS AND FUTURE

WORK

Whereas the MP3 algorithm compresses music data

by utilizing a psychoacoustic model that deﬁnes ac-

curacy of coefﬁcients of a discrete cosine transform,

we follow a similar but more elementary approach to

reduce triangle count: Visually less important areas

are lossy compressed by applying a continuous best

spline approximation. More important areas are low

pass ﬁltered with a wavelet transform. Thresholds for

wavelet coefﬁcients depend on rules that take ground

usage (cadastre data) and height differences into ac-

count. This is our “psycho visual” model. Results

are promising. However, it should be investigated if

an even better quality can be obtained by additional

smoothing of certain surfaces like streets. Future

work is required to explore the method in connection

with other sets of terrain data in order to see whether

experimentally determined thresholds still work well.

The mesh simpliﬁcation algorithm of Section 2 has

been chosen for simplicity. Further work is needed

to evaluate the ﬁltering technique as a pre-processing

step that generates input for established simpliﬁcation

methods.

REFERENCES

Bjørke, J. and Nilsen, S. (2003). Wavelets applied to simpli-

ﬁcation of digital terrain models. International Jour-

nal of Geographical Information Science,17:601-621.

Bruun, B. and Nilsen, S. (2003). Wavelet representation

of large digital terrain models. Computers & Geo-

sciences, 29:695–703.

Garland, M. and Heckbert, P. (1997). Surface simpliﬁca-

tion using quadric error metrics. In Proc. 24th An-

nual Conference on Computer Graphics and Interac-

tive Techniques, SIGGRAPH ’97, pages 209–216.

Hu, H., Ding, Y., Zhu, Q., Wu, B., Lin, H., Du, Z., Zhang,

Y., and Zhang, Y. (2014). An adaptive surface ﬁlter

for airborne laser scanning point clouds by means of

regularization and bending energy. ISPRS Journal of

Photogrammetry and Remote Sensing, 92:98–111.

Kalbermatten, M. (2010). Multiscale analysis of high res-

olution digital elevation models using the wavelet

transform. EPFL Lausanne, Lausanne.

Kopp, M. and Purgathofer, W. (1998). Interleaved dimen-

sion decomposition: A new decomposition method

for wavelets and its application to computer graph-

ics. In 2nd International Conference in Central Eu-

rope on Computer Graphics and Visualisation, pages

192–199.

Losasso, F. and Hoppe, H. (2004). Geometry clipmaps: Ter-

rain rendering using nested regular grids. ACM Trans.

Graph., 23(3):769–776.

Mahler, E. (2001). Scale-Dependent Filtering of High Reso-

lution Digital Terrain Models in the Wavelet Domain.

University of Zurich, Zurich.

Mocanu, B., Tapu, R., Petrescu, T., and Tapu, E. (2011). An

experimental evaluation of 3d mesh decimation tech-

niques. In International Symposium on Signals, Cir-

cuits and Systems, pages 1–4.

Mooney, D. and Swift, R. (1999). A Course in Mathe-

matical Modeling. The Mathematical Association of

America, Washington D.C.

Olanda, R., P´erez, M., Ordu˜na, J., and Rueda, S. (2014).

Terrain data compression using wavelet-tiled pyra-

mids for online 3d terrain visualization. Interna-

tional Journal of Geographical Information Science,

28:407–425.

Silva, F. and Gomes, A. (2004). Normal-based simpliﬁca-

tion algorithm for meshes. In Proc. Theory and Prac-

tice of Computer Graphics(TPCG 04), pages 211-218.

Strollnitz, E., DeRose, T., and Salesin, D. (1995). Wavelets

for computer graphics: a primer, part 1. IEEE Com-

puter Graphics and Applications, 15:77–84.

Vu, T. and Tokunaga, M. (2004). Filtering airborne laser

scanner data: A wavelet-based clustering method.

Photogrammetric Engineering & Remote Sensing,

70:1267–1274.

Wei, H. and Bartels, M. (2006). Unsupervised segmenta-

tion using gabor wavelets and statistical features in li-

dar data analysis. In 18th International Conference on

Pattern Recognition, pages 667–670.

Wiman, H. and Yuchu, Q. (2009). Fast compression and

access of lidar point clouds using wavelets. In Urban

Remote Sensing Event, pages 1–6.

W¨unsche, B. (1998). A Survey and Evaluation of Mesh Re-

duction Techniques. University of Auckland, Auck-

land.

Xu, L., Yang, Y., Jiang, B., and Li, J. (2007). Ground extrac-

tion from airborne laser data based on wavelet analy-

sis. In Proc. MIPPR/SPIE.

Context-sensitiveFilteringofTerrainDatabasedonMultiScaleAnalysis

113