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

Kernel density estimation Grasshopper Rhino 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.

Kernel density estimation Grasshopper

Create a list of random numbers in C#

Generating random numbers in C# isn’t the most intuitive process. I originally learnt to code back in the day using VB, where generating a random number was as simple as calling the Rnd() function.

In C#, there is slightly more work to do. We first have to create an instance of the Random class, which is non-static. Only then can we call a function to generate a random number for us.

The example below is written for the C# component in Grasshopper, but the gist of it can be used for any C# application.

C# code

  private void RunScript(ref object A)
  {

    var rand = new Random();
    var rtnlist = new List<double>();

    for (int i = 0; i < 100000; i++)
    {
      rtnlist.Add(rand.Next(1000));
    }
    A = rtnlist;

  }

What’s happening here?

We create an instantiation of the Random class. We call this ‘rand’. This is what contains the logic to generate our random numbers.

Since computers are inherently quite terrible at generating truly random numbers, programming languages tend to use pre-programmed lists of random numbers instead. The random class picks a place in the list to start from. Then, every time we call for a new random number, it is simply reading the next value in the list.

This is why the function to generate a new random number is called ‘Next’ – it is simply looking at the next value in its list of random numbers. The 1000 in the brackets specifies that the number returned will be scaled between 0 and 1000. (If we left the brackets blank, the random number would be between 0 and the highest possible integer.)

The script creates a list of 100000 random numbers between 0 and 1000. The output is sent to the variable ‘A’. This is a Grasshopper specific thing – if you were writing this as a method, you could change this to ‘return rtnlist’.

If you need many random numbers, you should still only create one instance of Random, and then use ‘Next()’ to build your collection of random numbers. If you don’t do this, and instead create an instance of ‘Random’ for every random number you need, the resulting numbers may not be random at all.

Triangulate a quad mesh in C#

How to take a mesh of quad faces and return a mesh with triangle faces.

The principle of the algorithm is to look at each quad and split them into two triangles. It chooses the split by the shortest diagonal across the quad.

The code is written in C# for Rhino/Grasshopper, using the Rhino mesh structure. However, it can be easily adapted to any mesh application.

  public Mesh Triangulate(Mesh x)
  {
    int facecount = x.Faces.Count;
    for (int i = 0; i < facecount; i++)
    {
      var mf = x.Faces[i];
      if(mf.IsQuad)
      {
        double dist1 = x.Vertices[mf.A].DistanceTo(x.Vertices[mf.C]);
        double dist2 = x.Vertices[mf.B].DistanceTo(x.Vertices[mf.D]);
        if (dist1 > dist2)
        {
          x.Faces.AddFace(mf.A, mf.B, mf.D);
          x.Faces.AddFace(mf.B, mf.C, mf.D);
        }
        else
        {
          x.Faces.AddFace(mf.A, mf.B, mf.C);
          x.Faces.AddFace(mf.A, mf.C, mf.D);
        }
      }
    }

    var newfaces = new List<MeshFace>();
    foreach (var mf in x.Faces)
    {
      if(mf.IsTriangle) newfaces.Add(mf);
    }

    x.Faces.Clear();
    x.Faces.AddFaces(newfaces);
    return x;
  }

Before and after:

triangulate mesh

Make buildings with heights in Elk for Grasshopper

I have recently been constructing whole cities in Rhino using the Elk plugin. This parses OpenStreetMap data into Grasshopper geometry, allowing you to visualise any city or area in the world in your Rhino viewport.

York in Grasshopper with Elk OpenStreetMap data

One of the most interesting ways of getting your cities to look like cities is getting the buildings right. It is possible to extract building data with Elk, but this isn’t intuitive if you aren’t already familiar with OpenStreetMap data. And even if you are, the output is a collection of points, which requires further work to turn it into actual buildings.

This article will show you how I created buildings by extracting the OSM data and converting them into 3D meshes.

OpenStreetMap: Understanding the data

A crash course into OpenStreetMap data: All (well, nearly all) objects in the map are either saved as points or as lists of points. Points are useful for small objects, such as benches, letterboxes or trees. For larger objects, such as roads and buildings, collections of points are used. For roads, these points join together to form lines along the centre of the road. For buildings, these points define the line around the building’s edge. These ‘point collections’ are essentially polylines.

Every object, whether a point or a polyline, has a list of data associated with it. Each item in this list has a key and a value. It is with these keys and values that we describe the object. The most common key is ‘name’, and we describe the name in the value field. Most objects have many tags. For example, if we go into the OpenStreetMap editor and have a look at a local branch of Costa Coffee:

