add more robust fix for edge cases

voxelise.py

@@ -16,6 +16,7 @@ import reebgraph
 import pymeshfix as pf
 import trimesh as tm
 from shapely.geometry import Point
+from scipy.spatial import KDTree
 
 
 def to_trimesh(mesh):
@@ -173,22 +174,25 @@ def vox_section_area(voxels, mesh_file_name, vox_dim=3):
             if poly.contains(Point(c[:2])):
                 SecArea[i] = poly.area
 
-                if abs(SecArea[i-1] - default_high) < 1:
-                    c1 = voxels.iloc[i-1,1:4].to_numpy()
-                    d1 = abs(sum(c-c1))/vox_dim
-                    print("Trying to fix voxel " + str(i-1) + ", dist: " + str(d1))
-                    if d1 <1.001:
-                        SecArea[i-1] = SecArea[i]
-                        print("voxel " + str(i-1) + " fixed")
-
         if SecArea[i] == 0:
             print("Uh OH unlinked voxel " + str(i))
             SecArea[i] = default_high
 
-        #if z>40:
-        #    print(set(SecArea))
-        #    break
     print("Unlinked Voxels: {}/{}".format(sum(SecArea==default_high), len(SecArea)))
+
+    fixed = 0
+    pts = voxels[['x', 'y', 'z']].to_numpy()
+    kd = KDTree(pts)
+    for i in np.where(SecArea==default_high)[0]:
+        c = pts[i,:]
+        inds = np.array(kd.query_ball_point(c, 1.1*vox_dim))
+        inds = inds[pts[inds,2]==c[2]]
+        SecArea[i] = max(SecArea[inds])
+        if SecArea[i] != default_high:
+            fixed += 1
+
+    print("Fixed {} of those".format(fixed))
+
     visualise_voxels(voxels, stresses=SecArea)