Search vertexes within a radius [SOLVED]
-
On 13/05/2016 at 09:22, xxxxxxxx wrote:
As noted by supergeordie, for general 3D distance queries, your best bet is some sort of hierarchical spatial partitioning structure (AABB tree, Octree, OOBB tree). The initial overhead is a bit steep and you must construct it carefully but the speed gains are phenomenally better than any other methods since these structures are basically nested divide-and-conquer machines. You can exclude large swaths of areas that aren't in the search very quickly.
-
On 13/05/2016 at 10:31, xxxxxxxx wrote:
I just had another idea to speed things up.
Will look into it and will report it here if it shows some good results.
Remember, I'm not a pro so all my solutions are very naive
-
On 13/05/2016 at 15:07, xxxxxxxx wrote:
If you are interested, here is a set of classes for constructing and using an AABB tree that I have used a couple times (C++!) :
http://www.kuroyumeszone.com/downloads/KDZAABBTree.zip
AABB is "Axis-Aligned Bounding Box" which means, in this case, that the boxes and sub-boxes are aligned along the World/Global axes. Instead of subdividing a space generically (ala Octree), it divides the object's space based on global splitting planes in all three dimensions. These are split iteratively to a particular grain-level forming a hierarchical tree of bounding boxes.
-
On 14/05/2016 at 03:37, xxxxxxxx wrote:
WOW!!! That is soooo complex.
Anyway, how would this new structure (if I manage to be able to implement it) will speedup the calculation of distances?
What I need is to store all the vertexes that are closer that a given distance, and the corresponding distance.
For speed purposes, I can use the squared distance without any problem.
So:If vertex A is closer to vertex B than a given distance , I must store (A, distance to B).
This for all the vertexes of a mesh.
But if A was stored, I don't need to store (B, distance to A).
However, I can't eliminate the B from the list as soon as I find distance match. Because, maybe vertex C or M or T would also need to be stored, because they are closer than the given distance to B.Is there any data structure that would be able to speed up this calculation? Would the AABB tree be of any use?
-
On 14/05/2016 at 08:08, xxxxxxxx wrote:
Hi Rui,
I don't know if this will help you or not. Because I'm not understanding exactly what you're attempting to do.
But when I wanted to rotate the selected polygons in an object. It was very slow because looping through the points array was very slow.
What I needed was a way to only process the selected polygons. And skip over the other ones that were not selected. Which is called "marking" in the SDK neighbor class.Robert showed me how to make my own marking scheme using a true/false array. And it worked great. It runs very fast even on dense meshes.
I know you're not selecting things. But you're still targeting a small subset of points in the object.
So perhaps using a marking array to skip over and not process most of the points is what you need?-ScottA
-
On 14/05/2016 at 08:45, xxxxxxxx wrote:
Here's the full code (minus serial number related classes) of my defunct CollisionDeformer plugin. It was written for R13 and earlier when people hadn't upgraded to R14+ yet (which has a built-in Collision Deformer object). It has the same AABBTree classes but shows how it was used to do quick deformations between objects. You would just need one tree for one object.
http://www.kuroyumeszone.com/downloads/CollisionDeformer.zip
-
On 14/05/2016 at 14:20, xxxxxxxx wrote:
Thank you all for the help that you are providing.
As I have it right now, it is working fine, but I wish I could make it faster. This is the result:
And, the difference between the old method (nested loops) and the new method (search in eight quadrant sub-lists) is as follows:
My current method, is as follows:
• create a list containing just the selected points or all the points if there is no selection.
• from this list, create eight sub-lists, for values of:
-x,-y,-z
-x,-y,+z
-x,+y,-z
-x,+y,+z
+x,-y,-z
+x,-y,+z
+x,+y,-z
+x,+y,+z• cycle through all the points that are in the initial list (p1).
• based on the signal of the coordinates of the point,calculate the index, to decide what sub-list to use.
• go through all the elements of the sub-list (p2).
• if a distance between p1 and p2 is less than a specific radius, store the points (p1, p2, distance) and leave the inner loop (I just need to find the first point p2 that is close enough to p1).
-----------
Is there any way to make this faster? I can't see how a data structure could improve this
-
On 14/05/2016 at 15:41, xxxxxxxx wrote:
You are basically doing a single-level Octree by dividing the points between octants. Tree-structures just go multiple levels which allow you to reduce the test space even faster. For instance, if all of the points of interest are only in one octant, the other seven octants are removed. If the remaining points are in only one section of the smaller 'octants', the others are removed (and so on). Think of it as having a bunch of nested objects with cumulative bounding boxes. If a higher-order bounding box (for a set of objects) is out of range, that object and all of its objects are ignored. You quickly end up only testing points that are (mostly) in range.
I wish that I had time to get an example using the AABBTree under your circumstances but I am working on the finishing touches for a new commercial plugin at the moment.
-
On 14/05/2016 at 16:44, xxxxxxxx wrote:
Mmmmmmm, so I ended up doing a simple (one level) octree
I will try to make it work with deeper levels, by subdividing each octant into another set of octants.
I came up to this method because I remembered about an article I read when I was 16 or 17 years old, about voxels.
Not bad, for someone who is not a professional programmer
Thank you so much, Robert. I will try to implement a deeper octree. -
On 17/05/2016 at 09:06, xxxxxxxx wrote:
I created my own implementation of octrees. I made them just 3 levels deep (I guess it is enough for what I need) and got huge speed gains
Here is a spreadsheet I created with the results from the old method (nested loops), a 1 level deep octree and 3 level deep octrees.
WOW!!!
-
On 17/05/2016 at 09:26, xxxxxxxx wrote:
I would love to see how you did that if possible Rui.
But only if it's not too much trouble.The theory of Octrees are covered a lot on the internet. But I'm not finding much code that works with C4D geometry so far.
If there is a way to show the process in C4D code it would help the rest of us.
But if it's too much code and too much trouble. Don't worry about it.-ScottA
-
On 17/05/2016 at 09:36, xxxxxxxx wrote:
I will prepare a commented snipped of the code I created.
I could have made it recursive, but I decided to make it linear.
Probably my code is not the cleanest or the most beautiful, but it is working
I have to leave now, But I will prepare a listing to publish here as soon as I get back. -
On 17/05/2016 at 10:17, xxxxxxxx wrote:
OK. Thanks.
It doesn't need to be pretty, fancy, or 100% correct. Whatever post will be good enough and appreciated.-ScottA
-
On 17/05/2016 at 11:16, xxxxxxxx wrote:
Good work, Rui! Now you see what I mean. Awesome stuff.
-
On 17/05/2016 at 12:56, xxxxxxxx wrote:
Here is the sample code.
Tell me what you guys think of it# the_selection is a variable that holds a Selection Set tag, # if just a sub-set of points is to be evaluated selection=None if the_selection!=None: selection=the_selection.GetBaseSelect() # now 'selection' is None or a selection bit list points=obj.GetAllPoints() if selection!=None: # create a list with only the selected points pts=[[i,p] for i,p in enumerate(points) if selection.IsSelected(i)] else: # no selection, so create a list with all the points pts=[[i,p] for i,p in enumerate(points)] # ******************************************************************************* # now pts contains a list with all the points coordinates and their index, # in the for [[index1,c4d.Vector()],...[indexN,c4d.Vector()]] # ******************************************************************************* # calculate the min and max values of the bounding box # defined by all points in the list min_x=min(a.x for a in pts) max_x=max(a.x for a in pts) min_y=min(a.y for a in pts) max_y=max(a.y for a in pts) min_z=min(a.z for a in pts) max_z=max(a.z for a in pts) # center coordinates of the cluster of points mid_coords=utils.MixVec(c4d.Vector(min_x,min_y,min_z),c4d.Vector(max_x,max_y,max_z),.5) # calculate all the octrees and their centers # three octrees will be created but it is easy # to create more levels. # # Octree 1 will contain 8 lists of points # Octree 2 will contain 64 lists of points (8*8) # Octree 3 will contain 512 lists of points (8*8*8) # # the oct_center lists contain the center coordinates of the # 2nd and 3rd octrees # # 'distance' is the radius around the points coordinates # in the case of the octrees it is useful to add extra points # that overlap, to create a bit of redundancy. # Otherwise, points that align with the limits of octrees clusters # could be ignored, if the resulting list is to be checked # for distances between points oct1,oct2,oct3,oct_center2,oct_center3 = Calc_All_Octress(pts,min_x,max_x,min_y,max_y,min_z,max_z,distance) # sample code to show how to check for points within a distance: for pts1 in pts: # get the index and the point coordinates i,pt=pts1[0],pts1[1] # get the index of the point in the octree 1 id1=GetQuadrant(pt,mid_coords) # if the point is found in octree #1... if oct1[id1]!=[[],[],[],[],[],[],[],[]]: # get the index of the point in the octree 2 id2=(id1*8)+GetQuadrant(pt,oct_center2[id1]) # if the point is found in octree #2... if oct2[id2]!=[[],[],[],[],[],[],[],[]]: id3=(id2*8)+GetQuadrant(pt,oct_center3[id2]) # if the point is found in octree #3... if oct3[id3]!=[[],[],[],[],[],[],[],[]]: # go through all the points in the relevant octree #3 subset for pts2 in oct3[id3]: # if we dont't need to check the point against itself # keep the following line if pts2[0]==i: continue # calculate the distance between the points # I used LengthSquared for speed reasons dist=(pts2[1]-pt).GetLengthSquared() # only go on if the distance is within the required radius if dist<=distance: # do something with the original point stored in 'pt' # or something with the newly found point, # whose index is stored in pts2[0] and # whose coordinates are stored in pts2[1] # # If just the first match is needed, we can # perform a 'break' and this will get back to the # main loop, increasing speed even further # ******************************************************************************* # ******************************************************************************* # ****************** Here are the octree creating routines ******************** # ******************************************************************************* # ******************************************************************************* def GetQuadrant(p,center) : return (4*(p.x<center.x))+(2*(p.y<center.y))+(p.z<center.z) # ******************************************************************************* def Create_Octs(ls,cx,cy,cz,dist) : # if the list is empty, return a list containing 8 empty lists # just to keep the structure intact if len(ls)==0: return [[],[],[],[],[],[],[],[]] # calculate the new limits with the distance padding added/subtracted cx1=cx-dist cx2=cx+dist cy1=cy-dist cy2=cy+dist cz1=cz-dist cz2=cz+dist # create the 8 new lists from the main volume list ls1=[a for a in ls if a[1].x<cx2 and a[1].y<cy2 and a[1].z<cz2] ls2=[a for a in ls if a[1].x<cx2 and a[1].y<cy2 and a[1].z>=cz1] ls3=[a for a in ls if a[1].x<cx2 and a[1].y>=cy1 and a[1].z<cz2] ls4=[a for a in ls if a[1].x<cx2 and a[1].y>=cy1 and a[1].z>=cz1] ls5=[a for a in ls if a[1].x>=cx1 and a[1].y<cy2 and a[1].z<cz2] ls6=[a for a in ls if a[1].x>=cx1 and a[1].y<cy2 and a[1].z>=cz1] ls7=[a for a in ls if a[1].x>=cx1 and a[1].y>=cy1 and a[1].z<cz2] ls8=[a for a in ls if a[1].x>=cx1 and a[1].y>=cy1 and a[1].z>=cz1] # return a list containing all the calculated sub-lists return [ls8,ls7,ls6,ls5,ls4,ls3,ls2,ls1] # ******************************************************************************* def Calc_All_Octress(pt_list,min_x,max_x,min_y,max_y,min_z,max_z,distance) : # calculate the center of the points volume cx,cy,cz=utils.MixNum(min_x,max_x,.5),utils.MixNum(min_y,max_y,.5),utils.MixNum(min_z,max_z,.5) # create the octree 1 # this first one is straightforward oct1=Create_Octs(pt_list,cx,cy,cz,distance) # create the octree 2 oct2=[] # octree 2 oct_center2=[] # centers 2 # store the previous min and max values because some lists my be empty and # python would complain if we tried to get a min or a max from an empty list xmin,ymin,zmin=min_x,min_y,min_z xmax,ymax,zmax=max_x,max_y,max_z # go through all the sub-lists inside oct1 for l in oct1: # retried the stored min and max values xmin,ymin,zmin=min_x,min_y,min_z xmax,ymax,zmax=max_x,max_y,max_z # if the sub-list is not empty, calculate new min and max values if len(l)>0: min_x=min(a[1].x for a in l ) max_x=max(a[1].x for a in l ) min_y=min(a[1].y for a in l ) max_y=max(a[1].y for a in l ) min_z=min(a[1].z for a in l ) max_z=max(a[1].z for a in l ) # calculate new center of the points volume in the sub-list # and add it to the list of centers of the new cluster for octree 2 cx,cy,cz=utils.MixNum(min_x,max_x,.5),utils.MixNum(min_y,max_y,.5),utils.MixNum(min_z,max_z,.5) oct_center2.append(c4d.Vector(cx,cy,cz)) # create 8 new cluters from this volume and add them to the octree 2 list oct2.extend(Create_Octs(l,cx,cy,cz,distance)) # create the octree 3 oct3=[] # octree 3 oct_center3=[] # centers 3 # go through all the sub-lists inside oct2 for l in oct2: # retried the stored min and max values xmin,ymin,zmin=min_x,min_y,min_z xmax,ymax,zmax=max_x,max_y,max_z # if the sub-list is not empty, calculate new min and max values if len(l)>0: min_x=min(a[1].x for a in l ) max_x=max(a[1].x for a in l ) min_y=min(a[1].y for a in l ) max_y=max(a[1].y for a in l ) min_z=min(a[1].z for a in l ) max_z=max(a[1].z for a in l ) # calculate new center of the points volume in the sub-list # and add it to the list of centers of the new cluster for octree 3 cx,cy,cz=utils.MixNum(min_x,max_x,.5),utils.MixNum(min_y,max_y,.5),utils.MixNum(min_z,max_z,.5) oct_center3.append(c4d.Vector(cx,cy,cz)) # create 8 new cluters from this volume and add them to the octree 3 list oct3.extend(Create_Octs(l,cx,cy,cz,distance)) # seeing how the code for calculating octree 2 and octree 3 is very similar, # it would be wuite easy to implement deeper octree calculations # for my purposes, three levels are enough # return the octrees and the octree centers lists return oct1,oct2,oct3,oct_center2,oct_center3 -
On 17/05/2016 at 15:49, xxxxxxxx wrote:
Thanks Rui,
I think you're missing [1] in one of your min max code blocks.
It would not run using a.x, etc...When I change the distance value. I don't see any difference in what gets grabbed.
Here is my entire code that I'm using in the script manager.
Create a cube with 8 subdivisions in x,y,z. Then make it editable. Then run this script.
At the end of the code I tested selecting the results from oct1,oct2,oct3,oct_center2, and oct_center3.
Seems to work as expected. But oct1 and oct2 select very few points before oct3 finally selects the rest of them. Not sure if this is normal?
But the main thing is the distance doesn't seem to do anything the way I'm using it.
Am I using it wrong?import c4d from c4d import utils # ******************************************************************************* # ****************** Here are the octree creating routines ******************** # ******************************************************************************* def GetQuadrant(p,center) : return (4*(p.x<center.x))+(2*(p.y<center.y))+(p.z<center.z) def Create_Octs(ls,cx,cy,cz,dist) : #If the list is empty, return a list containing 8 empty lists just to keep the structure intact if len(ls)==0: return [[],[],[],[],[],[],[],[]] #Calculate the new limits with the distance padding added/subtracted cx1 = cx-dist cx2 = cx+dist cy1 = cy-dist cy2 = cy+dist cz1 = cz-dist cz2 = cz+dist #Create the 8 new lists from the main volume list ls1=[a for a in ls if a[1].x<cx2 and a[1].y<cy2 and a[1].z<cz2] ls2=[a for a in ls if a[1].x<cx2 and a[1].y<cy2 and a[1].z>=cz1] ls3=[a for a in ls if a[1].x<cx2 and a[1].y>=cy1 and a[1].z<cz2] ls4=[a for a in ls if a[1].x<cx2 and a[1].y>=cy1 and a[1].z>=cz1] ls5=[a for a in ls if a[1].x>=cx1 and a[1].y<cy2 and a[1].z<cz2] ls6=[a for a in ls if a[1].x>=cx1 and a[1].y<cy2 and a[1].z>=cz1] ls7=[a for a in ls if a[1].x>=cx1 and a[1].y>=cy1 and a[1].z<cz2] ls8=[a for a in ls if a[1].x>=cx1 and a[1].y>=cy1 and a[1].z>=cz1] #Return a list containing all the calculated sub-lists return [ls8,ls7,ls6,ls5,ls4,ls3,ls2,ls1] def Calc_All_Octress(pt_list,min_x,max_x,min_y,max_y,min_z,max_z,distance) : #Calculate the center of the points volume cx,cy,cz = utils.MixNum(min_x,max_x,.5),utils.MixNum(min_y,max_y,.5),utils.MixNum(min_z,max_z,.5) #Create octree 1 #This first one is straightforward oct1 = Create_Octs(pt_list,cx,cy,cz,distance) #Create octree 2 oct2 = [] # octree 2 oct_center2 = [] # centers 2 #Store the previous min and max values because some lists my be empty and #python would complain if we tried to get a min or a max from an empty list xmin,ymin,zmin = min_x,min_y,min_z xmax,ymax,zmax = max_x,max_y,max_z #Go through all the sub-lists inside oct1 for l in oct1: #Retried the stored min and max values xmin,ymin,zmin=min_x,min_y,min_z xmax,ymax,zmax=max_x,max_y,max_z #If the sub-list is not empty, calculate new min and max values if len(l)>0: min_x=min(a[1].x for a in l ) max_x=max(a[1].x for a in l ) min_y=min(a[1].y for a in l ) max_y=max(a[1].y for a in l ) min_z=min(a[1].z for a in l ) max_z=max(a[1].z for a in l ) #Calculate new center of the points volume in the sub-list #Then add it to the list of centers of the new cluster for octree 2 cx,cy,cz = utils.MixNum(min_x,max_x,.5),utils.MixNum(min_y,max_y,.5),utils.MixNum(min_z,max_z,.5) oct_center2.append(c4d.Vector(cx,cy,cz)) #Create 8 new cluters from this volume and add them to the octree 2 list oct2.append(Create_Octs(l,cx,cy,cz,distance)) #The octree 2 list ended up with an extra set of [] around it #Remove the internal list from outer list to keep the structure if the octree 2 equal to the structure of octree 1 oct2=[val for sublist in oct2 for val in sublist] #Create the octree 3 oct3=[] #Octree 3 oct_center3=[] #Centers 3 #Go through all the sub-lists inside oct2 for l in oct2: #Retried the stored min and max values xmin,ymin,zmin = min_x,min_y,min_z xmax,ymax,zmax = max_x,max_y,max_z #If the sub-list is not empty, calculate new min and max values if len(l)>0: min_x = min(a[1].x for a in l ) max_x = max(a[1].x for a in l ) min_y = min(a[1].y for a in l ) max_y = max(a[1].y for a in l ) min_z = min(a[1].z for a in l ) max_z = max(a[1].z for a in l ) #Calculate new center of the points volume in the sub-list #Then add it to the list of centers of the new cluster for octree 3 cx,cy,cz = utils.MixNum(min_x,max_x,.5),utils.MixNum(min_y,max_y,.5),utils.MixNum(min_z,max_z,.5) oct_center3.append(c4d.Vector(cx,cy,cz)) #Create 8 new cluters from this volume and add them to the octree 3 list oct3.append(Create_Octs(l,cx,cy,cz,distance)) #The octree 3 list ended up with an extra set of [] around it #Remove the internal list from outer list to keep the structure if the octree 3 equal to the structure of octree 1 and octree 2 oct3 = [val for sublist in oct3 for val in sublist] #Seeing how the code for calculating octree 2 and octree 3 is very similar, #It would be easy to implement deeper octree calculations. But for my purposes three levels are enough #Return the octrees and the octree centers lists return oct1,oct2,oct3,oct_center2,oct_center3 def main() : #The selection variable holds a Selection Set tag if just a sub-set of points is to be evaluated #Selection is either None. Or a selection bit list selection = None if selection != None: selection = the_selection.GetBaseSelect() points = op.GetAllPoints() #Create a list named "pts" with only the selected points if selection != None: pts = [[i,p] for i,p in enumerate(points) if selection.IsSelected(i)] else: #No selection, so create a list using all the points in the object pts = [[i,p] for i,p in enumerate(points)] # ******************************************************************************* # pts now contains a list with all the points coordinates and their index, # in the for [[index1,c4d.Vector()],...[indexN,c4d.Vector()]] # **************************************************************************** #Calculate the min and max values of the object's bounding box defined by the points in the list min_x = min(a[1].x for a in pts) max_x = max(a[1].x for a in pts) min_y = min(a[1].y for a in pts) max_y = max(a[1].y for a in pts) min_z = min(a[1].z for a in pts) max_z = max(a[1].z for a in pts) #This gets the center coordinates of the cluster of points mid_coords = utils.MixVec(c4d.Vector(min_x,min_y,min_z), c4d.Vector(max_x,max_y,max_z), 0.5) # calculate all the octrees and their centers # three octrees will be created but it is easy to create more levels. # # Octree 1 will contain 8 lists of points # Octree 2 will contain 64 lists of points (8*8) # Octree 3 will contain 512 lists of points (8*8*8) # # the oct_center lists contain the center coordinates of the 2nd and 3rd octrees # # 'distance' is the radius around the points coordinates # In the case of the octrees it is useful to add extra points that overlap, to create a bit of redundancy # Otherwise, points that align with the limits of octrees clusters could be ignored, # if the resulting list is to be checked for distances between points distance = 1.0 oct1,oct2,oct3,oct_center2,oct_center3 = Calc_All_Octress(pts,min_x,max_x,min_y,max_y,min_z,max_z,distance) #Now let's select the points that the above code grabbed to visually see the results #Note: using oct3 seems to crash with array out of bounds error! PointS = c4d.BaseSelect() for p in xrange(len(oct1)) : #<----Also try oct2,oct3,oct_center2,oct_center3 PointS.Select(p) PointS.CopyTo(op.GetPointS()) c4d.EventAdd() if __name__=='__main__': main()-ScottA
-
On 18/05/2016 at 00:53, xxxxxxxx wrote:
The octrees don't store the points found within a distance.
They just subdivide the volume space in smaller chunks for faster verification. The distance provided for the creation of the octrees is just used to adjust the possible error that would occur, near the limits of the chunks.
Try the code below, that uses the octrees to go through all the points list and selects the points within a distance.
Oh, if you are using a cube to test things, since the mesh is regular, eventually, all points will find a match within the radius and everything will get selected.
Create more irregular shapes to test.distance =20.0 oct1,oct2,oct3,oct_center2,oct_center3 = Calc_All_Octress(pts,min_x,max_x,min_y,max_y,min_z,max_z,distance) #Now let's select the points that the above code grabbed to visually see the results #Note: using oct3 seems to crash with array out of bounds error! PointS = c4d.BaseSelect() PointS.DeselectAll() for pts1 in pts: # get the index and the point coordinates i,pt=pts1[0],pts1[1] # get the index of the point in the octree 1 id1=GetQuadrant(pt,mid_coords) # if the point is found in octree #1... if oct1[id1]!=[[],[],[],[],[],[],[],[]]: # get the index of the point in the octree 2 id2=(id1*8)+GetQuadrant(pt,oct_center2[id1]) # if the point is found in octree #2... if oct2[id2]!=[[],[],[],[],[],[],[],[]]: id3=(id2*8)+GetQuadrant(pt,oct_center3[id2]) # if the point is found in octree #3... if oct3[id3]!=[[],[],[],[],[],[],[],[]]: # go through all the points in the relevant octree #3 subset for pts2 in oct3[id3]: # don't check a point against itself if pts2[0]==i: continue # calculate the distance between the points dist=(pts2[1]-pt).GetLength() # only go on if the distance is within the required radius if dist<=distance: PointS.Select(i) PointS.CopyTo(op.GetPointS()) c4d.EventAdd() -
On 18/05/2016 at 02:02, xxxxxxxx wrote:
Thanks for posting guys. This looks like really useful stuff!
-
On 18/05/2016 at 07:50, xxxxxxxx wrote:
By the way:
I found a rather old link in my bookmarks. Maybe, it's interesting for you too:
http://thomasdiewald.com/blog/?p=1488 -
On 18/05/2016 at 08:47, xxxxxxxx wrote:
Thanks Rui.
I wasn't sure what to do with that for() loop in your original code (if it was optional or not).
But after inserting your second for() loop example. It works as expected.Thanks for posting the code,
-ScottA