Variational Subdivision
for Laplacian Splines

Technical Report
TR97-291

Henrik Weimer
Joe Warren

Rice University
Department of Computer Science
Houston, TX 77025

August 1997

The fundamental problem of geometric design is the representation of curved shapes. Traditionally such shapes are represented by parametric spline curves, e.g. NURBS, which are defined as the minimizers of variational problems. For example, cubic B-splines minimize the bending energy functional.

More recently subdivision curves and surfaces emerged as an alternative means for representing curved shapes. In this framework a curve or surface is defined as the limit of a repeated averaging process of control points. The subdivision scheme is determined by the subdivision mask S which specifies the weights used during the averaging.

A methodology for deriving subdivision schemes from homogeneous variational problems has been outlined in [Warren97]. This report shows the computational details of this approach and presents a method for deriving a subdivision scheme with predefined local support wich produces limit surfaces close to the minimizer of the variational problem. In particular this document contains the detailed derivation of Laplacian Splines discussed in section 4.2 of [Warren97].

The result will be a local stationary subdivision scheme for bounded rectangular grids which produces a limit surface that is close to the minimizer of the original variational problem. By increasing the support of the subdivision basis functions the resulting surface will come arbitrarily close to the real variationally optimal surface.

The Mathematica source code for this report can be found under URL http://www.cs.rice.edu/~jwarren/Laplace-tr.nb and an on-line HTML version under http://www.cs.rice.edu/~jwarren/Laplace-tr .

Laplacian Splines are defined as the minimizer to the variational problem derived from Laplace's equation

(min_f) 
      (∫_Ω)(
            ((((∂/(∂_((Subscript[t, 1])))) f (t)))^2) + 
              ((((∂/(∂_((Subscript[t, 2])))) f (t)))^2)) 
          dt

where the domain Ω is the real plane.The strong form of this equation is exactly Laplace's equation

((∂^2)/((∂_((Subscript[t, 1])))^2)) f (t) + 
        ((∂^2)/((∂_((Subscript[t, 2])))^2)) f (t) = 0.

We want to refer to functions that satisfy this equation as Laplacian Splines.

This report presents a subdivision scheme which produces the minimizer of the variational probelm as stated in (1) for a bounded rectangular grid. This subdivision scheme has global support. We will also show how  a local subdivision scheme can be constructed which produces a surface that is close to the minimizer of (1). By increasing the support of the local basis function of this approximating subdivision scheme, the resulting limit surface can lie arbitrarily close to the real Laplacian spline.

The energy matrix for Laplacian splines is derived from Laplace's equation. We use a generating function representation of the energy mask e(z).

d1[i_, n_, x_] := 
      Which[i == 0, 1 - x, i == n, -(1/x) + 1, True, 
        -(1/x) + 2 - x]d2[i_, n_, x_] := 
      Which[i == 0, 1 - x, i == n, -(1/x) + 1, True, 
        -(1/x) + 2 - x]Fact[i_, n_] := If[i == 0 || i == n, (1/2), 1]

e[i_, j_, n_] := 
      Expand[(Fact[j, n] d1[i, n, x] + Fact[i, n] d2[j, n, y]) (x^i) 
          (y^j)]ue[i_, j_, n_] := 
      Expand[(Fact[j, n] d1[i, n, (x^2)] + Fact[i, n] d2[j, n, (y^2)]) 
          (x^i) (y^j)]

GenRow[i_, j_, n_] := 
      Flatten[Transpose[
          Drop[Transpose[
              CoefficientList[
                e[i, j, n] + ((∑_(k = 0)^n) (x^k)) (y^(n + 1)), {
                  x, y}]], -1]]]GenRow2[i_, j_, n_] := 
      Flatten[Transpose[
          Drop[Transpose[
              CoefficientList[
                ue2[i, j, n] + 
                  ((∑_(k = 0)^n) (x^k)) (y^(n + 1)), {x, y}]], 
            -1]]]

GenRow[i,j,n] now generates the coefficients of the energy matrix for control point (i,j) in an (n+1)x(n+1) grid. You may want to re-evaluate GenRow below for different values of i and j. GenRow2[i,j,n] returns the coefficients of the energy matrix for the subdivided grid (only for control points corresponding to control points on the previous level grid, i.e. before subdivision).

MatrixForm[Partition[GenRow[1, 1, 4], 5]]

