Smallest Radius that includes x number of points


I’m trying to find the least expensive way to determine the smallest radius of a circle that circumscribes a X number of points. The goal is to find for each vertex in a mesh the smallest radius that circumscribes a given number of points from a list of points. The list of points, the mesh and the number of points to be circumscribed are inputs.

One method is to measure the distance from the mesh vertex to every point in the list and sort of distances. The Nth-1 number on the list will be the desired radius. This is a very slow method. Any other ideas?

A working example is:

private void RunScript(List<Point3d> pts, Point3d refPt, int numPt, ref object A, ref object B)
List<double> dist = new List<double>();

foreach(Point3d pt in pts)
  double d = Math.Sqrt((pt.X - refPt.X) * (pt.X - refPt.X) + (pt.Y - refPt.Y) * (pt.Y - refPt.Y));

A = (dist[numPt] + dist[numPt - 1]) / 2;

Another option would be to increment the radius of a circle until the desired number of point is inside. But I’m having problems implementing this and I’m not sure it will be a faster method.

The fastest way to perform a lot of cloud-closest-point searches is to use a spatial tree. Either an OcTree or a k-d tree. The proximity components in Grasshopper use QuadTrees and OcTrees (for the 2D, vs. 3D cases respectively). What you can do is find the N+1 closest points to your start point (+1 because you’ll need to include the point itself), then figure out which is the furthest. However that assumes the circle is centred at your start-point which need not be the case. If you want to really find the smallest possible sphere which contains a specific coordinate and at least N other points, then you will also need to work on minimum spanning spheres and possibly a lot of other algorithms that have no known efficient solutions.

It does not sound easy, or possibly even do-able with 100% accurate results.

You can access the GH octree code via the Grasshopper.Kernel.Geometry.SpatialTrees namespace.

Thanks @DavidRutten for the quick reply. It seems pretty complicated for my current knowledge of C# and the time I have available. The purpose of all this is implement an adaptative bandwidth Kernel Density Estimation for events in a city.
So I’m currently focusing on 2D and the circle is always centered on the start-point. I’m using the mesh vertices as start-points and the points being searched are events on a different list. The script I posted earlier takes about a second to execute for a single point. Which taking into account that the mesh may have 70000vertices and that there are 60000 points in the other list, I might be looking at many hours just to calculate a radius for each point. So I really need to find an alternative. Quadtrees are also in the same namespace? Or is there a different namespace for this type of data structure?

Actually I misremembered, this code has evolved since then. Just in Grasshopper.Kernel.Geometry you’ll find all sorts of classes for 2D and 3D spatial trees. You’ll want to start with Note2Tree or Node3Tree. Stay away from the SpatialTree namespace for now, I’m not 100% sure how well that code is tested.

Just a little detail as long as you use the algoritm above:
Since you don’t use the distances, skip the ‘expensive’ Math.Sqrt(...) method in the loop and collect only the square distances (the sorting will still do it’s thing) . Retrieve the Sqrt only from the final result before assigning A.

A for loop is a little faster than foreach.
And count down the loop instead of counting up.

    for (var i = pts.Count-1; i >= 0; i--)
        var pt = pts[i];
        dist.Add((pt.X - refPt.X) * (pt.X - refPt.X) + (pt.Y - refPt.Y) * (pt.Y - refPt.Y));
        // or 
        // dist.Add((pts[i].X - refPt.X) * (pts[i].X - refPt.X) + (pts[i].Y - refPt.Y) * (pts[i].Y - refPt.Y));

There should be a noticable difference in your case.

// Rolf

@RIL Thank you for you suggestions. I had already tried to remove the Math.Sqrt from the loop, but I didn’t see much difference. Now, I did it in another computer and I saw significant improvements. (I guess that I need to replug the data to the C# component before I can evaluate the performance of the script) Changing the foreach to a for loop isn’t making a bid difference. Overall I’m getting much better times.

@RIL suggested some micro-optimisations. They can certainly make a difference, but true improvement is achieved by lowering the big-O runtime of your algorithm. If you’re comparing every point to every other point, that means that you have to loop twice over your entire cloud. If n represents the size of your pointcloud, your runtime now is O(n^2). That in turn tells you that if you double the amount of points, you quadruple the amount of time needed to complete the algorithm. If you switch to a more economic closest-point search (for example spatial trees, or bucketing, or even just linear sorting along an axis), then the inner loop will no longer have to consider every single point in the cloud. It can be log(n) instead of plain old n. So your algorithm as a whole is now O(n \cdot log(n)). If you double the size of your cloud now, it’ll only take a little over twice as long.

Long time ago I wrote a series of blog posts about nearest point searches. It’s not state of the art stuff but hopefully understandable to students. First post.

The trees are def the fastest but I’ve found this method Point3dList.ClosestIndex to be really fast for this stuff and you won’t have to sort.(also requires very little lines of code) .

@DavidRutten Thanks for the link. Much clearer than Wikipedia or Google Searches on quadtrees.

I only find Grasshopper.Kernel.Geometry.SpatialTrees. Maybe I’m looking in the wrong place

Thanks I’ll have a look.


Here is c# multithreading code which can be used to write plugin.
With little modification it can be used in C# component for grasshopper.
I tested it on laptop with i7-7700HQ, Win 10.
It took 340 seconds to find minimum sphere containing 50 points for each of 70.000 vertices and 60.000 points.

    private List<int>[] DetermineMinimumSpheres(List<Point3d> vertices, List<Point3d> cloudPoints, int numOfPoints)
        Point3dList verticesPtL = new Point3dList(vertices);
        var cloudPointsCC = new Point3dList(cloudPoints);
        List<int>[] resultSpheres = new List<int>[vertices.Count];
        var tpBB = verticesPtL.BoundingBox;
        var cpBB = cloudPointsCC.BoundingBox;
        double diagonalSum = tpBB.Diagonal.Length + cpBB.Diagonal.Length + 1;
        Vector3d vectorAdd = new Vector3d(diagonalSum, diagonalSum, diagonalSum);
        Vector3d bbVector = cpBB.Center - tpBB.Center;
        if (!bbVector.IsTiny())
            vectorAdd = bbVector;
            if (!bbVector.Unitize())
                return null;
            vectorAdd = vectorAdd + bbVector * diagonalSum;
        ParallelLoopResult result = Parallel.For(0, vertices.Count, i =>
            resultSpheres[i] = PointIdsInsideSphere(vertices[i], cloudPointsCC, numOfPoints, vectorAdd);
        return resultSpheres;

    private static List<int> PointIdsInsideSphere(Point3d verticePt, Point3dList cloudPoints00, int numOfPoints, Vector3d vectorAdd)
        List<int> resultIDs = new List<int>();
        Point3dList cloudPoints = new Point3dList(cloudPoints00);
        if (cloudPoints.Count < numOfPoints)
            return resultIDs;
        for (int i = 0; i < numOfPoints; i++)
            var closestId = cloudPoints.ClosestIndex(verticePt);
            cloudPoints[closestId] = cloudPoints[closestId] + vectorAdd;
        return resultIDs;