[comp.lang.postscript] Bezier interpolation & BSplines

peterson-john@CS.Yale.EDU (John C. Peterson) (06/08/90)

Since not much practical information has been posted regarding Shamim's
original request, I thought I would add a little.  First, there is a
big difference between an interpolating curve and an approximating
curve.  An interpolating curve must pass through a set of given points,
while an approximating one should just pass near the control points.
For drawing purposes, approximating curves are more commonly used.  It
is possible to construct a good interpolant with cubic Bezier patches,
but the problem is underspecified.  Either initial and final slopes need
to be included or an initial slope and in initial curvature are needed.

While there are many different approximating curves, the BSpline is both
simple to construct (especially in terms of Bezier segments) and has
a locality property which minimizes the overall change in the curve when a
control point is moved.

Digging into my old postscript library, here is a routine to draw BSplines.
(This is written in PLisp, a Lisp front end to Postscript.  I can't abide
programming in RPN and screwing around with the stack).


;;  bspline.pl        Author: John Peterson (peterson-john@cs.yale.edu)

;;  The following routine constructs a BSpline from a set of control points.
;;  The input is an array of any length containing points.  Each point is
;;  an array of length 2.

;;  Source: foggy memories of Bob Barnhill's CAGD class.

;;  Make a triple knot at each end by repeating the first and last point
;;  two extra times.  A closed figure (periodic BSpline) would be obtained
;;  by instead repeating the last points before the first and first after
;;  the last.

(defun bspline (points) 0
   (let ((a (array (+ (length points) 4)))
	 (n (length points)))
     (setf (elt a 0) (elt points 0))
     (setf (elt a 1) (elt points 0))
     (dotimes (i n)
	 (setf (elt a (+ i 2)) (elt points i)))
     (setf (elt a (+ n 2)) (elt points (1- n)))
     (setf (elt a (+ n 3)) (elt points (1- n)))
     (moveto (elt (elt points 0) 0) (elt (elt points 0) 1))
     (internal-bspline a)))

;; Draw the bspline by picking up 6 points at a time from the control

(defun internal-bspline (pts) 0
   (let ((n (length pts)))
     (dotimes (i (- n 5))
	  (bspline-curve
		(elt pts i)
		(elt pts (+ i 1))
		(elt pts (+ i 2))
		(elt pts (+ i 3))
		(elt pts (+ i 4))
		(elt pts (+ i 5))))))

;; Generate a single bezier segment.  Lots of interpolation goes on here.

(defun bspline-curve (p1 p2 p3 p4 p5 p6) 0
  (let ((l1 (seg-length p1 p2))
	(l2 (seg-length p2 p3))
	(l3 (seg-length p3 p4))
	(l4 (seg-length p4 p5))
	(l5 (seg-length p5 p6)))
    (let* (
	   (tp1 (interp p3 p4 (/ l2 (+ l2 l3 l4))))
	   (tp2 (interp p3 p4 (/ (+ l2 l3) (+ l2 l3 l4))))
	   (tp3 (interp p4 p5 (/ l3 (+ l3 l4 l5))))
	   (tp4 (interp tp2 tp3 (/ l3 (+ l3 l4)))))
      (curveto (elt tp1 0) (elt tp1 1) (elt tp2 0) (elt tp2 1)
	       (elt tp4 0) (elt tp4 1)))))

;; Compute the distance between two points

(defun seg-length (p1 p2) 1
  (sqrt (+ (* (- (elt p1 0) (elt p2 0)) (- (elt p1 0) (elt p2 0)))
	   (* (- (elt p1 1) (elt p2 1)) (- (elt p1 1) (elt p2 1))))))

;; Interpolate between 2 points.  The interval (0,1) in alpha is
;; mapped onto (p1,p2).

(defun interp (p1 p2 alpha) 1
  (vector (+ (elt p1 0) (* alpha (- (elt p2 0) (elt p1 0))))
          (+ (elt p1 1) (* alpha (- (elt p2 1) (elt p1 1))))))


