Added chord length parameterisation.

ChangeLog

@@ -1,5 +1,11 @@
 2011-07-30  Rudy Neeser  <rudy.neeser@gmail.com>
 
+	* spline.lisp (initialize-instance): Added chord length
+	parameterisation, a non-uniform parameterisation.
+
+	* basis.lisp: Moved the b-spline basis code into basis.lisp from
+	utility.lisp.
+
 	* spline-interpolation.lisp (spline-interpolation): Added a method
 	that will take a series of points and return a b-spline that
 	interpolates all of the points.

basis.lisp

@@ -67,4 +67,289 @@
 					   -1/2  0    1/2  0
 					    0    2/2  0    0)))
 
-;;;-----------------------------------------------------------------------
+;;;-----------------------------------------------------------------------
+;;; Classes for b-spline basis functions
+
+(defclass b-spline-knots ()
+  ((knots :initform nil
+	  :initarg :knots
+	  :type (or null (simple-array double-float))
+	  :accessor knots
+	  :documentation "An ascending collection of knots, stored with
+	  repetition.")
+   (multiplicity :initform nil
+		 :initarg :multiplicity
+		 :type (or null (simple-array fixnum))
+		 :accessor multiplicity
+		 :documentation "The multiplicity of the above knots.")
+   (knot-count :initform 0
+	       :type fixnum
+	       :reader knot-count
+	       :documentation "The number of knots, taking in to account multiplicity.")
+   (no-multiplicity :initform nil
+		    :reader no-multiplicity-p
+		    :documentation "Indicates whether there is a
+		    multiplicity or not in the knots.")
+   (uniform :initform nil
+	    :reader uniform-p
+	    :documentation "A boolean indicating whether this is a
+	    uniform sequence of knots or not."))
+  (:documentation "Some docs. You will likely find it more convenient
+  to construct instance of this class using MAKE-KNOTS. Notice that
+  B-SPLINE-BASIS requires a phantom knot on either end of the knot
+  sequence. MAKE-KNOTS will automatically add those knots for you."))
+
+(defgeneric calculate-knot-count (knot-data)
+  (:documentation "Returns the number of knots in the data, taking in
+  to account multiplicity. Sets KNOT-COUNT accordingly.")
+  (:method ((knot-data b-spline-knots))
+    (setf (slot-value knot-data 'knot-count)
+	  (with-accessors ((multiplicity multiplicity)) knot-data
+	    (loop
+	       for m across multiplicity
+	       sum m)))))
+
+(defgeneric calculate-uniformity (knot-data)
+  (:documentation "Calculates whether the knot-data is uniform or
+  not. Sets UNIFORM-P accordingly.")
+  (:method ((knot-data b-spline-knots))
+    (with-accessors ((knots knots)
+		     (multiplicity multiplicity)) knot-data
+      (setf (slot-value knot-data 'no-multiplicity)
+	    (every #'(lambda (item)
+		       (= item 1))
+		   multiplicity))
+      (setf (slot-value knot-data 'uniform) 
+	    (and (no-multiplicity-p knot-data)
+		 (let ((distance (- (aref knots 1) (aref knots 0))))
+		   (loop
+		      for i from 0 below (length knots)
+		      for j from 1 below (length knots)
+		      always (equivalent distance (- (aref knots j)
+						     (aref knots i))))))))))
+
+(defmethod initialize-instance :after ((knot-data b-spline-knots) &key)
+  (with-accessors ((knots knots)
+		   (multiplicity multiplicity)) knot-data
+    (unless (= (length knots)
+	       (length multiplicity))
+      (error 'l-math-error
+	     :format-control "There are a different number of knots (~A) to specified multiplicity (~A)."
+	     :format-arguments (list (length knots) (length multiplicity))))
+    (calculate-knot-count knot-data)
+    (calculate-uniformity knot-data)))
+
+(defmethod (setf knots) :after (item (knot-data b-spline-knots))
+  (calculate-uniformity knot-data)
+  (calculate-knot-count knot-data))
+
+(defmethod (setf multiplicity) :after (item (knot-data b-spline-knots))
+  (calculate-uniformity knot-data)
+  (calculate-knot-count knot-data))
+
+(defun make-knots (knots &key
+		   (multiplicity (loop
+				    repeat (length knots)
+				    collect 1))
+		   (add-phantoms t))
+  "Creates a knot sequence. Notice that this will add some phantom
+knots to the beginning and end of the sequence, unless ADD-PHANTOMS is
+false. If MULTIPLICITY, a list of multiplicity values, is not given,
+each knot is assumed to be a simple (non-repeated) knot."
+  (declare (type list knots multiplicity))
+  (let* ((distance (- (second knots) (first knots)))
+	 (knots (if add-phantoms
+		    (append (cons (- (first knots) distance) knots)
+			    (list (+ (first (last knots)) distance)))
+		    knots))
+	 (multiplicity (if add-phantoms
+			   (append (cons 1 multiplicity)
+				   (list 1))
+			   multiplicity)))
+    (make-instance 'b-spline-knots
+		   :knots  (make-array (length knots) :element-type 'double-float
+				       :initial-contents (mapcar #'(lambda (knot)
+								     (coerce knot 'double-float))
+								 knots))
+		   :multiplicity (make-array (length multiplicity) :element-type 'fixnum
+					     :initial-contents (mapcar #'(lambda (multi)
+									   (coerce multi 'fixnum))
+								       multiplicity)))))
+
+(defgeneric all-knots (knots)
+  (:documentation "Given a knots data structure, this will return a
+  list of all knots, including multiplicities.")
+  (:method ((knots-data b-spline-knots))
+    (with-accessors ((knots knots)
+		     (multiplicity multiplicity)) knots-data
+      (let ((count (length knots)))
+	(labels ((handle-index (index count accum)
+		   (cond
+		     ((< count (aref multiplicity index))
+		      (handle-index index (1+ count) (cons (aref knots index) accum)))
+		     (t
+		      accum)))
+		 (determine-knots (index accum)
+		   (cond
+		     ((< index count)
+		      (determine-knots (1+ index) (append (handle-index index 0 nil)
+							  accum)))
+		     (t
+		      accum))))
+	  (nreverse (determine-knots 0 nil)))))))
+
+(defmethod print-object ((knots b-spline-knots) stream)
+  (print-unreadable-object (knots stream :identity t :type t)
+    (format stream "~{~A~^, ~}" (all-knots knots))))
+
+(defgeneric low-parameter (knot-data degree)
+  (:documentation "Given a knot sequence and the degree of the spline,
+  this returns the lowest possible parameter value the spline will
+  accept.")
+  (:method ((knots b-spline-knots) (degree integer))
+    (get-ith-knot knots degree)))
+
+(defgeneric high-parameter (knot-data degree)
+  (:documentation "Given a knot sequence and the degree of the spline,
+  this returns the highest possible parameter value the spline will
+  aceept.")
+  (:method ((knots b-spline-knots) (degree integer))
+    (get-ith-knot knots (- (knot-count knots) (1+ degree)))))
+
+(defun number-needed-knots (num-points degree)
+  "Given the number of points making up a spline, and the degree of
+  the spline, this will tell us the number of required knots."
+  (+ 1 degree num-points))
+
+
+(defgeneric get-ith-knot (knot-data i &optional offset)
+  (:documentation "Returns the ith knot, taking in to account
+  multiplicity. OFFSET should be positive number that array indices
+  are offset by.")
+  (:method ((knot-data b-spline-knots) (i integer) &optional (offset 0))
+    (with-accessors ((knots knots)
+		     (multiplicity multiplicity)) knot-data
+      (when (>= i (knot-count knot-data))
+	(error 'l-math-error :format-control "The index is larger than the number of knots."))
+      (cond
+	((no-multiplicity-p knot-data)
+	 ;; Can speed this up in the uniform case.
+	 (aref knots (+ i offset)))
+	(t
+	 (aref knots (loop
+			for index from 0 below (length multiplicity)
+			for m across multiplicity
+			sum m into count
+			while (<= count (+ offset i))
+			finally (return index))))))))
+
+(defun find-starting-knot-index (knot-data degree parameter)
+  "Given knot data, the degree of a spline, and a parameter, this
+  locates the first index of the knots which will be used to define
+  the point on the given spline. Returns first the b-spline-basis
+  family to be used, and then the offset index of the knot."
+  (check-type knot-data b-spline-knots)
+  (with-accessors ((knots knots)
+		   (multiplicity multiplicity)) knot-data
+    (let ((index (loop
+		    for knot across knots
+		    for mult across multiplicity
+		    for index = 0 then (+ index mult)
+		    while (<= knot parameter)
+		    finally (return (- index degree)))))
+      (values (- index degree)
+	      index))))
+
+(defun b-spline-basis (knot-data degree family parameter &key (offset degree) fast-spline-index)
+  "Ask for the value of a given basis function. This is based on
+Farin's 'Curves and Surfaces for Computer Aided Geometric Design. DEGREE is the degree of the curve.
+
+This has some pre-calculated basis functions for uniform knot
+secquences. You can access these by specifying an index using
+FAST-SPLINE-INDEX. This is an integer between 0 and DEGREE, and
+specifies which of the basis functions should be called."
+  (declare (type b-spline-knots knot-data)
+	   (type fixnum degree)
+	   (type number parameter)
+	   (type (or null fixnum) fast-spline-index))
+  (when (minusp degree)
+    (error 'l-math-error :format-control "The degree of a b-spline may not be negative."))
+  (let ((current (get-ith-knot knot-data family offset))
+	(before (get-ith-knot knot-data (1- family) offset))
+	(nth-after (get-ith-knot knot-data (+ family degree) offset))
+	(nth-after-1 (get-ith-knot knot-data (+ family (1- degree)) offset)))
+    (cond
+      ((zerop degree)
+       (if (and (<= before parameter)
+		(< parameter current))
+	   1
+	   0))
+      ((and fast-spline-index
+	    (uniform-p knot-data)
+	    (= degree 2)
+	    (equivalent (abs (- current before)) 1))
+       ;; Some code to speed up the quadratic case. First, we need
+       ;; to shift the parameter to begin at 0.
+       (let ((u (abs (second (multiple-value-list (truncate parameter))))))
+      	 (ecase fast-spline-index
+      	   (2 (* 1/2 (expt u 2)))
+      	   (1 (+ 1/2 u (- (expt u 2))))
+      	   (0 (+ 1/2 (- u) (* 1/2 (expt u 2)))))))
+      ((and fast-spline-index
+	    (uniform-p knot-data)
+      	    (= degree 3)
+	    (equivalent (abs (- current before)) 1))
+       ;; Some code to speed up the quadratic case. First, we need
+       ;; to shift the parameter to begin at 0.
+       (let ((u (abs (second (multiple-value-list (truncate parameter))))))
+      	 (ecase fast-spline-index
+      	   (3 (* 1/6 (expt u 3)))
+      	   (2 (+ 1/6 (* 1/2 u) (* 1/2 (expt u 2)) (* -1/2 (expt u 3))))
+      	   (1 (+ 4/6 (- (expt u 2)) (* 1/2 (expt u 3))))
+      	   (0 (+ 1/6 (* -1/2 u) (* 1/2 (expt u 2)) (* -1/6 (expt u 3)))))))
+      (t
+       (+ (if (equivalent (- nth-after-1 before) 0)
+	      0
+	      (* (/ (- parameter before)
+		    (- nth-after-1 before))
+		 (b-spline-basis knot-data (1- degree) family parameter :offset offset)))
+	  (if (equivalent (- nth-after current) 0)
+	      0
+	      (* (/ (- nth-after parameter)
+		    (- nth-after current))
+		 (b-spline-basis knot-data (1- degree) (1+ family) parameter :offset offset))))))))
+
+(defun chord-length-parameterisation (all-points degree)
+  "Produces a non-parameteric parameterisation for b-splines."
+  (cond
+    ((= (length all-points) (1+ degree))
+     (list -3 -1 0 1 2 3))
+    (t
+     (let* ((domain-points (rest (butlast all-points)))
+	    (result (loop
+		       for count from 0 below (- (length domain-points) 2)
+		       for i in domain-points
+		       for j in (rest domain-points)
+		       sum (lm:euclidean-distance i j) into sum
+		       collecting sum into collection
+		       finally (return collection)))
+	    (last (+ (first (last result))
+		     (euclidean-distance (first (last domain-points))
+					 (nth (- (length domain-points) 2) domain-points))))
+	    (result (append result (list last))))
+       (mapcar #'(lambda (item)
+		   (/ item last))
+	       (append (loop
+			  for i from (1- degree) above 0
+			  with distance = (- (second result)
+					     (first result))
+			  collecting (- (* distance i)))
+		       (cons 0 result)
+		       (loop
+			  for i from 1 below degree
+			  with distance = (- (first (last result))
+					     (first (last (butlast result))))
+			  collecting (+ last (* distance i)))))))))
+
+;;;-----------------------------------------------------------------------
+

spline.lisp

@@ -412,34 +412,41 @@ geometry points, or enought points."))
 
   (make-instance 'lm:b-spline :degree 3 :points *list-of-points* :uniform t)
 
+  A non-uniform parameterisation that attempts to equalise the
+  distance travelled over parameter sections is chord length parameterisation:
+  (make-instance 'lm:b-spline :degree 3 :points *list-of-points* :uniform t)
+
   Knots may be specified with the :knots argument, and should be
   constructed using lm:MAKE-KNOTS."))
 
-(defmethod initialize-instance :after ((spline b-spline) &key uniform)
+(defmethod initialize-instance :after ((spline b-spline) &key uniform chord-length)
   "UNIFORM: if true, will automatically create a set of uniform knots,
   based on the b-spline's degree and number of points in the b-spline
   polygon."
   (with-accessors ((points b-spline-points)
 		   (degree b-spline-degree)) spline
-    (when uniform
-      (let ((num-knots (number-needed-knots (length points) degree)))
+    (let ((num-knots (number-needed-knots (length points) degree)))
+      (when uniform
 	(setf (b-spline-knots spline)
 	      (make-knots (loop
 			     repeat num-knots
 			     for i from (- degree)
 			     collect i)
-			  (loop
-			     repeat num-knots
-			     collect 1)
-			  :add-phantoms nil))))
-    (when (null (b-spline-knots spline))
-      (error 'l-math-error
-	     :format-control "No knots have been provided, or requested to be generated (eg, see :uniform keyword)"))
-    (when (< (knot-count (b-spline-knots spline))
-	     (number-needed-knots (length points) degree))
-      (format t "Knot-count: ~A~%" (knot-count (b-spline-knots spline)))
-      (error 'l-math-error :format-control "This spline requires at least ~A knots."
-	     :format-arguments (list (number-needed-knots (length points) degree))))))
+			  :add-phantoms nil)))
+      (when chord-length
+	(setf (b-spline-knots spline)
+	      (make-knots (chord-length-parameterisation points degree))))
+      (when (and uniform chord-length)
+	(error 'l-math-error
+	       :format-control "Two different methods have been requested to produce knots."))
+      (when (null (b-spline-knots spline))
+	(error 'l-math-error
+	       :format-control "No knots have been provided, or requested to be generated (eg, see :uniform keyword)"))
+      (when (< (knot-count (b-spline-knots spline))
+	       num-knots)
+	(format t "Knot-count: ~A~%" (knot-count (b-spline-knots spline)))
+	(error 'l-math-error :format-control "This spline requires at least ~A knots."
+	       :format-arguments (list num-knots))))))
 
 (defmethod minimum-parameter ((spline b-spline))
     (low-parameter (b-spline-knots spline) (b-spline-degree spline)))