Quadratic B-Spline Surfaces with Free Parameters for the Interpolation of Curve Networks

Abstract

In this paper, we propose a method for constructing spline surfaces interpolating a B-spline curve network, allowing the presence of free parameters, in order to model the interpolating surface. We provide a constructive algorithm for its generation in the case of biquadratic tensor product B-spline surfaces and bivariate B-spline surfaces on criss-cross triangulations. Finally, we present graphical results.

1. Introduction

When modeling free form shapes, e.g., ship hull forms, it could be more convenient to generate smooth parametric surfaces, interpolating a curve network, instead of using control points. This approach has been studied in the literature from different view points and with different methods, such as blending-function methods [1,2,3,4,5,6], smooth regularly parametrized piecewise polynomial surfaces [7] or subdivision schemes [8,9].

Tensor product B-spline surfaces, despite their limits, due to the parametric domain rectangular topology, can be a good tool in CAGD since they are quite simple to use and not too tightly constrained for the object shape in some cases (see e.g., [10,11] (Chapter 10)). Possible unwanted oscillations can be produced by their high coordinate degree, but they can be successfully avoided by smooth parametric surfaces on triangulations, having a lower total degree [12,13,14].

In this paper, we compare biquadratic tensor product B-spline surfaces and bivariate B-spline surfaces on criss-cross triangulations in order to detect advantages and disadvantages of both constructions. We investigate the degree two because it is the first degree that guarantees continuity $C^1$ and it is computationally simpler with respect to higher degrees. We note that in [10] (Chapter 10) a method based on Gordon surfaces is proposed to construct a tensor product B-spline surface interpolating a curve network. In this case, it is necessary to define three surfaces belonging to different spline spaces to be merged together by the knot refinement algorithm, while the method proposed here is directly based on the knot vectors of the given B-spline curves.

Starting from a so called minimal configuration [15], where a $C^1$ quadratic B-spline curve network is interpolated by a quadratic B-spline surface, whose knots match the curve network knots, we consider a more general approach by introducing some free parameters that provide some degree of freedom in the surface shape and that are obtained by letting some surface knots not satisfying such a match. This means that some surface knots are not network knots. In this regard, we investigate some properties of the obtained surfaces, both biquadratic and quadratic.

The paper is organized as follows.

The B-spline curve network and its compatibility conditions are introduced in Section 2, while the B-spline surface and its existence conditions are presented in Section 3 with an algorithm for its construction.

In Section 4, we present some numerical and graphical results.

2. The Network of B-Spline Curves

Let us denote by ${\left\{P_j^{\left(r\right)}\right\}}_{j=0}^{n+2}$, $r=0,\dots ,R+1$ and ${\left\{Q_i^{\left(s\right)}\right\}}_{i=0}^{m+2}$, $s=0,\dots ,S+1$, $P_j^{\left(r\right)}$, $Q_i^{\left(s\right)}\in {\mathbb{R}}^3$, two given sets of control points, defining the $C^1$ quadratic B-spline curves

$$\begin{array}{c}{\displaystyle C_r\left(v\right)=\sum _{j=0}^{n+2}{P}_j^{\left(r\right)}{B}_j\left(v\right|V),r=0,\dots ,R+1,v\in [c,d]\subset \mathbb{R},}\hfill \\ {\displaystyle D_s\left(u\right)=\sum _{i=0}^{m+2}{Q}_i^{\left(s\right)}{B}_i\left(u\right|U),s=0,\dots ,S+1,u\in [a,b]\subset \mathbb{R},}\hfill \end{array}$$

(1)

with $m\ge R$ and $n\ge S$, respectively. As in [10] (Chapter 10) and in [15], here we suggest the following compatibility conditions to be satisfied by (1):

C1.

all the curves $D_s\left(u\right)$ are defined on

$$U=\{a=u_{-2}\equiv u_{-1}\equiv u_0<u_1<\cdots <u_m<u_{m+1}\equiv u_{m+2}\equiv u_{m+3}=b\}$$

and all the curves $C_r\left(v\right)$ are defined on

$$V=\{c=v_{-2}\equiv v_{-1}\equiv v_0<v_1<\cdots <v_n<v_{n+1}\equiv v_{n+2}\equiv v_{n+3}=d\}.$$

The quadratic B-spline functions $B_i\left(u\right|U)$ and $B_j\left(v\right|V)$, involved in (1) have supports $[u_{i-2},u_{i+1}]$ and $[v_{j-2},v_{j+1}]$, respectively [16,17]. We also consider the so called Greville abscissae, given by

$$s_i=\frac{u_{i-1}+u_i}{2},i=0,\dots ,m+2$$

and

$$t_j=\frac{v_{j-1}+v_j}{2},j=0,\dots ,n+2,$$

to define $(s_i,v_s)$, pre-image of $Q_i^{\left(s\right)}$ for all s, and $(u_r,t_j)$, pre-image of $P_j^{\left(r\right)}$ for all r (Figure 1);

C2.

there exist parameters

$$u_0=u_{i_0}<u_{i_1}<\cdots <u_{i_R}<u_{i_{R+1}}\equiv u_{m+1}$$

and

$$v_0=v_{j_0}<v_{j_1}<\cdots <v_{j_S}<v_{j_{S+1}}\equiv v_{n+1}$$

such that

$$I_{rs}=C_r\left(v_{j_s}\right)=D_s\left(u_{i_r}\right),r=0,\dots ,R+1,s=0,\dots ,S+1.$$

(2)

