sweep and prune OK

models/boneload.py

@@ -15,6 +15,7 @@
 import os
 import math
 import vtk
+import time
 colors = vtk.vtkNamedColors()
 # import numpy
 # from vtk.util.numpy_support import vtk_to_numpy
@@ -129,7 +130,7 @@ def load_off(meshfile):
     return nodes, tris
 
 
-def identify_nodes(coords, all_no, all_coords, eps=1e-3):
+def identify_nodes_brutal(coords, all_no, all_coords, eps=1e-3):
     """identify the nodes based on their coordinates "coords"
     among the list of coordinates "all_coords"
     all_no, given as argument, is an array of tags for the nodes defined in all_coords
@@ -141,6 +142,9 @@ def identify_nodes(coords, all_no, all_coords, eps=1e-3):
 
     returns ntags, a list of node tags of the identified nodes
     """
+    cpu_start = time.perf_counter()
+    print(f'identify_nodes_brutal...')
+
     notfound = []
     ntags = []
     for n in coords:  # for all nodes
@@ -161,11 +165,97 @@ def identify_nodes(coords, all_no, all_coords, eps=1e-3):
         if not found:
             notfound.append(n)
 
+    cpu_stop = time.perf_counter()
+
+    # print stats
+    print(f'identify_nodes_brutal done in {cpu_stop-cpu_start} seconds')
+
     if(len(notfound) != 0):
         print('WARNING: these nodes have not been identified:')
         for n in notfound:
             print(f'\tn={n}')
-        # raise Exception(f'ERROR: {len(notfound)} nodes have not been found!')
+        raise Exception(f'ERROR: {len(notfound)} nodes have not been found!')
+    return ntags
+
+
+def identify_nodes_SaP(coords, all_no, all_coords, eps=1e-3):
+    """identify the nodes based on their coordinates "coords"
+    among the list of coordinates "all_coords"
+    all_no, given as argument, is an array of tags for the nodes defined in all_coords
+
+    The algorithm used is "sweep & prune" along X.
+
+    Comparision with the brutal approach: (6567 nodes, 12696 triangles)
+        identify_nodes_brutal done in 26.600111499999997 seconds
+        sweep and prune done in 0.30412259999999947 seconds
+        (48500 tests instead of 818707890)
+
+    returns ntags, a list of node tags of the identified nodes
+    """
+    print('identify_nodes_SaP...')
+    cpu_start = time.perf_counter()
+
+    # define nod structure
+    class Nod:
+        def __init__(self, x, idx=-1):  # selected points have  a -1 index
+            self.x = x
+            self.idx = idx
+
+        def __str__(self):
+            return f'{self.idx}: {self.x}'
+
+    # build node list
+    nods = []
+    for i, p in enumerate(all_coords):
+        nods.append(Nod(p, all_no[i]))
+    for p in coords:
+        nods.append(Nod(p, -1))
+
+    # sort nodes according to their x coordinate
+    nods.sort(key=lambda nod: nod.x[0])
+
+    ntags = []
+
+    eps = 1e-3    # geometrical tolerance
+    ntests = 0
+    nnods = len(nods)
+    i = 0
+    while i < nnods:
+        n = nods[i]
+        x = n.x[0]
+        i2 = i
+        while True:
+            i2 += 1
+            if i2 == nnods:
+                break
+            n2 = nods[i2]
+            x2 = n2.x[0]
+            if abs(x2-x) > eps:  # test X
+                break
+            ntests += 1
+            if n.idx >= 0 and n2.idx >= 0:  # nodes from different sets
+                continue
+            if abs(n.x[1]-n2.x[1]) > eps:  # test Y
+                continue
+            if abs(n.x[2]-n2.x[2]) > eps:  # test Z
+                continue
+            no = max(n.idx, n2.idx)
+            ntags.append(no)
+            break # hyp: a max of 1 pair should be found per node
+            # print(f'{no} is selected')
+        i += 1
+
+    cpu_stop = time.perf_counter()
+
+    # print stats
+    print(f'sweep and prune done in {cpu_stop-cpu_start} seconds')
+    print(f'({ntests} tests instead of {len(coords)*len(all_coords)})')
+    print(f'{len(coords)} to be identified')
+    print(f'{len(ntags)} successfully identified')
+    if(len(coords) != len(ntags)):
+        print('WARNING: some nodes have not been identified!')
+        raise Exception(f'ERROR: in identify_nodes_SaP!')
+
     return ntags
 
 
