A short note on the generalized “smoothstep” function.

August 29, 2017

The Wikipedia page on Renderman’s “smoothstep” function,

\mathrm{smoothstep}(x)=3x^2-2x^3

features an entry on higher-order generalizations of this function, where each succeeding member has more of its higher order derivatives being set to zero at the endpoints 0 and 1.

Among other things, a formula for these generalized functions is presented, which I find rather unwieldy. I have thus written this short entry to present a much conceptually simpler formula.

One can use e.g. confluent divided differences to derive the following (IMNSHO cleaner) formula for n-th order smoothstep \mathrm{S}_n(x):

\displaystyle \mathrm{S}_n(x)=x^{n+1}\sum\limits_{k=0}^n\binom{n+k}{k}(1-x)^k

Here is a Mathematica demonstration of this identity for the first 21 members:

And @@ Table[InterpolatingPolynomial[{PadRight[{{0}, 0}, n + 2],
                                      PadRight[{{1}, 1}, n + 2]}, x] ==
             x^(n + 1) Sum[Binomial[n + k, k] (1 - x)^k, {k, 0, n}]
             // Simplify, {n, 0, 20}]

which should yield True.

In implementations for other languages (e.g. C++, JavaScript), one would want to use a suitable modification of Horner’s rule to evaluate the polynomial factor expressed in terms of 1-x, as well as exploiting the recurrence

\dbinom{n+k+1}{k+1}=\left(1+\dfrac{n}{k+1}\right)\dbinom{n+k}{k}

for generating the coefficients.

Advertisements

A compact Hilbert curve routine

January 30, 2017

Here is a relatively compact Mathematica routine for generating the b-th iterate of a Hilbert curve in [0,1]^n. The algorithm is due to Skilling:

HilbertCurve =
     Compile[{{x, _Integer}, {b, _Integer}, {n, _Integer}},
             Module[{t = BitXor[x, Quotient[x, 2]], p = 2, k, q, r, xx},
                    xx = Total[Table[BitAnd[Quotient[t, 2^(r (n - 1) + k)], 2^r],
                                     {r, b - 1, 0, -1}, {k, n - 1, 0, -1}]];
                    Do[q = p - 1;
                       Do[t = xx[[k]];
                          If[BitAnd[t, p] != 0,
                             xx[[1]] = xx[[1]] ~BitXor~ q,
                             t = BitAnd[t ~BitXor~ xx[[1]], q];
                             xx[[{1, k}]] = BitXor[xx[[{1, k}]], t]],
                          {k, n, 2, -1}];
                       If[BitAnd[xx[[1]], p] != 0,
                          xx[[1]] = xx[[1]] ~BitXor~ q];
                       p *= 2, {r, b - 1}];
                    2 xx/(p - 1) - 1], RuntimeAttributes -> {Listable}]

I’ve seen a variety of routines for generating the Hilbert curve before, but none were as compact or as general as this one. If you want a version that generates exact rational coordinates, it is easy to take out the requisite parts from the Compile[] function.

Here are the 2D and 3D versions of the curve:

With[{p = 6}, 
     Graphics[{ColorData[97, 1],
               Line[HilbertCurve[Range[0, 2^(2 p) - 1], p, 2]]}, Frame -> True]]

h2d

With[{p = 3}, 
     Graphics3D[Tube[HilbertCurve[Range[0, 2^(3 p) - 1], p, 3]], Axes -> True]]

h3d

Here’s a perspective projection of a four-dimensional Hilbert curve:

With[{p = 3, k = 2, f = 4, d = 3/2, r = 1/8}, 
     Graphics3D[{Directive[EdgeForm[], CapForm[None], GrayLevel[1/5], 
                           Glow[ColorData["Legacy", "SteelBlue"]],
                           Specularity[ColorData["Legacy", "Gold"], 16]], 
                 Tube[(f Delete[#, k])/(d - Extract[#, k]) & /@ 
                      HilbertCurve[Range[0, 2^(4 p) - 1], p, 4], r]}, 
                Boxed -> False, Lighting -> "Neutral"]]

h4d

(This visualization is admittedly a bit messy-looking; if you have suggestions for making a nice display of the four-dimensional version, tell me about it in the comments!)