Condition C2. is always satisfied by at least a curve network, even if it is supposed to be given. This can be proven following the same logical scheme given in [15]. Indeed, defining the following quantities:

$$\begin{array}{c}{\displaystyle {\sigma }_{r+1}=\frac{h_{r+1}}{h_r+h_{r+1}},{\sigma }_r^{\prime }=\frac{h_{r-1}}{h_{r-1}+h_r},}\\ {\displaystyle {\tau }_{s+1}=\frac{k_{s+1}}{k_s+k_{s+1}},{\tau }_s^{\prime }=\frac{k_{s-1}}{k_{s-1}+k_s},}\end{array}$$

(3)

for $0\le r\le m+2$ and $0\le s\le n+2$, with $h_{\mu }=u_{\mu }-u_{\mu -1}$, $k_{\nu }=v_{\nu }-v_{\nu -1}$ and $h_{-1}=k_{-1}=h_{m+2}=k_{n+2}=0$ (setting to zero a possible $\frac{0}{0}$ in (3)), since from [18] (Chapter X) and [19] (Section 1.5)

$$\begin{array}{cc}\hfill B_{i_r}\left(u_{i_r}\right|U)& ={\sigma }_{i_r+1},B_{i_r+1}\left(u_{i_r}\right|U)={\sigma }_{i_r+1}^{\prime },\hfill \\ \hfill B_{j_s}\left(v_{j_s}\right|V)& ={\tau }_{j_s+1},B_{j_s+1}\left(v_{j_s}\right|V)={\tau }_{j_s+1}^{\prime },\hfill \end{array}$$

and thanks to the local support property of the B-splines, condition C2. can be written as

$${\tau }_{j_s+1}{P}_{j_s}^{\left(r\right)}+{\tau }_{j_s+1}^{\prime }{P}_{j_s+1}^{\left(r\right)}={\sigma }_{i_r+1}{Q}_{i_r}^{\left(s\right)}+{\sigma }_{i_r+1}^{\prime }{Q}_{i_r+1}^{\left(s\right)},r=1,\dots ,R,s=1,\dots ,S.$$

Instead, on the boundary of the domain, we obtain

$$\begin{array}{cc}\hfill & {\sigma }_{i_r+1}{Q}_{i_r}^{\left(0\right)}+{\sigma }_{i_r+1}^{\prime }{Q}_{i_r+1}^{\left(0\right)}=P_0^{\left(r\right)},s=0,r=1,\dots ,R\hfill \\ \hfill & {\sigma }_{i_r+1}{Q}_{i_r}^{(S+1)}+{\sigma }_{i_r+1}^{\prime }{Q}_{i_r+1}^{(S+1)}=P_{n+2}^{\left(r\right)},s=S+1,r=1,\dots ,R\hfill \\ \hfill & {\tau }_{j_s+1}{P}_{j_s}^{\left(0\right)}+{\tau }_{j_s+1}^{\prime }{P}_{j_s+1}^{\left(0\right)}=Q_0^{\left(s\right)},r=0,s=1,\dots ,S\hfill \\ \hfill & {\tau }_{j_s+1}{P}_{j_s}^{(R+1)}+{\tau }_{j_s+1}^{\prime }{P}_{j_s+1}^{(R+1)}=Q_{m+2}^{\left(s\right)},r=R+1,s=1,\dots ,S.\hfill \end{array}$$

and

$$Q_0^{\left(0\right)}=P_0^{\left(0\right)},Q_{m+2}^{\left(0\right)}=P_0^{(R+1)},Q_0^{(S+1)}=P_{n+2}^{\left(0\right)},Q_{m+2}^{(S+1)}=P_{n+2}^{(R+1)}.$$

Therefore, taking into account that $m\ge R$, $n\ge S$, it is always possible to construct a curve network satisfying condition C2, because the number of its control points ${\overline{N}}_d=(n+3)(R+2)+(m+3)(S+2)$ is greater than the number of constraints $N_c=(R+2)(S+2)$. Moreover, we remark that the case of minimal configuration is given when $R=m$ and $S=n$.

3. Construction of the Interpolating ${\mathit{C}}^{\mathbf{1}}$ B-Spline Surface

Let us consider two possible subdivisions of the parametric domain $\mathsf{\Omega }=[a,b]\times [c,d]$: a rectangular partition $R_{mn}$ (Figure 1a), based on the knots $\{U\times V\}$ and a criss-cross triangulation $T_{mn}$ (Figure 1b), obtained by drawing the two diagonals in each rectangle of $R_{mn}$.

Here, we consider the two spline spaces it is possible to define on such partitions: $S_{2,2}^{1,1}\left(R_{mn}\right)$ of all $C^1$ biquadratic tensor product splines, being a bivariate polynomial of coordinate degree 2 in each rectangle of $R_{mn}$ and $S_2^1\left(T_{mn}\right)$ of all $C^1$ quadratic splines, being a bivariate polynomial of total degree 2 in each triangle of $T_{mn}$. The first space has a basis defined by the set of $(m+3)(n+3)$ B-splines $\{B_{i,j}(u,v)=B_i\left(u\right|U)·B_j\left(v\right|V),i=0,\dots ,m+2,j=0,\dots ,n+2\}$, while the second one is generated by the set of $(m+3)(n+3)$ B-splines $\{B_{i,j}(u,v),i=0,\dots ,m+2,j=0,\dots ,n+2\}$ with octagonal support (Figure 2) [16,17,20,21,22], being its dimension given by $(m+3)(n+3)-1$.

