integrated kdtree

git-svn-id: https://pycam.svn.sourceforge.net/svnroot/pycam/trunk@52 bbaffbd6-741e-11dd-a85d-61de82d9cad9

Makefile

-VER=0.1.8
+VER=0.1.9
 
 all:
 

Samples/STL/TestModel.cfg

@@ -9,9 +9,11 @@ MaxZ=3
 [config]
 ToolRadius=1
 TorusRadius=0.25
-Layers=4
-Samples=50
+Layers=1
+Samples=20
 Lines=20
 Cutter=SphericalCutter
 PathGenerator=PushCutter
-PathProcessor=ZigZagCutter
+PathProcessor=ContourCutter
+Direction=xy
+

pycam/Cutters/CylindricalCutter.py

@@ -115,26 +115,16 @@ class CylindricalCutter(BaseCutter):
 
     def intersect(self, direction, triangle):
         (cl_t,d_t) = self.intersect_circle_triangle(direction, triangle)
-        (cl_p1,d_p1) = self.intersect_circle_vertex(direction, triangle.p1)
-        (cl_p2,d_p2) = self.intersect_circle_vertex(direction, triangle.p2)
-        (cl_p3,d_p3) = self.intersect_circle_vertex(direction, triangle.p3)
-        (cl_e1,d_e1) = self.intersect_circle_edge(direction, Line(triangle.p1,triangle.p2))
-        (cl_e2,d_e2) = self.intersect_circle_edge(direction, Line(triangle.p2,triangle.p3))
-        (cl_e3,d_e3) = self.intersect_circle_edge(direction, Line(triangle.p3,triangle.p1))
         d = INFINITE
         cl = None
         if d_t < d:
             d = d_t
             cl = cl_t
-        if d_p1 < d:
-            d = d_p1
-            cl = cl_p1
-        if d_p2 < d:
-            d = d_p2
-            cl = cl_p2
-        if d_p3 < d:
-            d = d_p3
-            cl = cl_p3
+        if cl:
+            return (cl,d)
+        (cl_e1,d_e1) = self.intersect_circle_edge(direction, triangle.e1)
+        (cl_e2,d_e2) = self.intersect_circle_edge(direction, triangle.e2)
+        (cl_e3,d_e3) = self.intersect_circle_edge(direction, triangle.e3)
         if d_e1 < d:
             d = d_e1
             cl = cl_e1
@@ -144,13 +134,26 @@ class CylindricalCutter(BaseCutter):
         if d_e3 < d:
             d = d_e3
             cl = cl_e3
+        if cl:
+            return (cl,d)
+        (cl_p1,d_p1) = self.intersect_circle_vertex(direction, triangle.p1)
+        (cl_p2,d_p2) = self.intersect_circle_vertex(direction, triangle.p2)
+        (cl_p3,d_p3) = self.intersect_circle_vertex(direction, triangle.p3)
+        if d_p1 < d:
+            d = d_p1
+            cl = cl_p1
+        if d_p2 < d:
+            d = d_p2
+            cl = cl_p2
+        if d_p3 < d:
+            d = d_p3
+            cl = cl_p3
+        if cl:
+            return (cl,d)
         if direction.x != 0 or direction.y != 0:
-            (cl_p1,d_p1) = self.intersect_cylinder_vertex(direction, triangle.p1)
-            (cl_p2,d_p2) = self.intersect_cylinder_vertex(direction, triangle.p2)
-            (cl_p3,d_p3) = self.intersect_cylinder_vertex(direction, triangle.p3)
-            (cl_e1,d_e1) = self.intersect_cylinder_edge(direction, Line(triangle.p1,triangle.p2))
-            (cl_e2,d_e2) = self.intersect_cylinder_edge(direction, Line(triangle.p2,triangle.p3))
-            (cl_e3,d_e3) = self.intersect_cylinder_edge(direction, Line(triangle.p3,triangle.p1))
+            (cl_e1,d_e1) = self.intersect_cylinder_edge(direction, triangle.e1)
+            (cl_e2,d_e2) = self.intersect_cylinder_edge(direction, triangle.e2)
+            (cl_e3,d_e3) = self.intersect_cylinder_edge(direction, triangle.e3)
             if d_p1 < d:
                 d = d_p1
                 cl = cl_p1
@@ -160,6 +163,11 @@ class CylindricalCutter(BaseCutter):
             if d_p3 < d:
                 d = d_p3
                 cl = cl_p3
+            if cl:
+                return (cl,d)
+            (cl_p1,d_p1) = self.intersect_cylinder_vertex(direction, triangle.p1)
+            (cl_p2,d_p2) = self.intersect_cylinder_vertex(direction, triangle.p2)
+            (cl_p3,d_p3) = self.intersect_cylinder_vertex(direction, triangle.p3)
             if d_e1 < d:
                 d = d_e1
                 cl = cl_e1

pycam/Cutters/SphericalCutter.py

@@ -117,26 +117,16 @@ class SphericalCutter(BaseCutter):
 
     def intersect(self, direction, triangle):
         (cl_t,d_t) = self.intersect_sphere_triangle(direction, triangle)