OpenStreetMap tags values keys example

We have a list of keys and a list of values. The values are the data itself; the keys describe the type of data.

If you want to see this for yourself, go into the editor on OpenStreetMap. I personally prefer the Potlatch2 editor, which you can find in the drop-down next to the ‘edit’ button at the top. To see the tags, select an object on the map, then click ‘advanced’ in the bottom left.

Buildings in OpenStreetMap

Similarly, all buildings contain a collection of tags. The great thing about tags is that they are an open concept – you can use any tags that you like. With buildings, the most common tags include the building name and the building type. For example:

OpenStreetMap tags values keys example building

Using Elk to extract buildings

Elk can be used to extract the points that form the building outline. Like polylines in Rhino, the data is saved as a collection of points, and not a line itself, so it will be up to us to stitch the points together to create a building outline.

There is no built-in component that returns directly the building points, but we can use the GenericOSM to filter for data containing the ‘building’ key. The Polyline component joins the points together.

Grasshopper components Elk extract buildings

building outlines openstreetmap elk grasshopper

Make your buildings 3D

This is all well and good, but don’t they look a bit flat and boring? Let’s make them look a bit more realistic.

The easiest way to turn your polylines into 3D buildings is to extrude them, and then use patch to create a roof. But I can tell you through my own experience that this is a bad idea. For a few buildings it’s okay, but for the city scale, working with surfaces is very slow and you’ll quickly see your computer grind to a halt. To build our buildings at lightning speed, we have to use meshes.

There isn’t anything natively in Grasshopper that does this for us, but we can quickly write something in C#. Basically, we build up the walls with mesh faces, then use a Rhino command to cap the top for the roof.

This script takes in a list of points that describe the building outline, and a value for its height. (Let’s just send it a value of 10m or so for now.)

  private void RunScript(List<Point3d> pts, double ht, ref object A)
  {

    var building = new Mesh();

    //make walls
    for (int p = 0; p < pts.Count; p++)
    {
      building.Vertices.Add(pts[p]);
      building.Vertices.Add(new Point3d(pts[p].X, pts[p].Y, pts[p].Z + ht));
    }

    for (int p = 0; p < pts.Count - 1; p++)
    {
      int j = p * 2;
      building.Faces.AddFace(j, j + 1, j + 3, j + 2);
    }

    //make roof
    var roofpts = new List<Point3d>();
    for (int p = 0; p < pts.Count; p++)
    {
      roofpts.Add(building.Vertices[p * 2 + 1]);
    }
    var pline = new Polyline(roofpts);
    var roof = Mesh.CreateFromClosedPolyline(pline);
    building.Append(roof);

    A = building;

  }

This then gives us:

Elk Grasshopper example buildings

Elk Grasshopper example buildings

But not all buildings are 10 metres tall…

It made things a lot easier to assume a height for our buildings. But can we do any better?

Given that we have a height input, we should be able to put something a bit more sensible into it. Some buildings in OSM do come with height data, though most don’t. If they do, this information is under the height key. How can we extract the building’s height information and apply it to our buildings with Elk?

It currently isn’t easy as there isn’t a comprehensive tag filter within Elk. One solution is to use the GenericOSM with k=height. This will return everything with height data; we just have to assume that everything returned is a building (it usually is) and model it as such.

Set up your Grasshopper as below, using the same C# script as above to create the buildings:

Elk Grasshopper example buildings

For my example of York, barely a single building has any height data. But, very happily, one very patient mapper has recorded the height of all the different elements of York Minster, leading to a very pleasing model.

But I really want heights on ALL my buildings!

Then you’ll likely have to pay up.

At least in the UK, the only near-comprehensive source of building heights is via the Ordnance Survey Topography Layer. Limited amounts of data are available as free academic licences, but otherwise it ain’t cheap.

In the meantime, if you’ve made it this far, you’ve probably realised that there are a lot of gaps in the OSM database. These gaps are best filled by people with local knowledge where you live, so I would encourage you to register and to start mapping! 🙂

Modelling the city of York with Grasshopper, Elk and OpenStreetMap

Since my post last year on modelling cities in Grasshopper, I have learnt a few more things in how to get the most out of Grasshopper. I am now much more adept in using meshes, a much lighter data type than surfaces, and I have learnt a few more tricks in getting the buildings to look a little more building-like.