En[n] computes the energy matrix for a (n+1) x (n+1) grid, i.e. a (n+1)^2 x (n+1)^2 matrix computing the energy associated with each of the control points.

UEn[n] computes the energy matrix after upsampling, i.e. a matrix which computes the energy associated with the control points in a (n/2+1) x (n/2+1) grid and then upsamples them to the corresponding  (n+1) x (n+1) grid by inserting zeroes as th energy associated with new control points.

En[n_] := Flatten[Table[GenRow[i, j, n], {i, 0, n}, {j, 0, n}], 1]UEn[n_] := 
      Flatten[Table[
          If[EvenQ[i] && EvenQ[j], GenRow[(i/2), (j/2), (n/2)], 
            Table[0, {((((n/2) + 1))^2)}]], {i, 0, n}, {j, 0, n}], 1]

En[n] as defined above is singular, i.e. not invertible. Thus we can not solve En S = UEn by inverting En. To make En invertible we replace one of the (linearly dependent) rows in En by a constraint representing the preservation of moments. Call these modified matricec NewUn and NewUEn respectivey.

Moments[n_] := 
      Flatten[(1/(n^2)) 
          (Table[
              If[Mod[i, n] == 0 && Mod[j, n] == 0, 1, 
                If[Mod[i, n] == 0 || Mod[j, n] == 0, 2, 4]], {i, 0, n}, {j, 
                0, n}])]

NewEn[n_] := 
      ReplacePart[En[n], Moments[n], (1/2) ((((n + 1))^2) + 1)]NewUEn[n_] := 
      ReplacePart[UEn[n], Moments[(n/2)], (1/2) ((((n + 1))^2) + 1)]

Evaluate the following two lines to compute the exact subdivision matrix carrying control points from a 5x5 to a 9x9 control grid and from a 9x9 to a 17x17 grid. The resulting sequence of subdivision matrices will produce a limit surface that exactly solves Laplace's equation.

S5to9 = LinearSolve[N[NewEn[8]], NewUEn[8]];

S9to17 = LinearSolve[N[NewEn[16]], NewUEn[16]];

Unfortunately the resulting subdivision rules have global support (the exact solution must have global support as the variational problem is global). Here are some of the resulting exact subdivision rules rounded to three digits. Note the fast decay of the coefficients.

An internal vertex rule:

MatrixForm[
      Partition[
        Round[S9to17[[145]]*100.0]
            /100.0, 9]]

An internal horizontal edge rule:

MatrixForm[
      Partition[
        Round[S9to17[[146]]*100.0]
            /100.0, 9]]

An internal face rule:

MatrixForm[
      Partition[
        Round[S9to17[[163]] 100.0]
            /100.0, 9]]

Apply the subdivision matrices to some control points

ListPlot3D[
      Partition[
        S9to17 . 
          S5to9 . {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 
            1, 0, 0, 1, 1, 1}, 17], PlotRange -> All]

Note that the limit surface will be the exact solution to Laplace's equation.

The exact solution to the equation En S = U En as computed above leads to subdivision rules with global (but exponentially decaying) support.

In the sequel we want to find a subdivision scheme with local support that approximates the solution well. Our approach is to set up an unknown subdivision matrix with variable entries at locations that enforce the desired symmetries and sparsity. Call this subdivision matrix Shat. Now we can solve for the actual coefficients in Shat by minimizing the expression 
En Shat = U En in the infinity norm. Recall that the infinity norm of a matrix is the maximum of the sums of the absolute values of entries in any row.

Set up a sparse subdivision matrix with the desired symmetry and sparsity

Some helper functions