-        (cl_p1,d_p1) = self.intersect_sphere_vertex(direction, triangle.p1)
-        (cl_p2,d_p2) = self.intersect_sphere_vertex(direction, triangle.p2)
-        (cl_p3,d_p3) = self.intersect_sphere_vertex(direction, triangle.p3)
-        (cl_e1,d_e1) = self.intersect_sphere_edge(direction, Line(triangle.p1,triangle.p2))
-        (cl_e2,d_e2) = self.intersect_sphere_edge(direction, Line(triangle.p2,triangle.p3))
-        (cl_e3,d_e3) = self.intersect_sphere_edge(direction, Line(triangle.p3,triangle.p1))
         d = INFINITE
         cl = None
         if d_t < d:
             d = d_t
             cl = cl_t
-        if d_p1 < d:
-            d = d_p1
-            cl = cl_p1
-        if d_p2 < d:
-            d = d_p2
-            cl = cl_p2
-        if d_p3 < d:
-            d = d_p3
-            cl = cl_p3
+        if cl:
+            return (cl,d)
+        (cl_e1,d_e1) = self.intersect_sphere_edge(direction, triangle.e1)
+        (cl_e2,d_e2) = self.intersect_sphere_edge(direction, triangle.e2)
+        (cl_e3,d_e3) = self.intersect_sphere_edge(direction, triangle.e3)
         if d_e1 < d:
             d = d_e1
             cl = cl_e1
@@ -146,13 +136,26 @@ class SphericalCutter(BaseCutter):
         if d_e3 < d:
             d = d_e3
             cl = cl_e3
+        if cl:
+            return (cl,d)
+        (cl_p1,d_p1) = self.intersect_sphere_vertex(direction, triangle.p1)
+        (cl_p2,d_p2) = self.intersect_sphere_vertex(direction, triangle.p2)
+        (cl_p3,d_p3) = self.intersect_sphere_vertex(direction, triangle.p3)
+        if d_p1 < d:
+            d = d_p1
+            cl = cl_p1
+        if d_p2 < d:
+            d = d_p2
+            cl = cl_p2
+        if d_p3 < d:
+            d = d_p3
+            cl = cl_p3
+        if cl:
+            return (cl,d)
         if direction.x != 0 or direction.y != 0:
-            (cl_p1,d_p1) = self.intersect_cylinder_vertex(direction, triangle.p1)
-            (cl_p2,d_p2) = self.intersect_cylinder_vertex(direction, triangle.p2)
-            (cl_p3,d_p3) = self.intersect_cylinder_vertex(direction, triangle.p3)
-            (cl_e1,d_e1) = self.intersect_cylinder_edge(direction, Line(triangle.p1,triangle.p2))
-            (cl_e2,d_e2) = self.intersect_cylinder_edge(direction, Line(triangle.p2,triangle.p3))
-            (cl_e3,d_e3) = self.intersect_cylinder_edge(direction, Line(triangle.p1,triangle.p3))
+            (cl_e1,d_e1) = self.intersect_cylinder_edge(direction, triangle.e1)
+            (cl_e2,d_e2) = self.intersect_cylinder_edge(direction, triangle.e2)
+            (cl_e3,d_e3) = self.intersect_cylinder_edge(direction, triangle.e3)
             if d_p1 < d:
                 d = d_p1
                 cl = cl_p1
@@ -162,6 +165,11 @@ class SphericalCutter(BaseCutter):
             if d_p3 < d:
                 d = d_p3
                 cl = cl_p3
+            if cl:
+                return (cl,d)
+            (cl_p1,d_p1) = self.intersect_cylinder_vertex(direction, triangle.p1)
+            (cl_p2,d_p2) = self.intersect_cylinder_vertex(direction, triangle.p2)
+            (cl_p3,d_p3) = self.intersect_cylinder_vertex(direction, triangle.p3)
             if d_e1 < d:
                 d = d_e1
                 cl = cl_e1

pycam/Cutters/ToroidalCutter.py

@@ -181,26 +181,16 @@ class ToroidalCutter(BaseCutter):
 
     def intersect(self, direction, triangle):
         (cl_t,d_t) = self.intersect_torus_triangle(direction, triangle)
-        (cl_p1,d_p1) = self.intersect_torus_vertex(direction, triangle.p1)
-        (cl_p2,d_p2) = self.intersect_torus_vertex(direction, triangle.p2)
-        (cl_p3,d_p3) = self.intersect_torus_vertex(direction, triangle.p3)
-        (cl_e1,d_e1) = self.intersect_torus_edge(direction, Line(triangle.p1,triangle.p2))
-        (cl_e2,d_e2) = self.intersect_torus_edge(direction, Line(triangle.p2,triangle.p3))
-        (cl_e3,d_e3) = self.intersect_torus_edge(direction, Line(triangle.p3,triangle.p1))
         d = INFINITE
         cl = None
         if d_t < d:
             d = d_t
             cl = cl_t
-        if d_p1 < d:
-            d = d_p1
-            cl = cl_p1
-        if d_p2 < d:
-            d = d_p2
-            cl = cl_p2
-        if d_p3 < d:
-            d = d_p3
-            cl = cl_p3
+        if cl:
+            return (cl,d)
+        (cl_e1,d_e1) = self.intersect_torus_edge(direction, triangle.e1)
+        (cl_e2,d_e2) = self.intersect_torus_edge(direction, triangle.e2)
+        (cl_e3,d_e3) = self.intersect_torus_edge(direction, triangle.e3)
         if d_e1 < d:
             d = d_e1
             cl = cl_e1
