Cubic L1 Splines 170

1
Theory and Applications ofShape-preserving Cubic
L1 Splines
Shu-Cherng Fang North Carolina State
University Raleigh, NC, USA
September 20, 2006 at National Tsing Hua
University
Coauthors Hao Cheng, John E. Lavery, Yong Wang,
Wei Zhang, Yumin Lin, Yunbin Zhao
2
Outline
? Splines and shape-preservation
? Cubic L1 interpolating/smoothing splines

• Second-derivative based

• Univariate

• Bivariate

• Rectangular Sibson Elements
• Triangular Irregular Nets

• First-derivative based

• ? Mathematical models and computational results
• ? Domain decomposition for large scale
  applications

3
Spline Curves

• Spline functions provide a convenient way for
  drawing curves in 2- or
• 3-dimensional spaces.

• The term comes from drafting, where splines were
  flexible strips guided
• by points on a paper, used to draw curves.

• Spline functions possess many nice structural
  properties and
• excellent approximation powers.

• Commonly used for interpolation or approximation
  in real-life applications.

4
Classification

• Based on the number of variables, we have

• Univariate splines

• Bivariate splines

• Higher dimensional splines

• Based on the type of functions used, we have

• Polynomial-based splines

• B-splines

• Bezier curves

• L2 splines,

• Nonploynomial-based splines

• Exponential Splines
• Trigonometric splines
• Rational Splines (NURBS)

5
Shape-preservation

• Splines are useful for real-life applications,
  such as

• Terrain description

• Recognition of objects

• Virtual space simulation

• Fast 2-D and 3-D zooming

• Haptic devices design

• One fundamental requirement is that splines
  should be
• shape preserving.

• No universal, standard definition of
  shape-preservation,
• but no extraneous oscillation is preferred.

6
Illustration 1
7
Illustration 2
8
Illustration - 3
Mesh representing 128-point data set
9
Illustration 3 cont
10
Illustration - 4
Mesh representing 128-point data set
11
Illustration - 4 cont
12
Shape-preservation

• In general, shape-preservation means no
  nonphysical or extraneous
• oscillations.

• From geometric point of view, shape-preservation
  means the resulting
• curves retain geometric properties of the
  initial data, such as positivity,
• monotonicity, convexity, linear and planar
  section.

• Observation over the past 30 years indicates
  commonly used polynomial
• splines have inadequate shape-preserving
  capability.

• In particular, they cannot preserve shape well
  for arbitrary data
• with arbitrary changes in magnitude and in
  knot spacing.

13
Common Approaches

• Local interpolation procedure

• Piecewise linear interpolation

• Brodlie-Carlson-Fritsch-Butland interpolation

• Global minimization principle

• B splines

• L2 splines

• Hybrid fairing procedure

14
Univariate Cubic Lp Splines

• Given a strictly monotonic partition of a finite
  real interval x0, xI

?? lt x0 lt x1lt ??? lt xI?1 lt xI lt ?,
a cubic Lp spline z of the data
is a C1 smooth piecewise cubic
polynomial that minimizes
where 1? p lt ?.

• When p equals 1 or 2, z(x) is called a (second
  derivative
• based) univariate L1 or L2 spline,
  respectively.

• L2 splines are commonly used because they are
  easy to calculate.

• These are second derivative based.

15
Univariate Cubic L1 Spine

• z(x) is piecewise cubic such that on each
  subinterval xi, xi1, there is a cubic
• polynomial zi(x) satisfying z(x) zi(x),
  ?x?xi, xi1, i 0,,I-1.

• For each i, zi(x) can be uniquely expressed as

• Let hi xi1?xi and ?zi (zi1 ? zi) / hi, i
  0,,I ?1. Change variable from x to
• t (x ? xi ? hi/2) / hi, -½ ? t ? ½.

• Finding a 2nd derivative based univariate L1
  spline is equivalent to solving
• a nonsmooth optimization problem

s.t.
16
Geometric Programming Model
(Primal)
where
Nondifferentiable !
X cAc 0
17
Geometric Programming(GP)

• Primal geometric programming

(Primal)
where X is a cone and C is a convex subset in Rn.

