Jean-Laurent Mallet. Geomodeling
Подождите немного. Документ загружается.
148
CHAPTER
4.
DISCRETE
SMOOTH
INTERPOLATION
Among
all the
functions
</?
defined
on
fJ
and
interpolating
the
data
{<p(&}
:
I
6
I/},
we
propose calling "DSI solution"
the
function
(p
that
achieves
the
minimum
of the
generalized roughness criterion
R*
(</?)
where
F{<p(f,)
:
i
6
L} is the
family
of all the
functions
denned
on fi and
interpolating
the
given values
of
(p
on L.
Remark
Let us
consider
a
discrete model
A4
n
(£l,
N,
<p,C}.
If C
contains
no
cross con-
straint (see page
9),
then
it can be
verified
that
where
R*((f>
v
]
is the
individual roughness associated with
the
1-dimensional
model
A
/
(
1
(0,
JV,
(p
v
,C}.
In
this
case, minimizing
R*(tp)
is
equivalent
to
min-
imizing
separately each individual roughness
R*((p
v
\.
In
other
words,
from
the DSI
point
of
view,
if
there
is no
cross
constraint,
then
it is
equivalent
to
•
either globally interpolate
the
function
(p
associated
with
the
n-dimensional
model
M
n
($l,N,ip,C),
• or
individually interpolate each component
</?"
associated
with each 1-dimensional
model
M
l
(Sl,N,ip",C).
4.2.2
Genericity
At
first
glance,
the DSI
method seems
to
depend
on
arbitrary parameters
{v
v
(a,0}}
and
{(j,(a)}.
This
is
rather
confusing
for
someone
not
familiar
with
this method.
For the
sake
of
simplicity, these parameters could have been
fixed, but
this
would have hidden
the
genericity
of
this
method:
depending
on
the
actual choice
of
these parameters,
different
types
of
problems
can be
solved
and
different
families
of
solutions
can be
obtained.
There
is a
wide range
of
possible choices
for
these parameters. However,
we
will
see in
section (4.3)
that
the
uniqueness
of the DSI
solution
can be
guaranteed only when taking into account some slight
limitations
on the co-
efficients
{v"(a,
(3]}
and
/j,(a).
For
those
who are
curious
to
know
more about
these parameters, some examples
of
possible choices
are
presented below; other
examples
will
be
presented
on
pages 377,
and
366.
4.2.
THE
DSI
PROBLEM
149
Choosing
{v"(a,/?)}:
A
preliminary
remark
Most
of the
time,
for
simplicity's sake,
the
weightings
{v
v
(ot,fl}}
are
chosen
according
to the
following
model:
In
this expression,
the
coefficients
{r"}
are
assumed
to be
positive
and are
chosen
to
balance
the
different
orders
of
magnitude
of the
variations
of the
components
{(/^(a)},
while
the
weightings
{v(a,f3}}
define
the
"style"
of the
interpolation.
In
practice,
the
variations
of the
components
{(p
v
(ot)}
of
(p(a)
are
often
assumed
to
have approximately
the
same order
of
magnitude.
In
this
particular case,
it is
wise
to
choose
all the
coefficients
{r"}
equal
to 1.
The following
sections propose
a
couple
of
examples
of
coefficients
{v(a,
/3)}
honoring
the
conditions
of
theorem
3
page 157:
Choosing
{v(a,(3}}:
Harmonic
weightings
Let
|7V
0
(a:)|
be the
number
of
elements
of the
orbit
N°(a).
Harmonic weight-
ings
are
said
to be
used
if the
coefficients
{v"(a,0)
=
v(a,
/?)}
are
chosen
according
to the following
definition:
Harmonic
functions
have
the
characteristic property
of
being equal,
at any
point,
to
their
own
mean
on a
circle centered
on
this point
[87].
By
using
the
coefficients
{v(a,/3}}
defined
above,
we can
check
that
R((p\a]
— 0 if
(p(oi)
is
equal
to the
mean values
of the
</?(/?)
surrounding
the
node
a. For
this reason
and by
analogy with harmonic
functions,
we
propose
to
call these weighting
coefficients
harmonic.
Harmonic weightings
are
currently
the
most frequently used
in
practi-
cal
applications.
For
example,
the
interpolations presented
in figures
(1.11),
(1.12),
and
(1.13)
were obtained with such weightings. According
to the re-
mark
on
page
16, it can be
observed
that
the DSI
solutions
so
obtained
are
close
to
"Spline"
functions
(see also page
364).
Choosing
{v(a,/?)}:
Isotropic weightings
Most
of the
time,
the
design
of
natural
objects requires
the use of a d-
dimensional graph
Q([l,N)
with
the
following
characteristics:
• The
nodes
of
£/(f2,
N)
are
embedded
in the 3D
space
IR
3
.
• The
edges
of
G(£l,
N)
generate
a
mesh composed
of
adjacent
d-simplexes.
1
For
example:
—
if d — 1,
then
C/(J7,
N)
consists
of a
polygonal
curve,
—
if d — 2,
then
(?(fi,
N)
consists
of a
triangulated
surface,
and
1
A
1-simplex
is a
segment,
a
2-simplex
is a
triangle,
and a
3-simplex
is a
tetrahedron.
150
CHAPTER
4.
DISCRETE SMOOTH INTERPOLATION
Figure
4.2
Building
isotropic
weightings
in the
case
of
a
triangulated
surface.
—
if d = 3,
then
G(Q,
N)
consists
of a
tetrahedralized
volume.
If
the
distance
d(a,/3)
between
two
adjacent nodes
(a,/3)
is,
locally, approxi-
mately constant, then
(?(fi,
N) is
said
to be
isotropic
and the
harmonic weight-
ings
{^(ct,/3)}
yield excellent
results.
However,
as
suggested
in figure
(4.2),
if
(}(£l,N)
is
anisotropic, then
the
harmonic weightings
may
produce
a
DSI
solution inherited
from
the
same anisotropy.
For
this reason,
in
such cases
where
(7(f2,
N}
is
strongly anisotropic,
we
propose
to use the
isotropic weight-
ing
defined
by the
following
coefficients
{v(a,P}}:
In the
above expression,
t(a)
and
t(a,/3)
are
defined
as
follows
with
T(a)
= set of
simplexes sharing
the
node
a
T(ot.,
/3)
= set of
simplexes sharing
the
node
(a,
/3)
|T|
=
measure
2
of the
simplex
T
Choosing
{v(a,(3)}i
Anisotropic
weightings
Let us
assume
that
C/(Q,
N}
is a
given
ID,
2D or 3D
graph embedded
in the 3D
Euclidean space
JR
3
and let
^p
v
(a)
be a
physical parameter
to be
interpolated
with
DSI on
£(fi,
N}.
The
distribution
of
^(ot)
may be
anisotropic
in
7R
3
.
2
|T|
represents
the
length
of a
segment
T, the
area
of a
triangle
T, or the
volume
of a
tetrahedron
T.
4.2.
THE
DSI
PROBLEM
151
We
assume
that,
at
each node
a of
G(£l,
A
7
"),
the
location
x(a)
of the
node
is
given.
In
this
case,
we
propose using
the
anisotropic weightings defined
by
the
following
choice
of the
coefficients
{v
v
(a,0)}:
In
this expression,
•
N°(a)
is the
orbit
of the
node
a
(see
definition
1.3);
• e
G]0,1[
is
assumed
to be a
given
small,
positive
real
constant
3
chosen
to
prevent
v(a,/3)
from
vanishing;
•
(v(o;)
: a
G
£7} is a
given
field of
vectors
controlling
the
orientation
of the
anisotropy;
•
(j)(
ct,/?|x,v
) is
equal
to the
cosine
of the
angle
between
(x(/3)
—
x(a))
and
v(a)
in the
case
of an
"axial
anisotropy"
aligned
with
the
direction
v(a);
•
0(
a,
/?|x,
v ) is
equal
to the
sine
of the
angle
between
(x(/9)
—
x(a))
and
v(a)
in
the
case
of a
"planar
anisotropy"
orthogonal
to the
direction
v(a);
and
• p
>
0 is a
given
parameter
chosen
to
control
the
intensity
of the
anisotropy:
at
the
limit,
ifp
= 0, the
anisotropic
weightings
become
identical
to the
harmonic
weightings.
For
example,
in a
stack
of
geological layers,
if
(v(a)
: a €
17}
represents
the
field
of
vectors normal
to the
horizons, then these vectors
can be
used
to
define
planar
anisotropic weightings, which
in
turn
can be
used
to
interpolate
on a
3-grid
a
physical property
(e.g.,
a
seismic velocity) whose gradient
is
more
or
less orthogonal
to the
horizons.
An
example
of the
application
of
axial
anisotropic weightings
will
be
presented
on
page 377.
Choosing
{//(a)}
The
stiffness
function
{^(a)}
has
been introduced
to
modulate
the
smoothness
of
the DSI
solution locally (see equation
4.8
page 142).
The
only constraint
is
that
this
weighting function must
be
positive (see
section
4.3).
There
is a
wide
range
of
possible choices
and the
simplest
one
consists
in
using
a
uniform
weighting:
H(a)
= 1 V a e fi
A
more interesting choice
can be
proposed
in the
frequent
case where
the
coefficients
{v
v
(a,(3}
=
v(a,/3}}
are
independent
of v. In
this case,
let us
define
cr(a),
(f>(a)
and
(j,(a)
as
follows:
3
To
honor
the
second
constraint
(4.35),
the
value
of
e
must
be
chosen strictly positive.
For
example,
one can
choose
e
—
0.1.
152
CHAPTER
4.
DISCRETE SMOOTH INTERPOLATION
where
/t(o:)
is a
given positive weighting
function.
According
to the
comment
presented
on
page 143,
the
global roughness
R((f>]
appears
as the sum of the
squares
of the
discrepancies between
the
values
(p(a)
attached
to the
nodes
of
the
graph
£7(fJ,
N)
and the
mean values
(p(a)
of
their neighbors:
In the
above expression,
the
/}(o:)
coefficients
are
considered
"normalized
stiffness"
coefficients,
which must
be
positive.
In
practice,
the
simplest choice
consists
in
selecting
p,(a)
constant:
4.2.3
The
DSI
solution
Notations
The DSI
problem does
not
depend
on the
order
of the
values
{(p
v
(a)}
in the
column
matrix
(p.
Therefore,
to
simplify
notations,
we
will
assume
that
a
permutation
of the
elements
of the
matrix
<p
has
been performed
in
such
a
way
that
(p
can be
split into submatrices
4
<pj
and
(p
L
such
that
This
partition
of the
matrix
</?
induces
a
similar partition
of the
matrices
[W*
and
Q*
used
to
define
R*(p)
(see equation
(4.24)):
4
See
definition
of / and
L
on
page 147.
With these notations,
we can
express
R*((p]
as a
quadratic
function
of
<£/:
4.3.
UNIQUENESS
OF THE
DSI
SOLUTION
153
The
global
DSI
equation
As
suggested
in figure
(1.9),
the
graph
of
R*(<p)
as a
function
of
</?/
is a
paraboloid whose minimum corresponds
to the
solution
of our
problem.
This
minimum
is
reached when
the
vectorial derivative
dR
k
((p}/d(pi
is
equal
to the
null
vector. According
to
equation (4.32)
and to the
derivation rules
(4.26),
we
have
We
conclude
that
the
solutions
to our
problem
are
characterized
by the
fol-
lowing
matrix equation:
We
propose calling
this
equation
the
"global
DSI
equation."
As
will
be
seen
in
section (4.46),
another
formulation called
the
local
DSI
equation
allows
the
same solution
to be
obtained iteratively (provided
that
this solution
is
unique).
4.3
Uniqueness
of the DSI
solution
It
appears
that
looking
for a
solution
to the DSI
problem
is
nothing more
than minimizing
a
quadratic objective
function
-R*(v?)
that
leads
to the
linear
system (4.34).
The
rank
of the
square matrix
[Wfj]
defining
this
linear system
is
the
only essential
in
proving
the
uniqueness
of the
solution.
In
this section
we
show that, under
a set of
very large
sufficient
conditions
that
are
generally
honored
in
practice,
[Wj
f
]
is
invertible
and the DSI
solution
is
unique.
4.3.1
Consistency
Consistency
of the set L
The set L of
nodes where
ip
is
known
is
said
to be
DSI-consistent
or
more
simply
"consistent"
if
each connected component
of the
graph
(?(fi,
N)
con-
tains
at
least
one
node belonging
to L. In
other words,
L is
consistent
if, for
any
node
i
G
/,
there
is at
least
one
node
i G L
such
that
there
is a
path
\(i,t)
e
G(ft,N)
joining
i to L
154
CHAPTER
4.
DISCRETE SMOOTH INTERPOLATION
Consistency
of the
weightings
{v"(a,/3}}
and
/j,(a)
The
weightings
{v
v
(a,(3}}
and
n(a)
are
said
to be
DSI-consistent
or
more
simply
"consistent"
if, for any a G
0,
these
coefficients
are
such
that
Remark
Some
examples
of
consistent choices
for the
weighting
coefficients
{v
v
(a,fi}}
honoring conditions
(2) and (3) of the
above
definition
have already been
presented
in
section
(4.2.2).
4.3.2 Theorem 1
If
L is
consistent
and if the
weightings
{v
v
(a,fi}}
and
//(a)
are
consistent,
then each
of the
associated symmetrical matrices
{[VF
Z/
]
: v =
1,...,
n} of
size
(M
x M)
defined
by
equations (4.11),
(4.12),
and
(4.13)
has
only
one
eigen
value,
let's say,
AM,
equal
to
zero
and the
associated unit
M-dimensional
eigen vector
UM
has its M
components
all
equal
to
l/\/M:
Proof
For
simplicity's sake
and
without loss
of
generality,
it is
assumed
in the
fol-
lowing
proof
that
n = 1
and, therefore,
• no
distinction
will
be
made
between
</?"
and
ip,
• no
distinction
will
be
made
between
[W]
and
[W"],
and
• no
distinction
will
be
made
between
v(a,/3)
and
v
l
'(a,f3)
By
definition,
we
have
Let
A
Q
(a)
be the
coefficients,
such
that
Accounting
for the
fact
that
R(<p\oi)
is
nonnegative,
if the
fJ,(a)
are
consistent,
we
deduce that
[W
v
]
is
nonnegative also,
and we can
write
4.3.
UNIQUENESS
OF THE
DSI
SOLUTION
Taking into account
the
fact
that
the
weightings
{v(a,/3}}
are
also consistent,
we
are
sure
that
the
coefficients
{X@(a)}
have
the
following
property
for any
a
e
fh
This implies
that
we
have:
According
to the
definition
of the
notion
of
convex hull
and due to the
fact
that
the
coefficients
{X
/3
(a}}
are
strictly positive, (4.37) implies
that
(p(ci)
must honor
the
following
inequalities:
Let
a
e
0 be an
index corresponding
to the
minimum
of
</?
on
D:
Let
us now
consider
the
orbit
7V°(o!)
of a.
According
to
(4.38),
we can
thus
write
The
above property
and the
definition
(4.39)
of a
imply
that
Replacing
a by any
index
a
belonging
to
A
r
°(ct),
it can be
shown
in a
similar
way
that
If
the
graph
£?(0,
N)
has
only
one
connected component, then
all the
nodes
a
€
Q
are
connected,
and the
above proof
can be
recursively repeated
for
each
node
a of
this graph.
This
leads
to the
conclusion
that
According
to
equation (4.36),
this
proves
that,
if the
graph
C/(fi,
N)
has
only
one
component, then
[W]
=
\W"}
is
semipositive
definite
and has
only
one
eigen
value,
let's
say
AM,
equal
to
zero,
and the
associated
unit eigen vector
UM
is
such
that
5
5
Note
that,
according
to
(1.2),
M
is
equal
to the
number
|fi|
of
nodes
in
Q.
155
156
CHAPTER
4.
DISCRETE SMOOTH INTERPOLATION
Corollary
1
If
L is
consistent
and if the
weightings
{v
v
(a,(3}}
and
//(a)
are
consistent,
then
the
associated matrices
{[W"]
: v =
1,...,
n}
defined
by
equations
(4.11),
(4.12),
and
(4.13)
are
semipositive
definite,
and
each
of
them
has a
rank equal
to
(M-l).
4.3.3
Theorem
2
If
L is
consistent
and if the
weightings
{^^(a,/?)}
and
p,(a)
are
consistent,
then
any
square submatrix deduced
from
\W
V
\
by
removing symmetrically
one row
together with
the
corresponding column
is of
full
rank
and is
thus
invertible.
Proof
With
the
notations
of
theorem
1
above,
we
have:
Moreover,
we
know
that
[W]
=
\W
V
]
has
only
one
eigen value,
let's
say
AM,
equal
to
zero
and the
associated unit eigen vector
UM
is
such
that
If
we
consider
the
full
set of
eigen values
{Ai,...,
AM}
and the
associated unit
eigen
vectors
{t/i,...,
UM}
of
[W], then,
due to the
fact
that
[W] is
both
symmetrical
and
semipositive
definite,
we are
sure
that
According
to the
spectral decomposition
of
symmetrical semipositive
definite
matrices [15]
and
taking
into
account
the
fact
that
AM
=
0, we can
write
We
deduce
that
the
vectors
(Vi,...,
VM-I}
are
linearly independent
in
iR
M
,
which
implies
that
[V] is an
M
x
(M
—
1)
full
rank matrix:
Let us now
consider
the
columns
of
[V]*:
4.3. UNIQUENESS
OF THE
DSI
SOLUTION
Each
of
these columns
v
a
is a
vector
of
JR
-1
,
and it is
easy
to see
that
[W]
is
the
Gram matrix
of the set
{vi,...,
VM}-
Note
that,
under
the
consistency conditions
of
{v"(a,0}}
and
//(a)
(see equa-
tion
(4.35)),
this
equation implies
that
the
vectors
{v
a
}
are
non-null
because:
In
other words,
the following (M — 1) x
(M
—
1)
linear
system
has
only
one
solution:
This implies both
that
the
matrix
[V]*,
deduced
from
[V]*
by
removing
the
last
column
VM,
is
invertible
and
also
that
the
matrix
[W]
=
[W"],
deduced
from
[W]
=
\W
V
]
by
removing
the
last line
and the
last
column,
is
invertible
because
More generally,
we
conclude
that
any
submatrix
\W
V
]
deduced
from
[W"]
by
removing
symmetrically
one row and the
corresponding column
is
invertible.
4.3.4 Theorem
3
If
L is
consistent
and if the
weightings
{v"(a,0)}
and
(JL(O)
are
consistent,
then
the DSI
equation
has
only
one
solution.
Proof
Let us
consider
the
partition
of
[W*]
defined
by
equations
(4.30),
and
(4.31),
which induces
a
similar
partition
of W. It is
always possible
to
permute
Bearing
in
mind
that,
for i <
M,
the
vector
UM
is
orthogonal
to any
vector
V^,
we
can
easily
verify
the
following
relationship where
OM-I
is the
null vector
of-R^-
1
:
Taking into account
the
fact
that
[V]
has a
rank equal
to (M
—
1),
we
deduce
that
the
above relation
is the
only
linear dependency between
the
vectors
{vi,..
.,VM},
consequently,
the
non-null
vector
VM
can be
decomposed
in a
unique
way on
{t>i,...,
VM-I}'-
157