@@ -210,13 +200,29 @@ class ToroidalCutter(BaseCutter):
         if d_e3 < d:
             d = d_e3
             cl = cl_e3
+        if cl:
+            return (cl,d)
+        (cl_p1,d_p1) = self.intersect_torus_vertex(direction, triangle.p1)
+        (cl_p2,d_p2) = self.intersect_torus_vertex(direction, triangle.p2)
+        (cl_p3,d_p3) = self.intersect_torus_vertex(direction, triangle.p3)
+        if d_p1 < d:
+            d = d_p1
+            cl = cl_p1
+        if d_p2 < d:
+            d = d_p2
+            cl = cl_p2
+        if d_p3 < d:
+            d = d_p3
+            cl = cl_p3
+        if cl:
+            return (cl,d)
         (cl_t,d_t) = self.intersect_circle_triangle(direction, triangle)
         (cl_p1,d_p1) = self.intersect_circle_vertex(direction, triangle.p1)
         (cl_p2,d_p2) = self.intersect_circle_vertex(direction, triangle.p2)
         (cl_p3,d_p3) = self.intersect_circle_vertex(direction, triangle.p3)
-        (cl_e1,d_e1) = self.intersect_circle_edge(direction, Line(triangle.p1,triangle.p2))
-        (cl_e2,d_e2) = self.intersect_circle_edge(direction, Line(triangle.p2,triangle.p3))
-        (cl_e3,d_e3) = self.intersect_circle_edge(direction, Line(triangle.p3,triangle.p1))
+        (cl_e1,d_e1) = self.intersect_circle_edge(direction, triangle.e1)
+        (cl_e2,d_e2) = self.intersect_circle_edge(direction, triangle.e2)
+        (cl_e3,d_e3) = self.intersect_circle_edge(direction, triangle.e3)
         if d_t < d:
             d = d_t
             cl = cl_t
@@ -242,9 +248,9 @@ class ToroidalCutter(BaseCutter):
             (cl_p1,d_p1) = self.intersect_cylinder_vertex(direction, triangle.p1)
             (cl_p2,d_p2) = self.intersect_cylinder_vertex(direction, triangle.p2)
             (cl_p3,d_p3) = self.intersect_cylinder_vertex(direction, triangle.p3)
-            (cl_e1,d_e1) = self.intersect_cylinder_edge(direction, Line(triangle.p1,triangle.p2))
-            (cl_e2,d_e2) = self.intersect_cylinder_edge(direction, Line(triangle.p2,triangle.p3))
-            (cl_e3,d_e3) = self.intersect_cylinder_edge(direction, Line(triangle.p3,triangle.p1))
+            (cl_e1,d_e1) = self.intersect_cylinder_edge(direction, triangle.e1)
+            (cl_e2,d_e2) = self.intersect_cylinder_edge(direction, triangle.e2)
+            (cl_e3,d_e3) = self.intersect_cylinder_edge(direction, triangle.e3)
             if d_p1 < d:
                 d = d_p1
                 cl = cl_p1

pycam/Exporters/SVGExporter.py

-from pycam.Geometry import *
-
 class SVGExporter:
 
     def __init__(self, filename):
         self.file = file(filename,"w")
         self.file.write("""<?xml version='1.0'?>
 <svg xmlns='http://www.w3.org/2000/svg' width='640' height='800'>
-<g transform='translate(320,320) scale(30)' stroke-width='0.01' font-size='0.4'>
+<g transform='translate(320,320) scale(60)' stroke-width='0.01' font-size='0.2'>
 """)
         self._fill = 'none'
         self._stroke = 'black'
@@ -29,11 +27,11 @@ class SVGExporter:
             x = -7
         if y < -1000:
             y = -7
-        l = "<circle fill='" + self._fill +"'" + (" cx='%g'" % x) + (" cy='%g'" % -y) + " r='0.1'/>\n"
+        l = "<circle fill='" + self._fill +"'" + (" cx='%g'" % x) + (" cy='%g'" % -y) + " r='0.04'/>\n"
         self.file.write(l)
 
     def AddText(self, x, y, text):
-        l = "<text fill='" + self._fill +"'" + (" x='%g'" % x) + (" y='%g'" % -y) + ">" + text + "</text>\n"
+        l = "<text fill='" + self._fill +"'" + (" x='%g'" % x) + (" y='%g'" % -y) + " dx='0.07'>" + text + "</text>\n"
         self.file.write(l)
         
 

pycam/Exporters/__init__.py

-list = [ "SimpleGCodeExporter"]
+list = [ "SimpleGCodeExporter", "SVGExporter" ]
 __all__ = list

pycam/Geometry/Model.py

@@ -69,7 +69,11 @@ class Model:
             self._maxsize = max3(max(abs(self.maxx),abs(self.minx)),max(abs(self.maxy),abs(self.miny)),max(abs(self.maxz),abs(self.minz)))
         return self._maxsize
 
-    def triangles(self):
+    def triangles(self, minx=-INFINITE,miny=-INFINITE,minz=-INFINITE,maxx=+INFINITE,maxy=+INFINITE,maxz=+INFINITE):
+        if minx==-INFINITE and miny==-INFINITE and minz==-INFINITE and maxx==+INFINITE and maxy==+INFINITE and maxz==+INFINITE:
+            return self._triangles
+        if hasattr(self, "t_kdtree"):
+            return self.t_kdtree.Search(minx,maxx,miny,maxy)
         return self._triangles
 
     def subdivide(self, depth):

