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)
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)
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).
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
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.
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
Hii @5chmidt, I really do not have more experience in scripting However this script looks really good and even I can use it. I am currently writing my master thesis. I have to do some progress monitoring stuff. I have to compare BIM As-built model (i.e. point cloud) vs BIM As-planned model. I was also trying the same thing with âIspointinside_brepâ then it shoud turn out to color indicating that "brep have some points. but in progress monitoring its not somehow working.
could you please tell me what is this script is about? Also, I wants to know " Turn object into a Point cloud" which object we should turn out into a point cloud? My questions might seems to be silly. but it would be great if you could reply it.
Many Thanks
parimal
Is used to cast the point cloud which is stored as a more memory efficient object into an object that has enumerated points we can iterate through for the purposes of our test.
It will take a lot of work to make this script function to compare a BIM model with a point cloud. The script as it currently uses a hand drawn brep to section out a segment of points from within a larger cloud. To do this efficiently you would need a prep script to make rough boundaries around areas of the model where you expect to see shapes. If the BIM model specifies a pipe location you would want to grab only the points in the area where you expect a pipe to be present, so you would need a brep to test them.
Once you have the points you want to analyze you can run a best algorithm. Rhino includes rs.PlaneFitFromPoints which is helpful, but you will probably need to fit other types of shape. Here is some documentation on fitting a pipe: cylinder best fit. You would presumably need to develop other methods to fit I-beams and other architectural features.
@5chmidt Thankyou so much for your reply. I am testing on the road model and with Grashopper defination I have already merged and geomtry of the road. Its automatically merged as its Geo-referenced so. In this road there are 6 layers, I have control over each layer. I created this geometry using Grasshopper defination and also I have created the closed breps. Bascially its a layer detection process from Geomtry to PCL.
Looking forward to your reply. Thankyou so much
Could you please suggest me something I am also sharing pictures. Can I implement this âfit.planeâ approach here?
It looks like you already have your target road segments from your grasshopper model. You may not need to do a best fit plane. Instead you can create a bounding box around the segment of road you want to analyze and select a subset of points that are near the geometry you are matching to. Once you have a smaller subset of points you can find the distance from the point to the geometry and color the points accordingly. IE green = within tolerance, orange = slightly out, red = way out of tolerance.
As an additional note, Iâm not sure grasshopper is the best tool for this task. If you are creating your road models parametrically, grasshopper it maybe a good tool, but running large point cloud analysis in grasshopper will really weigh you down. In my previous example it was taking me 80 seconds to segment up 20 million points. Looking at your example, you will have dozens of road segments that would all need similar analysis. I generally try to keep grasshopper scripts light enough to run an iteration in a few seconds. If something is going to remain static, (like the as planned road) bake it out and so it is not part of the point cloud calculation script. Or if you can write your program in C#, python, or any of the programming languages Rhino supports it could help to prevent your grasshopper definition from becoming a giant spaghetti monster that take 15 minutes to redraw.
@parimalvala5022 what are your inputs to do this analysis? Is it a text file containing every XYZ describing the point cloud from your BIM model? And the meshes or breps you want to test for containment?
@rawitscher-torres@5chmidt
Thankyou for your reply.
I used âVolvoxâ plugin to import pointcloud in Grasshopper. Volvox plugin supports (.E57) point cloud files.
In comparsion I have to compare point cloud with Road model. I have created a Grasshopper defination to create road layers in Grasshopper (i.e. Asphalt layer, Base, subbase layer etc). In total there are 6 layers.
I have developed a script which is iterating the âpoints in brepâ
In total there are 6 road layers - In each layer there are currently 51 segments and the point cloud after subsampling process there are 21000 points. Now there is only one problem about " Threshold points".Threshold points are nothing but the " minimum number of points in each brep to iterate"
For example- I set threshold point values to 20.Then once the computation starts- if in one brep there are 19 points found then it will not change the color other wise it will change the color.
How to set this threshold points? (There should be a convincing reason to set the threshold point because if I set the threshold point to 25 then anyone can ask why 25?
Am I going on rightway. Kindly request you to suggest.
@5schmidt you are right computation times takes almost 15 minutes for one iteration of one road layer.
I am also sharing with you guys the sample GH defination. If you guys can help me out with this it would be a great help.
I would like to try this approach, could you please show me something visually by sharing some screen shots or drawing? Here you mentioned me that âfind the distance from point to Geometryâ here geomtry means âSegmentâ or the full road layer?
Lookinf forward to know more on this.
Thankyou verymuch in advance. It would be great if you could show me by drawing something or any post if you know related to this kind of task.
Thanks
This is a great thread, itâs great to know you can parse 20m point cloud in around 80s.
There are a lot of interesting features in Cloud Compare (free software) that are let down by lack of scriptability (which Gh has in spades). Have a look/read of their implementation of Octree (similar o RTree?) for point cloud rendering, I donât know if Rhino has anything comparable. The number of points in each box has implications for performance.
Below is a coloured comparison of the difference between two clouds, a very handy visual indicator of divergence.
Iâm not sure about fastest way of performing similar analysis in GH, something like the Fast Mesh Ray analysis in Sasquash might be an appropriate method, if this approach helps your problem.
I did some quick tests, on clipping a point cloud with given meshes.
The thing to clarify here is that I did not use the PointCloud class , but just a very large collection of points for a proof of concept.
The nice thing as well is that I am not outputting the actual geometry in the view port, so it has very fast display.
It would be great if you could tell me, for example, what do you do with the output from Volvox and why is your GH file taking 15 minutes to update, what are you trying to do?.
Perhaps soon In can code a new PointCloud class with some new implementations. But would be nice to hear your goals.
Below are some images of clipping 64 million points and only displaying those which are inside a geometry. For 64 million points the exact method takes 62 seconds and the whole code takes 1.1 minutes in total to execute.
@rawitscher-torres Thanks for your reply.
Basically I used volovx to import the point cloud (i.e .E57 file format point cloud) and also I have used that plugin for subsampling process and that plugin has much more capabilities. I am sharing with you a presentation in order to get some idea.
Its taking time because after subsampling process there are 21000 points and its checking in all Breps.
there are in total 51 Breps and the iteration process in script checking in each brep for points.
In my thesis, I have to do progress monitoring. In simple word, How much construction work is finished on that particular date (i.e. scan date)?
I got scans for 2 different dates. So, idea was to create Brep in road geometry and if the âpoints found in Brepâ considered as a âBuiltâ. In my first scan, there was an excavation work . So bascially in road there are 6 different layers. I have created Breps in each layer. By using my script,I have found that maximum number of âpoints found in 4th layerâ meaning "excavation is done till that 4thlayer.
problem is that I have run that script by setting âThresholdâ about 25 points. There is not a concrete reason that why 25? However, its giving correct result. there should be some relation with âpoints densityâ per square metre.
If you have any another idea to implement in Grashopper then I am looking forward to do it.
Looking forward to your reply. If you still didnât get in your head, I can provide more details.Point cloud detecttion_process.pdf (1.7 MB)
I am not exactly up-to-date on the latest point cloud processing plugins for grasshopper. But I suspect it will be rather difficult to implement the specific functionality your are looking without coding a few of your own grasshopper components.
Here is a great tutorial on creating a python component.
If you are not an experienced coder, I would recommend starting with python it is not strictly typed and the syntax is nicely readable.
