\listfiles \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{lmodern} \usepackage[scaled=.82]{luximono}%requires T1+textcomp \usepackage[T1]{fontenc} \newcommand\textSMC[1]{{\SMC #1}} \newcommand\acro[1]{\textSMC{#1}\@} \DeclareRobustCommand\cs[1]{\texttt{\char`\\#1}} \usepackage[cal=boondoxo]{mathalfa} \usepackage{textcomp} %\usepackage{url} %\def\url@ttstyle{% % \@ifundefined{selectfont}{\def\UrlFont{\tt}}{\def\UrlFont{\ttfamily\small}}} \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$ ($n \ge 2$) is a curve parametrized by the interval $[0,n]$, which is, except in degenerate cases, $C^2$-continuous (that is, has continuous curvature) and is on each interval $[k-1,k]$ ($0\}}] defines a clipping path 12{\tt pt} wide around the B-spline curve with control points {\tt B0}..{\tt B5}, then draws the {\tt } clipped to that path. \item[\cs{bspcurvenodes\{P\}\{Q\}}] creates a node sequence {\tt Q0} {\tt Q1},... from the position data in the arrays {\tt P.X}, {\tt P.Y}. \end{description} Details and examples are provided below. \section{Relaxed, Open B-spline} The algorithm to generate such a curve from a sequence of control points $B_0$, $\cdots$, $B_n$ is as follows: \begin{itemize} \item The curve starts at $B_0$ and ends at $B_n$. (Important: $n \ge 2$.) \item Divide each line segment $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).| \subsection{The \cs{bspcurvepoints} macros} There are two macros that provide for B-spline curves essentially the same functionality as the \verb|\pscurvepoints| macro from {\tt pstricks-add}. (That macro takes as input a parametric curve and constructs as output (at the PostScript level) arrays of data associated with the curve: the positions of points along the curve, the increment from the previous point and a normal vector to the curve. The principal uses for such data are (i) the \verb|\pspolylineticks| macro from {\tt pstricks-add}, which allows placement of ticks and other marks along a curve that has been approximated by a polyline; (ii) the \cs{polyIntersections} macro from \textsf{pst-node}, which allows you to find the points of intersection of the curve (approximated by a polyline) and an arbitrary line.) The macros \begin{verbatim} \bspcurvepoints{}{}{} \bspcurvepointsE{}{}{} \end{verbatim} work, in the first case, for a relaxed, uniform B-spline curve, and in the second, for such a curve with its initial and final segments removed, corresponding to the output from \verb|\psBsplineE| rather than \verb|\psBspline|. In both cases, you may set the keyword {\tt plotpoints} (default value: $50$) to change the number of sample points on each B\'ezier component. This will result in the construction of PostScript arrays with indices from $0$ to $n=$\textsf{num of segments}$\times$\textsf{(plotpoints-1)}. After running \begin{verbatim} \pnodes{B}(1,2)(3,-1)(4,1)(6,2)% define B0..B3 \bspcurvepoints[plotpoint=11]{B}{3}{P} \end{verbatim} the following PostScript arrays are created, each indexed from 0 to 30: \begin{verbatim} P.X, P.Y (position) PNormal.X, PNormal.Y (normal vector) PDelta.X, PDelta.Y (increment from previous position) \end{verbatim} and these may be used in the usual way to create nodes. For example, \begin{verbatim} \pnode(! P.X 8 get P.Y 8 get){Q} \pnode(! PNormal.X 8 get PNormal.Y get){Dir} \psrline(Q)(1cm;{(Dir)}) \end{verbatim} places {\tt Q} at the position on the curve with index 8, defines {\tt Dir} to be a normal vector at that point, then draws a line from {\tt Q} of length {\tt 1cm} in the direction of that normal. \subsection{Setting nodes on a B-spline curve} To set a node at parameter value $t$ on a B-spline curve after running \cs{bspcurvepoints}{\tt[E]}, call the macro \begin{verbatim} \bspNode{}{}{}{} \end{verbatim} For example, if I have constructed a B-spline curve using control points $B_0$,$\dots$,$B_5$, then \verb|\bspNode{B}{5}{2.1}{Q}| defines a node named $Q$ at $t=2.1$. The macro \cs{bspcurvenodes\{P\}\{R\}} creates a node sequence {\tt R0}..{\tt Rn} at the locations specified by the arrays {\tt P.X}, {\tt P.Y}. (Those arrays must first have been created with one of the \cs{bspcurvepoints} macros.) \subsection{B-spline function curves} By this we mean an open B-spline curve which is the graph of a function $y=f(x)$ and whose orientation is toward the right. It is not analytically simple to specify a formula for $f$ in most cases, and to compute $y$ from $x$ involves (a) finding the index of the B\'ezier segment containing $x$; (b) solving the cubic $x(t)=x$ for $t$; (c) substituting in $y(t)$. The package provides a macro to perform these calculations after generating the data using \cs{bspcurvepoints}{\tt[E]}: \begin{verbatim} \bspfnNode{}{}{}{} \end{verbatim} \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. If you used \verb|\pnodes| to do this, it would have constructed a macro \verb|\Snodecount| to store the value $100$. 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}{101} \end{verbatim} with any keywords options you wish. {\bf IMPORTANT:} The macro \verb|\psBsplineInterpC| modifies the node sequence {\tt S}, setting \verb|S101=S0|, and changing \verb|\Snodecount| to take the value $101$. This is convenient when you use the construction: \begin{verbatim} \pnodes{S}()()()()% sets \Snodecount to 3 \psBsplineInterpC{S}{\Snodecount}% constructs B-spline control pts SB0..SB4 \psBsplineNodesC{SB}{\Snodecount} \end{verbatim} 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} A B-spline curve can in many cases provide a good function interpolation mechanism, but the result is not guaranteed to be the graph of a function. \begin{verbatim} \begin{pspicture}(-.5,-.5)(6,4) \psdots(0,3.5)(1,.5)(3,2.5)(4,0)(5,2)(6,.5)% \pnodes{S}(0,3.5)(1,.5)(3,2.5)(4,0)(5,2)(6,.5)% S0..S5 \psBsplineInterp{S}{5}% construct SB0..SB5 \psBsplineNodes{SB}{5}% draw B-spline with control pts SB0..SB5 \bspcurvepoints[plotpoints=10]{SB}{5}{P} % construct the PS arrays \bspFnNode{SB}{5}{4.5}{QQ}% node QQ on curve at x=4.5 \psdot[linecolor=red](QQ)% \psaxes(0,0)(-.5,-.5)(6,4) \end{pspicture} \end{verbatim} \vspace{1pc} \begin{center} \begin{pspicture}(-.5,-.5)(6,4) \psdots(0,3.5)(1,.5)(3,2.5)(4,0)(5,2)(6,.5)% \pnodes{S}(0,3.5)(1,.5)(3,2.5)(4,0)(5,2)(6,.5)% S0..S5 \psBsplineInterp{S}{5}%SB0..SB5 \psBsplineNodes{SB}{5}% draw B-spline with control pts SB0..SB5 \bspcurvepoints[plotpoints=10]{SB}{5}{P} % construct the PS arrays \bspFnNode{SB}{5}{4.5}{QQ}% node QQ on curve at x=4.5 \psdot[linecolor=red](QQ)% \psaxes(0,0)(-.5,-.5)(6,4) \end{pspicture} \end{center} \vspace{12pt} \begin{verbatim} \documentclass{article} \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} \section{Thick B-spline curves} Inspired by the package {\tt pst-thick}, we provide a similar option for curves generated as B-spline interpolations. The new macro that accomplishes this is \begin{verbatim} \thickBspline#1#2#3#4 %#1=root | #2=nsegments | #3=thickness | #4=items to clip \end{verbatim} which expects the following data. \begin{itemize} \item A node sequence. This can be constructed with a command like \begin{verbatim} \pnodes{S}(0,0)(5,1)(4,4)(1,3)% \end{verbatim} which declares nodes S0..S3, and sets the macro \verb|\Snodecount| to 3. \item An interpolation command, such as \begin{verbatim} \psBsplineInterp{S}{\Snodecount}% \end{verbatim} which creates a framework of B-spline control points {\tt SB0..SB3}. \item Create the interpolating curve and the B\'ezier control points for its components, with names like {\tt SBR0..SBR2,SBL1..SBL3} etc, using \begin{verbatim} \psBsplineNodes[linestyle=none]{SB}{\Snodecount}% \end{verbatim} (The {\tt [linestyle=none]} may be omitted if you want the curve to show.) \item Create a clipping path of specified thickness around the interpolating curve and place graphics to be clipped: \begin{verbatim} \thickBspline[plotpoints=50,linestyle=none]{S}{3}{20pt}% {\psline[linecolor=red,linestyle=solid](0,0)(6,6)}% \end{verbatim} (The {\tt [linestyle=none]} controls whether the clipping path is rendered, and {\tt plotpoints} controls the number of subdivisions of each B\'ezier component. Its default value is 50.) \end{itemize} The clipping path is drawn by default positively oriented so that objects are clipped to its interior. By specifying the keyword {\tt reverseclip}, the clipping path will be reversed so that objects are clipped to the exterior. \begin{verbatim} \documentclass{article} \usepackage{pstricks} \usepackage{pst-bspline,pstricks-add} \begin{document} \begin{pspicture}[showgrid=true](-.5,-.5)(6,5) \pnodes{S}(1,0)(5,1)(4,4)(1,3)% \psdots(1,0)(5,1)(4,4)(1,3)% \psBsplineInterp{S}{\Snodecount}% \psBsplineNodes[linestyle=none]{SB}{\Snodecount}% \thickBspline[plotpoints=50,linestyle=none]{S}{3}{20pt}% {\psframe[fillstyle=crosshatch](-1,-1)(6,6)}% \end{pspicture} \end{document} \end{verbatim} \vspace{12pt} \begin{center} \begin{pspicture}[showgrid=true](-.5,-.5)(6,5) \pnodes{S}(1,0)(5,1)(4,4)(1,3)% \psdots(1,0)(5,1)(4,4)(1,3)% \psBsplineInterp{S}{\Snodecount}% \psBsplineNodes[linestyle=none]{SB}{\Snodecount}% \thickBspline[plotpoints=50,linestyle=none]{S}{\Snodecount}{20pt}{\psframe[fillstyle=crosshatch](-1,-1)(6,6)}% \end{pspicture} \end{center} \vspace{12pt} The \verb|\thickBspline| macro works as expected in the closed (periodic) case, taking advantage of automatic incrementing of the nodecount. Note that \verb|\thickBspline| interprets thickness as visual, unaffected by possible differences between {\tt xunit} and {\tt yunit}. \begin{verbatim} \documentclass{article} \usepackage{pstricks} \usepackage{pst-bspline,pstricks-add} \begin{document} \psset{yunit=1.5cm} \begin{pspicture}[showgrid=true](-.5,-.5)(6,5) \pnodes{S}(1.5,0)(5,1)(4,4)(1,3)% \psBsplineInterpC{S}{3}% % defines nodes SB0, SB1, SB2, SB3, SB4 --- the Bspline control points % increments \Snodecount by 1 for future macros % Don't use C form of \psBsplineNodes with this new \Snodecount \psBsplineNodes[linestyle=none,showpoints=false]{SB}{\Snodecount}% % Constructs the Bezier control points SBR0, SBL1, SBR1, etc \thickBspline[linestyle=none]{S}{\Snodecount}{22pt}% {\psframe[fillstyle=vlines](-1,-1)(6,6)}% \end{pspicture} \end{document} \end{verbatim} \vspace{12pt} \begin{center} \psset{yunit=1.5cm} \begin{pspicture}[showgrid=true](-.5,-.5)(6,5) \pnodes{S}(1.5,0)(5,1)(4,4)(1,3)% \psBsplineInterpC{S}{3}% % defines nodes SB0, SB1, SB2, SB3 --- the Bspline control points % increments \Snodecount by 1 for future macros \psBsplineNodes[linestyle=none,showpoints=false]{SB}{\Snodecount}% % Constructs the Bezier control points SBR0, SBL1, SBR1, etc \thickBspline[linestyle=none]{S}{\Snodecount}{22pt}% {\psframe[fillstyle=vlines](-1,-1)(6,6)}% %{\psframe[fillstyle=solid, fillcolor=lightgray](-1,-1)(6,6)}% \end{pspicture} \end{center} \end{document}