pycam/Geometry/Point.py

@@ -9,7 +9,7 @@ class Point:
         self.y = float(y)
         self.z = float(z)
     def __repr__(self):
-        return "Point<%g,%g,%g>" % (self.x,self.y,self.z)
+        return "Point%d<%g,%g,%g>" % (self.id,self.x,self.y,self.z)
     def mul(self, c):
         return Point(self.x*c,self.y*c,self.z*c)
     def div(self, c):

pycam/Geometry/PointKdtree.py

+
+import math
+
+from Point import *
+from Line import *
+from Triangle import *
+from kdtree import *
+
+
+class PointKdtree(kdtree):
+    def __init__(self, points=[],cutoff=5, cutoff_distance=0.5, tolerance=0.001):
+        self._n = None
+        self.tolerance=tolerance
+        nodes = []
+        for p in points:
+            n = Node();
+            n.point = p
+            n.bound = []
+            n.bound.append(p.x)
+            n.bound.append(p.y)
+            n.bound.append(p.z)
+            nodes.append(n)
+        kdtree.__init__(self, nodes, cutoff, cutoff_distance)
+
+    def dist(self, n1, n2):
+        dx = n1.bound[0]-n2.bound[0]
+        dy = n1.bound[1]-n2.bound[1]
+        dz = n1.bound[2]-n2.bound[2]
+        return dx*dx+dy*dy+dz*dz
+
+    def Point(self, x, y, z):
+        #return Point(x,y,z)
+        if self._n:
+            n = self._n
+        else:
+            n = Node();
+        n.bound = []
+        n.bound.append(x)
+        n.bound.append(y)
+        n.bound.append(z)
+        (nn,dist) = self.nearest_neighbor(n, self.dist)
+        if nn and dist<self.tolerance:
+            self._n = n
+            return nn.p
+        else:
+            n.p = Point(x,y,z)
+            self._n = None
+            self.insert(n)
+            return n.p
+
+

pycam/Geometry/Triangle.py

 from Point import *
 from Plane import *
 from utils import *
+from Line import *
 
 ORIENTATION_CCW = 2
 ORIENTATION_CW  = 3
@@ -15,15 +16,27 @@ except:
 class Triangle:
     id = 0
     # points are expected to be in ClockWise order
-    def __init__(self, p1=None, p2=None, p3=None):
+    def __init__(self, p1=None, p2=None, p3=None, e1=None, e2=None, e3=None):
         self.id = Triangle.id
         Triangle.id += 1
         self.p1 = p1
         self.p2 = p2
         self.p3 = p3
+        if (not e1) and p1 and p2:
+            self.e1 = Line(p1,p2)
+        else:
+            self.e1 = e1
+        if (not e2) and p2 and p3:
+            self.e2 = Line(p2,p3)
+        else:
+            self.e2 = e2
+        if (not e3) and p3 and p1:
+            self.e3 = Line(p3,p1)
+        else:
+            self.e3 = e3
 
     def __repr__(self):
-        return "Triangle<%s,%s,%s>" % (self.p1,self.p2,self.p3)
+        return "Triangle%d<%s,%s,%s>" % (self.id,self.p1,self.p2,self.p3)
 
     def name(self):
         return "triangle%d" % self.id
@@ -49,7 +62,7 @@ class Triangle:
                 self._sphere = gluNewQuadric()
             gluSphere(self._sphere, self._radius, 10, 10)
             glPopMatrix()
-        if False: # draw triangle id on triangle face
+        if True: # draw triangle id on triangle face
             glPushMatrix()
             cc = glGetFloatv(GL_CURRENT_COLOR)
             c = self.center()
@@ -102,9 +115,6 @@ class Triangle:
             # calculate normal, if p1-p2-pe are in clockwise order
             self._normal = self.p3.sub(self.p1).cross(self.p2.sub(self.p1))
             denom = self._normal.norm()
-#            # TODO: fix kludge: make surface normal point up for now
-#            if self._normal.z < 0:
-#                denom = -denom
             self._normal = self._normal.div(denom)
         return self._normal
 

pycam/Geometry/TriangleKdtree.py

@@ -6,74 +6,42 @@ from Line import *
 from Triangle import *
 from kdtree import *
 
-def BuildKdtree2d(triangles, cutoff=3, cutoff_distance=1.0):
-    nodes = []
-    for t in triangles:
-        n = Node();
-        n.triangle = t
-        n.bound = []
-        n.bound.append(min(min(t.p1.x,t.p2.x),t.p3.x))
-        n.bound.append(min(min(t.p1.y,t.p2.y),t.p3.y))
-        n.bound.append(max(max(t.p1.x,t.p2.x),t.p3.x))
-        n.bound.append(max(max(t.p1.y,t.p2.y),t.p3.y))
-        nodes.append(n)
-    return kd_tree(nodes, cutoff, cutoff_distance)
-
-tests = 0
-hits = 0
 overlaptest=True
 
