\documentclass[dvips,11pt]{article} \usepackage{amsmath} \usepackage{amsthm} \usepackage{graphicx} \usepackage{pstricks} \usepackage{multido,pst-node,pst-bspline,pstricks-add} \usepackage{amssymb} \usepackage[parfill]{parskip} \usepackage{hyperref} \hyphenation{Post-Script} %\date{} % Activate to display a given date or no date \begin{document} \begin{center}{\Large Cubic B-splines Using PSTricks\\[12pt] \large Michael Sharpe\\[10pt] msharpe@ucsd.edu} \end{center} A cubic uniform B-spline curve with control points $B_0 \ldots B_n$ is a curve parametrized by the interval $[0,n]$, which is $C^2$-continuous (that is, has continuous curvature) and is on each interval $[k-1,k]$ given by a cubic B\'ezier curve whose control points are derived from the $(B_k)$. These curves are discussed in any reasonably modern text on Numerical Analysis. One easily accessible source is the UCLA lecture notes of Kirby Baker: \noindent\url{http://www.math.ucla.edu/~baker/149.1.02w/handouts/dd_splines.pdf} I'll focus on two special cases: (i) relaxed, uniform B-splines; (ii) periodic, uniform B-splines. Uniform refers to the condition mentioned in the first paragraph: each B\'ezier sub-curve is parametrized by an interval of length~1. Relaxed means that the curvature at the endpoints $t=0, t=n$ is zero. Periodic means in effect that the $B_i$ repeat periodically, and the curve generated is a closed curve. \section{Relaxed, Open B-spline} The algorithm has the following steps. \begin{itemize} \item The curve starts at $B_0$ and ends at $B_n$. \item Divide each line $B_{k-1}B_k$ into equal thirds, with subdivision points labeled $R_{k-1}$, $L_k$ respectively, so that $B_k$ has $L_k$ as its immediate neighbor to the left, and $R_k$ as its immediate neighbor to the right. \item For $0}{} \psBsplineNodesC{}{} \psBsplineNodesE{}{} \end{verbatim} corresponding to the macros \verb|\psBspline|, \verb|\psBsplineC| and \verb|\psBsplineE|. The difference is that the macros with {\tt Nodes} in the name have as arguments the root node name and the last index, rather than the list of points. For example, with the above definition of {\tt P} in force, \verb|\psBsplineNodes{P}{2}| has exactly the same effect as \verb|\psBspline(2,1.5)(3,4)(5,1).| \section{B-spline Interpolation} This is the inverse problem. Being given points $(S_k)_{0\le k\le n}$, the goal is to produce the B-spline control points $B_k$ leading to the points $S_k$, so that the associated B-spline curve interpolates the $S_k$. \subsection{Open curve} We discuss first the case of an open, uniform B-spline curve with relaxed endpoints. According to the discussion above, we have to solve the equations \begin{align*} B_0&=S_0\\ B_0+4B_1+B_2&=6S_1\\ B_1+4B_2+B_3&=6S_2\\ \cdots&\\ B_{n-2}+4B_{n-1}+B_n&=6S_{n-1}\\ B_n&=S_n \end{align*} for the $B_k$. In matrix form, this becomes the tridiagonal system \[\begin{pmatrix}4&1\\ 1&4&1\\ &1&4&1\\ &&\cdots&&1\\ &&&1&4\end{pmatrix} \begin{pmatrix}B_1\\B_2\\B_3\\ \cdots\\ B_{n-1}\end{pmatrix}=§ \begin{pmatrix}6S_1-S_0\\6S_2\\6S_3\\ \cdots\\6S_{n-1}-S_{n}\end{pmatrix} \] The LU decomposition of the tridiagonal matrix may be seen to take the form \[ \begin{pmatrix}1\\ m_1&1\\ &m_2&1\\ &&\cdots\\ &&&m_{n-2}&1\end{pmatrix} \begin{pmatrix}m_1^{-1}&1\\ &m_2^{-1}&1\\ &&m_3^{-1}&1\\ &&&\cdots&1\\ &&&&m_{n-1}^{-1}\end{pmatrix} \] where $m_1=1/4$, $m_{k+1}=1/(4-m_k)$ for $k=1,\cdots,n-2$. The solution of the original system is therefore accomplished in two steps, introducing intermediate points $(R_k)$, by (in pseudo-code) \begin{verbatim} R_1=6*S_1-S_0 for i=2 to n-2 R_i=6*S_i-m_{i-1}* R_{i-1} R_{n-1}=(6*S_{n-1}-S_n)-m_{n-2}*R_{n-2} B_{n-1}=m_{n-1}*R_{n-1} for i=n-2 downto 1 B_i=m_i*(R_i-B_{i+1}) \end{verbatim} The code for the \verb|\psBsplineInterp| command uses this algorithm to solve for the $B_k$ as nodes, except that in order to save node memory, the $B$ nodes are substituted in place for the $R$ nodes, so that, for example, the first step becomes \verb|B_1=6*S_1-S_0|. Assuming you have previously defined nodes {\tt S0} $\cdots ${\tt S4}, \begin{verbatim} \psBsplineInterp{S}{4} \end{verbatim} will construct a sequence {\tt SB0} $\cdots ${\tt SB4} of nodes at the B-spline control points for the relaxed, uniform cubic B-spline interpolating the {\tt Sk}, and this curve may then be rendered with the command \begin{verbatim} \psBsplineNodes{SB}{4} \end{verbatim} If you don't care about keeping track of the internal operations and names for nodes, you may generate the curve directly with, for example, \begin{verbatim} \psbspline(0,0)(.5,.1)(1.5,.6)(2.5,1.4)(3.5,1.8)(4.5,1.7)% (5.8,1.0)(7.5,.25)(10,0) \end{verbatim} \subsection{Closed (periodic) case} We turn now to the periodic uniform B-spline curve interpolating $n$ points $S_0$,...,$S_{n-1}$. Extend the sequence periodically with period $n$, so that $S_n=S_0$, $S_{n+1}=S_1$, $S_{-1}=S_{n-1}$, and so on. In order to find the periodic control points $B_k$, we have to solve the $n$ equations \begin{align*} B_n+4B_1+B_2&=6S_1\\ B_1+4B_2+B_3&=6S_2\\ \cdots&\\ B_{n-2}+4B_{n-1}+B_n&=6S_{n-1}\\ B_{n-1}+4B_n+B_1&=6S_n \end{align*} for the $B_k$, $1\le k\le n$. In matrix form, this becomes the system \[\begin{pmatrix}4&1&&&1\\ 1&4&1\\ &1&4&1\\ &&\cdots&&1\\ 1&&&1&4\end{pmatrix} \begin{pmatrix}B_1\\B_2\\B_3\\ \cdots\\ B_{n}\end{pmatrix}=§ \begin{pmatrix}6S_1\\6S_2\\6S_3\\ \cdots\\6S_{n}\end{pmatrix} \] Let $(x_k,y_k)=6S_k$. We perform Gaussian elimination on the matrix \[\begin{pmatrix}4&1&&&1&x_1&y_1\\ 1&4&1&&&x_2&y_2\\ &1&4&1&&x_3&y_3\\ &&\cdots&&1\\ 1&&&1&4&x_n&y_n\end{pmatrix} \] As in the previous case, let $m_1=0.25$, $m_k=1/(4-m_{k-1})$ for $k\ge 2$. The factor $m_k$ will be the multiplier of row $k$ after the previous row operation, in order to normalize the row. These are the steps in the procedure. \begin{itemize} \item Initialize: multiply row 1 by $m_1$ so that its first entry (1,1) is 1. Replace $x_1$ by $m_1 x_1$ and $y_1$ by $m_1 y_1$. Entry $(1,n)$ is $m_1$. \item Subtract new row 1 from row 2 and multiply the resulting row by $m_2$. The leading entry (2,1) becomes $1$. Entry $(2,n)$ becomes $-m_1m_2$, and $x_2, y_2$ are updated to $m_2(x_2- x_1)$, $m_2(y_2-y_1)$. The superdiagonal entry (2,3) is the only other non-zero entry, and its new value is $m_2$. \item Subtract new row 1 from row $n$, so that its leading entry $(n,2)$ is $-m_1$. \item Subtract new row 2 from row 3 and multiply the result by $m_3$. The leading entry (3,3) becomes $1$ and the entry $(3,n)$ becomes $m_1m_2$, with $x_3, y_3$ updating to $m_3(x_3-x_2)$, $m_3(y_3-y_2)$. The superdiagonal entry (3,4) is now $m_3$. \item Subtract new row 2 times $-m_1$ from row $n$, whose leading entry $(n,3)$ is now $m_1m_2$. \item Continue in this way until row $n-2$ has been subtracted as above from row $n-1$, multiplying the result by $m_{n-1}$, and a suitable multiple has been subtracted from row $n$. The leading entry of row $n-1$ (column $n-1$) is $1$ and its $n^{\text{th}}$ entry is $1-(-1)^{n}m_1\cdots m_{n-2}$. Row $n$ has leading entry in column $n-1$, equal to $1$. \item Finally, subtract an appropriate multiple of row $n-1$ from row $n$ so that row $n$ has leading entry in column $n$. The resulting matrix is upper triangular, and we may now substitute back starting from the last row to give a complete reduction. \end{itemize} Here are the steps in pseudocode. We keep track of row $n$ with the array $b_k$, column $n$ with the array $c_k$. The indices for both run from 1 to $n$. \begin{verbatim} m(1)=0.25 for k=2 to n-1 m(k)=1/(4-m(k-1)) b(1)=1 b(n-1)=1 b(n)=4 c(n-1)=1% don't need c(n), =b(n) %multiply first row by m1 c(1)=m(1) x(1)=m(1)*x(1) y(1)=m(1)*y(1) for k=2 to n-1 %subtract normalized row k-1 from row k, renormalize row k c(k)=m(k)*(c(k)-c(k-1))%note that initially, c(k)=0 for 1}{} \end{verbatim} You must previously have defined a sequence, say {\tt S0} $\cdots$ {\tt S100} of \verb|\pnode|s that you plan to interpolate with a closed curve. Then \begin{verbatim} \psBsplineInterpC{S}{100} \end{verbatim} constructs the sequence {\tt SB0} $\cdots$ {\tt SB100} of B-spline control points (appending~{\tt B} to the root name) for a closed curve interpolating {\tt S0} $\cdots$ {\tt S100}, which may then be rendered with the command \begin{verbatim} \psBsplineNodesC{SB}{100} \end{verbatim} with any keywords options you wish. The following example illustrates that there is a difference between \verb|\psccurve| and B-spline interpolation, the former having a rounder appearance. Generally speaking, B-spline interpolation comes closer to minimizing the average curvature. \begin{verbatim} \documentclass{article} \usepackage{pstricks} \usepackage{pst-bspline,pstricks-add} \begin{document} \begin{pspicture}[showgrid=true](-.5,-.5)(6,5) \pnodes{P}(0,1)(2,0)(5,2)(6,4)(4,5)(2,4) \psBsplineInterpC{P}{5} \psBsplineNodesC*[linecolor=gray!40]{PB}{5} \psccurve[linecolor=red,showpoints=true](0,1)(2,0)(5,2)(6,4)(4,5)(2,4) \end{pspicture} \end{document} \end{verbatim} \vspace{1pc} \begin{center} Slight difference between psccurve and B-spline interpolation\\ \vspace*{2pc} \begin{pspicture}[showgrid=true](-.5,-.5)(6,5) \pnodes{P}(0,1)(2,0)(5,2)(6,4)(4,5)(2,4) \psBsplineInterpC{P}{5} \psBsplineNodesC*[linecolor=gray!40]{PB}{5} \psccurve[linecolor=red,showpoints=true](0,1)(2,0)(5,2)(6,4)(4,5)(2,4) \end{pspicture} \end{center} \newpage \begin{verbatim} \documentclass{article} \usepackage{graphicx} \usepackage{pstricks} \usepackage{pst-bspline,pstricks-add} \begin{document} \psset{unit=.25in} \begin{pspicture}[showgrid=true](-.5,-.5)(6,5) \pnodes{P}(0,1)(2,0)(5,2)(6,4)(4,5)(2,4) \pnode(3,3){C} \multido{\ra=0+.05,\rb=1+.05,\i=30+1}{40}{% \psBsplineC*[linecolor=blue!\i!brown]{B}% ([nodesep=\ra]{C}P0)([nodesep=\ra]{C}P1)% ([nodesep=\ra]{C}P2)([nodesep=\ra]{C}P3)% ([nodesep=\ra]{C}P4)([nodesep=\ra]{C}P5)} \end{pspicture} \end{document} \end{verbatim} \vspace{1pc} \begin{center} \psset{unit=.25in} \begin{pspicture}[showgrid=true](-.5,-.5)(6,5) \pnodes{P}(0,1)(2,0)(5,2)(6,4)(4,5)(2,4) \pnode(3,3){C} \multido{\ra=0+.05,\rb=1+.05,\i=30+1}{40}{% \psBsplineC*[linecolor=blue!\i!brown]{B}% ([nodesep=\ra]{C}P0)([nodesep=\ra]{C}P1)% ([nodesep=\ra]{C}P2)([nodesep=\ra]{C}P3)% ([nodesep=\ra]{C}P4)([nodesep=\ra]{C}P5)} \end{pspicture} \end{center} \end{document}