The combination of these means that I can now load in much larger areas into Grasshopper with much better performance. I decided to revisit the challenge of opening cities in Grasshopper and see where I could take it.

The city of York

In this post, I have modelled my historic home city of York. It is famous for its Minster and for its City Walls, and also has two rivers, numerous parks, and other historical buildings, which should make for some interesting renders.

Here is an actual photo of York:

Bird's eye view of York city centre

Image Wikipedia

And here is a rather awesome video of York:

Modelling in Grasshopper

As in my previous attempts, I am using OpenStreetMap map data. I download a map on the OSM export page, which provides an OSM file. I then use the Elk plugin for Grasshopper.

I have set up a standard Grasshopper file where I can connect it to any OSM file, and with no user intervention, it will generate a 3D model of the map data in this OSM file. The Grasshopper file is available here.

Results

Here is the above photo modelled in Grasshopper:

York in Grasshopper with Elk OpenStreetMap data

The buildings are mostly modelled using their outlines, which are available for most city centre locations and key landmarks. A very few buildings also have building height data, which allows the Minster to be modelled well. All buildings without height data are instead modelled by a random height of between 5 and 15 metres. (This gives a slightly more realistic appearance than, say, having all buildings to be 10m tall.)

York Minster

York Minster is a highly impressive structure, towering over the rest of the city.

_73551782_minster_wide

minster

Somebody has spent a lot of time and effort not only entering the height and outline of the Minster into the OSM database, but recording the heights and outlines of all the individual elements of the building. This means that we get a surprisingly intricate model of the Minster out of the box with no user intervention.

York in Grasshopper with Elk OpenStreetMap data - Minster

Foss Islands Road

Foss islands is an old industrial area near the city centre which is now mostly devoted to retail. A key landmark remains – an ominous Victorian chimney which stands slightly out of place next to a Morrisons supermarket.

Foss Islands York Morrisons and chimney

Image source

This also has height data, which produces some interesting views in its location next to the River Foss:

York in Grasshopper with Elk OpenStreetMap data - Foss Islands Road

Looking the other way, we can see the Minster in the distance:

York in Grasshopper with Elk OpenStreetMap data

York city walls

The city walls, around 800 years old, almost entirely circumnavigate the city centre.

York city walls, overlooking Minster

Image by Steve Nova: source

Unfortunately there is no automatic way to model the grass banks around the walls, but the walls are saved as polylines in OSM. These walls are made by converting the polylines into 2D meshes (the same as how the roads are generated) and then extruded 10m up. They are quite jagged though and could do with some more work.

York in Grasshopper with Elk OpenStreetMap data - walls and Minster

C# code and links

Some key code snippets used:

Roads

OSM roads are effectively polylines. In order to render them with colours, they need ‘expanding’ into 2D shapes that have a width. This is done by creating a mesh face for every polyline segment. The component for this is available here.

Buildings

Here is a more in-depth discussion on how to create buildings in 3D using OpenStreetMap data.

Grasshopper file

The whole Grasshopper file is available here.

Orient a Grasshopper mesh so that its normals face upwards

How to read the normal of a mesh in Grasshopper in C#, and flip the mesh if the normals are facing the wrong way.

Creating city buildings in Rhino

I recently revisited the task of generating cities in Rhino using Elk. One of the problems with this original solution was that I had no fast way of generating the roofs of the buildings. The most obvious solution to generate the building roofs was to use the ‘patch’ function of Rhino using the building outline as input. But surfaces are data- and processor-heavy, and files at the city scale would crash the computer. So I left the roofs out…

Rendering Park Street in Bristol in Grasshopper and Rhino using the Elk plugin and OpenStreetMap data.

Meshes are much lighter data structures than surfaces, and should be used where possible where performance is an issue. In this post, meshes were created by casting curves to meshes.

Mesh orientation

This works well, but the output meshes do have rendering issues. All mesh faces have a direction (a front and a back) and Rhino renders the mesh differently depending on what side it thinks it’s looking at. It is apparently quite random whether the above meshing trick generates a mesh which faces upwards or downwards.

Grasshopper mesh vector orientation problem on buildings

Flip backwards meshes

The solution is to ‘flip’ the offending meshes so that all meshes are facing the same way. The challenge lies in detecting the meshes which are facing the wrong way, and then to find a way to flip these meshes. The tidiest way is to do this with a C# component.

  private void RunScript(object x, object y, ref object A)
  {

    Mesh mesh = (Mesh) x;
    for(int i = 0; i < mesh.Faces.Count; i++)
    {
      var normal = mesh.Normals[0];
      if(normal.Z < 0)
      {
        mesh.Flip(true, true, true);
      }
    }

    A = mesh;

  }