-def ResetKdtree2dStats(overlap=True):
-    global tests, hits, overlaptest
-    hits = 0
-    tests = 0
-    overlaptest = overlap
-
-def GetKdtree2dStats():
-    global tests, hits
-    return (hits, tests)
-
-
 def SearchKdtree2d(kdtree, minx, maxx, miny, maxy):
     if kdtree.bucket:
         triangles = []
         for n in kdtree.nodes:
             global tests, hits, overlaptest
-            tests += 1
             if not overlaptest:
-#                print "not testing overlap"
                 triangles.append(n.triangle)
-                hits += 1
             else:
                 if not (n.bound[0]>maxx
-                        or n.bound[1]>maxy
-                        or n.bound[2]<minx
+                        or n.bound[1]<minx
+                        or n.bound[2]>maxy
                         or n.bound[3]<miny):
                     triangles.append(n.triangle)
-                    hits += 1
         return triangles
     else:
-#        return SearchKdtree2d(kdtree.lo, minx, maxx, miny, maxy)+SearchKdtree2d(kdtree.hi, minx, maxx, miny, maxy)
-
         if kdtree.cutdim==0:
             if maxx<kdtree.minval:
                 return []
-            elif maxx<=kdtree.cutval:
+            elif maxx<kdtree.cutval:
                 return SearchKdtree2d(kdtree.lo, minx, maxx, miny, maxy)
             else:
                 return SearchKdtree2d(kdtree.lo, minx, maxx, miny, maxy)+SearchKdtree2d(kdtree.hi, minx, maxx, miny, maxy)
         elif kdtree.cutdim==1:
-            if maxy<kdtree.minval:
+            if minx>kdtree.maxval:
                 return []
-            elif maxy<=kdtree.cutval:
-                return SearchKdtree2d(kdtree.lo, minx, maxx, miny, maxy)
+            elif minx>kdtree.cutval:
+                return SearchKdtree2d(kdtree.hi, minx, maxx, miny, maxy)
             else:
                 return SearchKdtree2d(kdtree.lo, minx, maxx, miny, maxy)+SearchKdtree2d(kdtree.hi, minx, maxx, miny, maxy)
         elif kdtree.cutdim==2:
-            if minx>kdtree.maxval:
+            if maxy<kdtree.minval:
                 return []
-            elif minx>kdtree.cutval:
-                return SearchKdtree2d(kdtree.hi, minx, maxx, miny, maxy)
+            elif maxy<kdtree.cutval:
+                return SearchKdtree2d(kdtree.lo, minx, maxx, miny, maxy)
             else:
                 return SearchKdtree2d(kdtree.lo, minx, maxx, miny, maxy)+SearchKdtree2d(kdtree.hi, minx, maxx, miny, maxy)
         elif kdtree.cutdim==3:
@@ -83,3 +51,21 @@ def SearchKdtree2d(kdtree, minx, maxx, miny, maxy):
                 return SearchKdtree2d(kdtree.hi, minx, maxx, miny, maxy)
             else:
                 return SearchKdtree2d(kdtree.lo, minx, maxx, miny, maxy)+SearchKdtree2d(kdtree.hi, minx, maxx, miny, maxy)
+
+
+class TriangleKdtree(kdtree):
+    def __init__(self, triangles, cutoff=3, cutoff_distance=1.0):
+        nodes = []
+        for t in triangles:
+            n = Node();
+            n.triangle = t
+            n.bound = []
+            n.bound.append(min(min(t.p1.x,t.p2.x),t.p3.x))
+            n.bound.append(max(max(t.p1.x,t.p2.x),t.p3.x))
+            n.bound.append(min(min(t.p1.y,t.p2.y),t.p3.y))
+            n.bound.append(max(max(t.p1.y,t.p2.y),t.p3.y))
+            nodes.append(n)
+        kdtree.__init__(self, nodes, cutoff, cutoff_distance)
+
+    def Search(self, minx, maxx, miny, maxy):
+        return SearchKdtree2d(self, minx, maxx, miny, maxy)

pycam/Geometry/kdtree.py

 #!/usr/bin/env python
 
 import math
+import pycam
+
+try:
+    from OpenGL.GL import *
+    from OpenGL.GLUT import *
+    from OpenGL.GLU import *
+except:
+    pass
 
 class Node:
     def __repr__(self):
@@ -31,18 +39,90 @@ def find_max_spread(nodes):
     return (maxspreaddim, maxspread)
 
 
-class kd_tree:
+class kdtree:
     id = 0
 
     def __repr__(self):
         if self.bucket:
+            if True:
                 return "(#%d)" % (len(self.nodes))
             else:
