MATLAB’s current default colormap, after years of (ab)using `jet`

, is called `parula`

. This colormap was named after a gaudy tropical bird.

## On emulating the “parula” colormap in Mathematica

March 5, 2018## A short note on the generalized “smoothstep” function.

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

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 and .

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 -th order smoothstep :

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 , as well as exploiting the recurrence

for generating the coefficients.

## A compact Hilbert curve routine

January 30, 2017Here is a relatively compact *Mathematica* routine for generating the -th iterate of a Hilbert curve in . 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]]

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

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"]]

(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!)

## A short note on Costa’s minimal surface

May 14, 2012The other day, I finally managed to simplify Alfred Gray’s parametric equations for Costa’s minimal surface. I might edit this post later, with details on how to manipulate the Weierstrass elliptic functions that show up in the equations, but enjoy these for now:

Here, and are respectively the Weierstrass elliptic function, its derivative, and the Weierstrass zeta function, with invariants and , and and are the usual beta and gamma functions. The invariants are the ones corresponding to the semi-periods and . The parameter ranges are and .

I was very delighted at my result, that I decided to make a stylized plot of the surface in *Mathematica*, using Perlin’s simplex noise for the coloring:

## Bunnies, teapots, and noise

March 16, 2012I should be able to have regular Internet access soon. For now, I’d like to share some of the experiments I did with coloring 3D objects with the help of Ken Perlin’s simplex noise in *Mathematica*. (Relatedly, I’ve found that the implementation of improved Perlin noise supplied in the *Mathematica* help file has a few subtle flaws, which I was able to fix. That would be for another blog entry, though.)

Here are some of the results of coloring the Utah teapot with simplex noise:

Here are some results of coloring the Stanford bunny with simplex noise:

I never would have thought that *Mathematica* would ever be able to do graphics as fancy as this. I’ll only note that none of these used the `Texture[]`

construct; both models in fact used simplex noise as solid textures (i.e. an `RGBColor[]`

associated with each point in the model).

I might write a (series of) more descriptive post(s) on these things much later; for now, these pictures, entirely generated in *Mathematica* (no post-processing at all!) will have to suffice.

## On making simulated landscapes

March 1, 2012This entry is based on an idea by Paul Bourke for generating simulated landscapes through the application of the two-dimensional discrete Fourier transform to an array of uniformly-distributed random numbers. I was also inspired by some of the color gradients featured in GRASS, and decided to see if I could make a somewhat realistic color gradient for depicting a simulated landscape. Here is the color gradient that resulted from my many experiments:

`topoFake[u_?NumericQ] :=`

Blend[{{0, RGBColor[0, 1/10, 1/4]}, {1/5, RGBColor[0, 3/5, 3/5]}, {1/5, RGBColor[206/255, 197/255, 12/17]}, {1/4, RGBColor[39/85, 9/26, 5/44]}, {11/40, RGBColor[14/25, 27/40, 2/9]}, {23/80, RGBColor[9/14, 19/25, 12/85]}, {1/3, RGBColor[7/15, 29/42, 1/255]}, {3/8, RGBColor[3/10, 3/5, 5/98]}, {17/40, RGBColor[0, 38/85, 0]}, {19/40, RGBColor[0, 6/17, 0]}, {21/40, RGBColor[0, 3/10, 0]}, {23/40, RGBColor[0, 1/4, 0]}, {5/8, RGBColor[0, 1/8, 0]}, {5/8, RGBColor[39/85, 9/26, 5/44]}, {43/60, RGBColor[3/5, 25/54, 1/6]}, {97/120, RGBColor[19/23, 7/10, 1/3]}, {151/180, RGBColor[7/8, 3/4, 13/36]}, {313/360, RGBColor[39/44, 15/17, 2/5]}, {9/10, RGBColor[13/17, 9/11, 16/51]}, {9/10, RGBColor[1/2, 1/2, 1/2]}, {91/100, RGBColor[1/2, 1/2, 1/2]}, {91/100, RGBColor[1, 50/51, 50/51]}, {19/20, RGBColor[1, 50/51, 50/51]}, {19/20, RGBColor[3/4, 3/4, 3/4]}, {24/25, RGBColor[3/4, 3/4, 3/4]}, {24/25, RGBColor[1, 50/51, 50/51]}, {1, RGBColor[1, 50/51, 50/51]}}, u] /; 0