Clear[basicRule]basicRule[i_, j_, y_, x_] := 
      Which[i == 0 && j == 0, 
        S1[0, 0] + x S1[0, 1] + y S1[0, 1] + x y S1[1, 1], 
        i == 1 && j == 0, 
        S2[0, 0] + x S2[0, 1] + y S2[1, 0] + x y S2[1, 1], 
        i == 1 && j == 1, 
        S3[0, 0] + x S3[0, 1] + y S3[0, 1] + x y S3[1, 1], 
        i == 2 && j == 0, 
        S4[0, 0] + x S4[0, 1] + y S4[1, 0] + x y S4[1, 1] + 
          (y^2) S4[0, 0] + x (y^2) S4[0, 1], i == 2 && j == 1, 
        S5[0, 0] + x S5[0, 1] + y S5[1, 0] + x y S5[1, 1] + 
          (y^2) S5[0, 0] + x (y^2) S5[0, 1], i == 2 && j == 2, 
        S6[0, 0] + x S6[0, 1] + (x^2) S6[0, 0] + y S6[0, 1] + 
          x y S6[1, 1] + (x^2) y S6[0, 1] + (y^2) S6[0, 0] + 
          x (y^2) S6[0, 1] + (x^2) (y^2) S6[0, 0], i == 3 && j == 0, 
        S7[0, 0] + x S7[0, 1] + y S7[0, 0] + x y S7[0, 1], 
        i == 3 && j == 1, 
        S8[0, 0] + x S8[0, 1] + y S8[0, 0] + x y S8[0, 1], 
        i == 3 && j == 2, 
        S9[0, 0] + x S9[0, 1] + (x^2) S9[0, 0] + y S9[0, 0] + 
          x y S9[0, 1] + (x^2) y S9[0, 0], i == 3 && j == 3, 
        S10[0, 0] + x S10[0, 0] + y S10[0, 0] + x y S10[0, 0], True, 
        Print[undefined]]

insideLeftBoundary[i_, j_, x_, y_] := 
      Which[EvenQ[i] && j == 0, basicRule[2, 0, x, y], EvenQ[i] && j == 1, 
        basicRule[2, 1, x, y], OddQ[i] && j == 0, basicRule[3, 0, x, y], 
        OddQ[i] && j == 1, basicRule[3, 1, x, y], True, Print[undefined]]insideRightBoundary[i_, j_, x_, y_] := 
      Which[EvenQ[i] && j == 0, basicRule[2, 0, x, (1/y)] y, 
        EvenQ[i] && j == 1, basicRule[2, 1, x, (1/y)] y, OddQ[i] && j == 0, 
        basicRule[3, 0, x, (1/y)] y, OddQ[i] && j == 1, 
        basicRule[3, 1, x, (1/y)] y, True, Print[undefined]]insideTopBoundary[i_, j_, x_, y_] := insideLeftBoundary[j, i, y, x]insideBottomBoundary[i_, j_, x_, y_] := 
      insideRightBoundary[j, i, y, x]

upperLeftCorner[i_, j_, x_, y_] := 
      Which[i ≥ j, basicRule[i, j, x, y], True, 
        basicRule[j, i, y, x]]upperRightCorner[i_, j_, x_, y_] := 
      upperLeftCorner[i, j, x, (1/y)] ylowerLeftCorner[i_, j_, x_, y_] := upperRightCorner[j, i, y, x]lowerRightCorner[i_, j_, x_, y_] := 
      upperRightCorner[i, j, (1/x), y] x

Clear[findRightRule]findRightRule[i_, j_, n_] := 
      With[{fromLeft = j, fromRight = 2 n - 2 - j, fromTop = i, 
          fromBottom = 2 n - 2 - i}, 
        Which[fromLeft < 2 && fromRight ≥ 2 && 
            fromTop ≥ 2 && fromBottom ≥ 2, 
          insideLeftBoundary[i, fromLeft, x, y], 
          fromLeft ≥ 2 && fromRight < 2 && 
            fromTop ≥ 2 && fromBottom ≥ 2, 
          insideRightBoundary[i, fromRight, x, y], 
          fromLeft ≥ 2 && fromRight ≥ 2 && 
            fromTop < 2 && fromBottom ≥ 2, 
          insideTopBoundary[i, j, x, y], 
          fromLeft ≥ 2 && fromRight ≥ 2 && 
            fromTop ≥ 2 && fromBottom < 2, 
          insideBottomBoundary[fromBottom, j, x, y], 
          fromLeft < 2 && fromTop < 2, upperLeftCorner[i, j, x, y], 
          fromLeft < 2 && fromBottom < 2, 
          lowerLeftCorner[fromBottom, fromLeft, x, y], 
          fromRight < 2 && fromTop < 2, 
          upperRightCorner[fromTop, fromRight, x, y], 
          fromRight < 2 && fromBottom < 2, 
          lowerRightCorner[fromBottom, fromRight, x, y], 
          fromLeft ≥ 2 && fromRight ≥ 2 && 
            fromTop ≥ 2 && fromBottom ≥ 2, 
          Which[EvenQ[i] && EvenQ[j], basicRule[2, 2, x, y], 
            EvenQ[i] && OddQ[j], basicRule[3, 2, y, x], OddQ[i] && EvenQ[j], 
            basicRule[3, 2, x, y], OddQ[i] && OddQ[j], 
            basicRule[3, 3, x, y]]]]