-            return "(%f<=%s,%d:%f,%s<=%f)" % (self.minval,self.lo, self.cutdim,self.cutval,self.hi,self.maxval)
+                s = "("
+                for n in self.nodes:
+                    if (len(s))>1:
+                        s += ","
+                        s += str(n.p.id)
+                        s += ")"
+                return s
+        else:
+            return "(%s,%d:%g,%s)" % (self.lo,self.cutdim,self.cutval,self.hi)
+
+    def to_OpenGL(self,minx,maxx,miny,maxy,minz,maxz):
+        if self.bucket:
+            glBegin(GL_LINES)
+            glVertex3d(minx,miny,minz)
+            glVertex3d(minx,miny,maxz)
+            glVertex3d(minx,maxy,minz)
+            glVertex3d(minx,maxy,maxz)
+            glVertex3d(maxx,miny,minz)
+            glVertex3d(maxx,miny,maxz)
+            glVertex3d(maxx,maxy,minz)
+            glVertex3d(maxx,maxy,maxz)
+
+            glVertex3d(minx,miny,minz)
+            glVertex3d(maxx,miny,minz)
+            glVertex3d(minx,maxy,minz)
+            glVertex3d(maxx,maxy,minz)
+            glVertex3d(minx,miny,maxz)
+            glVertex3d(maxx,miny,maxz)
+            glVertex3d(minx,maxy,maxz)
+            glVertex3d(maxx,maxy,maxz)
+
+            glVertex3d(minx,miny,minz)
+            glVertex3d(minx,maxy,minz)
+            glVertex3d(maxx,miny,minz)
+            glVertex3d(maxx,maxy,minz)
+            glVertex3d(minx,miny,maxz)
+            glVertex3d(minx,maxy,maxz)
+            glVertex3d(maxx,miny,maxz)
+            glVertex3d(maxx,maxy,maxz)
+            glEnd()
+        elif self.dim==6:
+            if self.cutdim == 0 or self.cutdim == 2:
+                self.lo.to_OpenGL(minx,self.cutval,miny,maxy,minz,maxz)
+                self.hi.to_OpenGL(self.cutval,maxx,miny,maxy,minz,maxz)
+            elif self.cutdim == 1 or self.cutdim == 3:
+                self.lo.to_OpenGL(minx,maxx,miny,self.cutval,minz,maxz)
+                self.hi.to_OpenGL(minx,maxx,self.cutval,maxy,minz,maxz)
+            elif self.cutdim == 4 or self.cutdim == 5:
+                self.lo.to_OpenGL(minx,maxx,miny,maxx,minz,self.cutval)
+                self.hi.to_OpenGL(minx,maxx,miny,maxy,self.cutval,maxz)
+        elif self.dim==4:
+            if self.cutdim == 0 or self.cutdim == 2:
+                self.lo.to_OpenGL(minx,self.cutval,miny,maxy,minz,maxz)
+                self.hi.to_OpenGL(self.cutval,maxx,miny,maxy,minz,maxz)
+            elif self.cutdim == 1 or self.cutdim == 3:
+                self.lo.to_OpenGL(minx,maxx,miny,self.cutval,minz,maxz)
+                self.hi.to_OpenGL(minx,maxx,self.cutval,maxy,minz,maxz)
+        elif self.dim==3:
+            if self.cutdim == 0:
+                self.lo.to_OpenGL(minx,self.cutval,miny,maxy,minz,maxz)
+                self.hi.to_OpenGL(self.cutval,maxx,miny,maxy,minz,maxz)
+            elif self.cutdim == 1:
+                self.lo.to_OpenGL(minx,maxx,miny,self.cutval,minz,maxz)
+                self.hi.to_OpenGL(minx,maxx,self.cutval,maxy,minz,maxz)
+            elif self.cutdim == 2:
+                self.lo.to_OpenGL(minx,maxx,miny,maxy,minz,self.cutval)
+                self.hi.to_OpenGL(minx,maxx,miny,maxy,self.cutval,maxz)
 
     def __init__(self, nodes, cutoff, cutoff_distance):
-        self.id = kd_tree.id
-        kd_tree.id += 1
+        self.id = kdtree.id
+        self.bucket = False
+        if nodes and len(nodes)>0:
+            self.dim = len(nodes[0].bound)
+        else:
+            self.dim = 0
+        kdtree.id += 1
         self.cutoff = cutoff
         self.cutoff_distance = cutoff_distance
 
@@ -62,9 +142,65 @@ class kd_tree:
                 self.minval = nodes[0].bound[cutdim]
                 self.maxval = nodes[-1].bound[cutdim]
                 self.cutval = nodes[median].bound[cutdim]
-                self.lo = kd_tree(nodes[0:median], cutoff, cutoff_distance)
-                self.hi = kd_tree(nodes[median:], cutoff, cutoff_distance)
+                self.lo = kdtree(nodes[0:median], cutoff, cutoff_distance)
+                self.hi = kdtree(nodes[median:], cutoff, cutoff_distance)
 
+    def dist(self, n1, n2):
+        dist = 0
+        for i in range(0,len(n1.bound)):
+            d = n1.bound[i]-n2.bound[i]
+            dist += d*d
+        return dist
 