Let ${\mathcal{B}}_{mn}=\{B_{i,j}(u,v),i=0,\dots ,m+2,j=0,\dots ,n+2\}$ be one of the two sets of spanning functions above mentioned. Then the $C^1$ B-spline parametric surface, interpolating the network (1), can be written as follows:

$$\mathcal{S}(u,v)=\sum _{i=0}^{m+2}\sum _{j=0}^{n+2}{C}_{i,j}{B}_{i,j}(u,v),(u,v)\in \mathsf{\Omega }.$$

(4)

If $B_{i,j}\in S_{2,2}^{1,1}\left(R_{mn}\right)$, then $\mathcal{S}$ is a tensor product surface, while, if $B_{i,j}\in S_2^1\left(R_{mn}\right)$, then $\mathcal{S}$ is a triangulated surface.

To obtain such a surface, it is enough to compute its control points $C_{i,j}\in {\mathbb{R}}^3$, $i=0,\dots ,m+2$, $j=0,\dots ,n+2$, having pre-images $(s_i,t_j)\in \mathsf{\Omega }$ and a way to do it is to consider the network as a set of isoparametric curves of $\mathcal{S}$:

$$S(u_{i_r},v)=C_r\left(v\right),r=0,\dots ,R+1,$$

(5)

$$S(u,v_{j_s})=D_s\left(u\right),s=0,\dots ,S+1.$$

(6)

In the next theorem, the isoparametric conditions (5) and (6) are translated into new conditions on the $\left\{C_{i,j}\right\}$’s.

Theorem 1.

The curve network (1) is a set of isoparametric curves of $\mathcal{S}$ if, and only if, the following conditions are satisfied:

$$\begin{array}{cc}\hfill & {\sigma }_{i_r+1}{C}_{i_r,j}+{\sigma }_{i_r+1}^{\prime }{C}_{i_r+1,j}=P_j^{\left(r\right)},r=1,\dots ,R,j=1,\dots ,n+1,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & {\tau }_{j_s+1}{C}_{i,j_s}+{\tau }_{j_s+1}^{\prime }{C}_{i,j_s+1}=Q_i^{\left(s\right)},s=1,\dots ,S,i=1,\dots ,m+1,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & C_{i,0}=Q_i^{\left(0\right)},i=0,\dots ,m+2,\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & C_{i,n+2}=Q_i^{(S+1)},i=0,\dots ,m+2,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & C_{0,j}=P_j^{\left(0\right)},j=0,\dots ,n+2,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & C_{m+2,j}=P_j^{(R+1)},j=0,\dots ,n+2,\hfill \end{array}$$

(12)

for all control points $\left\{C_{i,j}\right\}$ in (4).

Proof. 

Since the B-splines in ${\mathcal{B}}_{mn}$ have local support, $B_{i,j}(u_{\alpha },v)\equiv 0$ if $i<\alpha ,i>\alpha +1$, any j. Therefore, from (4) it results

$$\mathcal{S}(u_{i_r},v)=\sum _{i=0}^{m+2}\sum _{j=0}^{n+2}{C}_{i,j}{B}_{i,j}(u_{i_r},v)=\sum _{j=0}^{n+2}\left(C_{i_r,j}{B}_{i_r,j}(u_{i_r},v)+C_{i_r+1,j}{B}_{i_r+1,j}(u_{i_r},v)\right).$$

(13)

From the Bernstein–Bézier (BB-) coefficients of $B_{i_r,j}(u_{i_r},v)$, $B_{i_r+1,j}(u_{i_r},v)$, $B_j\left(v\right|V)$, $r=1,\dots ,R,j=0,\dots ,n+2$, [23,24], we get

$$B_{i_r,j}(u_{i_r},v)={\sigma }_{i_r+1}{B}_j\left(v\right|V),B_{i_r+1,j}(u_{i_r},v)={\sigma }_{i_r+1}^{\prime }{B}_j\left(v\right|V).$$

(14)

Then, from (13) and (14), we obtain

$$\mathcal{S}(u_{i_r},v)=\sum _{j=0}^{n+2}\left({\sigma }_{i_r+1}{C}_{i_r,j}+{\sigma }_{i_r+1}^{\prime }{C}_{i_r+1,j}\right)B_j\left(v\right|V).$$

(15)

Thus, from (1), (5) and (15), condition (7) holds.

Now, with $i_0=0$,

$$\mathcal{S}(u_0,v)=\sum _{i=0}^{m+2}\sum _{j=0}^{n+2}{C}_{i,j}{B}_{i,j}(u_0,v)=\sum _{j=0}^{n+2}{C}_{0,j}{B}_{0,j}(u_0,v)$$

(16)

and $B_{0,j}(u_0,v)=B_j\left(v\right|V)$ [23,24], from (1), (5) and (16), we get (11).

Similarly, by setting $u=u_{m+1}$ in (4), since $i_{R+1}=m+1$, (12) follows.

In the same way, from (6), we obtain (8)–(10). □

Remark 1.

From Theorem 1, we can immediately get that

$$\begin{array}{c}{C}_{0,0}=Q_0^{\left(0\right)}=P_0^{\left(0\right)},C_{m+2,0}=Q_{m+2}^{\left(0\right)}=P_0^{(R+1)},\hfill \\ C_{0,n+2}=Q_0^{(S+1)}=P_{n+2}^{\left(0\right)},C_{m+2,n+2}=Q_{m+2}^{(S+1)}=P_{n+2}^{(R+1)}.\hfill \end{array}$$

