Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

Faster de Boor #1

Merged
merged 1 commit into from

2 participants

@alang9

I noticed that your implementation of de Boor seemed to be quadratic in the number of control points rather than in the degree. This version is quadratic in the number of degrees instead, but it can still probably be made much faster.

Note also that this produces "natural" b-splines, which take the value 0 before and after the first and last knot resp. (I happen to prefer this, but it might not be what you want for your library).

Cheers,
Alex

@alang9 alang9 Faster version of deBoor. It assumes that the splines are `natural', …
…i.e. they take the value 0 before and after the last spline.
bb7fe7f
@mokus0
Owner

Thank you! I had intended to do this but never got around to it. The existing implementation is indeed quadratic in the number of control points. As I recall, laziness makes the cost proportional to (number of control points * number of knots less than x + degree^2), so not only is it worse than it should be, it's not even uniform across the spline.

I don't have a strong opinion about the behavior at the ends, but I do find the "natural" evaluation strategy inconvenient for clamped splines - they approach the last control point but then suddenly snap to zero when they reach it. I understand that it is mathematically elegant in some ways, but from a practical standpoint it really seems ugly to have to treat the last control point specially to avoid a bizarre (to the user) discontinuity.

My current inclination is to have 'evalBSpline' continue to have the same meaning it does now (so code that depends on the current behavior doesn't break) and add a new function 'evalNaturalBSpline'. That's a fairly long name though, so I'm open to other suggestions if you prefer something shorter.

@mokus0 mokus0 merged commit b86755d into from
@mokus0
Owner

I've restructured the code a bit - the old 'deBoor' output is used by the code to split b-splines, so I've separated your changes into a new function 'slice' that extracts the needed part of the spline before passing it to 'deBoor'. You may wish to proof-read the changes - my ad-hoc quickCheck tests indicate it behaves identically for all inputs it generated, but it's always possible that it missed some corner cases.

@alang9

Hi mokus0,

The changes look great to me, I won't have time to test them for another few days, but I'll tell you if I find anything.

Thanks a lot!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Feb 6, 2012
  1. @alang9

    Faster version of deBoor. It assumes that the splines are `natural', …

    alang9 authored
    …i.e. they take the value 0 before and after the last spline.
This page is out of date. Refresh to see the latest.
Showing with 30 additions and 8 deletions.
  1. +30 −8 src/Math/Spline/BSpline/Internal.hs
View
38 src/Math/Spline/BSpline/Internal.hs
@@ -5,6 +5,7 @@ module Math.Spline.BSpline.Internal
import Math.Spline.Knots
import Data.Monoid
+import Data.Ratio
import Data.Vector as V
import Data.VectorSpace
import Prelude as P
@@ -40,22 +41,27 @@ insertKnot spline x = spline
-- extend [1..5] -> [1,1,2,3,4,5,5]
extend vec
| V.null vec = V.empty
- | otherwise = V.cons (V.head vec) (V.snoc vec (V.last vec))
+ | otherwise = V.cons (V.head vec) (V.snoc vec (V.last vec))
-deBoor spline x = go us (controlPoints spline)
+deBoor spline x = go uLo . vtake (deg + 1) . vdrop (l - deg) $
+ controlPoints spline
where
+ l = maybe (-1) pred $ V.findIndex (> x) us
+ deg = degree spline
+ zero = fromRational (0 % 1)
+
us = knotsVector (knotVector spline)
-
+ uLo = stake (deg + 1) $ sdrop (l - deg) $ us
-- Upper endpoints of the intervals are the same for
-- each row in the table (they just line up differently
-- with the lower endpoints):
- uHi = V.drop (degree spline + 1) us
-
- -- On each pass, the lower endpoints of the
- -- interpolation intervals advance and the new
+ uHi = sdrop (l + 1) us
+
+ -- On each pass, the lower endpoints of the
+ -- interpolation intervals advance and the new
-- coefficients are given by linear interpolation
-- on the current intervals:
- go us ds
+ go us ds
| V.null ds = []
| otherwise = ds : go uLo ds'
where
@@ -63,6 +69,22 @@ deBoor spline x = go us (controlPoints spline)
ds' = V.zipWith4 (interp x) uLo uHi
ds (V.tail ds)
+ -- Try to take n, but if there's not enough, pan the rest with 0s
+ vtake n xs
+ | n <= V.length xs = V.take n xs
+ | otherwise = xs V.++ V.replicate (n - V.length xs) zeroV
+ stake n xs
+ | n <= V.length xs = V.take n xs
+ | otherwise = xs V.++ V.replicate (n - V.length xs) 0
+
+ -- Try to drop n, but if n is negative, pan the beginning with 0s
+ vdrop n xs
+ | n >= 0 = V.drop n xs
+ | otherwise = V.replicate (-n) zeroV V.++ xs
+ sdrop n xs
+ | n >= 0 = V.drop n xs
+ | otherwise = V.replicate (-n) 0 V.++ xs
+
interp x x0 x1 y0 y1
| x < x0 = y0
| x >= x1 = y1
Something went wrong with that request. Please try again.