+    def nearest_neighbor(self, node, dist=dist):
+        if self.bucket:
+            if len(self.nodes)==0:
+                return (None, 0)
+            best = self.nodes[0]
+            bestdist = self.dist(node, best)
+            for n in self.nodes:
+                d = self.dist(n,node)
+                if d<bestdist:
+                    best = n
+                    bestdist = d
+            return (best,bestdist)
+        else:
+            if node.bound[self.cutdim] <= self.cutval:
+                (best,bestdist) = self.lo.nearest_neighbor(node, dist)
+                if bestdist > self.cutval - best.bound[self.cutdim]:
+                    (best2,bestdist2) = self.hi.nearest_neighbor(node, dist)
+                    if bestdist2<bestdist:
+                        return (best2,bestdist2)
+                return (best,bestdist)
+            else:
+                (best,bestdist) = self.hi.nearest_neighbor(node, dist)
+                if bestdist > best.bound[self.cutdim] - self.cutval:
+                    (best2,bestdist2) = self.lo.nearest_neighbor(node, dist)
+                    if bestdist2<bestdist:
+                        return (best2,bestdist2)
+                return (best,bestdist)
+    
+    def insert(self, node):
+        if self.dim==0:
+            self.dim = len(node.bound)
+
+        if self.bucket:
+            self.nodes.append(node)
+            if (len(self.nodes)>self.cutoff):
+                self.bucket = False
+                (cutdim,spread) = find_max_spread(self.nodes)
+                self.cutdim = cutdim
+                self.nodes.sort(cmp=lambda x,y:cmp(x.bound[cutdim],y.bound[cutdim]))
+                median = len(self.nodes)/2
+                self.minval = self.nodes[0].bound[cutdim]
+                self.maxval = self.nodes[-1].bound[cutdim]
+                self.cutval = self.nodes[median].bound[cutdim]
+                self.lo = kdtree(self.nodes[0:median], self.cutoff, self.cutoff_distance)
+                self.hi = kdtree(self.nodes[median:], self.cutoff, self.cutoff_distance)
+        else:
+            if node.bound[self.cutdim] <= self.cutval:
+                self.lo.insert(node)
+            else:
+                self.hi.insert(node)
         
 

pycam/Importers/STLImporter.py

 import re
 
 from pycam.Geometry import *
+from pycam.Geometry.PointKdtree import PointKdtree
+from pycam.Geometry.TriangleKdtree import TriangleKdtree
+
+vertices = 0
+edges = 0
+epsilon = 0.0001
+kdtree = None
+
+def UniqueVertex(x,y,z):
+    global vertices
+    if kdtree:
+        last = Point.id
+        p = kdtree.Point(x,y,z)
+        if p.id == last:
+#            print p
+            vertices += 1
+        return p
+    else:
+        vertices += 1
+        return Point(x,y,z)
+
+def UniqueEdge(p1, p2):
+    global edges
+    if hasattr(p1,"edges"):
+        for e in p1.edges:
+            if e.p1 == p1 and e.p2 == p2:
+                return e
+            if e.p2 == p1 and e.p1 == p2:
+                return e
+    edges += 1
+    e = Line(p1,p2)
+    if not hasattr(p1,"edges"):
+        p1.edges = [e]
+    else:
+        p1.edges.append(e)
+    if not hasattr(p2,"edges"):
+        p2.edges = [e]
+    else:
+        p2.edges.append(e)
+    return e
+
 
 def ImportModel(filename):
-    model = Model()
+    global vertices, edges, kdtree
+    vertices = 0
+    edges = 0
+    kdtree = None
+
     f = file(filename,"r")
     solid = re.compile("\s*solid\s+(\w+)\s+.*")
     solid_AOI = re.compile("\s*solid\s+\"([\w\-]+)\"; Produced by Art of Illusion.*")
@@ -15,12 +60,15 @@ def ImportModel(filename):
     endloop = re.compile("\s*endloop\s+")
     vertex = re.compile("\s*vertex\s+(?P<x>[-+]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)?)\s+(?P<y>[-+]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)?)\s+(?P<z>[-+]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)?)\s+")
 
+    kdtree = PointKdtree([],3,1,epsilon)
     model = Model()
+
     t = None
     p1 = None
     p2 = None
     p3 = None
     AOI = False
+    
     for line in f:
         m = solid_AOI.match(line)
         if m:
@@ -31,6 +79,7 @@ def ImportModel(filename):
         if m:
             model.name = m.group(1)
             continue
+
         m = facet.match(line)
         if m:
             m = normal.match(line)
@@ -46,10 +95,9 @@ def ImportModel(filename):
         m = vertex.match(line)
         if m:
             if AOI:
-                p = Point(float(m.group('x')),float(m.group('z')),float(m.group('y')))
+                p = UniqueVertex(float(m.group('x')),float(m.group('z')),float(m.group('y')))
             else:
-                p = Point(float(m.group('x')),float(m.group('y')),float(m.group('z')))
-            # TODO: check for duplicate points (using kdtree?)
+                p = UniqueVertex(float(m.group('x')),float(m.group('y')),float(m.group('z')))
             if p1 == None:
                 p1 = p
             elif p2 == None:
@@ -66,17 +114,26 @@ def ImportModel(filename):
         if m:
             # make sure the points are in ClockWise order
             if n.dot(p3.sub(p1).cross(p2.sub(p1)))>0:
-                t = Triangle(p1, p2, p3)
+                t = Triangle(p1, p2, p3, UniqueEdge(p1,p2), UniqueEdge(p2,p3), UniqueEdge(p3,p1))
             else:
-                t = Triangle(p1, p3, p2)
+                t = Triangle(p1, p3, p2, UniqueEdge(p1,p3), UniqueEdge(p3,p2), UniqueEdge(p2,p1))
             t._normal = n
             n=p1=p2=p3=None
-#            if t.normal().z < 0:
-#                continue
+            if t.normal().z < 0:
+                continue
             model.append(t)
             continue
         m = endsolid.match(line)
         if m:
             continue
+
+    model.p_kdtree = kdtree
+    model.t_kdtree = TriangleKdtree(model.triangles())
+
+    print "Imported STL model: ", vertices, "vertices,", edges, "edges,", len(model.triangles()),"triangles"
+    vertices = 0
+    edges = 0
+    kdtree = None
+
     return model
     

