Point related operations execution speed

I’ve been doing some research into mimicking the logic found here and [here] (http://research.microsoft.com/pubs/115797/paper.pdf) to create stippled distribution via grasshopper+python. I’ve been able to script the logic, but am finding that the execution time to be way too long. I’m not sure what kind of times the authors of the research were getting, but I am hoping to make mine faster.

A lot of the operations in the script are related to solving for the distance between two points. Currently I am using rhinocommon’s Rhino.Geometry.Point3d.DistanceTo( ). Would it be faster to solve it mathematically? Previously, I’ve been told that you can add and subtract points like numbers. I’m just wondering if going through rhinocommon is adding anything to the time. I have found that ghpythonlib.component operations seem to take longer to execute relative to rhinocommon, which is why I avoid using them and why I am posing this question to bypass rhinocommon.

Out of curiosity, has anyone else created a good stippling script?

Example from research authors:


I’d test this to make your own method and compare speed.

Options to consider: to get the distance you need to square root. However in some cases the squared distance could suffice, at least to filter out distant points.

Another option is to filter out distant points by first comparing individual values of X ,Y and Z. Without the need to calculate exact distance if the comparison already renders them out of range.

Also are you creating 2D or 3D dots? If so, you want to use the Rhino.Geometry.Point2d as a 3d point will needlessly calculate with the z values as well.


The calculation behind this function is pretty simple…

Not to mention very robust when dealing with very small distance. If you feel you can get away with it you can switch to a custom distance metric which only calculates the squares, that will be somewhat faster.

However the true speed-up here is to figure out which distances not to compute. It sounds like what you need is a faster nearest-neighbour lookup. I have a series of 6 blog posts about nearest neighbour search. You can find the first one here.

Thanks for pointing me to your blog posts. I read them and think it’s interesting how I’ve pondered similar shortcuts.

After my initial post, I’ve discovered that nearest neighbor lookup is what I was coding manually. The equivalent would be Rhino.Geometry.PointCloud.ClosestPoint(), correct? Does this Rhino command use the faster search methods that you blogged about?

I am coding this script in ghcomp as part of a grasshopper script. I saw somewhere in one of the papers that their script was executing in 30-40 seconds to stipple something simple. My script, when iterating over a similar point count, is taking about 6 hours (I have not implementing multithreading yet). Pretend that I am not using any rhino related commands and just doing mathematical calculations- do python scripts execute slower through grasshopper than it would as a standalone program?

This class is really nothing more than an array of points. Thus, it’s ClosestPoint function is not super efficient…

Is there any nearest-neighbor lookup function in rhinocommon that I can use then? The PointCloud method is the only one I’ve ever used (and it’s working faster than my previous script version since implementing this afternoon). Is there a function that uses the methods in @DavidRutten 's blog posts? Or do I have to script it myself?

You might take a look at Rhino.Geometry.RTree.


I´ve tried to write some kind of nearest neighbour search function some time ago to generate normals for pointclouds with


method. Searching the cloud for neighbours in a cloud with 50K Points and searching for 6 neighbours it runs with approx. 2000pts/sec. How large are your pointclouds ?


Just as an update, doing dist=(p1.X-p2.X)^2+(p1.Y-p2.Y)^2 seems to take longer than dist=Rhino.Geometry.Point3d.DistanceTo(). I expected it to be faster too.

you might speed this up further by importing the method at the head of your script eg:

import Rhino
from Rhino.Geometry.Point3d import DistanceTo as Dist2Pt
# and later use
dist = p1.Dist2Pt(p2)


To get a good stippled image I need to compare at least 10,000+ stipple points (called sites) with 100,000+ underlying points.

I am unfamiliar with RTree. I see that you can set up a RTree by inputting a PointCloud, but how do you get the closest point/points in that tree to a predetermined point? And why is this faster than the regular PointCloud.ClosestPoint() method?

Side question: does your Dist2Pt suggestion make the script run faster or just make typing the code faster?

I guess this might be much faster than my case because i am constantly adding and removing points to/from the RTree which might not be required if only one single neighbour per point is sufficient. Here is a good example using RTree, note that it searches mesh vertices but the overal routine is quite similar when searching closest points in a pointcloud:


It makes the code run slightly faster because the lookup time for the method falls away. You might timeit to see the difference. This only makes sense if you call a method thousands of times.



for your neighbour search you might find these two links interesting. The first is using numpy which i never got to work under x64. The second goes a bit deeper but does not require numpy.


Use an RTree for this. It is exactly what they were designed for.


can you provide a python example which shows how to find eg. 5 nearest neighbours of a pointcloud point without adding/removing points from the RTree during search ?


Looking a little closer, I can see how with what I’ve exposed that this is not very easy or obvious on how to solve closest point testing with an RTree. I’ll see what I can come up with for an example.


@stevebaer, A concise example would be appreciated.

@clement , Overnight I ran the script with PointCloud.ClosestPoint implemented and at half the point count as my weekend test (29 hours) and it only took 2.4 hours, which is a really big improvement. But after reviewing the results closely, I identified a few thing: the site positions are barely distributed well after a few iterations (more are necessary), the iterations are taking way to long, and my sample points look better than the sites do as a stipple pattern.

This is the logic I’ve been following (from http://graphics.uni-konstanz.de/publikationen/2009/capacityconstrainedpointdistributions/website/):

I discovered this on Friday to use as my weighted sample point generator and it seems to give pretty good results that could pass as a stipple pattern (section 6.0 is where the logic is detailed): http://www.eng.uwaterloo.ca/~alopyrev/cs791/stippling_report.html

Here’s an example reusing an RTree to find the 5 closest points to any given point in space. Try creating something like a MeshTorus with 300x300 rows/columns and run this script against that object. A lot of the time the closest 5 points can be found in a very small number of samplings.

Hopefully we can come up with easier to use functions in V6 that would significantly reduce the size of this script (and potentially give a slight performance bump.)

import rhinoscriptsyntax as rs
import Rhino.Geometry as rg
import random

def search_callback(obj, args):
    cloud, points, search_point, call_count = args.Tag
    args.Tag[3] = call_count + 1
    location = cloud[args.Id].Location
    index_to_replace = -1
    distance = location.DistanceTo(search_point)
    # find furthest point from test point
    for i, point in enumerate(points):
        d = point.DistanceTo(search_point)
        if d > distance:
            distance = d
            index_to_replace = i
    if index_to_replace >= 0:
        points[index_to_replace] = location
        args.SearchSphere = rg.Sphere(search_point, distance)
def main():
    mesh_id = rs.GetObject("Select Mesh", rs.filter.mesh)
    if mesh_id is None: return
    cloud = rg.PointCloud(rs.MeshVertices(mesh_id))
    tree = rg.RTree.CreatePointCloudTree(cloud);
    search_count = 5
    while True:
        search_point = rs.GetPoint("Search Point")
        if search_point is None: break
        points = []
        radius = -1
        #seed point list with 'search_count' items for starters
        sample = random.sample(xrange(cloud.Count), search_count)
        for i, sample_index in enumerate(sample):
            r = points[i].DistanceTo(search_point)
            if r > radius: radius = r
        tag = [cloud, points, search_point, 0]
        tree.Search(rg.Sphere(search_point, radius), search_callback, tag)
        print "search called {0} times for mesh with {1} points".format(tag[3], cloud.Count)


You should see large performance improvements using my RTree sample over the PointCloud.ClosestPoint function when there are a significant number of points involved.