Jean-Laurent Mallet. Geomodeling
Подождите немного. Документ загружается.
158
CHAPTER
4.
DISCRETE SMOOTH INTERPOLATION
symmetrically both rows
and
columns
of
W//
to
write this matrix
in the
following form:
where
[Wjj]
is
deduced
from
\W
V
]
by
removing symmetrically
the
rows
and
columns
I
6
L.
Such
a
permutation does
not
change
the
rank
of
[W//].
Taking into account
theorem
2
above,
we can
observe
that
where
|/| is the
number
of
elements
of the set
I.
Taking equation
(4.24)
into
account,
we can
write
Observing
that
[W//]
is an
(n-1/|
x
n-1/|)
symmetrical matrix,
from
equation
(4.40),
it can be
concluded
that
[W/
7
]
is a
full
rank matrix
and is
thus invert-
ible.
As a
direct consequence,
the
global
DSI
equation
(4.34)
has
only
one
solution: this guarantees
the
uniqueness
of the
solution
of the DSI
problem
whatever
the
constraints
may be.
4.3.5 Factorization
of the DSI
solution
Let us
consider
a
1-dimensional
model
M.
l
(Q,
TV,
</?,
C]
and let us
assume
that
the two
following
conditions
are
satisfied:
• The set L of
Control-Nodes
is
consistent (see page 153).
• C
contains
only
the
hard
constraints
C
L
(see
page
14).
Under
these
conditions,
it is
possible
to
show
that
the DSI
solution
</?/
can be
decomposed
on the
basis
of \L\
blending functions
{B\}
independent
of the
function
</r.
Blending
functions
B\
According
to the
global
DSI
equation
(4.34),
it is
straightforward
to
verify
that,
under
the
conditions
(4.41),
the
solution
tpi
can
always
be
written
in
the
following
way:
4.3.
UNIQUENESS
OF THE
DSI
SOLUTION
159
The
matrix
B^
of
size
(|/|
x
\L\)
can be
decomposed into
a
series
of \L\
column
matrices Bj, each of them having |/| rows:
This allows
us to
write
the DSI
solution
as
The
B\
functions play
a
role similar
to the
"Blending Functions" used
in
Bezier
or
Spline
functions
(e.g.,
[77]),
and it
will
be
shown
that
they have
the
following
similar property
if
there
is no
constraint other
than
C
L
and if L and
the
weightings
{v(a,0}}
and
n(oi)
are
consistent:
However,
contrary
to the
Bezier
and
Spline functions,
the DSI
solution
is not
an
approximation,
but an
interpolation method implying
that
the
B
i
(i]
values
does
not
necessarily belong
to
[0,1}.
Computing
the
blending
functions
B\
Let
A(
be the
"unit
impulse," defined
as
follows:
With
these
notations,
and
under
the
conditions
(4.41),
it is
obvious
that
we
have:
In
other words,
B\
is the
response
of the
"DSI
filter" to the
"unit impulse"
A£.
The
property
(4.42)
provides
an
easy
way to
compute
the
Bj
blending
functions
thanks
to
DSI.
Property
of the
blending
functions
B\
If
L and the
weightings
{v"(a,fi)}
and
//(a)
are
consistent (see page 153),
then, under
the
conditions (4.41),
we
have:
160
CHAPTER
4.
DISCRETE SMOOTH INTERPOLATION
Proof
If
there
is no
constraint other than
C
L
and if L and the
weightings
{v(a,
/3)}
and
/A(a)
are
consistent, then
the
conditions
of
theorem
1 on
page
154 are
satisfied
and we
have:
We
deduce
from
this relation
that
Moreover,
according
to
theorem
2
page
156,
we are
sure
that
the
inverse
W
n
l
exists because
[W] is of
rank
(M
—
1), and we can
write
4.4
The
local
DSI
equation
In
this
section,
instead
of
solving
the DSI
matrix
equation (4.34) directly,
we
propose
an
iterative approach making
it
unnecessary
to
compute
and
store
[W*].
The
major
advantage
of
such
an
approach
is
that
it
provides
an
easy
way
to
account
for
hard
and
dynamic
constraints
4.4.1
Computing
the
variations
of
R*((p)
Variation
of
R((p)
The
following
can be
deduced
from
the
definition
of
R((p\k):
This implies
that
4.4.
THE
LOCAL
DSI
EQUATION
161
We
conclude
from
this
relation
that
we
have:
Variation
of
p((p)
For any
soft
constraint
of c €
C~,
we can
observe
that
Using
the
vectorial derivation rules
(4.26),
we
obtain
This
implies
that
d\A
t
c
•
</?
—
b
c
\
2
1d(p
v
(a)
can be
written
as
From this relation
and
equation
(4.16),
we
deduce
that
with
:
6
Equations
(4.34),
(4.46),
(4.49),
and
(4.50)
are
considered
to be the
backbone
of
this
book
and,
for
this reason, have been
put
into boxes. Dozens
of
DSI
constraints
are
pre-
sented
throughout this
book,
and,
for
each
of
them,
the aim is to
determine
the
coefficients
{A^(Q)},6
c>
7?(a),andr^a|^).
162
CHAPTER
4.
DISCRETE SMOOTH INTERPOLATION
According
to
equation
(4.17),
we
conclude
that
Variation
of
R*(p}
Taking into account
the
definition
of
R*(y>)
(see page 145)
and
according
to
equations (4.43)
and
(4.44),
it is
easy
to see
that
4.4.2
The
local
DSI
equation
The DSI
solution
<p
is
such
that
dR*((f>)/d(f>"(a)
= 0;
hence
the
component
(p
v
(oL)
of
(f>
must satisfy
the
following
equation:
6
where
the
terms
and
are
defined
by
4.4.
THE
LOCAL
DSI
EQUATION
163
with:
We
propose calling equation
(4.46)
the
"local
DSI
equation"
at
node
a.
Cross
constraints
It is
interesting
to
note
that,
in
equation
(4.49),
the
term
x^(a\ip}
is
different
from
0 if, and
only
if, c is a
cross
constraint
(see page
9).
Therefore,
we
will
say
that
x^.(a\(p)
is the
"cross constraint component"
of c on
(p
at
node
a.
Remark
1
We
can
observe
that
the
term
T^(a\(p)
can
also
be
written according
to the
following
equation:
In the
case where
A^.
(a) is not
equal
to
zero,
this equation
can
also
be
written
as
Remark
2
Taking
into
account
the
constraints
normalization (see equation
(4.4))
as
well
as the
values
g
v
(oi)
and
7^
(a)
denned
above,
it is
easy
to
check
that
the
square
matrices
[W] and
[A
2
]
denned
by
(4.12),
and
(4.18)
are
such
that
164
CHAPTER
4.
DISCRETE SMOOTH INTERPOLATION
We
conclude
that
the
balancing ratio
w
defined
by
equation
(4.25)
can
also
be
computed according
to the
following
formula:
Taking into account equation
(4.4),
we
also have:
Example: Harmonic weightings
By
developing
the
local
DSI
equation (4.46) with
the
harmonic weightings
v(a,
/#),
as
denned
by
equation
(4.29),
the
following
expressions
for the
terms
7
G
v
(ct\(p)
and
g
v
(a)
are
obtained:
Example:
Fuzzy
Control-Node constraint
Let us
consider
a
discrete model
Ai
n
(fi,
AT,
</?,
C) and let a
G
/ be a
given node
for
which
the
unknown component
(p
v
(ot)
of
(p(oi)
is
known
to be
close
to a
given
value
</>£:
Such
a
piece
of
information
can
easily
be
turned into
a DSI
constraint
c
called
"fuzzy
Control-Node" belonging
to
C~
and
denned
by
{c
honored}
with
:
7
In
these
expressions,
|A
r
°(o:)|
represents
the
number
of
elements
of the
orbit
TV
0
(a) of
the
node
a
(see
figure
(1.3)).
Note
that,
in the
case
of
harmonic
weightings,
the
term
g
v
(a)
does
not,
in
fact,
depend
on v.
4.4.
THE
LOCAL
DSI
EQUATION
165
This
DSI
constraint
is
certainly
the
simplest
imaginable,
and it is
easy
to
build
the
associated terms
T^.(o\(p)
and
7^(0:)
as
defined
by
equation
(4.49):
4.4.3
Proposal
for an
iterative algorithm
The
local
DSI
equation presented
in the
previous section suggests
the
following
iterative algorithm
in
estimating
the
solution
(p:
II-
DSI
algorithm
let / be the set of
nodes where
<p(a)
is
unknown
let
(p
=
<£>[o]
be a
given initial approximated solution
let
w
be an
estimated value
of the
balancing
factor
while
(
more iterations
are
needed
) {
} // end :
for_all(
a
G
/ )
} // end :
while(
more
iterations
are
needed
)
As
you can
see, this very simple algorithm does
not
explicitly
use the
matrix
\W*\
occurring
in the DSI
equation. However,
it
must continually recompute
products
{p,((3}
•
v
v
(alphai,0)
•
v"(a2,/3)},
which
are
precisely
the
ones used
to
build
[W*]
(see equation
4.12).
If
the
initial solution
y?[
0
]
is
close
to the
actual
DSI
solution, then
few
itera-
tions
are
needed
and the
computation cost becomes acceptable.
For
example,
such
a
situation occurs
in an
interactive context where
the
initial solution
is
taken
as
equal
to the
previous solution
before
any
modifications
of the
constraints
and/or
the
values
{<£>"(•£)
:
i
G L} are
performed
by the
user.
To
show
that
the
proposed algorithm actually converges, notice
that
R*(<p)
can be
expressed
as a
function
of
ip
v
(a)
according
to the
following formula:
The
coefficients
A,
B and C are
independent
of
(p
v
(oi)
and
The
minimum
of the
global criterion
R*(y>)
as a
function
of
<p"(a)
is
achieved
for
the
value
(p
v
(oi)
=
—B/A,
which,
by the
way, corresponds
to the
value
given
by the
equation
(4.46).
Thus,
we are
able
to
conclude
that
our
algorithm
converges since,
at
each
step
of the
iterative
process,
the
value
of the
positive
or
null
function
R*(tp]
decreases.
166
CHAPTER
4.
DISCRETE SMOOTH INTERPOLATION
It
must
be
noted
that,
if the
conditions
of the
theorem
3
presented
on
page
157 are
honored, then
the final
solution
(p
is
unique:
the
initial solution
(/?[QJ,
while
leaving
the
result unchanged, influences only
the
speed
of
convergence.
If
the
conditions
of
this theorem
are not
satisfied, then
the
iterative algorithm
will
yield
a final
solution
that
depends
on the
initial
solution.
Remark
If we consider a node a with no constraint attached, then the terms Tv(a\(p)
an d
7"
(a)
are
both equal
to
zero
and the
local
DSI
equation
(4.46)
becomes:
{ no
constraint
attached
to a }
In
such
a
case,
the
shape
of
this equation, which tends
to
suggest
that
(p"(a}
is
free
from
any
constraint,
is
quite
confusing. However,
this
is not
true.
It
is
important
to
understand
that,
recursively through
the
graph
Q(£t,
N),
this
value
y>
v
(oi)
depends
on
values
{<£>"(/?)}
where
constraints
are
applied,
so
that
(p
v
(a)
is
also dependent
on
these constraints.
4.4.4 Dynamic constraints
Prom
a
theoretical point
of
view,
the
constraints
c 6 C
should always remain
unchanged
during
the
iterative process presented
in
section
(4.4.3).
How-
ever,
in
practice, dynamically updating these constraints
at
each step
of the
iterative process
can be
extremely
useful.
For
example:
• If a
constraint
c is not
linear,
then
it is
always
possible
to
linearize
it in
the
neighborhood
of the
current
solution
(see
section
(4.4.5)).
This
current
solution
evolves
at
each
iteration
making
it
necessary
to
update
such
a
lin-
earization.
•
Consider
the
example
in figure
(1.11),
where
(f
represents
the
geometry
of a
surface
H
that
has to fit a set of
points
scattered
in the 3D
space.
Each
data
point
is
assumed
to
attract
the
impact
point
corresponding
to its
projection
onto
the
surface
in a
given
direction.
It is
clear
that,
keeping
the
directions
of
attraction
constant,
the
impact
points
change
during
the
iteration
process,
thus
it is
wise
to
update
them:
this,
in
turn,
implies
that
the
constraints
modeling
these "attractions"
need
to be
updated.
To
permit
the
updating
of
such "dynamic constraints,"
we
propose
modifying
the
iterative algorithm presented
in the
previous section
as
follows:
//—
DSI
algorithm
with
DSI
constraints
let
/
be the set of
nodes
where
<p(ot)
is
unknown
let
<p
=
tp[Q]
be a
given
initial
approximated
solution
let
w
be an
estimated
value
of the
balancing
factor
while
(
more
iterations
are
needed
) {
for_all(
a
€
/ ) {
for_all(
i/
) {
4.4.
THE
LOCAL
DSI
EQUATION
167
Figure
4.3
Smoothing
a 3D
surface
with DSI:
the
initial
surface
is a
sphere
whose
geometry
has
been
randomly
perturbed.
}
II end :
for_all(
a
€
/ )
} // end :
while(
more
iterations
are
needed
)
It
should
be
noted
that,
from
a
theoretical point
of
view,
the
uniqueness
of the
solution
can no
longer
be
guaranteed when such dynamic constraints
are
used.
In
general,
the final
solution depends
both
on the
initial solution
</?[o]
and the
ordering
of the
nodes
of
O
used
in the
main loop
of the DSI
algorithm:
However,
in
practice,
if the
modifications
of the
constraints between
two
suc-
cessive
iterations
are
small, then
the
convergence
is
still observed,
and the
final
solution
is
weakly dependent
on the
initial solution.
Dynamic
constraints appear
to
provide
a
very
efficient
way of
solving prob-
lems
with complex constraints.
The
global approach consisting
of
solving
the
matrix equations
(4.34)
does
not
allow such dynamic constraints
to be
taken
into account: accordingly,
we
consider
the
iterative approach based
on the
above algorithm
as the
only
one to be
used
in
practical applications.
Example
1:
Smoothing
Let us
assume
that
we
want
to
"smooth"
the
component
(p
v
of a
function
(p
belonging
to a
discrete model
^"(OjTV,
</9,C).
More precisely,
we
would like
to
replace
the
current values
{^^(a)
: a 6 fi} by new
values
{ip
v
(a)
:
a
<E
$1}
in
such
a way
that
}
for_all(
c 6 C ) {
if(
c is
dynamic
) •
update
c
}
initial
solution
c^[
ordering of
one particular solution