pycam/Importers/TestModel.py

@@ -5,27 +5,47 @@ def TestModel():
     points.append(Point(-2,1,4))
     points.append(Point(2,1,4))
     points.append(Point(0,-2,4))
-
     points.append(Point(-5,2,2))
     points.append(Point(-1,3,2))
     points.append(Point(5,2,2))
-
     points.append(Point(4,-1,2))
     points.append(Point(2,-4,2))
     points.append(Point(-2,-4,2))
     points.append(Point(-3,-2,2))
 
+    lines = []
+    lines.append(Line(points[0],points[1]))
+    lines.append(Line(points[1],points[2]))
+    lines.append(Line(points[2],points[0]))
+    lines.append(Line(points[0],points[3]))
+    lines.append(Line(points[3],points[4]))
+    lines.append(Line(points[4],points[0]))
+    lines.append(Line(points[4],points[1]))
+    lines.append(Line(points[4],points[5]))
+    lines.append(Line(points[5],points[1]))
+    lines.append(Line(points[5],points[6]))
+    lines.append(Line(points[6],points[1]))
+    lines.append(Line(points[6],points[2]))
+    lines.append(Line(points[6],points[7]))
+    lines.append(Line(points[7],points[2]))
+    lines.append(Line(points[7],points[8]))
+    lines.append(Line(points[8],points[2]))
+    lines.append(Line(points[8],points[9]))
+    lines.append(Line(points[9],points[2]))
+    lines.append(Line(points[9],points[0]))
+    lines.append(Line(points[9],points[3]))
+
     model = Model()
-    model.append(Triangle(points[0],points[1],points[2]))
-    model.append(Triangle(points[0],points[3],points[4]))
-    model.append(Triangle(points[0],points[4],points[1]))
-    model.append(Triangle(points[1],points[4],points[5]))
-    model.append(Triangle(points[1],points[5],points[6]))
-    model.append(Triangle(points[1],points[6],points[2]))
-    model.append(Triangle(points[2],points[6],points[7]))
-    model.append(Triangle(points[2],points[7],points[8]))
-    model.append(Triangle(points[2],points[8],points[9]))
-    model.append(Triangle(points[2],points[9],points[0]))
-    model.append(Triangle(points[0],points[9],points[3]))
+    model.append(Triangle(points[0],points[1],points[2],lines[0],lines[1],lines[2]))
+    model.append(Triangle(points[0],points[3],points[4],lines[3],lines[4],lines[5]))
+    model.append(Triangle(points[0],points[4],points[1],lines[5],lines[6],lines[0]))
+    model.append(Triangle(points[1],points[4],points[5],lines[6],lines[7],lines[8]))
+    model.append(Triangle(points[1],points[5],points[6],lines[8],lines[9],lines[10]))
+    model.append(Triangle(points[1],points[6],points[2],lines[10],lines[11],lines[1]))
+    model.append(Triangle(points[2],points[6],points[7],lines[11],lines[12],lines[13]))
+    model.append(Triangle(points[2],points[7],points[8],lines[13],lines[14],lines[15]))
+    model.append(Triangle(points[2],points[8],points[9],lines[15],lines[16],lines[17]))
+    model.append(Triangle(points[2],points[9],points[0],lines[17],lines[18],lines[2]))
+    model.append(Triangle(points[0],points[9],points[3],lines[18],lines[19],lines[3]))
     return model
 

pycam/PathGenerators/DropCutter.py

@@ -30,7 +30,8 @@ class DropCutter:
                     t_max = None
                     cl_last = None
                     self.cutter.moveto(p)
-                    for t in self.model.triangles():
+                    triangles = self.model.triangles(p.x-self.cutter.radius,p.y-self.cutter.radius,z0,p.x+self.cutter.radius,p.y+self.cutter.radius,+INFINITE)
+                    for t in triangles:
                         if t.normal().z < 0: continue;
                         cl = self.cutter.drop(t)
                         if cl and (cl.z > z_max or cl_max is None):
@@ -85,7 +86,8 @@ class DropCutter:
                     t_max = None
                     cl_last = None
                     self.cutter.moveto(p)
-                    for t in self.model.triangles():
+                    triangles = self.model.triangles(p.x-self.cutter.radius,p.y-self.cutter.radius,z0,p.x+self.cutter.radius,p.y+self.cutter.radius,+INFINITE)
+                    for t in triangles:
                         if t.normal().z < 0: continue;
                         cl = self.cutter.drop(t)
                         if cl and (cl.z > z_max or cl_max is None):

pycam/PathGenerators/PushCutter.py

@@ -138,7 +138,13 @@ class PushCutter:
             prev = Point(x,y,z)
             hits.append(Hit(prev, None, 0, None))
 
-            for t in model.triangles():
+            triangles = None
+            if dx==0:
+                triangles = model.triangles(minx-self.cutter.radius,y-self.cutter.radius,z,maxx+self.cutter.radius,y+self.cutter.radius,INFINITE)
+            else:
+                triangles = model.triangles(x-self.cutter.radius,miny-self.cutter.radius,z,x+self.cutter.radius,maxy+self.cutter.radius,INFINITE)
+
+            for t in triangles:
                 #if t.normal().z < 0: continue;
                 # normals point outward... and we want to approach the model from the outside!
                 n = t.normal().dot(forward)