----------------------------------------------------------------
For those without a Plisp compiler, here is the same function in
(absolutely unreadable) real Postscript.

/BSPLINE--FR-1 2 dict def /**MAIN**--FR-2
1 dict def /INTERNAL-BSPLINE--FR-3 1 dict def
/**MAIN**--FR-4 1 dict def /BSPLINE--frame 1
dict def /BSPLINE-CURVE--FR-5 5 dict def
/BSPLINE-CURVE--FR-6 1 dict def /BSPLINE-CURVE--FR-7
1 dict def /BSPLINE-CURVE--frame 6 dict def
/BSPLINE-CURVE--FR-8 1 dict def /INTERP--frame 3
dict def /BSPLINE-CURVE--FR-9 1 dict def
/INTERNAL-BSPLINE--frame 1 dict def /SEG-LENGTH--frame
2 dict def /a 0 def /points 0 def /pts
0 def /bspline { BSPLINE--frame begin /points
exch def /points load length 4 add array
/points load length BSPLINE--FR-1 begin /n
exch def /a exch def /a load 0 /points
load 0 get put /a load 1 /points load
0 get put 0 1 /n load 1 sub {
**MAIN**--FR-2 begin /i exch def /a load /i
load 2 add /points load /i load get put
end } for /a load /n load 2 add /points
load /n load 1 sub get put /a load /n
load 3 add /points load /n load 1 sub
get put /points load 0 get 0 get /points
load 0 get 1 get moveto /a load
internal-bspline end end } def /internal-bspline {
INTERNAL-BSPLINE--frame begin /pts exch def /pts
load length INTERNAL-BSPLINE--FR-3 begin /n exch def
0 1 /n load 5 sub 1 sub {
**MAIN**--FR-4 begin /i exch def /pts load
/i load get /pts load /i load 1 add get
/pts load /i load 2 add get /pts load /i
load 3 add get /pts load /i load 4 add
get /pts load /i load 5 add get
bspline-curve end } for end end } def
/bspline-curve { BSPLINE-CURVE--frame begin /p6 exch
def /p5 exch def /p4 exch def /p3 exch def
/p2 exch def /p1 exch def /p1 load /p2 load
seg-length /p2 load /p3 load seg-length /p3 load
/p4 load seg-length /p4 load /p5 load seg-length
/p5 load /p6 load seg-length BSPLINE-CURVE--FR-5
begin /l5 exch def /l4 exch def /l3 exch def
/l2 exch def /l1 exch def /p3 load /p4 load
/l2 load /l2 load /l3 load add /l4 load add
div interp BSPLINE-CURVE--FR-6 begin /tp1 exch def
/p3 load /p4 load /l2 load /l3 load add /l2
load /l3 load add /l4 load add div interp
BSPLINE-CURVE--FR-7 begin /tp2 exch def /p4 load
/p5 load /l3 load /l3 load /l4 load add /l5
load add div interp BSPLINE-CURVE--FR-8 begin /tp3
exch def /tp2 load /tp3 load /l3 load /l3
load /l4 load add div interp BSPLINE-CURVE--FR-9
begin /tp4 exch def /tp1 load 0 get /tp1
load 1 get /tp2 load 0 get /tp2 load
1 get /tp4 load 0 get /tp4 load 1 get
curveto end end end end end end } def
/seg-length { SEG-LENGTH--frame begin /p2 exch def
/p1 exch def /p1 load 0 get /p2 load 0
get sub /p1 load 0 get /p2 load 0 get
sub mul /p1 load 1 get /p2 load 1 get
sub /p1 load 1 get /p2 load 1 get sub
mul add sqrt end } def /interp { INTERP--frame
begin /alpha exch def /p2 exch def /p1 exch
def [ /p1 load 0 get /alpha load /p2 load
0 get /p1 load 0 get sub mul add /p1
load 1 get /alpha load /p2 load 1 get
/p1 load 1 get sub mul add ] end }
def

----------------------------------------------------------------

Hope this code is of use.

    John Peterson
    peterson-john@cs.yale.edu

    Think Granite!