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:
privatevoid RunScript(List<Point3d> pts, Mesh mesh, refobject 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> publicdouble 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.
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
privatevoid RunScript(refobject 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.
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;
}
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.
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:
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:
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.
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.)
privatevoid RunScript(List<Point3d> pts, double ht, refobject A)
{
var building = new Mesh();
//make wallsfor (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 roofvar 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:
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:
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! 🙂
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.
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:
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.
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.
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.
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.
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.
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…
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.
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.
privatevoid RunScript(object x, object y, refobject 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.
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.
This method essentially replicates the Grasshopper components below:
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:
publicstatic 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:
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
privatevoid RunScript(bool x, List<object> y, refobject A)
{
if(x)
{
//instantiate new value listvar vallist = new Grasshopper.Kernel.Special.GH_ValueList();
vallist.CreateAttributes();
//customise value list positionint 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 componentthis.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.
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>privatereadonly RhinoDoc RhinoDocument;
/// <summary>Gets the Grasshopper document that owns this script.</summary>privatereadonly GH_Document GrasshopperDocument;
/// <summary>Gets the Grasshopper script component that owns this script.</summary>privatereadonly 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>privatereadonlyint 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:
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!