• Dual geometric programming

(Dual)
where
18
Conjugate Transform

• Given a function w(z) over a domain W?Rn, the
  conjugate
• transform of w(z) is a function ?(?) over a
  domain ? ?Rn, where

and
19
Optimality Conditions

• x and y are optimal solutions to the primal
  problem and dual
• problem, respectively, if and only if

(I)
(II)
(III)
20
Dual GP

• The dual problem is a convex program with a
  linear objective function and
• quadratic constraints.

(Dual)
where
Y is the row space of A.
is a convex set
where
21
Dual to Primal Transformation

• Optimality conditions indicate the dual optimal
  solution can be
• transformed to a primal optimal solution by
  solving an LP

(T)
Where ,
if the dual solution
is on the boundary of ?i,,
otherwise,
if
is on the boundary
and
of ?i, otherwise,
22
Univariate Results - Theory

• Theoretical proof of shape preserving
  properties.
• GP approach allows us to prove that under
  various circumstances, cubic L1
• splines preserve linearity and
  convexity/concavity of data. In particular,
• (i) when four consecutive data points lie on
  a straight line, the cubic L1 spline is
• linear in the interval between the second and
  third data points
• (ii) cubic L1 splines of convex/concave data
  preserve convexity/concavity if the
• divided differences of the data do not
  increase/decrease too rapidly
• (iii) when cubic L1 splines do not preserve
  convexity/concavity, they do not
• cross the piecewise linear interpolant and,
  therefore, they do not have extraneous
• oscillation.
• References Cheng, H., Fang, S.-C., and Lavery,
  J. (2002, 2005a, 2005b)

23
Univariate Results Theory cont
24
Univariate Results - Algorithm

• Development of a computationally efficient
  active set algorithm.
• The algorithm requires only simple
  algebraic operations to obtain
• an exact optimal solution in a finite number
  of iterations.
• Takes milliseconds to solve problems with a
  thousand knots on Pentium-III
• 1GHz PC.
• In both stability and efficiency, it
  outperforms a currently used
• discretization-based algorithm.
• Reference Cheng, H., Fang, S.-C., and Lavery,
  J. (2004)

25
Illustration - 1
2nd derivative based L1 univariate Spline
2nd derivative based L2 univariate Spline
26
Illustration 2
L1 univariate Spline
Brodlie-Fritsch-Butland interpolant
27
Univariate Results - Applications

• Application to term structure analysis of bond
  market.
• Reference Chiu, N.-C., Fang, S.-C., Lin, J.-Y.,
  Wang, Y., and Lavery, J. (2005)

28
Bivariate Cubic L1 Splinesover Tensor Product
Grid

• The ordered sets

and
,
form a partition (tensor-product grid) ? over
x0, xIy0, yJ.

• zij is given at each knot (xi, yj).

• Find a smooth piecewise cubic function z(x, y)
  to interpolate
• the data set (xi, yj, zij), i 0,, I, j
  0,, J.

29
Sibson Elements

• Given a rectangle xi, xi1yj, yj1, which
  is divided into
• four triangles by drawing the two diagonals.

30
Sibson Elements cont

• A Sibson element (Han and Schumaker 1997) is a
  bivariate
• piecewise polynomial z(x, y) that is cubic in
  each triangle, and

• is C1 at the lines separating the four triangles,

• is C1 with the Sibson elements in the adjacent
  rectangles,

• has derivative ?z/ ?x that is linear along the
  edges
• x xi and x xi1,

• has derivative ?z/ ?y that is linear along the
  edges
• y yj and y yj1.

31
Sibson Element Cubic L1 Splines

• A piecewise cubic polynomial (x, y) is called
  a Sibson (2nd
• derivative based) L1 spline, if

where z(x, y) is a Sibson element, Tij is the
rectangular area xi, xi1yj, yj1.
32
Sibson Element Cubic L1 Splines cont
z(x, y) can be expressed on each triangle Tijk, k
1, , 4, of rectangle Tij as
33
Sibson Element Cubic L1 Splines cont
34
C1 Smooth Constraints
Over segments (A, C), (B, C), (C, D) and (C, E)
35
C1 Smooth Constraints cont
Over segments (A, B), (B, E), (D, E) and (A, D)
All these constraints can be written as a cone
constraint Ac 0.
36
Math Model of Sibson Element Cubic L1 Spline