@@ -973,82 +1063,38 @@ if __name__ == "__main__":
     # print(help(vtk.vtkCellArray))
     # import sys; sys.exit()
 
-    # basic implementation of "sweep & prune" algorithm.
+    # test of the "sweep & prune" algorithm.
+
+    npts = 10
 
     # build a list of points with random coordinates
     import random
-    all_pts = []
-    for i in range(10):
-        all_pts.append(Pt([random.random(), random.random(), random.random()]))
+    all_coords = []
+    all_no = []
+    for i in range(npts):
+        all_coords.append([random.random(), random.random(), random.random()])
+        all_no.append(i)
 
     # build a list of selected points
-    pts = []
-    pts.append(all_pts[1])
-    pts.append(all_pts[2])
-    pts.append(all_pts[5])
-    pts.append(all_pts[8])
+    coords = []
+    coords.append(all_coords[1])
+    coords.append(all_coords[2])
+    coords.append(all_coords[5])
+    coords.append(all_coords[8])
 
     # add some variants of the selected pts
-    for p in pts[0:3]:
-        all_pts.append(Pt([p.x[0], 0.0, 1.0]))
-        all_pts.append(Pt([p.x[0], p.x[1], 1.0]))
-
+    for i, p in enumerate(all_coords[0:3]):
+        all_coords.append([p[0], 2.0, 2.0]); all_no.append(npts); npts+=1
+        all_coords.append([p[0], p[1], 2.0]); all_no.append(npts); npts+=1
+        
+    # display all points
     print(f'all points:')
-    for no, p in enumerate(all_pts):
+    for no, p in zip(all_no, all_coords):
         print(f'\t{no}: {p}')
 
-    print(f'selected:')
-    for p in pts:
+    # display selected coordinates
+    print(f'selected coordinates:')
+    for p in coords:
         print(f'\t{p}')
 
-    class Nod:
-        def __init__(self, x, idx=-1):  # selected points have  a -1 index
-            self.x = x
-            self.idx = idx
-
-        def __str__(self):
-            return f'{self.idx}: {self.x}'
-
-    # build node list
-    nods = []
-    for no, p in enumerate(all_pts):
-        nods.append(Nod(p.x, no))
-    for p in pts:
-        nods.append(Nod(p.x, -1))
-
-    # sort nodes according to their x coordinate
-    nods.sort(key=lambda nod: nod.x[0])
-
-    print('sorted list:')
-    for n in nods:
-        print(f'\t{n}')
-
-    print('sweep and prune:')
-    eps = 1e-3    # geometrical tolerance
-    ntests = 0
-    nnods = len(nods)
-    i = 0
-    while i < nnods:
-        n = nods[i]
-        x = n.x[0]
-        i2 = i
-        while True:
-            i2 += 1
-            if i2 == nnods:
-                break
-            n2 = nods[i2]
-            x2 = n2.x[0]
-            if abs(x2-x) > eps:  # test X
-                break
-            ntests += 1
-            if n.idx >= 0 and n2.idx >= 0:  # nodes from different sets
-                continue
-            if abs(n.x[1]-n2.x[1]) > eps:  # test Y
-                continue
-            if abs(n.x[2]-n2.x[2]) > eps:  # test Z
-                continue
-            no = max(n.idx, n2.idx)
-            print(f'{no} is selected')
-        i += 1
-    print(
-        f'sweep and prune done ({ntests} tests instead of {len(pts)*len(all_pts)})')
+    identify_nodes_SaP(coords, all_no, all_coords, eps=1e-3)

models/mandiblemodel.py

@@ -314,7 +314,8 @@ def create_group(mesh_file_or_coord, all_no, all_coords, domain, groups, gname=N
     else:
         raise Exception(f'Data is not understood: {mesh_file_or_coord}')
     # identify nodes of the mesh among the given vertices of the main mesh
-    ntags = boneload.identify_nodes(nodes, all_no, all_coords)
+    # ntags = boneload.identify_nodes_brutal(nodes, all_no, all_coords)
+    ntags = boneload.identify_nodes_SaP(nodes, all_no, all_coords)
 
     # create a metafor group
     groupset = domain.getGeometry().getGroupSet()