oneRow[i_, j_, n_] := 
      Which[i < 2 && j < 2, findRightRule[i, j, n], 
        i ≥ 2 && j < 2, 
        findRightRule[i, j, n] (x^Floor)[((i - 1)/2)], 
        i < 2 && j ≥ 2, 
        findRightRule[i, j, n] (y^Floor)[((j - 1)/2)], 
        i ≥ 2 && j ≥ 2, 
        findRightRule[i, j, n] (x^Floor)[((i - 1)/2)] 
          (y^Floor)[((j - 1)/2)]]

oneFlatRow[i_, j_, n_] := 
      Flatten[Transpose[
          Drop[Transpose[
              CoefficientList[
                oneRow[i, j, n] + 
                  ((∑_(k = 0)^(n - 1))(x^k)) (y^n), {x, y}]], 
            -1]]]

Definition of the Sparse Matrix

Shat[n_] := 
      Flatten[Table[oneFlatRow[j, i, n], {j, 0, 2 n - 2}, {i, 0, 2 n - 2}], 
        1]

Now Shat[n] computes a sparse (variable) subdivision matrix carrying the control points from a bounded nxn grid to a (2n-1)x(2n-1) grid. 

For example Shat[3] subdivides a 3x3 grid into a 5x5 grid. Note that we enforced sparsity and symmetry in the structure of this matrix.

MatrixForm[Shat[3]]

Set up a Linear Programming Problem to Solve for the Actual Coefficients in the Subdivision Matrix

Starting from the equation E S = U E we set up a linear programming problem to solve for S. The linear program will find an S such that the norm || E S - U E || is minimal, where ||.|| denotes the L_infinity norm, i.e. the sum of the absolut values of the elements in any one row should be minimal.

One crucial observation is that the symmetry and sparsity structure of S as defined above laeds only to a finite and fixed number of different rows in E S - U E, i.e. it is sufficient to minimize this small and fixed number of rows once to find a minimal solution for any size of S.

We want to solve here separately for the purely interior subdivisin rules as solving the whole system is beyond the abilities of Mathematica. We solved for the complete system using the CPLEX linear programming package.

Some helper functions to set up the linear program

GetPosNegSubst[vars_, basename_] := 
      Table[vars[[i]] -> 
          ToExpression[basename<>"P"<>ToString[i]] - 
            ToExpression[basename<>"N"<>ToString[i]], {i, 1, 
          Dimensions[vars][[1]]}]

GetSubst[vars_, basename_] := 
      Table[vars[[i]] -> 
          ToExpression[basename<>ToString[i]], {i, 1, 
          Dimensions[vars][[1]]}]

ExtractVars[vectorOfRelations_] := 
      Union[Flatten[
          Table[Variables[
                vectorOfRelations[[i, 1
                  ]]] ⋃ 
              Variables[
                vectorOfRelations[[i, 2
                  ]]], {i, 1, 
              Dimensions[vectorOfRelations][[1
                ]]}]]]

FindExpr[expr_, list_] := 
      Which[First[list] === expr, 1, True, 1 + FindExpr[expr, Rest[list]]]

RowSumConstraint[oneRow_, allExpressions_, objective_] := 
      (∑_(i = 1)^(Dimensions[oneRow][[1
              ]]))With[{
              thisExpr = oneRow[[i]]}, 
            With[{thisIndex = FindExpr[thisExpr, allExpressions]}, 
              ToExpression["SlackPos"<>ToString[thisIndex]] + 
                ToExpression["SlackNeg"<>ToString[thisIndex]]]] 
        ≤ objective

SlackEqualities[allExpressions_] := 
      Table[With[{
            thisExpr = 
              allExpressions[[i]]}, 
          thisExpr == 
            ToExpression["SlackPos"<>ToString[i]] - 
              ToExpression["SlackNeg"<>ToString[i]]], {i, 1, 
          Dimensions[allExpressions][[1
            ]]}]

