PointCloud Analysis


#1

Greetings,
I am attempting to do some analysis of PointCloud objects within Rhino.
I have a pretty simple script in which I am attempting to build a best fit plane from a subset of points within a point cloud. The inputs are a cubic polysurface volume, which is used to limit the points used to create a best fit plane and a point cloud.

def PointTest(Brep,Cloud):
    if rs.IsPolysurfaceClosed(Brep)!=True: return
    brep=rs.coercebrep(Brep)
    Cloud=rs.coercegeometry(Cloud)
    Result=[]
    tol=scriptcontext.doc.ModelAbsoluteTolerance
    for pt in Cloud.GetEnumerator():
        point=rs.coerce3dpoint(pt.Location)
        if brep.IsPointInside(point,tol,True)==True:
            Result.append(point)
    Plane=rs.PlaneFitFromPoints(Result)   
    Plane=rs.coerceplane(Plane)
    rc,Curves,points=Rhino.Geometry.Intersect.Intersection.BrepPlane(brep,Plane,tol)
    if rc==True:
        PlanarSrf=rs.AddPlanarSrf(Curves)
    return PlanarSrf

Currently I am dealing with a point cloud that is around half a million points. Processing a cloud this size takes around 3.2 seconds. The lower limit for size of point cloud that I am dealing with is around 20 million points. With this in mind I attempted to create a multi-threaded solution for parsing through the enumerated point cloud, code below. I have not been able to get past the “The type arguments for method ‘ForEach’ cannot be inferred from the usage. Try specifying the type arguments explicitly.” error message. Should it possible to parse through a single point cloud using multipul threads? If so what is the proper syntax for doing so?

def Threaded_PointTest(Brep,Cloud):
    ## Make Sure Test Polysurface is Closed ##
    if rs.IsPolysurfaceClosed(Brep)!=True:
        print "Test brep is Not Closed"
        return
    
    ## Get a Private Copy for Multi-Threading ##
    brep=rs.coercebrep(Brep)
    brep.EnsurePrivateCopy()
    
    ## Turn Object into a PointCloud ##
    Cloud=rs.coercegeometry(Cloud)
    
    ## Holder for Points inside Surface ##
    Result=[]
    tol=scriptcontext.doc.ModelAbsoluteTolerance
    
    def TestPoints(Pt):
        try:
            print Pt
        except:
            pass
    
    CloudPts=Cloud.GetEnumerator()
    tasks.Parallel.ForEach(CloudPts,TestPoints)

Thanks,
5chmidt


[RhinoCommon] RTree is too limited
(Menno Deij - van Rijswijk) #2

You may want to look into ghpython parallel routines, as mentioned below because, to quote Steve Baer, “it is pretty hard to get the syntax correct in python for using the generic Parallel.For loop”. (which is exactly where you are getting problems)


#3

Hi Peter,

I’m not an expert but I think you’ll have to call this parallel task with a range of integers, for example:
tasks.Parallel.ForEach(range(Cloud.Count),TestPoints)
And you’ll have to create a list for each result, like:
Result = range(Cloud.Count)

You cannot append the results to a list as you will not get a meaningful printing inside a paralleled task.

If I understand correctly what you are trying to do, i.e. filtering points which are inside a box-like volume I would probably just check the coordinates numerically. You can transform the points to fit to the box orientation. This may be faster than brep.IsPointInside (if it does not have a special case for boxes).


#4

Thanks for the replies. I’ll test code up both options and do some testing with larger point cloud samples and see how it goes.


#5

I have assembled the following using script using ghpython.

def Threaded_PointTest(Brep,Cloud):
    ## Make Sure Test Polysurface is Closed ##
    if rs.IsPolysurfaceClosed(Brep)!=True:
        print "Test brep is Not Closed"
        return
    
    ## Get a Private Copy for Multi-Threading ##
    brep=rs.coercebrep(Brep)
    brep.EnsurePrivateCopy()
    
    ## Turn Object into a PointCloud ##
    Cloud=rs.coercegeometry(Cloud)
    
    ## Holder for Points inside Surface ##
    Result=[]
    tol=scriptcontext.doc.ModelAbsoluteTolerance
    
    def TestPoints(Pt):
        try:
            point=rs.coerce3dpoint(Pt.Location)
            if brep.IsPointInside(point,tol,True)==True:
                Result.append(point)
        except:
            pass
    
    CloudPts=Cloud.GetEnumerator()
    ghpythonlib.parallel.run(TestPoints,CloudPts)

    Plane=rs.PlaneFitFromPoints(Result)
    Plane=rs.coerceplane(Plane)
    r,Curves,points=Rhino.Geometry.Intersect.Intersection.BrepPlane(brep,Plane,tol)
    if r==True:
        PlanarSrf=rs.AddPlanarSrf(Curves)

The script runs fine, I am able to scan through and evaluate a set of 20 million points in around 80 seconds. I am going to explore some alternative methods for accelerating this calculation. It looks to me the Cloud.GetEnumerator() for a large cloud is a rather lengthy way to pass the points into the parallel.run().

Also I think I am having issues with garbage collection. In the script I am creating 20 million 3dpoint objects in memory, some of which are used for calculating a best fit plane and never added to the document. Other Rhino objects have .Dispose() function that address this but I don’t see that option for a 3dpoint. After running the script I’m seeing around 4 GB of RAM being used. After panning around in Rhino a bit the garbage collection seems to kick in, realize the 3dpoints are not being used and returns the expected levels. I have also tried the ClearUndo command to clear the history and seen no drop in RAM consumption.

If anyone can point me the right direction that would be great.
Thanks,
5chmidt


(Menno Deij - van Rijswijk) #6

The reason that Point3d does not have dispose, is that it is a value-type, meaning that it is allocated on the stack, rather than the heap. When these go out of scope, the memory will be cleared. Any Point3d in a collection (list, dictionary, …) will need be garbage collected when the collection is collected. You can force garbage collection using System.GC.Collect()

The way I would approach this is with an RTree data structure of the point cloud. Unfortunately, the RTree in RhinoCommon is too limited to support searching on different levels in the tree, but there are python packages that have RTree functionality so you could use these I guess.

The idea is that an RTree contains hierarchical data on bounding boxes that contain points. So, you would traverse the tree from the root on up and test if the bounding box at each level is contained within the Brep. If at one level it is contained within the Brep, you don’t need to test the higher levels of the tree. This way, you can more quickly decide if a (possibly large) set of points is inside the Brep. The overhead of creating the RTree will likely outweigh the is-inside-brep test of 20M points.


[RhinoCommon] RTree is too limited
(Matt Gaydon) #7

I realise this is an old thread but is it possible to edit the code to return just the points inside the Brep, i have tried to edit the code and commented out the last portion regarding the plane but not had much luck.

I have been doing this is a Grasshopper Python Component mind. So maybe this is why its failing

Thanks Matt


#8

Hi Matt,

please note my 2 scripts posted in this thread. Both use multi-processing, the second is a bit more complex.

c.