Now, we want to count the number of free parameters, useful to model the surface $\mathcal{S}$, i.e., the number of control points in the surface definition that can be arbitrarily chosen. Since $m\ge R$ and $n\ge S$, we suppose $m=R+k$ and $n=S+\ell$, with $k,\ell \ge 0$.

Theorem 2.

In order to model the surface $\mathcal{S}$, the number of free control points in $\mathcal{S}$ is

$$N=\left\{\begin{array}{cc}k\ell +\ell +k+1\hfill & in case S_{2,2}^{1,1}\left(R_{mn}\right)\hfill \\ k\ell +\ell +k\hfill & in case S_2^1\left(T_{mn}\right)\hfill \end{array}\right..$$

Proof. 

The number of free parameters (denoted by N) is given by the number of parameters we have to compute (denoted by $N_1$) minus the number of constraints given by the interpolation conditions (denoted by $N_2$).

$N_1$ is the dimension of the spline space used for the construction of the surface:

$$N_1=\left\{\begin{array}{cc}(m+3)(n+3)\hfill & in case S_{2,2}^{1,1}\left(R_{mn}\right)\hfill \\ (m+3)(n+3)-1\hfill & in case S_2^1\left(T_{mn}\right)\hfill \end{array}\right..$$

From Theorem 1, the number of interpolation conditions is $R(n+1)+S(m+1)+2(m+3)+2(n+3)$, but, taking into account the redundant conditions at the four corners of the domain and the compatibility conditions C2. for each inner point $(u_{i_r},v_{j_s})$, for $r=1,\dots ,R$ and $s=1,\dots ,S$, then the number of interpolation conditions is

$$\begin{array}{cc}\hfill N_2& =R(n+1)+S(m+1)+2(m+3)+2(n+3)-4-RS\hfill \\ \hfill & =mn+3m+3n+8-k\ell -k-\ell .\hfill \end{array}$$

Therefore, we obtain

$$N=N_1-N_2=\left\{\begin{array}{cc}k\ell +\ell +k+1\hfill & in case S_{2,2}^{1,1}\left(R_{mn}\right)\hfill \\ k\ell +\ell +k\hfill & in case S_2^1\left(T_{mn}\right)\hfill \end{array}\right..$$

 □

The aim of this section is the construction of the interpolating surface, where the constrained coefficients depend on the control points ${\left\{P_j^{\left(r\right)}\right\}}_{j=0}^{n+1}$ and ${\left\{Q_i^{\left(s\right)}\right\}}_{i=0}^{m+1}$ of the curve network. From (7) and (8), we can observe that

$${\displaystyle C_{i_{r+1},j}=\frac{P_j^{\left(r\right)}}{{\sigma }_{i_r+1}^{\prime }}-C_{i_r,j}\frac{{\sigma }_{i_r+1}}{{\sigma }_{i_r+1}^{\prime }},C_{i,j_{s+1}}=\frac{Q_i^{\left(s\right)}}{{\tau }_{j_s+1}^{\prime }}-C_{i,j_s}\frac{{\tau }_{j_s+1}}{{\tau }_{j_s+1}^{\prime }}}$$

Therefore, the $C_{i,j}$’s that are excluded from (7) and (8), i.e., for

$$i\ne i_r,i_{r+1}\mathrm{and}j\ne j_s,j_{s+1},\mathrm{with}i=1,\dots ,m+1,j=1,\dots ,n+1,r=1,\dots ,R,s=1,\dots ,S$$

(17)

are free parameters. The other parameters can be computed depending on them. Moreover, we can notice that (see Figure 3)

-

the knots in the partitions U and V not associated to a curve of the network are barriers identifying subdomains in $\mathsf{\Omega }$. In each subdomain, we have a minimal configuration;

-

the number of subdomains in $\mathsf{\Omega }$ is equal to $(k+1)(\ell +1)$;

-

the control point in the left-down corner of each subdomain is identified as a free parameter;

-

the control point $C_{1,1}$ is always a free parameter.

The following lemma provides the necessary relations for the computations of the constrained control points.

Lemma 1.

If curves (1) are isoparametric curves of the surface (4), then the constrained control points have the following expression:

$$C_{i,j}={\mathsf{\Gamma }}_{i,j}+{(-1)}^{i+j+{\overline{r}}_{\gamma }+{\overline{s}}_{\delta }}\frac{k_j}{k_{{\overline{s}}_{\delta }+1}}\frac{h_i}{h_{{\overline{r}}_{\gamma }+1}}{C}_{1+{\overline{r}}_{\gamma },1+{\overline{s}}_{\delta }}$$

(18)

for $i=1,\dots ,m+1$, $j=1,\dots ,n+1$ (excluding the free parameters in (17)), with

$$\begin{array}{cc}{\mathsf{\Gamma }}_{i,j}=\hfill & {\displaystyle \sum _{r=1}^{i-{\overline{r}}_{\gamma }-1}{(-1)}^{r+1}\frac{(h_{i-r}+h_{i-r+1})h_i}{h_{i-r}{h}_{i-r+1}}{P}_j^{(i-r-\alpha )}}\hfill \\ \hfill & {\displaystyle +{(-1)}^{i+{\overline{r}}_{\gamma }}\frac{h_i}{h_{{\overline{r}}_{\gamma }+1}}\sum _{s=1}^{j-{\overline{s}}_{\delta }-1}{(-1)}^s\frac{(k_{j-s}+k_{j-s+1})k_j}{k_{j-s}{k}_{j-s+1}}{Q}_{1+{\overline{r}}_{\gamma }}^{(j-s-\beta )}}\hfill \end{array}$$

(19)

and ${\sum }_{\ell =1}^0·=0$, where (see Figure 3)

-

${\overline{r}}_{\gamma }$ is the index of the last barrier before the control point $C_{i,j}$ in the direction U, i.e., ${\overline{r}}_{\gamma }<i<{\overline{r}}_{\gamma +1}$;

-

${\overline{s}}_{\delta }$ is the index of the last barrier before the control point $C_{i,j}$ in the direction V, i.e., ${\overline{s}}_{\delta }<j<{\overline{s}}_{\delta +1}$;

-

α is the number of barriers in the direction U before the control point $C_{i,j}$;

-

β is the number of barriers in the direction V before the control point $C_{i,j}$.

Proof. 

Let us consider a domain, as in Figure 3, and the first subdomain, where minimal configuration is present. From (8), we have ${\tau }_{s+1}{C}_{i,s}+{\tau }_{s+1}^{\prime }{C}_{i,s+1}=Q_i^{\left(s\right)}$, $s=1,\dots ,{\overline{s}}_1$, $i=1,\dots ,{\overline{r}}_1$. For $i=1$, applying it repeatedly, we can write (see Figure 4)

$$C_{1,j}=\sum _{s=1}^{j-1}{(-1)}^{s+1}\frac{Q_1^{(j-s)}}{{\tau }_{j-s+1}^{\prime }}\prod _{k=1}^{s-1}\frac{{\tau }_{j-k+1}}{{\tau }_{j-k+1}^{\prime }}+{(-1)}^{j+1}{C}_{1,1}\prod _{s=1}^{j-1}\frac{{\tau }_{j-s+1}}{{\tau }_{j-s+1}^{\prime }}$$

(20)

where $j=1,\dots ,{\overline{s}}_1$ and ${\prod }_{l=1}^0·=1$.

Then, from (7) it holds ${\sigma }_{r+1}{C}_{r,j}+{\sigma }_{r+1}^{\prime }{C}_{r+1,j}=P_j^{\left(r\right)}$ for $r=1,\dots ,{\overline{r}}_1$, $j=1,\dots ,{\overline{s}}_1$ and similarly, for a general j such that $j\le {\overline{s}}_1$ we have

$$C_{i,j}=\sum _{r=1}^{i-1}{(-1)}^{r+1}\frac{P_j^{(i-r)}}{{\sigma }_{i-r+1}^{\prime }}\prod _{k=1}^{r-1}\frac{{\sigma }_{i-k+1}}{{\sigma }_{i-k+1}^{\prime }}+{(-1)}^{i+1}{C}_{1,j}\prod _{r=1}^{i-1}\frac{{\sigma }_{i-r+1}}{{\sigma }_{i-r+1}^{\prime }}$$

(21)

with $i=1,\dots ,{\overline{r}}_1$. Then, substituting (20) in (21) we have:

$$C_{i,j}={\mathsf{\Gamma }}_{i,j}+{(-1)}^{i+j}{C}_{1,1}\prod _{s=1}^{j-1}\frac{{\tau }_{j-s+1}}{{\tau }_{j-s+1}^{\prime }}\prod _{r=1}^{i-1}\frac{{\sigma }_{i-r+1}}{{\sigma }_{i-r+1}^{\prime }},i=1,\dots ,{\overline{r}}_1,j=1,\dots ,{\overline{s}}_1,$$

with

$${\mathsf{\Gamma }}_{i,j}=\sum _{r=1}^{i-1}{(-1)}^{r+1}\frac{P_j^{(i-r)}}{{\sigma }_{i-r+1}^{\prime }}\prod _{k=1}^{r-1}\frac{{\sigma }_{i-k+1}}{{\sigma }_{i-k+1}^{\prime }}+{(-1)}^{i+1}\prod _{r=1}^{i-1}\frac{{\sigma }_{i-r+1}}{{\sigma }_{i-r+1}^{\prime }}\sum _{s=1}^{j-1}{(-1)}^{s+1}\frac{Q_1^{(j-s)}}{{\tau }_{j-s+1}^{\prime }}\prod _{k=1}^{s-1}\frac{{\tau }_{j-k+1}}{{\tau }_{j-k+1}^{\prime }}.$$

Now, we follow the same scheme in the second subdomain. So, from (8) ${\tau }_{s+1}{C}_{i,s}+{\tau }_{s+1}^{\prime }{C}_{i,s+1}=Q_i^{\left(s\right)}$, $s=1,\dots ,{\overline{s}}_1$, $i={\overline{r}}_1+1,\dots ,{\overline{r}}_2$, we are able to obtain

$$C_{i,j}={\mathsf{\Gamma }}_{i,j}+{(-1)}^{i+j+{\overline{r}}_{\gamma }}{C}_{{\overline{r}}_{\gamma }+1,1}\prod _{s=1}^{j-1}\frac{{\tau }_{j-s+1}}{{\tau }_{j-s+1}^{\prime }}\prod _{r=1}^{i-{\overline{r}}_{\gamma }-1}\frac{{\sigma }_{i-r+1}}{{\sigma }_{i-r+1}^{\prime }}$$

for $i={\overline{r}}_{\gamma }+1,\dots ,{\overline{r}}_{\gamma +1}$ e $j=1,\dots ,{\overline{s}}_1$ with

$$\begin{array}{cc}\hfill {\mathsf{\Gamma }}_{i,j}=& \sum _{r=1}^{i-{\overline{r}}_{\gamma }-1}{(-1)}^{r+1}\frac{P_j^{(i-r-1)}}{{\sigma }_{i-r+1}^{\prime }}\prod _{k=1}^{r-1}\frac{{\sigma }_{i-k+1}}{{\sigma }_{i-k+1}^{\prime }}\hfill \\ & +{(-1)}^{i+{\overline{r}}_{\gamma }+1}\prod _{r=1}^{i-{\overline{r}}_{\gamma }-1}\frac{{\sigma }_{i-r+1}}{{\sigma }_{i-r+1}^{\prime }}\sum _{s=1}^{j-1}{(-1)}^{s+1}\frac{Q_{{\overline{r}}_{\gamma }+1}^{(j-s)}}{{\tau }_{j-s+1}^{\prime }}\prod _{k=1}^{s-1}\frac{{\tau }_{j-k+1}}{{\tau }_{j-k+1}^{\prime }}.\hfill \end{array}$$

In the third subdomain, from (8), we obtain

$$C_{i,j}={\mathsf{\Gamma }}_{i,j}+{(-1)}^{i+j+{\overline{s}}_{\delta }}{C}_{1,{\overline{s}}_{\delta }+1}\prod _{s=1}^{j-{\overline{s}}_{\delta }-1}\frac{{\tau }_{j-s+1}}{{\tau }_{j-s+1}^{\prime }}\prod _{r=1}^{i-1}\frac{{\sigma }_{i-r+1}}{{\sigma }_{i-r+1}^{\prime }}$$

for $i=1,\dots ,{\overline{r}}_1$, $j={\overline{s}}_{\delta }+1,\dots ,{\overline{s}}_{\delta +1}$, with

$$\begin{array}{cc}\hfill {\mathsf{\Gamma }}_{i,j}=& \sum _{r=1}^{i-1}{(-1)}^{r+1}\frac{P_j^{(i-r)}}{{\sigma }_{i-r+1}^{\prime }}\prod _{k=1}^{r-1}\frac{{\sigma }_{i-k+1}}{{\sigma }_{i-k+1}^{\prime }}\hfill \\ & +{(-1)}^{i+{\overline{s}}_{\delta }+1}\prod _{r=1}^{i-1}\frac{{\sigma }_{i-r+1}}{{\sigma }_{i-r+1}^{\prime }}\sum _{s=1}^{j-{\overline{s}}_{\delta }-1}{(-1)}^{s+1}\frac{Q_1^{(j-s-1)}}{{\tau }_{j-s+1}^{\prime }}\prod _{k=1}^{s-1}\frac{{\tau }_{j-k+1}}{{\tau }_{j-k+1}^{\prime }}.\hfill \end{array}$$

Lastly, in the fourth subdomain, from (8), we are able to obtain

$$C_{i,j}={\mathsf{\Gamma }}_{i,j}+{(-1)}^{i+j+{\overline{r}}_1+{\overline{s}}_1}{C}_{{\overline{r}}_1+1,{\overline{s}}_1+1}\prod _{s=1}^{j-{\overline{s}}_1-1}\frac{{\tau }_{j-s+1}}{{\tau }_{j-s+1}^{\prime }}\prod _{r=1}^{i-{\overline{r}}_1-1}\frac{{\sigma }_{i-r+1}}{{\sigma }_{i-r+1}^{\prime }}$$

for $i={\overline{r}}_1+1,\dots ,{\overline{r}}_2$, $j={\overline{s}}_1+1,\dots ,{\overline{s}}_2,$ with

$$\begin{array}{cc}\hfill {\mathsf{\Gamma }}_{i,j}=& \sum _{r=1}^{i-{\overline{r}}_1-1}{(-1)}^{r+1}\frac{P_j^{(i-r-1)}}{{\sigma }_{i-r+1}^{\prime }}\prod _{k=1}^{r-1}\frac{{\sigma }_{i-k+1}}{{\sigma }_{i-k+1}^{\prime }}\hfill \\ & +{(-1)}^{i+{\overline{r}}_1+1}\prod _{r=1}^{i-{\overline{r}}_1-1}\frac{{\sigma }_{i-r+1}}{{\sigma }_{i-r+1}^{\prime }}\sum _{s=1}^{j-{\overline{s}}_1-1}{(-1)}^{s+1}\frac{Q_{1+{\overline{r}}_1}^{(j-s-1)}}{{\tau }_{j-s+1}^{\prime }}\prod _{k=1}^{s-1}\frac{{\tau }_{j-k+1}}{{\tau }_{j-k+1}^{\prime }}.\hfill \end{array}$$

Generalizing, taking into account the definitions (3), and since

$$\frac{{\sigma }_{\ell }}{{\sigma }_{\ell }^{\prime }}=\frac{h_{\ell }}{h_{\ell -1}},\frac{{\tau }_{\ell }}{{\tau }_{\ell }^{\prime }}=\frac{k_{\ell }}{k_{\ell -1}},$$

we obtain (18).  □

For the special case of minimal configuration, the expression is simpler and we get lemma 3.1 of [15].

Thus, the surface $\mathcal{S}$ can be written as follows.

Theorem 3.

Let the B-spline curve network (1) satisfy the compatibility conditions C1. and C2. Then the surface (4) with isoparametric curves (1) is:

$$\mathcal{S}(u,v)={\mathcal{S}}_b(u,v)+{\mathcal{S}}_{\mathsf{\Gamma }}(u,v)+{\mathcal{S}}_{fp}(u,v),$$

(22)

where

$$\begin{array}{cc}{\mathcal{S}}_b(u,v)=\hfill & {\displaystyle \sum _{j=0}^{n+2}\left(P_j^{\left(0\right)}{B}_{0,j}(u,v)+P_j^{\left(m\right)}{B}_{m+2,j}(u,v)\right)+\sum _{i=1}^{m+1}\left(Q_i^{\left(0\right)}{B}_{i,0}(u,v)+Q_i^{\left(n\right)}{B}_{i,n+2}(u,v)\right),}\hfill \\ {\mathcal{S}}_{\mathsf{\Gamma }}(u,v)=\hfill & \underset{\mathit{excluding}\mathit{the}\mathit{free}\mathit{parameters}\mathit{in}\left(17\right)}{\underbrace{{\displaystyle \sum _{i=1}^{m+1}\sum _{j=1}^{n+1}}}}{\mathsf{\Gamma }}_{i,j}{B}_{i,j}(u,v),\hfill \\ & \\ {\mathcal{S}}_{fp}(u,v)=\hfill & {\displaystyle \sum _{\mathit{in}\mathit{each}\mathit{subdomain}\mathit{of}\mathit{type}[u_{{\overline{r}}_{\gamma }},u_{{\overline{r}}_{\gamma +1}}]\times [v_{{\overline{s}}_{\delta }},v_{{\overline{s}}_{\delta +1}}]}{C}_{1+{\overline{r}}_{\gamma },1+{\overline{s}}_{\delta }}\left(B_{1+{\overline{r}}_{\gamma },1+{\overline{s}}_{\delta }}(u,v)\right.}\hfill \\ & \left.{\displaystyle +\sum _{i,j\mathit{involved}\mathit{in}\mathit{the}\mathit{considered}\mathit{subdomain}}{(-1)}^{i+j+{\overline{r}}_{\gamma }+{\overline{s}}_{\delta }}\frac{k_j}{k_{{\overline{s}}_{\delta }+1}}\frac{h_i}{h_{{\overline{r}}_{\gamma }+1}}{B}_{i,j}(u,v)}\right)\hfill \end{array}$$

with ${\mathsf{\Gamma }}_{i,j}$ defined in (19).

Proof. 

The expression (22) follows directly from Theorem 1 and Lemma 1. □

Concerning the behaviour of the $C_{i,j}$ generation process, given in (18), we can state the following theorem related to the round-off error growth, showing that such a growth is linear.

Theorem 4.

If the sequence of partitions $\{U\times V\} of \Omega$ is A-quasi uniform, i.e., there exists a constant $A\ge 1$ such that $0<{max}_{i,j}\{h_i,k_j\}/{min}_{i,j}\{h_i,k_j\}\le A$, then the round-off error growth is linear.

We omit the proof of the above theorem because it can be obtained by following the same reasoning given in Theorem 3.3 of [15]. If the sequence of partitions $\{U\times V\}$ is uniform, i.e., $A=1$, then Theorem 4 holds.

Remark 2.

Although the surface $\mathcal{S}$ can be written as in Theorem 3, both for the case $S_{2,2}^{1,1}\left(R_{mn}\right)$ and for the case $S_2^1\left(T_{mn}\right)$, we can note that in the case $S_2^1\left(T_{mn}\right)$ the number of free parameters is one less than in the case $S_{2,2}^{1,1}\left(R_{mn}\right)$, since in $S_2^1\left(T_{mn}\right)$ the $B_{i,j}$’s are linearly dependent, the dependence relationship being [22,25]:

$$\sum _{i=1}^{m+1}\sum _{j=1}^{n+1}{(-1)}^{i+j}{h}_i{k}_j{B}_{i,j}(u,v)=0.$$

(23)

It is possible to write one of the free parameters in (22) depending on the others. In the special case of minimal configuration, we recall that (see [15]) this fact leads to the uniqueness of the interpolating surface $\mathcal{S}$. Moreover, if in (22) we consider all equal free parameters, we can use the dependence relationship (23) to obtain a unique surface interpolating the curve network also in the case of non minimal configuration.

4. Numerical Results

In this section we present several applications. We propose comparisons between surfaces interpolating the curve network both in $S_{2,2}^{1,1}\left(R_{mn}\right)$ and $S_2^1\left(T_{mn}\right)$. Moreover, we highlight, with an example, the behaviour of interpolating surfaces in $S_2^1\left(T_{mn}\right)$ as pointed out in Remark 2.

4.1. Test 1

In this first example, we show a curve network to be interpolated (Figure 5) and compare tensor product and triangulated surfaces (Figure 6). Here, we consider two subdomains of the parametric domain $\mathsf{\Omega }$, being the barrier at knot $v_2$. As remarked in Section 1, the higher degree of tensor product appears in the higher oscillations.

4.2. Test 2

Here again, we compare tensor product and triangulated surfaces in three different cases (Figure 7, Figure 8 and Figure 9), for given curve networks.

4.3. Test 3

In the following example, we show what was noticed in Remark 2: in the case of all equal free control points, we obtain again the uniqueness of the triangulated surface (Figure 10b and Figure 11b), described in [15] for a minimal configuration, while under the same hypothesis tensor product surfaces change their shape (Figure 10a and Figure 11a).

4.4. Test 4

Regarding the choice of the free parameters, it is possible to choose them such that the surface is “smooth”, i.e., the energy of a thin plate described by the surface is minimal, e.g., see [27]. Since in a generic subdomain $[u_{{\overline{r}}_{\gamma }},u_{{\overline{r}}_{\gamma +1}}]\times [v_{{\overline{s}}_{\delta }},v_{{\overline{s}}_{\delta +1}}]$ there is one free parameter, we can extend the method proposed in [15] (Section 4.2) and we can consider the following approximation

$$\mathsf{\Theta }={\int }_{u_{{\overline{r}}_{\gamma }}}^{u_{{\overline{r}}_{\gamma +1}}}{\int }_{v_{{\overline{s}}_{\delta }}}^{v_{{\overline{s}}_{\delta +1}}}[{\left({\mathcal{S}}^{uu}(u,v)\right)}^2+2{\left({\mathcal{S}}^{uv}(u,v)\right)}^2+{\left({\mathcal{S}}^{vv}(u,v)\right)}^2]\mathrm{d}u\mathrm{d}v,$$

of the thin plate energy, with

$${\mathcal{S}}^{uu}(u,v)=\frac{{\partial }^2\mathcal{S}(u,v)}{\partial u^2},{\mathcal{S}}^{uv}(u,v)=\frac{{\partial }^2\mathcal{S}(u,v)}{\partial u\partial v},{\mathcal{S}}^{vv}(u,v)=\frac{{\partial }^2\mathcal{S}(u,v)}{\partial v^2}.$$

Then the free parameter $C_{1+{\overline{r}}_{\gamma },1+{\overline{s}}_{\delta }}$ can be chosen by minimizing $\mathsf{\Theta }$ with respect to it. By setting $\frac{\mathrm{d}\mathsf{\Theta }}{\mathrm{d}{C}_{1+{\overline{r}}_{\gamma },1+{\overline{s}}_{\delta }}}=0$, we obtain

$$\begin{array}{c}{C}_{1+{\overline{r}}_{\gamma },1+{\overline{s}}_{\delta }}\hfill \\ =-h_{1+{\overline{r}}_{\gamma }}{k}_{1+{\overline{s}}_{\delta }}\frac{{\int }_{u_{{\overline{r}}_{\gamma }}}^{u_{{\overline{r}}_{\gamma +1}}}{\int }_{v_{{\overline{s}}_{\delta }}}^{v_{{\overline{s}}_{\delta +1}}}[({\mathcal{S}}_b^{uu}+{\mathcal{S}}_{\mathsf{\Gamma }}^{uu}){\mathcal{S}}_{fp}^{uu}+2({\mathcal{S}}_b^{uv}+{\mathcal{S}}_{\mathsf{\Gamma }}^{uv}){\mathcal{S}}_{fp}^{uv}+({\mathcal{S}}_b^{vv}+{\mathcal{S}}_{\mathsf{\Gamma }}^{vv}){\mathcal{S}}_{fp}^{vv}]}{{\int }_{u_{{\overline{r}}_{\gamma }}}^{u_{{\overline{r}}_{\gamma +1}}}{\int }_{v_{{\overline{s}}_{\delta }}}^{v_{{\overline{s}}_{\delta +1}}}[{\left({\mathcal{S}}_{fp}^{uu}\right)}^2+2{\left({\mathcal{S}}_{fp}^{uv}\right)}^2+{\left({\mathcal{S}}_{fp}^{vv}\right)}^2]}.\hfill \end{array}$$

(24)

Such a point is the one minimizing $\mathsf{\Theta }$. The integrals in (24) can be exactly evaluated by a composite tensor-product Gauss–Legendre quadrature formula with $3\times 3$ nodes, in the case $S_{2,2}^{1,1}\left(R_{mn}\right)$, and a composite Gaussian cubature on triangular domains (see [28]), implemented by the Matlab function triquad [29], in the case $S_2^1\left(T_{mn}\right)$.

An example of the above choice of free parameters is proposed in Figure 12, where, on the left, it is immediately clear its effect on the tensor product surface with respect to the corresponding one in Figure 13.

We note that, in this example, the surfaces in Figure 12b and Figure 13b are the same, being a case of minimal configuration.

4.5. Test 5

In this example, we compare the methods proposed here with another one that produces a Gordon-type $C^1$ biquadratic B-spline surface, [10] (Chapter 10).

In such a method, the surface is written as the Boolean sum of two skinned surfaces, interpolating the curves $C_r\left(v\right)$ and $D_s\left(u\right)$ given in (1), respectively, and one tensor product surface interpolating the points $I_{rs}$, given in (2) (for details of the construction see [10] (Chapter 10)). In order to obtain a standard B-spline representation of the resulting surface, it is necessary to make the three surfaces compatible in the B-spline sense (defined on the same knot vectors) by applying the knot refinement algorithm [10].

Given the curve network of Figure 14a, with $m=n=R=S=2$, $U=V=\{0,0,0,\frac{1}{3},\frac{2}{3},1,1,1\}$, we construct the three surfaces interpolating the network: in Figure 14b, the biquadratic Gordon-type B-spline surface; in Figure 15a, the biquadratic tensor product surface; and in Figure 15b, the quadratic triangulated surface. The free control point for the tensor product surface is determined by the minimization process (24) and we use it also for the triangulated one.

The shape of the three surfaces is comparable from the graphical point of view, but we note that the Gordon-type surface is not defined on the rectangular partition defined by the curve network knot vectors, due to the knot refinement.

5. Conclusions

In this paper, we generalize the method proposed in [15] for the definition of a B-spline surface interpolating a B-spline curve network by the introduction of free parameters. The cases of biquadratic tensor product and quadratic triangulated surfaces are studied, compared and confirmed by graphical examples. In order to obtain more flexibility in the shape modeling of such surfaces, the possible introduction of double knots in their definition is an interesting topic that we intend to address in a forthcoming paper.