DiscardZeros[sortedRow_] := 
      If[ToString[First[sortedRow]] == "0", 
        DiscardZeros[Rest[sortedRow]], sortedRow]

UniqueRows[matrix_] := 
      Union[Table[
          DiscardZeros[
            Sort[matrix[[i]]]], {i, 1, 
            Dimensions[matrix][[1
              ]]}]]

grepRules[rules_, interestingVars_] := 
      Complement[
        Table[If[
            Complement[
                Variables[rules[[i]]], 
                interestingVars] == {}, 
            rules[[i]], 0], {i, 1, 
            Dimensions[rules][[1
              ]]}], {0}]

AllMatrixConstraints[matrix_, objective_] := 
      With[{allExprs = Union[Flatten[matrix]]}, 
        SlackEqualities[allExprs] ⋃ 
          Table[RowSumConstraint[
              matrix[[i]], allExprs, 
              objective], {i, 1, 
              Dimensions[matrix][[1
                ]]}]]

DumpConstraints[constraints_, filename_] := 
      With[{stream = OpenWrite[filename, PageWidth -> 1000000]}, 
        For[i = 1, 
          i ≤ 
            Dimensions[constraints][[1
              ]], i++, 
          Write[stream, 
            constraints[[i]]]]; 
        Close[stream]; ]

Get the distinct rows in E S - U E

DSmUD[n_] := En[2 n - 2] . Shat[n] - UEn[2 n - 2]

Evaluate the following line to find all the different rows in E S - U D.

uniqueRows = UniqueRows[DSmUD[5]];

Note that using a larger system leads to the same distinct rows in E S - U D:

moreUniqueRows = UniqueRows[DSmUD[7]];

Dimensions[uniqueRows]Dimensions[moreUniqueRows]

Partition these rows of E S - U E  into three categories: 
     purely interior
     interior and boundary
     interior, boundary and corner

cornerVars = {S1[0, 0], S1[0, 1], S1[1, 1], S2[0, 0], S2[1, 0], 
        S2[0, 1], S2[1, 1], S3[0, 0], S3[0, 1], S3[1, 1]}; boundaryVars = {S4[0, 0], S4[0, 1], S4[1, 0], S4[1, 1], S5[0, 0], 
        S5[0, 1], S5[1, 0], S5[1, 1], S7[0, 0], S7[0, 1], S8[0, 0], 
        S8[0, 1]}; interiorVars = {S6[0, 0], S6[0, 1], S6[1, 1], S9[0, 0], S9[0, 1], 
        S9[0, 2], S10[0, 0]};

interiorRows = grepRules[uniqueRows, interiorVars]; boundaryRows = 
      grepRules[uniqueRows, interiorVars ⋃ boundaryVars]; cornerRows = 
      grepRules[uniqueRows, 
        interiorVars ⋃ boundaryVars ⋃ cornerVars];

Dimensions[interiorRows]Dimensions[boundaryRows]Dimensions[cornerRows]

Set up the systems of constraints

Off[General::"spell1"]interiorConstraints = AllMatrixConstraints[interiorRows, objective]; 
    boundaryConstraints = AllMatrixConstraints[boundaryRows, objective]; 
    cornerConstraints = AllMatrixConstraints[cornerRows, objective]; 
    On[General::"spell1"]

Set up constraints to enforce constant and linear precision

One[n_] := Flatten[Table[1, {x, 1, n}, {y, 1, n}]]

tmpOne = One[5]; tmpOneSubdiv = One[9]; tmpS5 = Shat[5];

eigenConstraints[relevant_] := 
      Flatten[Table[{
            tmpS5[[relevant[[i
                    ]]]] . tmpOne == 
              tmpOneSubdiv
                [[relevant[[i
                  ]]]]}, {i, 1, 
            Dimensions[relevant][[1
              ]]}]]

interiorEigenConstr = eigenConstraints[{21, 30, 31}]; boundaryEigenConstr = eigenConstraints[{19, 20, 21, 28, 29, 30, 31}]; 
    cornerEigenConstr = 
      eigenConstraints[{1, 10, 11, 19, 20, 21, 28, 29, 30, 31}];

Finally solve for the coefficients for the interior

Remember that all variables cann be positive or negative but the solver handles only positive variables. Thus all variables X in the system have to replaced by terms of the form (Xpositive - Xnegative).