• The objective function is nondifferentiable!

• Variables are determined by the two gradients at
  each node.

37
Primal GP
(Primal)
where
Nondifferentiable !
X cAc 0
C zij?0? R9 ?zi,j1 ?0 ? R9 ?
zi1,j1 ? 0 ? zi1,j ? R9
38
Dual GP
(Dual)
where
Y is the row space of A.
39
Dual GP cont

• The dual problem is a convex program with a
  linear objective function and
• convex cubic constraints.

is a convex set
where
where
40
Dual to Primal Transformation

• According to the optimality conditions, the dual
  optimal solution can be
• transformed to a primal optimal solution by
  solving the following LP

(T)
where P?R(4IJ)?(48IJ)
41
Bivariate Results

• Theoretical proof of linearity preservation
  property.
• Under some dual interior assumption, if four
  data points lie on a plane,
  
• then the bivariate cubic L1 spline preserves
  linearity over the rectangular area.

• Development of a computationally efficient
  compressed primal-dual algorithm.
• The algorithm generates discretized bivariate
  cubic L1 splines. It takes less than
• one hour to solve problems with 100100 knots
  on Pentium-IV 2.8GHz PC.

• Real application to terrain description.
• Reference Wang, Y., Fang S.-C., and Lavery, J.
  (2005a, 2005b)

42
Linearity Preservation
43
Terrain Application
44
L1 vs. L2 Sibson Element Splines
Sibson 2nd derivative based L1 Spline
Sibson 2nd derivative based L2 Spline (conventiona
l spline)
45
Bivariate L1 Splines over TIN

• A triangulated irregular network (TIN) is a
  triangulation
• of data locations (xi,yi), i 1,,N, formed
  by a set of triangles
• such that
• whose vertices are the data locations,
• whose interiors do not meet, and
• whose union covers the convex hull of the data
  locations.

• zi is given at each knot (xi, yi).

• Find a smooth piecewise cubic function z(x, y)
  to interpolate
• the data set (xi, yi, zi), i 0,, N.

46
rHCT Elements

• rHCT(reduced Hsieh-Clough
• -Tocher) element is defined over
• a triangular area with vertices
• (x1,y1), (x2, y2) and (x3,y3).

• It is determined
• by the interpolation value and
• the first partial derivatives with
• respective to x and y at each
• vertex, i.e., (z1, z1x, z1y), (z2, z2x,
• z2y), (z3, z3x, z3y).

• According to the barycenter of
• the footprint triangle, it is
• partitioned into three
• subtriangles, B1, B2, B3. On each
• subtriangle, a bivariate cubic
• polynomial is defined.

47
rHCT Elements cont

• An rHCT element is a bivariate piecewise
  polynomial z(x, y)
• that is cubic in each triangle.

• The rHCT element has following properties

• It is C1 at the lines separating the three
  subtriangles.

• It is C1 with the rHCT elements in the adjacent
  triangles.

• Along each edge of the footprint triangle, the
  derivative taken
• in the direction perpendicular to the edge
  varies linearly.

48
rHCT Element Cubic L1 Spline

• A TIN partitions domain ? into K triangles B1,
  , BK.
• Triangle Bi has vertices
  

and is further partitioned into , ,
according to the barycenter of Bi.

• The math model becomes

when and are common vertices.
s.t.
49
rHCT Element Cubic L1 Spline

• GP models have been developed.
• Compressed primal-dual algorithm works well.
• More flexible than Sibson elements, but
  triangulation presents a new twist
• for bi-level optimization.

• Reference Zhang, W., Wang, Y., Fang, S.-C., and
  Lavery, J. (2005)

50
rHCT Elements Illustration
grid (a)
grid (b)
grid (c)
51
rHCT Elements Illustration
rHCT based L1 spline on grid (a)
rHCT based L1 spline on grid (b)
rHCT based L1 spline on grid (c)
Sibson based L1 spline
52
1st Derivative Based Univariate Cubic L1 Spline