The result is that the roofs now are all the same colour. The red lines, representing the normals of the mesh faces, are all now facing the same direction.

Grasshopper vectors on mesh direction problem

Convert a Brep to a Mesh in Rhino/Grasshopper C#

How to convert a Rhino Brep to a Rhino Mesh using RhinoCommon. This method is suitable for use in Grasshopper development and the Grasshopper C# component.

Create meshes from Brep

The code below produces a list of customised meshes based upon an input Brep and some mesh settings. The examples are written as static extension classes for compiled components (so you can access it directly as a Brep method), though can be easily adapted for the C# component too.

        public static List<Mesh> BrepToMeshes(this Brep brep, double maxEdge)
        {
            Mesh[] mesh;
            MeshingParameters mp = new MeshingParameters();
            mp.MaximumEdgeLength = maxEdge;
            mesh = Mesh.CreateFromBrep(brep, mp);
            return mesh.ToList<Mesh>();
        }

This method essentially replicates the Grasshopper components below:

Grasshopper component mesh to Brep

The MeshingParameters class replicates the ‘Mesh Settings’ component. You can assign settings to your mesh by creating an instance of MeshingParameters and accessing its properties, much as I have done with MaximumEdgeLength.

Joining the meshes

The output is a list of meshes. If you want to truly replicate the Mesh Brep component, you will also need to join all meshes in the list into a single mesh. This can be done with the Append method:

        public static Mesh JoinMeshes(this List<Mesh> meshes)
        {
            var rtnmesh = new Mesh();
            foreach (Mesh mesh in meshes)
            {
                rtnmesh.Append(mesh);
            }
            return rtnmesh;
        }

And for completeness, you can call the methods in a single line:

        public static Mesh BrepToMesh(this Brep brep, double maxEdge)
        {
            return JoinMeshes(BrepToMeshes(brep, maxEdge));
        }

Get points from polyline in RhinoCommon/Grasshopper C# component

Polylines are highly lightweight forms of curves that have many useful applications in Rhino and Grasshopper.

What is a polyline?

A polyline is a list of points, nothing more. When Rhino reads in a polyline, it creates the curve by doing a dot-to-dot of the points.

It can also be thought of as an interpolation of points with degree = 1. The curve will always pass through the inputted points.

polyline

Extracting the points from a polyline

Let’s say you have a polyline as input in Grasshopper. You want to extract the list of points used to create that polyline within a C# script.

The most intuitive method is to look for a property of the polyline called ‘Points’ or something like that. But if we try and find it…

polyline_grasshopper_code

…there’s nothing there! So if there’s no property containing the points, where is this data?

Method 1: Get the whole list of points

As shown in the RhinoCommon SDK documentation, PolyLine inherits from the Point3dList, which itself inherits from the RhinoList class.

We can use the ToList() method to get the list of points:

  private void RunScript(Polyline pline, object y, ref object A)
  {
    var pts = new List<Point3d>();

    pts = pline.ToList();

    A = pts;
  }

Method 2: Get a single point

We can extract a single point from the ToList() method:

  private void RunScript(Polyline pline, object y, ref object A)
  {
    var pt = new Point3d();

    pt = pline.ToList()[0];

    A = pt;
  }

Or, even shorter, we can treat the polyline like a list itself and access the points directly:

  private void RunScript(Polyline pline, object y, ref object A)
  {
    var pt = new Point3d();

    pt = pline[0];

    A = pt;
  }

Thanks to Fraser Greenroyd for our joint effort in unpicking this deceptively unintuitive challenge – he also wrote about it here.

Instantiate a Value List Grasshopper component with C#

In a previous post I instantiated a slider automatically using only C#. This means it is possible for a component to add components to your canvas automatically, potentially saving a lot of time in many repetitive workflows.

Here is a similar piece of code, but for the Value List:

C# code

  private void RunScript(bool x, List<object> y, ref object A)
  {
    if(x)
    {
      //instantiate  new value list
      var vallist = new Grasshopper.Kernel.Special.GH_ValueList();
      vallist.CreateAttributes();

      //customise value list position
      int inputcount = this.Component.Params.Input[1].SourceCount;
      vallist.Attributes.Pivot = new PointF((float) this.Component.Attributes.DocObject.Attributes.Bounds.Left - vallist.Attributes.Bounds.Width - 30, (float) this.Component.Params.Input[1].Attributes.Bounds.Y + inputcount * 30);

      //populate value list with our own data
      vallist.ListItems.Clear();
      var item1 = new Grasshopper.Kernel.Special.GH_ValueListItem("Red", "r");
      var item2 = new Grasshopper.Kernel.Special.GH_ValueListItem("Green", "g");
      var item3 = new Grasshopper.Kernel.Special.GH_ValueListItem("Yellow", "y");
      vallist.ListItems.Add(item1);
      vallist.ListItems.Add(item2);
      vallist.ListItems.Add(item3);

      //Until now, the slider is a hypothetical object.
      // This command makes it 'real' and adds it to the canvas.
      GrasshopperDocument.AddObject(vallist, false);

      //Connect the new slider to this component
      this.Component.Params.Input[1].AddSource(vallist);
    }
  }

This code is written for the C# component in Grasshopper. If you are writing components within Visual Studio, there is a little extra work to do in defining GrasshopperDocument and Component, otherwise it’s pretty much a copy-and-paste job too.

GrasshopperDocument and Component variables in Grasshopper from C# component to Visual Studio

In the C# component for Grasshopper, it is easy to perform functions on either itself or on the document you are in. For example:

Component.Params.Input[0].SourceCount;

…will give you the number of things connected to the first input of the particular C# component containing that line of code.

This is made possible by some variables that are pre-defined within all C# components. Start a new C# component and scroll to lines 49-61 and you should see:

#region Members
  /// <summary>Gets the current Rhino document.</summary>
  private readonly RhinoDoc RhinoDocument;
  /// <summary>Gets the Grasshopper document that owns this script.</summary>
  private readonly GH_Document GrasshopperDocument;
  /// <summary>Gets the Grasshopper script component that owns this script.</summary>
  private readonly IGH_Component Component;
  /// <summary>
  /// Gets the current iteration count. The first call to RunScript() is associated with Iteration==0.
  /// Any subsequent call within the same solution will increment the Iteration count.
  /// </summary>
  private readonly int Iteration;
#endregion

This is all well and good if you want to use the C# component, but what if you’re developing something a bit more meaty, and want to move over to Visual Studio?

I often use the GrasshopperDocument and Component properties in particular, and wanted to be able to use these in VS.

Can’t I just call GrasshopperDocument and Component in Visual Studio?

Within the C# component, some magic goes on that gives these variables their abilities. So copy and paste of the above variables isn’t going to work.

Instead, what we need to do is create our own variables with the same names as above, and manually give them their functionality.

How to set up GrasshopperDocument and Component

Like pretty much everything on this blog, I write about stuff as I learn it, and there may be a better way to do things! But the method below works for me.

Open or start a new Grasshopper component project as usual. (See here for setting up.)

At the very top of your class, create some class variables:

GH_Document GrasshopperDocument;
IGH_Component Component;

Then, within Solve_Instance, we can give meaning to these variables:

Component = this;
GrasshopperDocument = this.OnPingDocument();

Or alternatively,

GrasshopperDocument = Instances.ActiveCanvas.Document;

So, your whole component class should look like:

    public class GHtestComponent : GH_Component
    {
        GH_Document GrasshopperDocument;
        IGH_Component Component;

        public GHtestComponent()
            : base("GH test", "Nickname",
                "Description",
                "Category", "Subcategory")
        {
        }

        protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager)
        {
        }

        protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager)
        {
        }

        protected override void SolveInstance(IGH_DataAccess DA)
        {
            Component = this;
            GrasshopperDocument = this.OnPingDocument();
        }

        protected override System.Drawing.Bitmap Icon
        { //..
        }

        public override Guid ComponentGuid
        { //..
        }
    }

Notes

According to this GH forum post, it is possible for the OnPingDocument() to return null. I haven’t explored further on when this might happen, but you may wish to wrap this line in an ‘if’ or a ‘try catch’.

It seems unnecessary to me to have to assign values to the variables in the Solve_Instance, but this is a way I have working. I couldn’t get it to work when I put it in the constructor.

I still haven’t found the equivalent line for RhinoDoc, mostly because I haven’t needed to. When the time comes for me to figure it out (or if you happen to know yourself and would like to share it with the world!) I’ll update this page.

Update: RhinoDoc

Accessing the RhinoDoc variable is actually really easy – at the top of your C# file, just add a ‘using’:

using Rhino;

Then you can use the RhinoDoc variable like you would use it in the C# component. Simple!