Off[General::"spell1"]posNegSubstitution = 
      GetPosNegSubst[interiorVars ⋃ boundaryVars ⋃ cornerVars, 
        "foo"]; On[General::"spell1"]

allInteriorConstraints = 
      interiorConstraints ⋃ interiorEigenConstr /. posNegSubstitution; 
    allInteriorVars = ExtractVars[allInteriorConstraints];

interiorCoeffs = 
      ConstrainedMin[objective, N[allInteriorConstraints], 
          N[allInteriorVars]][[2]]

Here are the resulting subdivision rules for the interior:

InteriorS5 = N[tmpS5] /. N[posNegSubstitution] /. interiorCoeffs;

MatrixForm[
      Partition[InteriorS5[[31]], 5]]

MatrixForm[
      Partition[
        Round[InteriorS5[[40]]*100.0]
            /100.0, 5]]

MatrixForm[
      Partition[
        Round[InteriorS5[[41]]*100.0]
            /100.0, 5]]

Sparse Solution for the Complete Grid

Mathematica's linear optimizer "ConstrainedMin" simply seems to run forever on the linear program arising from the complete set of constraints. We used an external solver, CPLEX, to solve for these coefficients instead. 

Note that by construction these subdivision rules have constant precision (enforced in the constraints).

Here are the results:

cornerCoeffs = {S1[0, 0] -> 1.430966, 
        S1[0, 1] -> -0.227169, S1[1, 1] -> 0.02337218, 
        S10[0, 0] -> 0.25, S2[0, 0] -> 0.4309657, 
        S2[0, 1] -> 0.09808379, S2[1, 0] -> 0.4475783, 
        S2[1, 1] -> 0.02337218, S3[0, 0] -> 0.2178049, 
        S3[0, 1] -> 0.2834315, S3[1, 1] -> 0.2153321, 
        S4[0, 0] -> -0.1427126, S4[0, 1] -> 0.05264114, 
        S4[1, 0] -> 1.450619, S4[1, 1] -> -0.2704763, 
        S5[0, 0] -> 0.004644042000000001, 
        S5[0, 1] -> 0.05624038000000001, S5[1, 0] -> 0.4709391, 
        S5[1, 1] -> 0.4072921, S6[0, 0] -> 0.002322021, 
        S6[0, 1] -> -0.1111111, S6[1, 1] -> 1.435156, 
        S7[0, 0] -> 0.4198616000000001, S7[0, 1] -> 0.0801384, 
        S8[0, 0] -> 0.2308056, S8[0, 1] -> 0.2691944, 
        S9[0, 0] -> 0.02777778, S9[0, 1] -> 0.4444444};

Sparse5 = tmpS5 /. cornerCoeffs;

Sparse9 = Shat[9] /. cornerCoeffs;

MatrixForm[
      Partition[Sparse5[[11]], 5]]

MatrixForm[
      Partition[Sparse5[[37]], 5]]

MatrixForm[
      Partition[Sparse5[[40]], 5]]

MatrixForm[
      Partition[Sparse5[[41]], 5]]

Now apply the sparse S to some control points and compare this visually to the exact solution

ListPlot3D[
      Partition[{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 1, 0, 0, 1, 1, 1, 
          0, 0, 1, 1, 1}, 5], PlotRange -> All, Boxed -> False, 
      Axes -> False, BoxRatios -> {1, 1, 1}, 
      HiddenSurface -> False]

The sparse approximation based on 5x5 basis functions:

ListPlot3D[
      Partition[
        Sparse9 . 
          Sparse5 . {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 1, 0, 0, 1, 1, 
            1, 0, 0, 1, 1, 1}, 17], PlotRange -> All, 
      Boxed -> False, Axes -> False, BoxRatios -> {1, 1, 1}]

The exact subdivision process with globally supported basis functions:

ListPlot3D[
      Partition[
        S9to17 . 
          S5to9 . {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 1, 0, 0, 1, 1, 
            1, 0, 0, 1, 1, 1}, 17], PlotRange -> All, 
      Boxed -> False, Axes -> False, BoxRatios -> {1, 1, 1}]

References

[Warren97] Joe Warren: Subdivision Schemes for Variational Splines. July 1997, submitted for publication to CAGD.

Acknowledgements

This work has been supported in part by NSF grant CCR-9500572 and TATP grant 003604-010.