• The problem is equivalent to solving the
  following optimization problem.

s.t.
53
Dual GP

• Dual problem is a convex program with a linear
  objective and convex cubic
• constraints.

54
1st Derivative Based rHCT L1 Spline

• The math model becomes

s.t.
where
55
Comparisons
First derivative based L1 spline
First derivative based L2 spline
Second derivative based L1 spline
Second derivative based L2 spline
56
Comparisons cont
First derivative based L1 spline
First derivative based L2 spline
Second derivative based L1 spline
Second derivative based L2 spline
57
On-Going Research

• Exact solution method to the 2nd derivative
  based bivariate cubic
• L1 splines.

• Optimal triangulation for rHCT elements.

• First derivative based bivariate cubic L1splines.

• Error propagation and domain decomposition for
  large scale
• applications.

• Haptic device trivariate applications.

58
Error Propagation
59
Error Propagation cont
60
Domain Decomposition
Overlap 1
Overlap 2
61
Domain Decomposition cont
62
L1 Splines for Urban Terrain
2nd derivative L1 Spline Representation of
Baltimore Football Stadium
63
Baltimore City 1
64
Baltimore City 2
65
Baltimore City 3
66
Baltimore City 4
67
Conclusion

• Changing from L2 norm to L1 norm has huge
  positive
• influence on shape preservation.

• First derivative based L1 spline preserves shape
  better than
• second derivative based L1 spline.

• No constraints, penalties, a posteriori
  filtering or user
• interaction required (but can be used if
  desired)

• Same geometric programming framework for
  univariate and
• multivariate.

• Useful for terrain (line of sight, view shed),
  geophysical features,
• geographical data, biological objects,
  mechanical objects, image
• compression / analysis, etc.

68
References

• Cheng, H., Fang, S.-C., and Lavery, J.,
  Univariate cubic L1 splines a geometric
• programming approach, Math Methods of
  Operations Research (2002) 56, 197-229.

• Cheng, H., Fang, S.-C., and Lavery, J., An
  efficient algorithm for generating univariate
• cubic L1 splines, Computational Optimization
  and Applications (2004) 29, 219-253.

• Cheng, H., Fang, S.-C., and Lavery, J., A
  geometric programming framework for
• univariate cubic L1 smoothing splines, Annals
  of Operations Research (2005) 133,
• 229-248.

• Cheng, H., Fang, S.-C., and Lavery, J.,
  Shape-preserving properties of univariate cubic
• L1 splines, Computational and Applied
  Mathematics (2005) 174, 361-382.

• Chiu, N.-C., Fang, S.-C., Lin, J.-Y., Wang, Y.,
  and Lavery, J., Approximating term
• structure of interest rate using cubic L1
  splines, submitted to European Journal of
• Operational Research.

• Lavery, J., Univariate cubic Lp splines and
  shape-preserving, multiscale interpolation
• by univariate cubic L1 splines, Computer
  Aided Geometric Design. (2000) 17, 319-336.

69
References cont

• Lavery, J., Shape-preserving, multiscale
  interpolation by bi- and multivariate cubic L1
• splines, Computer Aided Geometric Design.
  (2000) 18, 321-343.

• Wang, Y., Fang, S.-C., and Lavery, J., A
  geometric programming approach for bivariate
• cubic L1 splines, Computers and Mathematics
  with Applications (2005), 49, 481-514.

• Wang, Y., Fang S.-C., and Lavery, J., A
  compressed primal-dual algorithm for
• bivariate cubic L1 splines, Journal of
  Computational and Applied Mathematics (2006),
• to appear.

• Zhang, W., Wang, Y., Fang, S.-C., and Lavery,
  J., Cubic L1 splines on triangulated
• irregular networks, Pacific Journal of
  Optimization (2006), 2, 289-317.

• Zhao, Y., Fang S.-C., and Lavery, J., Geometric
  dual formulation of the first
• derivative based univariate cubic L1 spline
  functions, Journal of Global
• Optimization (2006), to appear.

70
Thank You !fang_at_eos.ncsu.eduwww.ie.ncsu.edu/fan
group