I wanted to try it out on a simpler function first, and decided to use the `peaks()`

function of MATLAB as a test case:

`Plot3D[3 (1 - x)^2 E^(-x^2 - (y + 1)^2) - 10 (x/5 - x^3 - y^5) E^(-x^2 - y^2) - 1/3 E^(-(x + 1)^2 - y^2), {x, -3, 3}, {y, -3, 3},`

Axes -> None, BoundaryStyle -> None, Boxed -> False,

ColorFunction -> (topoFake[#3] &), Mesh -> False,

PerformanceGoal -> "Quality", PlotPoints -> 85, PlotRange -> All,

ViewPoint -> {-1.3, -2.4, 2.}]

Not too shabby, eh? Here’s what results if we try out the gradient on a modification of Bourke’s idea:

`With[{n = 256},`

BlockRandom[SeedRandom[20, Method -> "Legacy"];

ListPlot3D[Abs[InverseFourier[

Fourier[RandomReal[{0, 1}, {n, n}]] * Table[

Sech[(i/n - 0.5)/0.025] * Sech[(j/n - 0.5)/0.025],

{i, n}, {j, n}]]],

Axes -> None, Boxed -> False, BoundaryStyle -> None,

BoxRatios -> {1, 1, 1/10}, ColorFunction -> (topoFake[#3] &),

Mesh -> False, PerformanceGoal -> "Quality"]]]

The resulting archipelago looks a tad chaotic, but it’s a start. I’m still looking for a way to modify Bourke’s technique to produce bigger “islands”.

## Faking marbles

February 29, 2012Still on the matter of Perlin noise, I finally managed to figure out how to wrap Perlin noise on a sphere (that is, as a `Texture[]`

). It was a bit inconvenient that it seems (at least, to me) that not one of the usual references on Perlin noise mentioned the explicit formulae needed for a proper spherical embedding of Perlin noise. Making use of the functions already defined in the doc page for `Compile[]`

, here are two *Mathematica* functions you can use for generating textures that can be wrapped on a sphere:

`sphericalPerlin1 = With[{octaves = 4},`

Compile[{{theta, _Real}, {phi, _Real}, {amplitude, _Real},

{frequency, _Real}, {gain, _Real}, {lacunarity, _Real}},

Module[{noiseVal = 0.0, x, y, z, sf, freq = frequency,

amp = amplitude},

sf = Sin[phi];

x = sf*Cos[theta] + 1; y = sf*Sin[theta] + 1; z = Cos[phi] + 1;

Do[

noiseVal += amp*signedNoise[x*freq, y*freq, z*freq];

freq *= lacunarity; amp *= gain,

{octaves}

];

noiseVal

], CompilationOptions -> {"InlineExternalDefinitions" -> True}

]]

`sphericalPerlin2 = With[{octaves = 4},`

Compile[{{theta, _Real}, {phi, _Real}, {amplitude, _Real},

{frequency, _Real}, {gain, _Real}, {lacunarity, _Real}},

Module[{noiseVal = 0.0, x, y, z, sf, freq = frequency,

amp = amplitude},

sf = Sin[phi];

x = sf*Cos[theta] + 1; y = sf*Sin[theta] + 1; z = Cos[phi] + 1;

Do[

noiseVal += amp*Abs[signedNoise[x*freq, y*freq, z*freq]];

freq *= lacunarity; amp *= gain,

{octaves}

];

noiseVal

], CompilationOptions -> {"InlineExternalDefinitions" -> True}

]]

As an example, here is how to use `sphericalPerlin2[]`

along with `DensityPlot[]`

to generate a texture:

`DensityPlot[sphericalPerlin2[th, ph, 0.5, 2.0, 0.5, 2.0],`

`{th, 0, 2 Pi}, {ph, 0, Pi},`

AspectRatio -> Automatic, Frame -> None, PerformanceGoal -> "Quality", PlotPoints -> 200, PlotRange -> All]

and here is the result, after adding appropriate coloring (with thanks to Rob Collyer for the color gradient):

In fact, this is the texture I used for upper leftmost “marble” in the first image.