summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/source/metapost/splines/splines.dtx
blob: da89dda74d69e2e9a67df96921f06e0200fd36f2 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
% \iffalse
% -------------------------------------------------------------------
%
% Copyright 2002--2004, Daniel H. Luecking
%
% Splines.mp may be distributed and/or modified under the conditions of the
% LaTeX Project Public License, either version 1.3 of this license or (at
% your option) any later version. The latest version of this license is in
%    <http://www.latex-project.org/lppl.txt>
% and version 1.3 or later is part of all distributions of LaTeX version
% 2003/12/01 or later.
%
% Splines has maintenance status "author-maintained". The Current Maintainer
% is Daniel H. Luecking. The Base Interpreter is MetaPost (or Metafont).
%
%<*driver>
\ProvidesFile{splines.dtx}
  [2005/02/05 v0.2 Metafont/post macros to compute splines.]%
\documentclass[draft]{ltxdoc}
\usepackage{docmfp}

\addtolength{\textwidth}{.5878pt}

\def\mytt{\upshape\mdseries\ttfamily}
\renewcommand\marg[1]{{\mytt \{#1\}}}
\renewcommand\oarg[1]{{\mytt  [#1]}}
\renewcommand\parg[1]{{\mytt  (#1)}}
\renewcommand{\meta}[1]{{$\langle$\rmfamily\itshape#1$\rangle$}}
\DeclareRobustCommand\cs[1]{{\mytt\char`\\#1}}
\def\prog#1{{\mdseries\scshape #1}}




\def\MF{\prog{meta\-font}}
\def\MP{\prog{meta\-post}}

\def\CMF{\prog{Meta\-font}}
\def\CMP{\prog{Meta\-post}}
\def\opt#1{{\sffamily\upshape#1}}
\def\mfc#1{{\mytt#1}}
\let\env\mfc
\let\file\mfc
\let\gbc\mfc
\renewcommand\{{{\mytt\char`\{}}
\renewcommand\}{{\mytt\char`\}}}
\renewcommand\|{${}\mathrel{|}{}$}

\makeatletter
\newcommand\bsl{{\mytt\@backslashchar}}
% Stupid lists!
\def\@listi{\leftmargin\leftmargini
  \parsep \z@ \@plus\p@ \@minus\z@
  \topsep 4\p@ \@plus\p@ \@minus2\p@
  \itemsep\parsep}
\let\@listI\@listi \@listi
\renewcommand\labelitemi{\normalfont\bfseries \textendash}
\renewcommand\labelitemii{\textasteriskcentered}
\renewcommand\labelitemiii{\textperiodcentered}
\leftmargini\parindent
% Stupid index!
\def\usage#1{\textrm{#1}}
\def\index@prologue{\section*{Index}\markboth{Index}{Index}Numbers
refer to the page where the corresponding entry is described.}
\def\IndexParms{%
  \parindent \z@ \columnsep 15pt
  \parskip 0pt plus 1pt
  \rightskip 5pt plus2em \mathsurround \z@
  \parfillskip=-5pt \small
  % less hanging:
  \def\@idxitem{\par\hangindent 20pt}%
  \def\subitem{\@idxitem\hspace*{15pt}}%
  \def\subsubitem{\@idxitem\hspace*{25pt}}%
  \def\indexspace{\par\vspace{10pt plus 2pt minus 3pt}}}
\renewcommand\routinestring{}
\renewcommand\variablestring{\space(var.)}
% Why does every command have to be indexed twice?
\renewcommand\SpecialMfpIndex[3]{\@bsphack
  \index{%
    \string#1\actualchar
    \string\verb\quotechar*\verbatimchar\string#1\verbatimchar
    #2 \encapchar usage}%
  \@esphack}
\def\close@crossref{\SpecialEscapechar{:}}
\makeatother
\def\VariableIndex#1{\SpecialMfpIndex{#1}{\variablestring}{}}
\def\RoutineIndex #1{\SpecialMfpIndex{#1}{}{}}

\title{Macros to Compute Splines\thanks{This file has version number
        \fileversion, last revised \filedate.}}
\author{Dan Luecking}
\date{\filedate}
\SpecialEscapechar{:}
\def\bslash{:}
\DisableCrossrefs
\CodelineIndex
\AlsoImplementation

\begin{document}
  \DeleteShortVerb{\|}
  \DocInput{splines.dtx}
\end{document}
%</driver>
%\fi
%
% \CheckSum{25}
% \CharacterTable
%  {Upper-case    \A\B\C\D\E\F\G\H\I\J\K\L\M\N\O\P\Q\R\S\T\U\V\W\X\Y\Z
%   Lower-case    \a\b\c\d\e\f\g\h\i\j\k\l\m\n\o\p\q\r\s\t\u\v\w\x\y\z
%   Digits        \0\1\2\3\4\5\6\7\8\9
%   Exclamation   \!     Double quote  \"     Hash (number) \#
%   Dollar        \$     Percent       \%     Ampersand     \&
%   Acute accent  \'     Left paren    \(     Right paren   \)
%   Asterisk      \*     Plus          \+     Comma         \,
%   Minus         \-     Point         \.     Solidus       \/
%   Colon         \:     Semicolon     \;     Less than     \<
%   Equals        \=     Greater than  \>     Question mark \?
%   Commercial at \@     Left bracket  \[     Backslash     \\
%   Right bracket \]     Circumflex    \^     Underscore    \_
%   Grave accent  \`     Left brace    \{     Vertical bar  \|
%   Right brace   \}     Tilde         \~}
%
% \catcode`\_=12
% \GetFileInfo{splines.dtx}
% \maketitle
%
% \StopEventually{\PrintIndex}
% A cubic spline through a set of points is a curve obtained by joining
% each point to the next with a cubic parametrized curve, where adjoining
% cubics must have matching first and second derivative at their common
% point.
%
% It is possible for \MP{} to compute the necessary controls.
% Unfortunately, the controls are not uniquely determined unless the curve
% is required to be closed. For open curves, there is need for two
% additional equations at the end points. A `relaxed spline' is produced
% if we require that the second derivative is $(0,0)$ at those points.
% For a closed curve, the equality of the first and second derivatives at
% the common beginning/ending point gives the needed additional equations.
%
% Note that this equates \emph{time} derivatives, so this works best when
% points are relatively evenly spaced and so the speed is relatively
% uniform. If points are differently spaced then the relatively slower
% speed between closely spaced points allows sharper turns without large
% second derivatives. Curves produced tend to have a more natural look,
% and relaxed splines are most suitable for smoothing data that is
% obtained by taking observations at evenly spaced times. Still, the
% technique is somewhat unstable when points are closely spaced, for
% example when a small change in the position of one point can produce a
% large change in its direction when viewed from another point.
%
% Start with version control information.
%    \begin{macrocode}
%<*package>
if known splines_fileversion: endinput fi;
string splines_fileversion;
splines_fileversion := "2006/09/25, v0.2a";
message "Loading splines.mp " & splines_fileversion;

%    \end{macrocode}
%
% \DescribeRoutine{list_to_array}
% Now for a command that takes a variable name \gbc{(suffix arr)} and
% copies a list of pairs to \mfc{arr1}, \mfc{arr2}, etc..
% The suffix must be declared by the calling program so that \mfc{arr}
% is numeric but \mfc{arr[\,]} are pairs, for example by ``\mfc{save arr; pair
% arr[\,];}''.
%    \begin{macrocode}
def list_to_array (suffix arr) (text list) =
  arr := 0;
  for _itm = list :
    arr[incr arr] := _itm;
  endfor
enddef;

%    \end{macrocode}
%
% \DescribeRoutine{compute_spline}
% In this command we generate the equations common to all cubic
% splines: the equality of derivatives at all interior points.
% This command accepts three suffixes: \gbc{points}, \gbc{pr}, and
% \gbc{po} which should represent previously declared arrays of pairs.
% \gbc{points} is the array of points to be connected, and must be known.
% Arrays \gbc{pr} and \gbc{po} \emph{must} be unknown and will hold the
% computed control points. See the code of \mfc{mkrelaxedspline} for
% an example of how this was arranged for the variables \mfc{rs_pr} and
% \mfc{rs_po}.
%
% Contrary to the previous (unreleased) version, \mfc{compute_spline}
% takes a boolean argument and appends the same equations at the first (=
% last) point.
%
% \DescribeRoutine{mkrelaxedspline}
% The first of these macros appends the necessary additional equations to
% get zero second derivatives at the ends.
% \DescribeRoutine{mkclosedspline}
% The second simply calls \mfc{compute_spline} with the boolean set to
% true. Both return the computed path. In theory the knowledgeable user
% can call \gbc{compute_spline (false)}, append a choice of equations for
% the ends, and then call \gbc{mksplinepath (false)}.
%
% \DescribeRoutine{dospline}
% This version accepts a list of pairs and produces a spline through
% them. It simply stores the list in an array and calls the appropriate
% version that operates on an array.
%    \begin{macrocode}
def compute_spline (expr closed) (suffix points, pr, po) =
  % interior equations:
  for j= 2 upto points - 1 :
    % equate first derivatives:
    po[j] + pr[j] = 2 points[j];
    %   and second derivatives:
    pr[j+1] + 2 pr[j] = 2 po[j] + po[j-1];
  endfor
  % for a closed curve, the first and last are also interior:
  if closed:
   po 1 + pr 1 = 2 points 1;
   po[points] + pr[points] = 2 points[points];
   pr 2 + 2 pr 1 = 2 po 1 + po[points];
   pr 1 + 2 pr[points] = 2 po[points] + po[points-1];
  fi
enddef;

vardef mksplinepath (expr closed) (suffix points, pr, po) =
  points1..controls po1 and
  for j = 2 upto points if not closed: -1 fi:
    pr[j]..points[j]..controls po[j] and
  endfor
    if closed: pr 1..cycle else: pr[points]..points[points] fi
enddef;

vardef mkrelaxedspline (suffix pnts) =
  save rs_pr, rs_po;
  pair rs_po[], rs_pr[];
  % Equate second derivative to zero at both end points
  rs_pr 2 + pnts 1 = 2 rs_po 1 ;
  pnts[pnts] + rs_po[pnts-1] = 2 rs_pr[pnts];
  compute_spline (false) (pnts, rs_pr, rs_po);
  mksplinepath   (false) (pnts, rs_pr, rs_po)
enddef;

vardef mkclosedspline (suffix pnts) =
  save cs_pr, cs_po;
  pair cs_pr[], cs_po[];
  compute_spline (true) (pnts, cs_pr, cs_po);
  mksplinepath   (true) (pnts, cs_pr, cs_po)
enddef;

vardef dospline (expr closed) (text the_list) =
  save _sp; pair _sp[];
  list_to_array (_sp) (the_list);
  if closed :
    mkclosedspline (_sp)
  else:
    mkrelaxedspline (_sp)
  fi
enddef;

%    \end{macrocode}
% The above computations produce a $2$-dimensional spline. A $1$-dimensional
% cubic spline would be a function $f(t)$ with numeric values rather
% than pair values. Such are often used to interpolate functions. That is,
% given pairs $(x\sb j,y\sb{j})$, and assuming they lie on the graph of
% some function (generally unknown), fill in the graph with $y = f(x)$
% where $f$ is a cubic function of $x$ in each interval $x\sb j < x < x\sb
% {j+1}$, making sure that the resulting graph is as smooth as possible at
% the points $x\sb j$.
%
% The requirements on our $2$-dimensional path are the following:
% \begin{enumerate}
%   \item The $j$th link should connect $(x\sb{j},y\sb{j})$ to $(x\sb{j+1},
%       y\sb{j+1})$.
%   \item The $x$-part of that link should increase linearly from $x\sb{j}$ to
%     $x\sb{j+1}$ as $t$ goes from $0$ to $1$.
%   \item The $y$-part should be a cubic $y = f(x)$.
%   \item The $x$-derivatives $df/dx$ and $d^2f/dx^2$ should match at the
% connecting points.
% \end{enumerate}
%
% Two necessary equations for converting between $x$ and $t$ coordinates
% are:
% \begin{equation}\label{first}
%         x = x\sb{j} + t \Delta x\sb{j}
% \end{equation}
% (where $\Delta x\sb{j} = x\sb{j+1} - x\sb{j}$) and
% \begin{equation}\label{second}
%       \frac{df}{dt} = \frac{dx}{dt}\frac{df}{dx} = (x\sb{j+1} -
%       x\sb{j}) \frac{df}{dx}.
% \end{equation}
% Thus we want to choose controls so that (\ref{first}) is maintained and
% so that $x$-derivatives match. It turns out that this requires controls
% at
% \begin{equation}
% \begin{array}{c}
%   (x\sb{j}, y\sb{j}) - (\Delta x\sb{j-1}, s\sb{j} \Delta x\sb{j-1})/3\\
%   (x\sb{j}, y\sb{j}) + (\Delta x\sb{j}  , s\sb{j} \Delta x\sb{j}  )/3
% \end{array}
% \end{equation}
% where $s\sb{j}$ is the slope (derivative) at $x\sb{j}$. These control
% points will produce matching slopes regardless of the values
% chosen for the $s\sb j$. To get matching second derivatives we need the
% same conditions as in parametric splines. But those equations simplify
% to the form:
% \begin{displaymath}
%  s\sb{j+1} dx\sb{j} - 2s\sb{j} (dx\sb{j} + dx\sb{j-1}) +
%       s\sb{j-1}dx\sb{j-1}
%           =  3y\sb{j+1} - 3y\sb{j-1}.
% \end{displaymath}
% As with 2-D splines there can be almost any equations at the end points.
% For a relaxed spline we equate the second derivatives to 0. To get a
% periodic function, we equate the slope and second derivative at
% beginning to those at the end. This makes it possible to put a shifted
% copy of the graph with starting point at the end of the original and
% have the same smoothness at that connection as at the other points.
%
% \DescribeRoutine{compute_fcnspline}
% This issues the equation for the slopes (array \gbc{sl} of
% \emph{unknown} numerics). The array \gbc{points} contains the $(x,y)$
% values and \gbc{dx} is a temporary numeric array which will be
% overwritten if known.
%
% \DescribeRoutine{mkfcnsplinepath}
% This simply assembles the path from the information computed by the
% above equations (and the extra equations given in the calling command).
%
% \DescribeRoutine{mkrelaxedfcnspline}
% This sets up arrays for the  \gbc{dx} and \gbc{sl} parameters of
% \mfc{compute_fcnspline}, emit the necessary endpoint equations (zero
% second derivatives) and calls the previous two routines.
%
% \DescribeRoutine{mkperiodicfcnspline}
% This does the same as the previous command, but the endpoint equations
% make the first and second derivatives at the start equal to those at
% the end.
%
% \DescribeRoutine{fcnspline}
% Finally, this command copies a list of pairs into an array and calls the
% appropriate command to process them.
%    \begin{macrocode}
def compute_fcnspline (suffix points, dx, sl) =
  % Get delta_x:
  for j = 1 upto points - 1: dx[j] := xpart (points[j+1]-points[j]);
  endfor
  for j=2 upto points - 1:
    sl[j + 1] * dx[j] + 2sl[j]*(dx[j] + dx[j-1]) + sl[j-1]*dx[j-1]
        = 3*ypart(points[j+1] - points[j-1]);
  endfor
enddef;

vardef mkfcnsplinepath (suffix points, dx, sl) =
  points1..controls (points1 + (1, sl1)*dx1/3) and
    for j = 2 upto points - 1:
      (points[j] - (1, sl[j])*dx[j-1]/3) ..points[j]..
        controls (points[j] + (1,sl[j])*dx[j]/3) and
    endfor
    (points[points] - (1,sl[points])*dx[points-1]/3)..points[points]
enddef;

vardef mkperiodicfcnspline (suffix pnts) =
  save _sl, _dx; numeric _dx[], _sl[];
  compute_fcnspline (pnts, _dx, _sl);
  % periodicity equations:
  _sl 1 = _sl[pnts];
  _sl 2 * _dx 1 + 2 _sl 1 * _dx 1 + 2 _sl[pnts] * _dx[pnts-1]
      + _sl[pnts-1] * _dx[pnts-1]
          = 3 * ypart(pnts[2] - pnts[pnts-1]);
  mkfcnsplinepath  (pnts, _dx, _sl)
enddef;

vardef mkrelaxedfcnspline (suffix pnts) =
  save _sl, _dx; numeric _dx[], _sl[];
  compute_fcnspline (pnts, _dx, _sl);
  % relaxation equations.
  _sl 2 * _dx 1 + 2 _sl1 * _dx 1 = 3 * ypart(pnts2 - pnts1);
  _sl[pnts-1] * _dx[pnts-1] + 2 _sl[pnts] * _dx[pnts-1]
        = 3 * ypart(pnts[pnts] - pnts[pnts-1]);
  mkfcnsplinepath  (pnts, _dx, _sl)
enddef;

vardef fcnspline (expr periodic) (text the_list) =
  save _fs; pair _fs[];
  list_to_array (_fs) (the_list);
  if periodic:
    mkperiodicfcnspline (_fs)
  else:
    mkrelaxedfcnspline (_fs)
  fi
enddef;

%</package>
%    \end{macrocode}
% \clearpage
%\Finale