Mathematics Archives - James Ramsden

A kernel density estimation implementation for Grasshopper

Kernel density estimation is a way of generating a smoothed function of density of a collection of data.

Randal Olson used Kernel density estimation to perform a statistical likelihood on the location of Waldo in the Where’s Waldo (Where’s Wally) books. Perhaps surprisingly, this didn’t reveal a random distribution, but showed that Waldo was more likely to appear in certain areas of the page.

The mathematics of kernel density estimation are closely linked to kernel smoothing, which I use a lot for my interpolation development. And since pretty much everything I do is in Grasshopper, I wanted to see if I could quickly implement KDE and make a few pretty pictures.

Output

Here, we can input and move a collection of points. In regions with high point density, the colour becomes more red.

And of course, making the most of the parametric capability of Grasshopper, we can move the points and get an updated Kernel Density Estimation in real time.

Grasshopper script

We input a mesh (or a Brep to convert to a mesh), which will provide our analysis surface for colouring later. We also input a list of data points.

A bit of C# then calculates the ‘score’ for each mesh point – or an estimation for the density over that mesh point. This returns a list of numbers. We can then colour the mesh according to these values.

The C# script – the ‘density’ component – contains the following code:

  private void RunScript(List<Point3d> pts, Mesh mesh, ref object A)
  {

    var rtnlist = new List<double>();

    foreach(var v in mesh.Vertices)
    {
      double tempval = 0;
      foreach (var pt in pts)
      {
        tempval += kernel(pt.DistanceTo(v), 0.2); //0.2 is a smoothing value - you can change this 
      }
      rtnlist.Add(tempval);
    }

    A = rtnlist;
  }

  // <Custom additional code> 

  public double kernel(double input, double factor)
  {
    return Math.Exp(-(factor * input));
  }

Gravity wells

The kernel density estimation can sometimes be thought of as a smoothed histogram. Another way to visualise regions of high density is to take advantage of the Z component and give our mesh some height.

To do this, I reconstructed each mesh point using the vertex score as an input for the mesh point Z component. So that the input mesh points are still easily visible, I multiplied these scores by -1. This gave a visualisation quite akin to a collection of gravity wells, the points perhaps being tiny planets floating in a 2D space.

The mesh can be modified by deconstructing the mesh and then deconstructing the points within this mesh. We then supply a modified Z value, then re-construct the mesh.

Calculate the cross product: C# code

A method written in C# for calculating the cross product between two vectors:

  Vector3d CrossProduct(Vector3d v1, Vector3d v2)
  {
    double x, y, z;
    x = v1.Y * v2.Z - v2.Y * v1.Z;
    y = (v1.X * v2.Z - v2.X * v1.Z) * -1;
    z = v1.X * v2.Y - v2.X * v1.Y;

    var rtnvector = new Vector3d(x, y, z);
    rtnvector.Unitize(); //optional
    return rtnvector;
  }

The cross product enables you to find the vector that is ‘perpendicular’ to two other vectors in 3D space. The magnitude of the resultant vector is a function of the ‘perpendicularness’ of the input vectors.

Angle between two vectors

Something rather simple but I keep forgetting how to do it: what is the angle between any two vectors? It’s simple when you know how (as well as being an opportunity to try LaTeX for WordPress!):

$cos\theta=\frac{\sum{uv}}{\left\{(\sum u^2)(\sum v^2)\right\}^{0.5}}=\frac{\textbf{u'v}}{\left\{(\textbf{u'u})(\textbf{v'v})\right\}^{0.5}}$

Example

Find the angle between two vectors in 3D space:

$\textbf{u}=\begin{bmatrix}5&0&0\end{bmatrix}$
$\textbf{v}=\begin{bmatrix}10&2&5\end{bmatrix}$

$cos\theta=\frac{\sum{uv}}{\left\{(\sum u^2)(\sum v^2)\right\}^{0.5}}=\frac{5\times10+0\times2+0\times5}{\left\{(5^2+0^2+0^2)(10^2+2^2+5^2)\right\}^{0.5}}=\frac{50}{\left\{25\times129\right\}^{0.5}}=0.880$ $\therefore\theta=28.3\textdegree$

This technique can be used for any number of dimensions. For 2D space (e.g. vectors on a graph on a piece of paper) u and v will each contain two values instead of three, and the calculation is then done in the same way.

C# code example

  //returns angle between two vectors
  //input two vectors u and v
  //for 'returndegrees' enter true for an answer in degrees, false for radians
  double AngleBetween(Vector3d u, Vector3d v, bool returndegrees)
  {
    double toppart = 0;
    for (int d = 0; d < 3; d++) toppart += u[d] * v[d];

    double u2 = 0; //u squared
    double v2 = 0; //v squared
    for (int d = 0; d < 3; d++)
    {
      u2 += u[d] * u[d];
      v2 += v[d] * v[d];
    }

    double bottompart = 0;
    bottompart = Math.Sqrt(u2 * v2);


    double rtnval = Math.Acos(toppart / bottompart);
    if(returndegrees) rtnval *= 360.0 / (2 * Math.PI);
    return rtnval;
  }